This is a test version of Biostars. For the public version, visit https://www.biostars.org.
comparing two fastq files

I have 2 fastq files F1.fastq and F2.fastq. F2.fastq is a smaller file which is a subset of reads from F1.fastq. I want reads in F1.fastq which ARE NOT in F2.fastq. The following python code does not seem to work. Can you suggest edits?

needed_reads = []
reads_array = []
chosen_array = []

for x in Bio.SeqIO.parse("F1.fastq","fastq"):
        reads_array.append(x)

for y in Bio.SeqIO.parse("F2.fastq","fastq"):
        chosen_array.append(y)

for y in chosen_array:
        for x in reads_array:

                if str(x.seq) != str(y.seq) : needed_reads.append(x)

output_handle = open("DIFF.fastq","w")
SeqIO.write(needed_reads,output_handle,"fastq")
output_handle.close()
next-gen sequencing
• 10,682 views • link
updated by Ram Ram • written by sumithrasank75 sumithrasank75

Instead of reinventing the wheel, why not use existing programs? cmpfastq is one such program.

• 0 views • link
written by komal.rathi komal.rathi

Try it with tiny input files (like 2 records each, with only one in common) and you should be able to track down the bug.

• 0 views • link
written by matted matted

as for the reason it does not work: your program checks each read in file 1 against each read in file 2 and adds it to output every it these two do not match, basically adding the same read many many times to the output

• 0 views • link
written by Istvan Albert Istvan Albert

Yes, this is what is happening, but I am unable to get a way to fix this.

• 0 views • link
written by sumithrasank75 sumithrasank75

Frankly you should look at the results others gave you since this code has many other problems as well regarding its scalability - may not work all that well on large files.

As for your code you need to store the data first in a set and match against it something like this (not tested just typing it in):

results = []
seen = set()
for s in Bio.SeqIO.parse("F1.fastq","fastq"):
    seen.add(str(s.seq))

for s in Bio.SeqIO.parse("F2.fastq","fastq"):
    if str(s.seq) not in seen:
        results.append(x)

output_handle = open("DIFF.fastq","w")
SeqIO.write(results,output_handle,"fastq")
output_handle.close()
• 0 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

Instead of reinventing the wheel, why not use something like:

join -t '\t' -v1 -1 2 -2 2 \
<(gunzip -c f1.gz | paste - - - - | sort -t '\t' -k2,2) \
<(gunzip -c f2.gz | paste - - - - | sort -t '\t' -k2,2)
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

1 answer

With BBMap:

filterbyname.sh in=F1.fastq names=F2.fastq out=filtered.fastq include=f

You can also use dedupe for this with the "uniqueonly" flag if the names are different and you want the unique reads by sequence.

• 0 views • link
updated by Ram Ram • written by Brian Bushnell Brian Bushnell

I need the solution to be within the python environment

• 0 views • link
written by sumithrasank75 sumithrasank75

Log in to answer this question.

More posts like this