Python：构建共识序列

Question

我想从 python 中的几个序列构建一个共识序列，我正在寻找最有效/最 pythonic 的方法来实现这一点。

我有一个这样的字符串列表：

sequences = ["ACTAG", "-TTCG", "CTTAG"]

我还有这样一个字母表：

alphabet = ["A", "C", "G", "T"]

和这样的位置频率矩阵：

如果一个字符在某个位置出现次数最多，则此字符为一致序列。

此外，当 2 个或更多字符在同一位置出现相同的情况时，会出现其他字符（在此示例中，位置 0 => A 或 C = M，请参阅 IUPAC Codes）

因此，我的示例的预期共识序列是“MTTAG”。

编辑：

根据给定的字母表和位置频率矩阵，获得此一致序列的最有效/最pythonic 方法是什么？

Answer 1

如果您已经有了位置频率矩阵，您可以将其作为 pandas DataFrame 处理。我选择了它的方向，使字母表成为索引（注意最后的 transpose 调用）：

freq = pd.DataFrame([[1, 1, 0, 0], [0, 1, 0, 2], [0, 0, 0, 3], [2, 1, 0, 0], [0, 0, 3, 0]], columns=['A', 'C', 'G', 'T']).transpose()

给予

   0  1  2  3  4
A  1  0  0  2  0
C  1  1  0  1  0
G  0  0  0  0  3
T  0  2  3  0  0

您只想查看最常见的核苷酸：

most_common = freq[freq == freq.max(axis=0)]

给予

     0    1    2    3    4
A  1.0  NaN  NaN  2.0  NaN
C  1.0  NaN  NaN  NaN  NaN
G  NaN  NaN  NaN  NaN  3.0
T  NaN  2.0  3.0  NaN  NaN

然后创建一个函数，根据 IUPAC 代码从上述矩阵的单列中确定共识：

codes = {
    'A': 'A', 'C': 'C', 'G': 'G', 'T': 'T', 
    'AG': 'R', 'CT': 'Y', 'CG': 'S', 'AT': 'W', 'GT': 'K', 'AC': 'M', 
    'CGT': 'B', 'AGT': 'D', 'ACT': 'H', 'ACG': 'V', 
    'ACGT': 'N'
}
def freq_to_code(pos):
    top_nucs = pos.dropna().index
    key = ''.join(sorted(top_nucs))
    return codes[key]

将该函数应用于每一列并形成一个字符串以获得最终结果：

consensus = most_common.apply(freq_to_code, axis=0)
print(''.join(consensus))

给出 MTTAG

Python：构建共识序列

Python: build consensus sequence

python

sequence

dna-sequence