2
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ But the number of reads in the FastQ sample file that was fed into STAR was 59,194,421 (I stated it was much fewer than was reported). I believe that the 185M and 206M (391M total) referred to by FastQC as "reads" are actually "mappings" since I fed FastQC a bam file instead of a FastQ file, which is what this question is about. Perhaps this is a misunderstanding stemming from the recent edit. I will fix it. $\endgroup$ – hepcat72 Jul 14 '20 at 17:03
  • $\begingroup$ I understand the general caveats about running FastQC on RNA-Seq FastQ data. I regularly check for rRNA contamination and check over-represented sequences for rRNA segments. My question however pertains to why you can run FastQC on a bam file, what its utility is, how to interpret the results in that context (e.g. a "read" is a "mapping"), and what the pitfalls are WRT using a bam as input. Also, I knew it was a bad sequencing run, but this question is not about identifying bad sequencing runs. $\endgroup$ – hepcat72 Jul 14 '20 at 17:14
  • $\begingroup$ Hi I think @swbarnes has answered the question very nicely, they are concerned with two different measurements and swb proivide the underlying rationale ... one was for DNA genomes not for RNAseq. It would be cool to accept the answer. $\endgroup$ – M__ Jul 14 '20 at 17:31
  • $\begingroup$ If you think that FASTQC should have a built-in filter for primary alignments only when reading a bam, tell the FASTQC developers, not me. $\endgroup$ – swbarnes2 Jul 14 '20 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.