Sanity check · Bulk RNA-seq
How to Check Library Strandedness in Bulk RNA-seq
One minute with infer_experiment.py or salmon -l A tells you which strand flag to use, before a wrong guess quietly eats your counts.
By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed October 2026 · 5 min read
You got BAM files back from the core, picked a counter, copied a command from an old pipeline, and got a count matrix. Nothing errored. The strand flag (featureCounts -s, htseq-count --stranded) has a default that may not match your library, and a wrong value does not fail. It just counts the wrong reads.
The cost is large. In a published demonstration of the htseq-count --stranded parameter, using the wrong setting on a dUTP library gave about 700k counted reads instead of 21 million. Every downstream step, PCA, DESeq2, pathway analysis, then runs on a matrix that looks fine and is mostly empty or antisense.
This page gives you a check you can finish in the next hour: infer the strandedness from the data, confirm it with a second tool, set the flag, and verify the counting summary matches what you expected.
What it looks like when it's happening
- featureCounts .summary shows a very large Unassigned_NoFeatures or Unassigned_Ambiguity row, and Assigned is far below what the mapping rate suggests
- htseq-count output ends with __no_feature or __ambiguous totals that rival or exceed the sum of all gene counts
- Total counted reads per sample are a small fraction of uniquely mapped reads from the STAR or HISAT2 log
- Highly expressed genes you know should be there (housekeeping genes, a gene you knocked out or overexpressed) have strangely low counts
- Antisense or overlapping genes on the opposite strand show unexpectedly high counts next to a well-known sense gene
- Salmon logs a library type different from the one you passed with -l, or reports a low percentage of fragments compatible with the library type
- Switching -s 1 to -s 2 changes the assigned read count dramatically instead of slightly
Why it happens
Standard poly-A RNA-seq throws away which DNA strand the RNA came from. Stranded protocols keep it. The dUTP method, used in Illumina TruSeq Stranded Total RNA, marks the second cDNA strand with dUTP and degrades it before PCR, so only first-strand cDNA is amplified. Other protocols mark the template strand with adapters instead. The two families produce opposite read orientations relative to the gene, even though both are called "stranded".
A counter then has to decide whether a read's strand must match the gene's strand. featureCounts -s 0 ignores strand, -s 1 requires the read to be on the same strand as the feature, and -s 2 requires the opposite strand. HTSeq maps the same ideas to --stranded=no, yes and reverse. For a first-strand dUTP library you need -s 2 or --stranded=reverse. HTSeq's default is yes, which is the wrong answer for dUTP data.
When the setting is wrong, most reads land on the strand the counter refuses to accept. They are dropped as unassigned or called ambiguous, and the few that survive come from antisense transcription or overlapping genes on the other strand. Using -s 0 on a stranded library does less damage: you keep the reads, but reads from overlapping genes on opposite strands now collide and get counted for the wrong gene or flagged as ambiguous.
The failure is silent because nothing in the BAM says what protocol made it. The strandedness lives in the lab's library prep kit, which the core may not report, and a sample sheet from a public dataset often omits it. You have to infer it from the reads.
The checks
Run them in order. Each one tells you what healthy looks like and what the problem looks like.
0/8 checked · saved in this browser
Look up the kit in the core's submission form or the GEO/SRA metadata. Illumina TruSeq Stranded Total RNA is a first-strand dUTP protocol, so it needs
featureCounts -s 2orhtseq-count --stranded=reverse. Treat this as a hypothesis, not an answer, because kit names get copied wrongly between projects.- Healthy
- The kit name is recorded and maps to a specific strand setting that you write down before you run anything.
- Red flag
- The kit is unknown, listed as 'RNA-seq' only, or the sample sheet mixes kits across samples. Go to the data-based checks and do not trust any default.
Run it against a gene model in BED format. It samples reads, checks which gene strand each read overlaps, and reports the fraction of reads in each orientation pattern. Run it per sample, not once per project, so a mix-up at the core shows up.
bashinfer_experiment.py -r annotation.bed -i sample1.bam > sample1.strand.txt cat sample1.strand.txt- Healthy
- One pattern dominates by a wide margin and the other is near zero. For paired-end data, the two lines labelled with the '++,--' and '+-,-+' patterns tell you the answer: a large fraction for one and a small fraction for the other means stranded, with the dominant one deciding between `-s 1` and `-s 2`.
- Red flag
- The two fractions are both close to half, which means unstranded (`-s 0`), or the file also has a large 'failed to determine' fraction. The packet does not give a formal cutoff for 'clear', so judge by the gap between the two fractions and compare across samples.
Run Salmon with
-l Aon one or two samples and read the library type it reports in the log. Salmon's codes have three parts: orientation (I, O, M), strandedness (S or U) and direction (F or R), for example ISF, ISR, IU, OSR. Automatic detection cannot work if the upstream aligner was told to report only strand-aware alignments, so do not run it on a BAM that was already filtered that way.bashsalmon quant -l A -t transcripts.fasta -a sample1.bam -o salmon_sample1- Healthy
- The detected type agrees with infer_experiment.py. A dUTP library appears as an `R` direction, such as ISR, for paired-end data, and an unstranded library appears as `IU`.
- Red flag
- Salmon and RSeQC disagree, or Salmon picks IU for a library you believe is stranded. Stop and find out which tool saw the real data, usually by checking that the BAM is the one you think it is.
After a run with an explicit library type (for example
-l ISR), look at the percentage of fragments compatible with it in the Salmon output. With the correct setting, stranded protocols show about 92 to 100% of fragments matching the specified type.- Healthy
- A very high compatible fraction for the type you set, similar across all samples in the experiment.
- Red flag
- A noticeably lower percentage points to a misspecified type or a failed stranded protocol for that sample. One low sample among good ones is often a different kit or a swap.
Collect the infer_experiment.py outputs from every BAM into one table and sort. Samples from a single experiment should share a protocol, but samples pooled from public sources or from multiple batches may not.
bashfor f in *.strand.txt; do echo "$f"; grep -E "Fraction|\+\+|\+-" "$f"; done- Healthy
- Every sample reports the same dominant pattern.
- Red flag
- One sample or one batch flips pattern or looks unstranded. Count it with its own flag, or check that it belongs in the experiment at all. The packet does not cover partially stranded libraries, so a gradual mix of fractions needs investigation, not a guess.
Pick one BAM and run featureCounts three times with
-s 0,-s 1and-s 2, keeping every other option the same. Then compare the Assigned row in each.summaryfile. Use-pfor paired-end data.bashfor s in 0 1 2; do featureCounts -p -t exon -g gene_id -a annotation.gtf -s $s -o test_s$s.txt sample1.bam done grep -H Assigned test_s*.txt.summary- Healthy
- Two of the three runs clearly lose reads compared with the correct one, and the winner matches what infer_experiment.py said. For an unstranded library, `-s 0` wins and the other two lose roughly half.
- Red flag
- All three runs assign similar numbers, which suggests the library is unstranded or the annotation is wrong, or `-s 1` and `-s 2` are both low, which suggests the wrong GTF or a chromosome naming mismatch rather than a strand problem.
For each sample, divide the featureCounts Assigned count (or the sum of gene counts from HTSeq) by the uniquely mapped read count from the STAR or HISAT2 log. Do this for every sample and plot it as a bar chart or look at the sorted table.
- Healthy
- A consistently high fraction across samples, with only small variation from sample to sample.
- Red flag
- A fraction far lower than expected, or one sample far from the rest. The packet's example of about 700k assigned against 21 million is the extreme case, but even a halving is enough to change which genes pass filtering.
Pick a highly expressed gene, or a gene you edited, that has a neighbor on the opposite strand. Load the BAM in IGV with strand-colored reads and see which strand the reads cover. Then confirm the count matrix gives the sense gene the counts.
- Healthy
- Reads lie on the strand predicted by your protocol, and the count matrix puts them on the sense gene.
- Red flag
- Counts sit on the antisense neighbor, or a well-known gene has almost nothing in the matrix while IGV shows clear coverage.
What to do about it
Recount with the strand flag the data supports
When: infer_experiment.py and Salmon agree on a stranded library and your earlier count used a different setting.
Rerun the counter on the original BAMs with the right flag and discard the old matrix. For first-strand dUTP libraries use featureCounts -p -t exon -g gene_id -a annotation.gtf -s 2 -o counts.txt input.bam, or htseq-count --stranded=reverse -f bam annotation.gtf input.bam > counts.txt. For the other stranded family use -s 1 or --stranded=yes. For unstranded data use -s 0 or --stranded=no.
Caveat: Anything built on the old matrix (normalization, filtering, DE results, figures) must be rerun. Do not patch around it.
Let Salmon detect the library type at quantification
When: You are using Salmon for quantification and do not know the protocol, or you are processing a mix of public datasets.
Pass -l A so Salmon infers orientation and strandedness from the reads while it quantifies, then record the type it chose in your pipeline log for each sample. Afterward, set the type explicitly in the final production run so the choice is reproducible and written down.
Caveat: Automatic detection does not work if the aligner was told to report only strand-aware alignments. The packet also does not say whether it is equally reliable for single-end and paired-end data, so confirm single-end results with infer_experiment.py.
Set strandedness per sample, not per project
When: The checks show samples from different kits or batches within one experiment.
Keep a sample sheet column for the strand setting, fill it from the infer_experiment.py output, and have your pipeline read the flag from that column when it builds each counting command.
Caveat: If strandedness is confounded with condition (all treated samples from one kit), that is a design problem on top of a counting problem. Counting each sample correctly does not remove the batch difference.
Fix the pipeline default
When: You inherited a script or workflow that hardcodes a strand flag.
Remove the hardcoded value and make the strand setting a required parameter that must be passed explicitly. Add a step that runs infer_experiment.py and fails the run if its result disagrees with the declared setting.
Caveat: This takes an hour to build and saves every later project from the same silent error, but it adds a failure mode for datasets with unusual protocols that you will need to handle by hand.
When not to "fix" it
Do not change the flag just because a few genes look off. If infer_experiment.py shows a near 50/50 split, the library is unstranded and -s 0 or --stranded=no is correct, even if the protocol was sold as stranded; count it as unstranded and note the discrepancy to the core. Also leave antisense signal alone when it is real biology: genuine antisense transcripts and overlapping genes produce a small fraction of opposite-strand reads even in a perfectly stranded library, so a small minority pattern is expected and is not a reason to switch to -s 0. And if a project was finished and published with the correct setting, there is nothing to recount.
Five things experienced analysts do here
- Run infer_experiment.py on every BAM before the first count, and save the output next to the BAM. It takes about a minute and it is the cheapest check in the whole pipeline.
- Write the strand setting and the evidence for it (tool, output, date) into the project README, so a reviewer or a future you does not have to rederive it.
- Never trust the default of a counting tool. HTSeq's default `--stranded=yes` is wrong for dUTP libraries, and many featureCounts examples online use `-s 0` without saying so.
- Compare Assigned reads to uniquely mapped reads for every sample as a routine QC table. A strand error shows up there long before it shows up in a PCA.
- Treat disagreement between two tools as information. If RSeQC and Salmon disagree, look at the inputs (which BAM, which annotation, whether alignments were filtered) before choosing a winner.
Questions people ask
- How do I check if RNA-seq data is stranded?
Run RSeQC
infer_experiment.py -r annotation.bed -i input.bamon the aligned BAM and read the fractions it reports for each orientation pattern. A strong imbalance means stranded, and a near even split means unstranded. You can confirm withsalmon quant -l A, which detects the library type automatically.- What does featureCounts -s 0, 1 and 2 mean?
-s 0ignores strand,-s 1counts reads on the same strand as the gene, and-s 2counts reads on the opposite strand. First-strand dUTP libraries such as Illumina TruSeq Stranded Total RNA need-s 2.- What is reverse stranded RNA-seq?
It is a library where read 1 comes from the strand opposite the transcript, which is what dUTP first-strand protocols produce. In HTSeq this is
--stranded=reverse, in featureCounts it is-s 2, and in Salmon it appears as a library type with an R direction, such as ISR.- What happens if I use the wrong strandedness setting?
You do not get an error. In a published htseq-count example, the wrong setting on a dUTP library gave about 700k counted reads instead of 21 million. The reads are dropped as unassigned or ambiguous, and the surviving counts include antisense signal.
- Can Salmon figure out strandedness by itself?
Yes,
-l A(or--libType auto) infers it during quantification. It cannot do this if the aligner was told to report only strand-aware alignments, so verify with infer_experiment.py when you are working from a BAM.
Related pages
Related reading on the blog
Sources
- Strandedness during cDNA synthesis, the stranded parameter in htseq-count and analysis of RNA-Seq data — dUTP vs second-strand protocols, htseq-count --stranded options, and the 700k vs 21 million read example
- Salmon library type documentation — Three-part library type codes, -l A auto-detection, and the 92 to 100% matching fragments guide
- featureCounts manual - strand settings — featureCounts -s usage for strand-specific counting
- HTSeq 2.0 documentation: htseq-count — Official --stranded parameter options
- Batch Effect: To Correct or Not for Bulk RNA-seq Data — Related pre-DE check: look at PCA before differential expression
Part of the Strandedness series.