我有这种格式的数据:
> ABC12
ATCGACAG
> def34
ACCGACG
等等
我已存储的每个基因与具有>开始作为值的线的字典。 所以,字典是一样的东西{“ABC12”:“ATCGACAG”等}
现在,我希望能够将每个基因进行比较,从而计算在每个站点的,T的,C的,或G的数量。
我能想出的唯一的事情就是要打破字典成每个核苷酸网站列表,并使用拉链()与抗衡。 这是最好的方式,如果是这样,我怎么破解字典到每个站点列表?
我有这种格式的数据:
> ABC12
ATCGACAG
> def34
ACCGACG
等等
我已存储的每个基因与具有>开始作为值的线的字典。 所以,字典是一样的东西{“ABC12”:“ATCGACAG”等}
现在,我希望能够将每个基因进行比较,从而计算在每个站点的,T的,C的,或G的数量。
我能想出的唯一的事情就是要打破字典成每个核苷酸网站列表,并使用拉链()与抗衡。 这是最好的方式,如果是这样,我怎么破解字典到每个站点列表?
使用collections.Counter
:
>>> from collections import Counter
>>> Counter('ATCGACAG')
Counter({'A': 3, 'C': 2, 'G': 2, 'T': 1})
>>> Counter('ACCGACG')
Counter({'C': 3, 'A': 2, 'G': 2})
有没有理由不使用Biopython?
from Bio import AlignIO
alignment =AlignIO.read("alignment.fas", "fasta")
n=0
while n<len(alignment[0]):
A=alignment[:,n].count('A')
C=alignment[:,n].count('C')
G=alignment[:,n].count('G')
T=alignment[:,n].count('T')
gap=alignment[:,n].count('-')
print "at position %s there are %s A's, %s C's, %s G's, %s T's and %s gaps" % (n, A, C, G, T, gap)
n=n+1
确保你有一个真正的比对(即序列长度相同)。
PS我道歉print语句的丑陋格式...
这将返回
at position 0 there are 2 A's, 0 C's, 0 G's, 0 T's and 0 gaps
at position 1 there are 0 A's, 1 C's, 0 G's, 1 T's and 0 gaps
at position 2 there are 0 A's, 2 C's, 0 G's, 0 T's and 0 gaps
at position 3 there are 0 A's, 0 C's, 2 G's, 0 T's and 0 gaps
at position 4 there are 2 A's, 0 C's, 0 G's, 0 T's and 0 gaps
at position 5 there are 0 A's, 2 C's, 0 G's, 0 T's and 0 gaps
at position 6 there are 1 A's, 0 C's, 0 G's, 0 T's and 1 gaps
at position 7 there are 0 A's, 0 C's, 2 G's, 0 T's and 0 gaps
s1 = "ATCGACAG"
s2 = "ACCGACG"
alignment = [s1[i] if s1[i] == s2[i] else "-" for i in range(len(min([s1,s2],key=len)))]
print "".join(alignment)
A-CGAC-