Comparison · read counting
featureCounts vs HTSeq-count: Which One Should You Use?
Same BAM, same GTF, but a tenfold-plus speed gap and different defaults for overlapping and multi-mapped reads mean the two tools rarely produce identical count matrices.
By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 4 min read
The verdict
Default to featureCounts for anything bulk RNA-seq, ChIP-seq, or ATAC-seq scale. On a benchmark with roughly 4.4M RNA-seq reads, featureCounts finished in 1.0 minute using 16 MB of memory against HTSeq-count's 22.7 minutes and 101 MB; on roughly 5.4M ChIP-seq fragments the gap was 0.9 minutes versus 36.0 minutes. That's not a rounding difference, it's the gap between counting a cohort of a hundred BAMs over lunch versus queuing an overnight cluster job. featureCounts also multithreads a single job natively with -T, ships as an R package (Rsubread) that drops straight into a DESeq2/edgeR workflow, and accepts a plain SAF table when you want to count over custom regions like peaks without building a GTF.
Reach for HTSeq-count when you need its named overlap modes. intersection-strict or intersection-nonempty give you an explicit, documented answer for reads that straddle exon boundaries or hit overlapping gene annotations, and --nonunique fraction/random give you a stated, chosen policy for multi-mapped reads instead of featureCounts' fragment-voting heuristic. It's also the right call when you're reproducing a published analysis that names htseq-count and a version, since swapping counters changes the exact count matrix even under "equivalent" settings, for the reasons in the table below. For a handful of samples, the speed penalty simply doesn't matter enough to justify breaking that reproducibility.
Both tools do the same job on paper: take an aligned BAM and a gene model, and produce one count per feature per sample. The difference is what happens when a read touches more than one feature, or aligns more than once, and how fast the arithmetic runs, because both of those choices are baked into the algorithm, not just the command-line interface.
featureCounts is a single C program from the Subread package. It streams through the BAM once, assigns each read or read pair (with -p) to the feature it overlaps, and by default silently drops a read that overlaps more than one feature unless you turn that on with -O. It multithreads a single counting job across cores with -T, and it reads either a full GTF/GFF or a bare-bones SAF table (GeneID, Chr, Start, End, Strand) you can build by hand for custom regions.
HTSeq-count is a Python program with its own GFF parser and its own vocabulary for ambiguity. -m switches between union (the default), intersection-strict, and intersection-nonempty for how it resolves a read spanning an intron or two adjacent exons, and --nonunique switches between none (default, drop), all, fraction, and random for a read that lands in more than one feature. That flexibility costs wall-clock time, because it's pure Python doing per-read interval arithmetic instead of C.
Head to head
| Criterion | featureCounts | HTSeq-count | Edge |
|---|---|---|---|
| Speed on RNA-seq data (~4.4M reads, SEQC dataset) | Completes in about 1.0 minute using roughly 16 MB of memory. | Takes about 22.7 minutes using roughly 101 MB of memory. | featureCounts |
| Speed on ChIP-seq data (~5.4M fragments, H3K27me3) | Finishes in about 0.9 minutes using roughly 4 MB of memory. | Takes about 36.0 minutes using roughly 31 MB of memory. | featureCounts |
| Multithreading | Native multithreading via `-T` splits one BAM's counting across cores; it's the only major counter with this. | Single-threaded per invocation; you parallelize by launching one process per sample yourself. | featureCounts |
| Implementation and where it runs | Written in C, shipped as the Subread command-line tool and the Rsubread R/Bioconductor package, runs on Windows, macOS, and Linux. | Written in Python, distributed as the HTSeq package, primarily a Unix command-line tool. | featureCounts |
| Overlap-resolution control | One core rule: by default a read overlapping more than one feature is dropped; `-O` counts it for every feature it overlaps, with optional fractional assignment via `fraction=TRUE`. | Three named modes via `-m`: union (default), intersection-strict, and intersection-nonempty, each with a distinct, documented rule for ambiguous overlaps. | HTSeq-count |
| Ambiguous / multi-mapped read handling | Excludes multi-overlap reads by default; fragment voting can assign a fragment to the gene whose exons more of its mate reads overlap. | `--nonunique` defaults to `none` (drop ambiguous reads), with `all`, `fraction`, or `random` as explicit alternatives you pick per run. | Tie Both let you keep or split ambiguous reads; featureCounts votes by evidence, HTSeq-count lets you choose a global policy. |
| Feature boundary (interval) convention | Treats the rightmost base of a feature as included (closed interval), matching the GFF specification. | Treats the rightmost base of a feature as excluded (open interval), which can shift counts by a read at feature edges. | featureCounts |
| Annotation input formats | Reads GTF/GFF and also SAF, a plain 5-column format (GeneID, Chr, Start, End, Strand) you can hand-build for custom regions without touching GTF syntax. | Reads GTF/GFF only, through HTSeq's own Python GFF parser. | featureCounts |
| Documentation depth | User's Guide lives inside the Subread/Rsubread documentation plus the original 2014 Bioinformatics paper. | A dedicated, example-heavy readthedocs page walks through each overlap mode and each `--nonunique` setting with worked examples. | HTSeq-count |
Use featureCounts when
- You're counting dozens to hundreds of BAMs for a bulk RNA-seq or ChIP-seq cohort and wall-clock time or memory on a shared cluster actually matters.
- Your downstream analysis is already in R/Bioconductor (DESeq2, edgeR, limma) and you want to call Rsubread's featureCounts() directly instead of shelling out.
- You need to count reads over custom regions like ChIP-seq peaks or enhancers and want a simple SAF table instead of a full GTF.
- You have multiple cores available and want a single counting job to use them via -T instead of running samples one at a time.
- You're on Windows or macOS without a Python/Unix environment already set up for HTSeq.
Use HTSeq-count when
- You're reproducing a published analysis whose methods section names htseq-count specifically, and matching its exact numbers matters more than speed.
- You need one of the named overlap modes, intersection-strict or intersection-nonempty, to resolve reads over heavily overlapping gene annotations in a documented, citable way.
- You're counting a small number of samples, where a 22.7-minute run instead of a 1-minute run has no practical cost.
- You want explicit, chosen control over ambiguous multi-mapped reads via --nonunique fraction or random rather than featureCounts' fragment-voting heuristic.
- Your existing validated pipeline already uses htseq-count and switching would perturb count-based thresholds you've already published or locked in.
Switching between them
Merging output differs in shape: featureCounts writes one table with Geneid/Chr/Start/End/Strand/Length plus one column per BAM, so you must drop the five meta columns before joining samples (the merge approach in lesson 85 does exactly this); HTSeq-count writes one two-column file per sample and appends summary rows (__no_feature, __ambiguous, __too_low_aQual, __not_aligned, __alignment_not_unique) that must be filtered out before you build a count matrix. Defaults aren't equivalent even when they sound like it: featureCounts drops reads overlapping more than one feature unless you add -O, and HTSeq-count's default union mode plus --nonunique none also drops ambiguous reads, but the two algorithms don't discard the same reads, so total counted library size shifts when you switch tools even with "default" settings on both sides. The interval convention differs too: HTSeq-count excludes a feature's rightmost base, featureCounts includes it, which matters most for short features or reads that just barely clear a boundary. Re-derive your strand setting on the new tool rather than assuming a flag value carries over, and re-check any count or CPM threshold you tuned on one tool's output before trusting it on the other's.
Pitfalls with either
- Running featureCounts with the default -s 0 on a stranded library inflates counts for antisense genes; check the library prep protocol and set -s 1 or -s 2 before trusting the numbers.
- Treating featureCounts' default exclude-multi-overlap behavior as equivalent to HTSeq-count's default union mode; read what each algorithm actually does before assuming 'default' means the same thing across tools.
- Feeding featureCounts' raw output straight into a merge script without dropping the Chr/Start/End/Strand/Length columns breaks a simple join across samples; strip those five columns first, as in lesson 85's merge approach.
- Loading HTSeq-count files into DESeq2/edgeR without removing the trailing __no_feature/__ambiguous/__too_low_aQual/__not_aligned/__alignment_not_unique rows counts them as fake genes; filter out any row whose ID starts with __ before building the matrix.
- Running htseq-count serially across dozens of samples wastes time since it won't multithread itself; launch one process per sample in parallel with GNU parallel or a job array instead of waiting on a serial loop.
- Assuming a pipeline swap from HTSeq-count to featureCounts (or back) reproduces the exact same count matrix on the same BAM and GTF; expect small shifts from the interval convention and overlap-resolution differences and re-validate any gene near a significance cutoff.
Questions people ask
- Is featureCounts more accurate than HTSeq-count?
Neither is more accurate in a statistical sense; both are deterministic algorithms that differ in which reads get counted at feature boundaries and overlaps. featureCounts includes a feature's rightmost base (closed interval, matching the GFF spec) while HTSeq-count excludes it (open interval), so counts at feature edges can differ by a read or two on identical input.
- Why is featureCounts so much faster than HTSeq-count?
featureCounts is a C program and HTSeq-count is Python; on roughly 4.4M RNA-seq reads featureCounts finished in about 1.0 minute versus HTSeq-count's 22.7 minutes, and featureCounts also multithreads natively with -T while HTSeq-count runs one BAM per process.
- Can I mix featureCounts and HTSeq-count outputs in the same project?
Avoid it for a single differential expression comparison. The two tools resolve overlapping and multi-mapped reads differently by default, so counts for the same BAM and GTF won't match exactly, and mixing them into one count matrix introduces a technical batch effect that has nothing to do with biology.
- Does the strand setting mean the same thing in both tools?
No, treat it as a separate decision for each tool. featureCounts uses a numeric -s flag (0 unstranded, 1 stranded, 2 reversely stranded), and getting it wrong on a stranded library silently changes counts for genes with antisense transcription, so check the library prep protocol before copying a flag value between tools.
- Should I switch an existing published pipeline from HTSeq-count to featureCounts?
Only if you're prepared to re-validate downstream results. Because of the interval convention and overlap-resolution differences, the two tools won't produce identical count matrices even with matched settings, so a mid-project switch can shift low-count genes across your significance threshold.
Related pages
- Glossary · Count matrix
Related reading on the blog
Sources
- featureCounts: an efficient general purpose program for assigning sequence reads to genomic features — Speed/memory benchmark, interval convention, fragment voting, multithreading.
- HTSeq documentation: htseq-count counting reads within features — Default union mode and --nonunique settings.
- featureCounts documentation and user guide — Default overlap exclusion, -O flag, SAF format.
- Rsubread Bioconductor package documentation — Rsubread/Subread version numbers.