我有一个fastq文件,比如reads.fastq
。我有一个7-mer
字符串的列表。对于每次读入reads.fastq
,我想检查它是否至少包含列表中的一个7-mer
字符串。条件是,如果找到匹配项(hamming distance ==0
),则读操作将写入一个数组chosen_reads
,并且从fastq文件中读取的下一个匹配。如果未找到匹配项,则循环将继续,直到找到匹配项为止。输出数组由唯一的读取组成,因为一旦找到第一个匹配项,匹配循环就会终止。我写了下面的代码,但是输出数组中的读取不是唯一的,因为所有与hamming距离为零的匹配都会被报告。请建议编辑:
def hamming(s1, s2):
#Return the Hamming distance between equal-length sequences
if len(s1) != len(s2):
raise ValueError("Undefined for sequences of unequal length")
return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))
for x in Bio.SeqIO.parse("reads.fastq","fastq"):
reads_array.append(x)
nmer = 7
l_chosen = ['gttattt','attattt','tgctagt']
chosen_reads = []
for x in reads_array:
s2 = str(x.seq)
for s in [s2[i:i+nmer] for i in range(len(s2)-nmer-1)]:
for ds in l_chosen:
dist = hamming(ds,s)
if dist == 0:
print s2, s,ds,dist
chosen_reads.append(x)
当前代码不会跳出循环从
reads.fastq
读取下一个read
当它找到一个汉明距离为0的字符串时,应该使用标志来决定何时中断,并在需要中断时为该标志指定真值-您确定要将
^{pr2}$x
附加到chosen_reads
中,这似乎是错误的,为了获得唯一的匹配,您应该添加s2
字符串和匹配的ds
对吗?如果这是您想要的,您可以像下面这样将一个元组附加到chosen_reads
中,而不是当前的附加逻辑-如果我明白你在问什么,汉明距离正在试图找到至少一个“选择”字符串准确。你正在做的迭代是缓慢的,尝试突破可能是丑陋的。在
我可能会建议a regex是什么会有帮助。您可以自动创建匹配字符串:
你将很难击败正则表达式引擎的速度
相关问题 更多 >
编程相关推荐