How are multimapping reads handled right now by th...
# rnaseq
m
How are multimapping reads handled right now by the pipeline by default? Are they used for count matrix generation? One of my conditions has far more multimapping reads than the other. Can I just use the output from the default settings and take it into differential abundance?
j
By default (
--aligner star_salmon
), STAR aligns to genome + transcriptome, and Salmon quantifies from the transcriptome BAM in alignment mode. Multimapping reads aren't simply discarded: if a read maps to 2-20 loci, all alignments are kept and Salmon's EM algorithm probabilistically assigns it across those transcripts based on their relative abundance, then that rolls up to gene counts via tximport. So most multimappers do get used. But there's a hard ceiling: we set
--outFilterMultimapNmax 20
, so a read mapping to more than 20 loci is dropped entirely, STAR marks it unmapped, and it never reaches Salmon. That shows up as "% mapped to too many loci" in STAR's
Log.final.out
, separate from "% mapped to multiple loci" (the 2-20 loci reads that do get used). So before taking the default output into DA: check that STAR log split for both conditions. • Mostly in "multiple loci" -> EM handles it, counts should be fine. • Meaningful chunk in "too many loci" -> those reads are just gone, worth chasing down first, before you do a lot of analysis. That pattern usually means rRNA/repeat contamination or gDNA rather than ordinary paralog ambiguity. If you do want to raise the ceiling rather than just diagnose, the threshold is overridable via
--extra_star_align_args
, e.g.
--extra_star_align_args '--outFilterMultimapNmax 50'
. The pipeline merges this with its own STAR args and lets your value win on any flag it repeats, so you don't need to touch the pipeline config. Worth testing as a rerun on the affected samples to see how many "too many loci" reads get rescued into the multimapping bucket, rather than applying it blind. There's some info here about which files to use with DA: nf-co.re/rnaseq/3.26.0/docs/output#…
m
I should have clarified that I am using the prokaryotic pipeline so am not using STAR, just Bowtie2/Salmon. I did look in those Bowtie2 logs but it only provides a number of reads aligned concordantly >1 times. Is there any guidance specific to the prokaryotic pipeline?
j
I'm not a bacterial genomics expert I'm afraid, so take the following with a pinch of salt and seek better advice. But I don't think Bowtie2 does any exclusion of multi-mappers, so the story would be similar. You'll want to understand what the reason for that multi-mapping is. It's possible it's to do with rRNAs, try
--remove_ribo_rna --ribo_removal_tool bowtie2
and see if that does anything. If that doesn't change matters much, then you're on to reasons derived from annotation and I'm hesitant to advise. We derive the transcriptome for the prokaryotic pathway based on the CDS from the GTF, there's no clustering or anything going on. It may be that you have a multi-copy family of some kind that's more important in one of your conditions than the other. Try the ribo-removal first 🙂