我有一个 fastq 文件,例如 reads.fastq
。我有一个 7-mer
字符串列表。对于 reads.fastq
中的每次读取,我想检查它是否至少包含列表中的 7-mer
字符串之一。条件是,如果找到匹配项(汉明距离 ==0
),则读取内容将写入数组 chosen_reads
,并且来自 fastq 文件的下一个读取内容会匹配。如果未找到匹配项,则循环将继续,直到找到匹配项。输出数组由唯一的读取组成,因为一旦找到第一个匹配项,匹配循环就会终止。我编写了以下代码,但输出数组中的读取不是唯一的,因为报告了汉明距离为零的所有匹配项。请提出修改建议:
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)
最佳答案
当您当前的代码找到汉明距离为 0 的字符串时,您当前的代码不会跳出循环从 reads.fastq
读取下一个 read
,您应该使用标志来决定何时突破,并在需要突破时为该标志指定 True 值 -
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)
breakFlag = False
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)
breakFlag = True
break;
if breakFlag:
break;
您确定要将 x
附加到 chosen_reads
中,这似乎是错误的,为了获得唯一的匹配项,也许您应该附加 s2
字符串和匹配的 ds
代替,对吗?如果这就是您想要的,您可以将一个元组附加到 chosen_reads
中,如下所示,而不是当前的附加逻辑 -
chosen_reads.append((ds, s2))
关于python - 选择汉明距离为零的读数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30929484/