Получение прочитанных имен из файла fastq для определенного seq

Question

У меня есть файл fastq.Чтения имеют длину 75 б.п. Вот так выглядит файл fastq

@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=??;

У меня небольшая нуклеотидная последовательность около 30 п.н.

Допустим, это мой нуклеотид

CTGTTGGGGGCCCGTGC

Что я хочу сделать, это найти эту последовательность в файле fastq и извлечь соответствующее имя чтения, в котором существует последовательность. Таким образом, имя для чтения в этом случае будет

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

Также я хочу разрешить частоту несовпадения 1 в последовательности.

Любой сценарий или командная строка?

Спасибо

С уважением

Answer 1

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

Если по "коэффициенту несоответствия 1" вы хотите иметь возможность сопоставить любые 29 из этих 30 (например, CTG.TGGGGGCCCGTGC , это немного сложнее.

...

Эх, не намного сложнее

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