Question: Searching for the common sequences in multiple (>2) fastq files
1
gravatar for BPors
2.3 years ago by
BPors40
BPors40 wrote:

Hi,

I am trying to see that if my fastq files ( 15 of them) coming from different sources are sharing the same or very similar reads (allowing 1-2 mismatches), and then I would like to grab each whole sequence that they were sharing ( Even though they are sharing 50 bp, I want to grab the whole read, which is 150 bp).

I checked Biopieces and Bowtie2, but I could not find multiple comparisons.

I would appreciate if you can share your thoughts about it.

Thank you!

common align sequence fastq • 999 views
ADD COMMENTlink modified 2.3 years ago by guillaume.rbt750 • written 2.3 years ago by BPors40

Hi,

To do that I would align the reads from all the fastq you have on your reference genome.

Then I would detect regions where a minimum coverage is present on all of 15 bams, then extract the reads from those regions.

If they align all to the same regions, they should be similar.

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by guillaume.rbt750

Thank you for your reply! I would like to try this approach, however can you explain why should I look for minimum coverage regions? Because they(the common sites) would be more rare than the others? I couldnt understand, sorry

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by BPors40

The regions where you will have a minimum coverage on all bams will have matching reads from all the fastq.

For example if you take a threshold of 5X, for a given position on the genome, at least 5 reads on each fastq is mapped on this position.

As all those reads will be matching the same regions, you can deduct that those reads are similar between each others.

ADD REPLYlink written 2.3 years ago by guillaume.rbt750

True..Many thanks for the explanation. Would you suggest the usage of samtools ?

ADD REPLYlink written 2.3 years ago by BPors40

To get the coverage I personally use bedtools :

bedtools genomecov -d -ibam your_bam_1.bam -g your_genome.fa > coverage_1.txt &

ADD REPLYlink written 2.3 years ago by guillaume.rbt750

I doubt it'd be fast, but I don't know how big your fastqs are, but you could parse all the reads in through BioPython. You can then just use python's string operons such as:

if mystring in fastq.seq:
    print(fastq.seq)

Not sure off the top of my head how to deal with allowing a couple of mismatches, but I'm sure some sort of module would allow it (maybe pairwise2), but that would be a LOT of aligning which would probably make it even slower.

I definitely think you should consider assembling the reads first though and then just find common contigs or something, since you may get a huge, huge dataset from comparing mutiple fastqs.

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by Joe15k

I am trying to see that if my fastq files ( 15 of them) coming from different sources are sharing the same or very similar reads (allowing 1-2 mismatches),

and

and then I would like to grab each whole sequence that they were sharing ( Even though they are sharing 50 bp, I want to grab the whole read, which is 150 bp).

Those seem to contradict each other, as if there is only a 50bp overlap, the reads are not shared. Do you actually want reads that share kmers? Or reads that come from the same organism? Perhaps you could explain exactly what you are trying to do, and why; otherwise "similar reads" is a pretty vague requirement.

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by Brian Bushnell17k

Sure thing. What I meant is, my reads are 150 bp each in all fastq files, and I want to get the reads which are shared by other fastq files. The ultimate goal would be to grab the ones which are matching in more than 140 bps. However, I can put this threshold to 50 bp. So, if the read1 from fileA.fastq is sharing 50 bp in read12 of fileB.fastq and read56 and read78 of fileC.fastq, I want to retrieve the read1.

ADD REPLYlink written 2.3 years ago by BPors40

You could do this approximately with BBDuk:

bbduk.sh in=fileA.fq ref=fileB.fq,fileC.fq outm=matched.fq k=50 mm=f

That will give you all of the reads in fileA that contain at least 50bp contiguous regions totally covered by 31-mers in fileV or fileC. It's not exactly the same as sharing 50-mers, but very similar.

ADD REPLYlink written 2.2 years ago by Brian Bushnell17k

Just to throw some fresh perspective:

I would produce a consensus FASTA separately for both and then align one to the other to check the alignment percent. I believe that SAMtools allows you to produce a consensus FASTA via pileup or mpileup.

You could also go back to the first bioinformatics tools like BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and compare them there. BLAST allows multiple comparisons. Just go to Nucleotide BLAST and check the 'Align two or more sequences' checkbox - it allows you to upload FASTA.

Kevin

ADD REPLYlink written 2.2 years ago by Kevin Blighe52k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 697 users visited in the last hour