comparing two fastq files

Question

[sumithrasank75]

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()

[komal.rathi]

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

[matted]

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.

[Istvan Albert]

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

[sumithrasank75]

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

[Istvan Albert]

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()

[Pierre Lindenbaum]

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)

[Pierre Lindenbaum]

Hello sumithrasank75!

It appears that your post has been cross-posted to another site: http://stackoverflow.com/questions/30942734/getting-records-which-are-different-from-two-fastq-files

This is typically not recommended as it runs the risk of annoying people in both communities.

Answer 1

[Brian Bushnell]

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.

[sumithrasank75]

I need the solution to be within the python environment