我目前正在处理一些DNA序列数据，并且需要为每个位点创建一个频率矩阵。例如，如下所示：
A T G C
0.2 0.3 0.3 0.2
0.3 0.4 0.1 0.2
0.7 0.1 0.1 0.1
输入是许多DNA序列的列表，例如：
te_seqs = ["ATCTACTGATG", "ATACAGTACATAGA", "ATAGACAGTTGTGCG", "GTCGATACGT", ...]
每个序列长数千个字符，并且有数千个序列。输出是一个numpy矩阵，其中包含每个站点的频率计数。例如，在上面的数据中，第一个站点具有3个As和1个G，因此A的频率为3/4 = 0.75，G的频率为1/4 = 0.25。T和C的频率均为0。
我目前拥有所有序列的列表，并且通过将1加到Nx4矩阵来获得频率。问题是我正在处理大量序列，并且嵌套for循环在时间上并不理想：
for seq in te_seqs:
for i,nuc in enumerate(seq):
if nuc == "A":
te_pwm[i, 0] = te_pwm[i, 0] + 1
elif nuc == "T":
te_pwm[i, 1] = te_pwm[i, 1] + 1
elif nuc == "G":
te_pwm[i, 2] = te_pwm[i, 2] + 1
elif nuc == "C":
te_pwm[i, 3] = te_pwm[i, 3] + 1
for seq in gene_seqs:
for i,nuc in enumerate(seq):
if nuc == "A":
gene_pwm[i, 0] = gene_pwm[i, 0] + 1
elif nuc == "T":
gene_pwm[i, 1] = gene_pwm[i, 1] + 1
elif nuc == "G":
gene_pwm[i, 2] = gene_pwm[i, 2] + 1
elif nuc == "C":
gene_pwm[i, 3] = gene_pwm[i, 3] + 1
我的问题是1）检查字符串列表中的每个字符串是否还有更Python的方式？2）是否有更好的方法来创建基频矩阵？
谢谢！
您可以用于itertools.izip_longest()遍历站点，然后用于collections.Counter进行计数：
import collections
import itertools
te_seqs = ["ATCTACTGATG", "ATACAGTACATAGA", "ATAGACAGTTGTGCG", "GTCGATACGT"]
sites = map(collections.Counter, itertools.izip_longest(*te_seqs))
for site in sites:
A = site.get("A", 0)
T = site.get("T", 0)
G = site.get("G", 0)
C = site.get("C", 0)
total = float(A + T + G + C)
print A / total, T / total, G / total, C / total
这产生
0.75 0.0 0.25 0.0
0.0 1.0 0.0 0.0
0.5 0.0 0.0 0.5
0.0 0.25 0.5 0.25
1.0 0.0 0.0 0.0
0.0 0.25 0.25 0.5
0.5 0.5 0.0 0.0
...
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句