This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How can I find the position of a sequence in a fastq file

Hello, I have a fastq file with 4 lines per sequence. This is an example of one of the entries:

@m54284U_200831_121835/3/ccs
ATGGGTGGGCTGGCAGTAGCCAGGGACGATGGGCTCTTCTCTGGGGATCCCAACTGGTTTCCTAAGAAATCCAAGGAGAACCCTCGGAACTTCTCGGACAACCAGTTGCAAGAGGGCAAGAACGTGATTGGGTTGCAGATGGGCACCAACCGTGGAGCATCTCAGGCCGGCATGACCGGCTATGGGATGCCACGGCAGATCCTCTGATCATACTCTCTCTCCTTCCCCTGCCCTCCATGAATGGTTAATATATATGTATATATATGTTTTAGCAGACATTCCCTGAGAGCCCCTGGATTGCTGAACCCCCCTCTGCCAGGGTCCAGGCCAGCCTATCTTGTCACCACTGGCAGGGCCTGATAATTGCCTCTCTCTCTCTCTCTCTTTCTCTCTCTCTCTCTCTCTGGGCTTACTAATGCATTCCTCCCCCCACATTATTCCCACAGTCTCAAGCACGTGGATTCTGCTGTAGTCGTACGCCGATGCGAAACATCGGCCACGTCGCTATTGCAGCGAGTAGATCGGAAGAGCACACGTCTGAACTCCA
+
5555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555

I have the following sequence: CAGACGTGTGCTCTTCCGATCT. I would like to know:

1) In which entries this sequence is present (both the number of entries and which entries).

2) Where in the sequence does this string match (which nucleotide position within the entry).

For example, I would like something like: The sequence is present in all entries. In entry 1, the sequence starts at nucleotide 10. etc.

I have tried with grep and awk but I certainly lack the knowledge to do this. Any help would be so appreciated!!!

fastq sequence position find
• 2,125 views • link
updated by GenoMax GenoMax • written by pomodoro_sinensis • pomodoro_sinensis

1 answer

You can use bbduk.sh in filter mode for question 1.

bbduk.sh -Xmx2g in=file.fq literal=CAGACGTGTGCTCTTCCGATCT k=10 outm=file.fq

Capture STDERR (where logs are written), they will show you how many sequences had that sequence. A guide for bbduk is available.

• 0 views • link
written by GenoMax GenoMax

Thank you! I used your command but I am getting the following error:

!/content/bbmap/bbduk.sh -Xmx2g in=small.fastq literal=CAGACGTGTGCTCTTCCGATCT k=10 outm=small.fastq

java -ea -Xmx2g -Xms2g -cp /content/bbmap/current/ jgi.BBDuk -Xmx2g in=small.fastq literal=CAGACGTGTGCTCTTCCGATCT k=10 outm=small.fastq
Error: Unable to initialize main class jgi.BBDuk
Caused by: java.lang.NoClassDefFoundError: stream/ConcurrentReadStreamInterface
• 0 views • link
written by pomodoro_sinensis • pomodoro_sinensis

Do you have java (1.8 and above) installed? You just unzipped the bbmap software bundle and did not move any files/folders inside it anywhere. You can use conda to install bbmap otherwise.

• 0 views • link
written by GenoMax GenoMax

That worked thank you!! I am getting the following result in a reduced version with only 250 reads:

0.028 seconds.
Initial:
Memory: max=2058m, total=2058m, free=2014m, used=44m

Added 13 kmers; time:   0.005 seconds.
Memory: max=2058m, total=2058m, free=2003m, used=55m

Input is being processed as unpaired
Started output streams: 0.015 seconds.
Processing time:        0.148 seconds.

Input:                      250 reads       250791 bases.
Contaminants:               134 reads (53.60%)  127677 bases (50.91%)
Total Removed:              134 reads (53.60%)  127677 bases (50.91%)
Result:                     116 reads (46.40%)  123114 bases (49.09%)

Time:                           0.170 seconds.
Reads Processed:         250    1.47k reads/sec
Bases Processed:        250k    1.47m bases/sec

I am assuming the "contaminants" are my sequence of interest?

• 0 views • link
written by pomodoro_sinensis • pomodoro_sinensis

You can easily check that. Add rcomp=f if you only want to search in the reads in forward direction.

• 0 views • link
written by GenoMax GenoMax

This is great thanks so much!! Do you know if this tool allows for checking where in the read this sequence was found? Or do you know of any other tool? I played around with grep but couldn't make it work...

• 0 views • link
written by pomodoro_sinensis • pomodoro_sinensis

Not easily I would think. You will need to parse alignments (CIGAR strings) or simply convert the sequence you find to fasta and do a multiple-sequence alignment.

• 0 views • link
written by GenoMax GenoMax

Log in to answer this question.

More posts like this