Hello,
In my current job, I am dealing with perturbation rna-seq data from cancer biopsy. The sequencing library is prepared using heat lysis and polyA tail enrichment to only select mRNA data.
Here is my problem, after trimming and assessing sequencing quality (which is good), I notice a very low mapping ratio on my data (between 10 to 30%). I use kallisto for mapping with kmer=31. I check for my reads size distribution, and they are most of the time over the required 31 n lenght.
I checked for contamination using fastq_screen, and most of the reads (>90%) mapped to human genome on a subset of 2M reads.
Do you have any idea where these non-mapping reads comes from? Is it possible that my reads are mostly from intronic region (even with polyA purification, explaining why we see high mapping ratio with bowtie2 from fastq_screen)?
•
written
by
meetmet •
1 answer
Hey,
Your polyA enrichment should primarily capture mature mRNA, but it can also include some pre-mRNA transcripts that retain introns, especially in cancer samples where splicing may be dysregulated. This could explain why FastQ-Screen with Bowtie2 shows high mapping to the human genome (>90%), as Bowtie2 aligns to the full genome and tolerates intronic sequences, while Kallisto pseudoaligns specifically to the transcriptome and may discard reads with intronic content.
The low Kallisto mapping (10-30%) and moderate STAR alignment (49%) suggest many reads are not fully exonic. Your BLAST results on unmapped reads hitting BAC clones (human genomic fragments) further indicate these may derive from intronic or intergenic regions, possibly due to genomic DNA contamination during library preparation via heat lysis, which may not fully eliminate DNA despite RNA-focused protocols.
To investigate:
- Analyze your STAR BAM files for read distribution across exons, introns, and intergenic regions using RSeQC:
infer_experiment.py -i your_sample.bam -r hg38.bed
This will quantify strandedness and intron/exon overlap.
- Or use featureCounts from Subread for similar categorization:
featureCounts -a hg38.gtf -o counts.txt -t exon -g gene_id your_sample.bam
Adjust for intron counting by modifying the GTF.
- Check for ribosomal RNA contamination by mapping reads to rRNA sequences with Bowtie2:
bowtie2 -x rRNA_index -1 reads_R1.fastq -2 reads_R2.fastq -S rRNA_mapped.sam
High rRNA hits could explain non-mapping to transcriptome.
- Retry pseudoalignment with Salmon using a decoy-aware index (include full genome as decoys) to handle potential DNA contamination:
Follow Salmon documentation to build the index with --decoy option.
If decoy inclusion boosts mapping, it confirms genomic DNA presence. Sequence deeper if on-target counts remain low, but address library issues upstream.
Kevin
Log in to answer this question.
-
Nanopore DRS data
written by baibhu1234 •
-
ONT Improved RCA Barcoded Plasmid Sequencing
written by thewiggler •
-
Estimate large insertions size in WGS data.
written by meetmet •
-
Low exonic mapping of RNASeq data to draft genome. Should I suspect poor annotation?
written by CephBirk •
-
Matching read IDs to mapped transcripts/genes
written by joshcylee •
-
Mapping both PE reads and Single-end reads into transcriptome
written by pbigbig
-
error while creating kover dataset from fastq
written by ekaterina.avershina •
-
Illumina Rnaseq: Adaptor Sequence Followed By Polya?
written by Biomonika (Noolean)
-
Simulation Of Rnaseq Read Mapping With Reads Having Different Polya Tail Lengths.
written by Gaima Kushwaha •
-
Why Are There Many Rna-Seq Hits To Intronic Regions?
written by lh3
Could be genomic DNA contamination. Enrichment is just that, enrichment, not perfect selection without noise. Poor RNA quality usually increases noise. In the end it does not really matter, since in silico magic cannot save library issues. If on-target counts are too low you have to sequence deeper.
Kallisto is a pseudoaligner, so, as ATpoint suggested, it could be DNA contamination. You can confirm this using an aligner and a reference genome.
Looks like OP has done that, assuming that the samples are human.
I did contamination analysis, and I found nothing. I also run another aligner, STAR, and I could align "only" 49% of my reads with default parameters. So I am quiet lost. When I blast the unmapped reads, they ofter map to BAC clone of different chromosomes for human genoem (my interest) with very high confidence....
It’s possibly due to ribosomal RNA contaminating — those are regions that are difficult to assemble.
I recommend using bowtie2 or kallisto to map solely against ribosomal RNA and see the percentage of reads that map.
It's worth noting that the k-mer size is different from your insert / read size. The value of k _should_ be smaller than your read size, since the k-mer spectrum is used to identify the gene/transcript a read came from.
If you have a lot of reads mapping to the human _genome_ but not the _transcriptome_, this suggests that you have some kind of DNA contamination, which might come from an issue in the polyA enrichment step, or even reverse transcription.
What reference transcriptome did you use for kallisto? Does it contain decoy sequences? I'm wondering what would happen if you add a set of decoy sequences like the salmon docs suggest and try again. If you do this decoy-aware quantification with the human genome as the set of decoys and you find a much higher percentage of reads successfully mapping, then you likely have some human DNA that has made its way through the RNA protocol.
But your follow-up comment about BLASTing the unmapped reads suggests you might have bacterial contamination in your samples, not just DNA contamination in your RNA samples. It might be worth talking with your local wet lab biologist or lab tech to get some advice on this.
Those aren’t bacteria; they are BAC clones — often used for storing fragments of larger genomes.
There’s no reason to re-run pseudoalignment; the STAR BAM results can be used to infer how many reads originate from exon, intron, or intergenic regions (RSeqc or subread should be able to do this).