If you can ignore the abysmal mapping results and large amount of PCR duplication, I am finding something funny going on with rmdup, single-end Torrent PGM data, and mapQ scores:
Let's say I have 554109 mapped reads.
After rmdup -s I am left with 16761 nodup mapped reads.
If I check mapQ >= 10, I have 122 reads.
If I check mapQ >= 20, I have 1 read.
Now let's say I repeat, only this time I remove all the unmapped reads and all the mapped reads of mapQ < 10.
I start with 554109 mapped reads. Of these, 171375 have mapQ >= 10.
I isolate them with samtools view –F 4 –q 10 –b –h .
After rmdup -s I am left with 15688 nodup mapped reads
If I check mapQ >= 10, I have 15688 reads.
If I check mapQ >= 20, I have 29 reads.
Now let's say I repeat again, only this time I am going to remove all reads with mapQ less than 20.
I start with 554109 mapped reads. Of these, 4138 have mapQ >= 20.
After isolation and rmdup -s, I am left with 3154 nodup mapped reads.
If I check mapQ >= 10, I have 3154 reads.
If I check mapQ >= 20, I have 3154 reads.
rmdup -s is performing differently based on the mapped quality scores of the reads that I give it, when I was under the impression that for single end reads it takes the 5' position and only retains the one with the highest mapping quality. It is obviously not doing that here, as it seems to be preferentially choosing lower mapQ reads over reads with higher mapQ scores. At least, that's the only explanation I can think of - I shouldn't be able to magically pull out 3153 more mapQ >20 from the same original 554109 mapped reads.
Thoughts?
Let's say I have 554109 mapped reads.
After rmdup -s I am left with 16761 nodup mapped reads.
If I check mapQ >= 10, I have 122 reads.
If I check mapQ >= 20, I have 1 read.
Now let's say I repeat, only this time I remove all the unmapped reads and all the mapped reads of mapQ < 10.
I start with 554109 mapped reads. Of these, 171375 have mapQ >= 10.
I isolate them with samtools view –F 4 –q 10 –b –h .
After rmdup -s I am left with 15688 nodup mapped reads
If I check mapQ >= 10, I have 15688 reads.
If I check mapQ >= 20, I have 29 reads.
Now let's say I repeat again, only this time I am going to remove all reads with mapQ less than 20.
I start with 554109 mapped reads. Of these, 4138 have mapQ >= 20.
After isolation and rmdup -s, I am left with 3154 nodup mapped reads.
If I check mapQ >= 10, I have 3154 reads.
If I check mapQ >= 20, I have 3154 reads.
rmdup -s is performing differently based on the mapped quality scores of the reads that I give it, when I was under the impression that for single end reads it takes the 5' position and only retains the one with the highest mapping quality. It is obviously not doing that here, as it seems to be preferentially choosing lower mapQ reads over reads with higher mapQ scores. At least, that's the only explanation I can think of - I shouldn't be able to magically pull out 3153 more mapQ >20 from the same original 554109 mapped reads.
Thoughts?