Retrival of read names from fastq file for a particular seq

Question

I have a fastq file.The reads are of length 75bp. This is how a fastq file looks like

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

I have a small nucleotide sequence of around 30bp.

Lets say this is my nucleotide seq

CTGTTGGGGGCCCGTGC

What I want to do is search for this sequence in the fastq file and extract the corresponding read name in which the sequence exists. So the read name in this case would be

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

Also I want to allow a mismatch rate of 1 in the sequence.

Any script or command line ?

Thanks

Regards

For those of us unfamiliar with fastq, please provide a small sample of input data, and what you want the output to look like. — glenn jackman, Jan 07 '14 at 20:59
@glenn : i added the format, let me know if you still don't get it Thanks. — user3138373, Jan 07 '14 at 21:39

glenn jackman · Accepted Answer · 2014-01-07T22:22:57.807

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

If, by "mismatch rate of 1", you want to be able to match any 29 out of those 30 (such as CTG.TGGGGGCCCGTGC, that's quite a bit more complex.

...

Eh, not that much more complex:

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

This works i might ask you further development in the code once i understand it. Thanks — user3138373, Jan 07 '14 at 22:43

Retrival of read names from fastq file for a particular seq

1 Answers1