найти штрих-код ДНК с несоответствиями в последовательности

Question

У меня есть 36-nt, как это: atcttgttcaatggccgatcXXXXgtcgacaatcaa в файле fastq с XXXX - это разные штрих-коды. Я хочу искать штрих-код в файле в точном положении (от 21 до 24) и печатать последовательности с точностью до 3 несоответствий в последовательности, а не штрих-код.

Например: У меня есть штрих-код: aacg поиска, что штрих-код между положением 21 до 24 в fastq файл с разрешением 3 несовпадений в последовательности, как:

atcttgttcaatggccgatcaacggtcgacaatcaC# it has 1 mismatch 
ttcttgttcaatggccgatcaacggtcgacaatcaC# it has 2 mismatch 
tccttgttcaatggccgatcaacggtcgacaatcaC# it has 3 mismatch 

Я пытался найти уникальные строки первого использования awk и искать несоответствия, но мне очень скучно искать и находить их.

awk 'NR%4==2' 1.fq |sort|uniq -c|awk '{print $1"\t"$2}' > out1.txt 

Есть ли какой-либо быстрый способ найти?

спасибо.

Answer 1

Использование Python:

strs = "atcttgttcaatggccgatcaacggtcgacaatcaa" 

with open("1.fq") as f: 
    for line in f: 
     if line[20:24] == "aacg": 
      line = line.strip() 
      mismatches = sum(x!=y for x, y in zip(strs, line)) 
      if mismatches <= 3: 
       print line, mismatches 

atcttgttcaatggccgatcaacggtcgacaatcac 1 
ttcttgttcaatggccgatcaacggtcgacaatcac 2 
tccttgttcaatggccgatcaacggtcgacaatcac 3 

Answer 2

Использование Python:

import re 
seq="atcttgttcaatggccgatcaacggtcgacaatcaa" 
D = [ c for c in seq ] 
with open("input") as f: 
    for line in f: 
     line=line.rstrip('\n') 
     if re.match(".{20}aacg", line): 
      cnt = sum([ 1 for c,d in zip(line,D) if c != d]) 
      if cnt < 4: 
       print cnt, line 

Answer 3

Использование модуля Python регулярное выражение позволяет указать количество несовпадений

import regex #intended as a replacement for re 
from Bio import SeqIO 
import collections 

d = collections.defaultdict(list) 
motif = r'((atcttgttcaatggccgatc)(....)(gtcgacaatcaa)){e<4}' #e<4 = less than 4 errors 
records = list(SeqIO.parse(open(infile), "fastq")) 
for record in records: 
    seq = str(record.seq) 
    match = regex.search(motif, seq, regex.BESTMATCH) 
    barcode = match.group(3) 
    sequence = match.group(0) 
    d[barcode].append(sequence) # store as a dictionary key = barcode, value = list of sequences 
for k, v in d.items(): 
    print("barcode = %s" % (k)) 
    for i in v: 
     print("sequence = %s" % (i)) 

с использованием групп захвата, четвертая группа (3), будет штрих-код