Python:多个共识序列
Python: Multiple Consensus sequences
从 dna 序列列表开始,我必须在 return 中拥有所有可能的共识(结果
每个位置核苷酸频率最高的序列)序列。如果在某些位置核苷酸有
相同的最高频率,我必须获得最高频率的所有可能组合。
我还必须在 return 中包含配置文件矩阵(每个序列的每个核苷酸的频率矩阵)。
到目前为止,这是我的代码(但它 return 只有一个共识序列):
seqList = ['TTCAAGCT','TGGCAACT','TTGGATCT','TAGCAACC','TTGGAACT','ATGCCATT','ATGGCACT']
n = len(seqList[0])
profile = { 'T':[0]*n,'G':[0]*n ,'C':[0]*n,'A':[0]*n }
for seq in seqList:
for i, char in enumerate(seq):
profile[char][i] += 1
consensus = ""
for i in range(n):
max_count = 0
max_nt = 'x'
for nt in "ACGT":
if profile[nt][i] > max_count:
max_count = profile[nt][i]
max_nt = nt
consensus += max_nt
print(consensus)
for key, value in profile.items():
print(key,':', " ".join([str(x) for x in value] ))
TTGCAACT
C : 0 0 1 3 2 0 6 1
A : 2 1 0 1 5 5 0 0
G : 0 1 6 3 0 1 0 0
T : 5 5 0 0 0 1 1 6
(可以看到,在第四个位置,C和G得分相同,这意味着我必须获得两个共识序列)
是否可以修改此代码以获得
所有可能的序列,或者你能解释一下如何获得正确结果的逻辑(伪代码)吗?
非常感谢您!
我相信有更好的方法,但这是一个简单的方法:
bestseqs = [[]]
for i in range(n):
d = {N:profile[N][i] for N in ['T','G','C','A']}
m = max(d.values())
l = [N for N in ['T','G','C','A'] if d[N] == m]
bestseqs = [ s+[N] for N in l for s in bestseqs ]
for s in bestseqs:
print(''.join(s))
# output:
ATGGAACT
ATGCAACT
从 dna 序列列表开始,我必须在 return 中拥有所有可能的共识(结果 每个位置核苷酸频率最高的序列)序列。如果在某些位置核苷酸有 相同的最高频率,我必须获得最高频率的所有可能组合。 我还必须在 return 中包含配置文件矩阵(每个序列的每个核苷酸的频率矩阵)。
到目前为止,这是我的代码(但它 return 只有一个共识序列):
seqList = ['TTCAAGCT','TGGCAACT','TTGGATCT','TAGCAACC','TTGGAACT','ATGCCATT','ATGGCACT']
n = len(seqList[0])
profile = { 'T':[0]*n,'G':[0]*n ,'C':[0]*n,'A':[0]*n }
for seq in seqList:
for i, char in enumerate(seq):
profile[char][i] += 1
consensus = ""
for i in range(n):
max_count = 0
max_nt = 'x'
for nt in "ACGT":
if profile[nt][i] > max_count:
max_count = profile[nt][i]
max_nt = nt
consensus += max_nt
print(consensus)
for key, value in profile.items():
print(key,':', " ".join([str(x) for x in value] ))
TTGCAACT
C : 0 0 1 3 2 0 6 1
A : 2 1 0 1 5 5 0 0
G : 0 1 6 3 0 1 0 0
T : 5 5 0 0 0 1 1 6
(可以看到,在第四个位置,C和G得分相同,这意味着我必须获得两个共识序列)
是否可以修改此代码以获得 所有可能的序列,或者你能解释一下如何获得正确结果的逻辑(伪代码)吗?
非常感谢您!
我相信有更好的方法,但这是一个简单的方法:
bestseqs = [[]]
for i in range(n):
d = {N:profile[N][i] for N in ['T','G','C','A']}
m = max(d.values())
l = [N for N in ['T','G','C','A'] if d[N] == m]
bestseqs = [ s+[N] for N in l for s in bestseqs ]
for s in bestseqs:
print(''.join(s))
# output:
ATGGAACT
ATGCAACT