What pitfalls exist with running FastQC on a bam file?

Question

I was analysing a bad sequencing run of some RNA data using FastQC, I supplied it RNA-STAR-aligned bam files (and the human reference). The output of MultiQC had many more "unique reads" in the sequence counts table than there were reads in the FastQ files that were fed to RNA-STAR.

At first, I realized that all the FastQC stats were about "mappings", not reads (as they were labeled), but "unique mappings" doesn't seem to make sense either, because even if I consider the possibility that segments of reads were mapping in different places, the numbers don't make sense. Here is the worst case:

• Max read length after trimming: 47 bases
• Number of reads in the original FastQ file that went into STAR: 59M

Bam Stat Table Numbers:

• Total records: 391 million
• Unique mappings: 3.1 million
• Non-primary Hits: 335 million

FastQC Sequence Counts Bar graph:

• "Unique reads": 185 million
• "Duplicate reads": 206 million

While I found documentation stating that you can run FastQC on bam files, there were no caveats I could find about it needing to be unmapped reads.

Can anyone explain precisely why there can be 185 million "unique reads" when bam stats says there were 3 million unique mappings and the number of reads is 59M?

Answer 1

I don't see why you are trying to compare two totally different measurements.

When FASTQC says a read is unique, it means that no other read shares its sequence.

When STAR says a read is uniquely mapped, that means it maps to one and only place. That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.

You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.

FASTQC was designed for DNA sequencing, not RNASeq. Some of its warnings make no sense in the context of RNASeq. If I was doing RNASeq of a filtered set of pancreatic alpha cells, I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.

You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes. The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.

That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong. Find out where all those reads are aligning, you might have rRNA contamination. If you are doing an organism with a high level of duplication in its genome that might be normal.