Lower alignment rate when using collapsed reads

Question

[blz]

Hello everyone!

I made a RNA-Seq mapping two ways: 1) using all reads, and 2) using collapsed reads.

I got a lower percentage of mapped reads in the second case and I can't understand why: collapsed reads represent unique reads, should they not map to same locations? I mean, the sequence is the same, only the quantity is varying (many per read type in the first case and only one per read type in the second case). Is there a simple explanation to this? (I used fastx_collapser to collapse reads)

bowtie2 all reads mapping

bowtie2 -p 6 -x genome_index -U seq_noAdapt.fastq.gz -S seq_noAdapt.sam 2>seq_noAdapt.log

27557790 reads; of these:
  27557790 (100.00%) were unpaired; of these:
    3488459 (12.66%) aligned 0 times
    4281314 (15.54%) aligned exactly 1 time
    19788017 (71.81%) aligned >1 times
87.34% overall alignment rate

bowtie2 reads collapsed mapping

bowtie2 -p 6 -x genome_index -f seq_noAdapt_collapsed.fasta -S seq_noAdapt_collapsed.sam 2>seq_noAdapt_collapsed.log

1166266 reads; of these:
  1166266 (100.00%) were unpaired; of these:
    504804 (43.28%) aligned 0 times
    192571 (16.51%) aligned exactly 1 time
    468891 (40.20%) aligned >1 times
56.72% overall alignment rate

Some notes about the dataset

• It's a old one, so it's single-end, 36bp reads.
• It's intend to identify short RNAs.
• As you can see above, there is a high number of multi-mapping reads and this is expected: the loci giving rise to them are repetitive (and the genome itself is very repetitive).

Thank you!

Answer 1

[LChart]

Sequences are collapsed if they're sufficiently identical. Garbage (random) reads won't collapse. If you collapse the "good" reads and don't collapse the "garbage", then the proportion of "garbage" must go up.

Answer 2

[Brian Bushnell]

It looks like fastx_collapser does not support mismatches. As per LChart's answer, you are enriching for low-quality reads by getting rid of exact duplicates. There are two ways of dealing with this. One is to use BBTool's Clumpify which does allow mismatches when removing duplicates, although with such short reads you'd have to use non-default parameters like shorter kmers (here, s is the number of mismatches allowed):

clumpify.sh in=reads.fq out=deduped.fq dedupe k=19 s=2 passes=6

Alternatively, you can deduplicate after mapping:

dedupebymapping.sh in=mapped.sam out=deduped.sam

This is more robust to errors than sequence-based deduplication, although it will still enrich for garbage because the unmapped reads don't get deduplicated. You can also do both!

However, in this case, I don't recommend any deduplication. The reason is that with 36bp SE reads in RNA-seq data you expect a lot of duplicates and there is no way to differentiate PCR duplicates from coincidental duplicates caused by high coverage. Duplicate removal in RNA-seq is questionable in general if you are doing anything quantitative (such as peak-calling), but especially so if you only have SE reads. You can still remove optical duplicates but not PCR duplicates.