选择与汉明距离零读取

我有一个fastq文件，说reads.fastq。我有一个7-mer字符串的列表。对于reads.fastq中的每个读取，我想检查它是否至少包含列表中的7-mer字符串中的一个。条件是，如果发现匹配（hamming distance ==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)

来源

2015-06-19 Ssank

您目前的代码不会跳出循环读取下一个read从reads.fastq时，它已经发现了汉明距离0的字符串，则应该使用标志来决定何时爆发，并指定该标志的真值，当你需要打破 -

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))

来源

2015-06-19 04:01:47

我删除的spead最后，如果breakFlag：打破;因为这不是通过读取迭代，其他建议工作。我希望将整个记录写入chosen_reads以备后用 – Ssank

哦，对不起，不，你不需要最后一个'如果breakFlag'，将会修复答案。 –

如果我理解你的要求，海明距离为试图找到的至少一个3“精确选择”字符串。迭代，因为你正在做的很慢，试图爆发可能是丑陋的。

我也许会建议a regex是什么在这里会有所帮助。您可以自动创建匹配的字符串：

import re 
chosen_re = re.compile('|'.join(l_chosen)) 

chosen_reads = [x for x in reads_array if chosen_re.search(str(s.seq))]

你会很很难击败正则表达式引擎

来源

2015-06-19 04:21:08

选择与汉明距离零读取

回答

相关问题