Hi all,
I'm trying to screen paired Illumina RNAseq reads against a database of putative contaminants prior to assembly. BWA (version 7.0.3a, in case it matters) ran beautifully. Now I'm trying to extract the names of all the reads that mapped to something in my contaminant database. Basically, if a read OR its mate matched a putative contaminant, I want to toss it out. A colleague told me that the output is easy to parse, given that all the unmapped reads will have a "*" in the third column of the output line. I wrote up a quick perl script to parse the output based on the "*".
Later, I learned that its possible to convert BWA's sam output to bam (via samtools) and then count the mapped reads with the following command:
samtools view -c -F4 <bam file>
(I found that here: http://seqanswers.com/forums/showthread.php?t=16743)
The issue is that I get two different answers, and I can't figure out why. What's wrong with my approach?
I'm trying to screen paired Illumina RNAseq reads against a database of putative contaminants prior to assembly. BWA (version 7.0.3a, in case it matters) ran beautifully. Now I'm trying to extract the names of all the reads that mapped to something in my contaminant database. Basically, if a read OR its mate matched a putative contaminant, I want to toss it out. A colleague told me that the output is easy to parse, given that all the unmapped reads will have a "*" in the third column of the output line. I wrote up a quick perl script to parse the output based on the "*".
Later, I learned that its possible to convert BWA's sam output to bam (via samtools) and then count the mapped reads with the following command:
samtools view -c -F4 <bam file>
(I found that here: http://seqanswers.com/forums/showthread.php?t=16743)
The issue is that I get two different answers, and I can't figure out why. What's wrong with my approach?
Comment