Retratação dos nomes lidos do arquivo fastq para um determinado seq

Question

Eu tenho um arquivo fastq. As leituras são de comprimento 75bp. É assim que um arquivo fastq se parece com

@ERR030899.2391 HWI-BRUNOP16X_0001:4:1:16426:3738#0/1
NNCGTGTCAGTGGCCAGCAGCCCACACTGCGCATGGCTGATACTGTGGACCCCCTGGACTGGCTTTTTGGGGAGT
+
###########################################################################
@ERR030899.2392 HWI-BRUNOP16X_0001:4:1:16582:3744#0/1
NNAAAGTCCTGCGCTGCGGAGGACAGGAAGCACCCCCTCAATAGCCAGCACCCACAGTGAGCTGAGCACTTACAG
+
##'(''((((5/333**+)(10-11+1,,,,1/1/F<<<<FF:FFFFFF=FFFFFFFFFFFFFFFFFFFFFFFF#
@ERR030899.2393 HWI-BRUNOP16X_0001:4:1:16911:3745#0/1
NNGGATTTAGCGGGGTGATGCCTGTTGGGGGCCCGTGCCCTCCTACTTGGGGGGCAGGGGCTAGGCTGCAGAGGT
+
###########################################################################
@ERR030899.2394 HWI-BRUNOP16X_0001:4:1:18194:3739#0/1
NNACAAGCAATTTAGTGATAATGTCCAGAGGGCCAAGGATGCGGACCACCTTTTGCAGAACTCATATCTCGAGCA
+
##*+*)'+++FFFFFFFFFF58588==?:?FFFFFFFFFFFFFF<FFFFFFFFFF=FFFFFFFFFFFFFF6=??;

Eu tenho uma pequena sequência nucleotídica de cerca de 30 pb.

Vamos dizer que este é o meu nucleotídeo seq

CTGTTGGGGGCCCGTGC

O que eu quero fazer é procurar por essa seqüência no arquivo fastq e extrair o nome de leitura correspondente no qual a seqüência existe. Então, o nome de leitura, neste caso, seria

@ERR030899.2393 HWI-BRUNOP16X_0001:4:1:16911:3745#0/1

Além disso, desejo permitir uma taxa de incompatibilidade de 1 na sequência.

Qualquer script ou linha de comando?

Obrigado

Atenciosamente

    

Answer 1

awk -v seq="CTGTTGGGGGCCCGTGC" '
  NR%4 == 1 {name = $0}
  NR%4 == 2 && index($0, seq) {print name}
' filename

Se, por "taxa de incompatibilidade de 1", você quiser ser capaz de corresponder 29 dessas 30 (como CTG.TGGGGGCCCGTGC , isso é um pouco mais complexo.

...

Eh, não muito mais complexo:

awk -v seq="CTGTTGGGGGCCCGTGC" '
  NR%4 == 1 {name=$0}
  NR%4 == 2 {
    if (index($0, seq))
      print "found seq \"" seq "\" in " name
    else
      for (i=1; i<=length(seq); i++) {
        patt = substr(seq, 1, i-1) "." substr(seq, i+1)
        if (match($0, patt)) {
          print "found pattern \"" patt "\" in " name
          break
        }
      }
  }
' filename