Install Subread (featureCounts) and HTSeq for RNA-seq Read Counting (conda)

Once STAR or HISAT2 has aligned your RNA-seq reads and handed you a BAM file, you still do not have an expression matrix — you have millions of individual read coordinates. The step that turns those alignments into a gene-by-sample counts table is called read summarization, and the two tools nearly everyone reaches for are featureCounts (from the Subread package) and htseq-count (from HTSeq). This guide installs both from bioconda into one isolated conda environment, verifies each one works, and walks through counting a real BAM against a GTF two ways.

Why install two counters?

featureCounts and htseq-count solve the exact same problem — assign aligned reads to genes using a GTF annotation — but they come from different traditions and disagree on a few defaults that trip people up constantly. Installing both side by side is genuinely useful: you can cross-check a suspicious result, or match whichever tool a published protocol used.

Neither tool is "wrong" — they are both counting the same overlaps, just with different default assumptions about strand and multi-mapping. The strandedness default is the single most common source of a counts table that looks broken when it is actually just mismatched to your library prep.

Prerequisites

You need a working conda installation. If you have not set one up yet, follow Install Miniconda on Apple Silicon first, then come back here.

As with every bioconda install on this site, we install only from conda-forge and bioconda and pass --override-channels so conda never falls back to the Anaconda defaults channel, which is gated behind a Terms-of-Service prompt on recent conda releases.

TL;DR: copy-paste install

conda create -n bu-counts --override-channels -c conda-forge -c bioconda subread htseq
conda activate bu-counts
featureCounts -v
htseq-count --version

The rest of this guide walks through each step and what to expect.

Step by step

1. Create the environment

We call the environment bu-counts and ask for both packages explicitly. conda-forge is listed first so it takes channel priority; bioconda is where subread and htseq actually live.

conda create -n bu-counts --override-channels -c conda-forge -c bioconda subread htseq

Conda will solve the environment and print a plan. You should see something like this (exact build strings and versions vary depending on when you run it, so they are omitted below):

Channels:
 - conda-forge
 - bioconda
Platform: osx-arm64
Collecting package metadata (repodata.json): done
Solving environment: done

## Package Plan ##

  added / updated specs:
    - htseq
    - subread

The following NEW packages will be INSTALLED:

  htseq              bioconda/osx-arm64::htseq-...
  subread             bioconda/osx-arm64::subread-...
  numpy               conda-forge/osx-arm64::numpy-...
  pysam               bioconda/osx-arm64::pysam-...
  python              conda-forge/osx-arm64::python-...
  ... (htslib, zlib, openssl, and other shared libraries)

Two things worth noticing. First, subread itself has almost no dependencies — it is a compiled C toolkit, so featureCounts shows up alongside just a handful of system libraries. Second, htseq pulls in a real Python stack: numpy for array math and pysam for BAM parsing, since htseq-count is a Python program built on pysam/htslib under the hood.

2. Activate the environment

conda activate bu-counts

Both featureCounts and htseq-count are now on your PATH, and only while this environment is active.

Apple Silicon note

Bioconda has published native osx-arm64 builds for subread and htseq for some time now, so on an Apple Silicon Mac the command above should resolve directly to arm64 packages with no emulation involved. You can confirm this after activating:

conda list -n bu-counts subread htseq
file $(which featureCounts)

The channel line for each package should read bioconda/osx-arm64, and file should report an arm64 Mach-O binary. If you ever hit a PackagesNotFoundError for either package on Apple Silicon — for example if you pin an unusually old version that predates arm64 support — fall back to an x86_64 environment run under Rosetta 2:

CONDA_SUBDIR=osx-64 conda create -n bu-counts_x86 --override-channels -c conda-forge -c bioconda subread htseq
conda activate bu-counts_x86
conda config --env --set subdir osx-64

For current subread and htseq releases you should not need this — the native install above is the right default.

Verify the install

Check that both binaries resolve inside the environment and report a version:

which featureCounts htseq-count
featureCounts -v
htseq-count --version

featureCounts -v prints its version banner to stderr rather than stdout — a quirk shared with several other Subread-suite tools — so do not be surprised if you need to redirect it (featureCounts -v 2>&1) when scripting.

Smoke test: count a tiny BAM against a GTF

The most convincing check is to actually count something. Build a one-contig reference and a single read, align it, write a one-gene GTF, then run both counters against the same BAM.

cd $(mktemp -d)
printf '>chr1\nACGTACGTACGTAGCTAGCTAGCTAGGCTAGCTAGCATCGATCGATCGTACGATCGATCGA\n' > ref.fa
printf '@r1\nAGCTAGCTAGCATCGATCG\n+\nIIIIIIIIIIIIIIIIIII\n' > reads.fq

bwa index ref.fa
bwa mem ref.fa reads.fq | samtools sort -o aligned.bam -
samtools index aligned.bam

cat > genes.gtf << 'EOF'
chr1	toy	exon	1	62	.	+	.	gene_id "geneA"; transcript_id "geneA.1";
EOF

featureCounts -a genes.gtf -o fc_counts.txt aligned.bam
htseq-count --format=bam --order=pos aligned.bam genes.gtf > htseq_counts.txt

Success looks like a geneA row with a count of 1 in both fc_counts.txt and htseq_counts.txt. If you instead get 0, the most common culprit on a toy example like this is the GTF's exon coordinates not overlapping the read at all, or a strand mismatch — see the errors table below.

Reading the output

featureCounts writes a tab-delimited file with a # comment line, then a header row of Geneid, Chr, Start, End, Strand, Length, followed by one column per BAM. It also writes a companion .summary file that tabulates how many reads were assigned versus unassigned (no feature, ambiguous, multi-mapping, and so on).

htseq-count instead writes plain gene_id-tab-count lines to stdout, with a handful of special rows appended at the bottom: __no_feature, __ambiguous, __too_low_aQual, __not_aligned, and __alignment_not_unique. Both files carry the same information; you just have to know which rows to strip before feeding either one into R.

From BAM to a real counts table

A realistic multi-sample featureCounts run for paired-end, reverse-stranded RNA-seq (a common setup for kits like Illumina's TruSeq Stranded) looks like this:

featureCounts \
  -a annotation.gtf \
  -o gene_counts.txt \
  -T 4 \
  -p --countReadPairs \
  -s 2 \
  sample1.bam sample2.bam sample3.bam

Here -a points at the GTF, -T 4 uses four threads, -p --countReadPairs counts paired-end fragments rather than individual mates (both flags are required together in current Subread releases), and -s 2 tells featureCounts the library is reverse-stranded — use -s 0 for unstranded and -s 1 for forward-stranded.

The htseq-count equivalent, using the same strandedness convention, counts all three samples in one call and runs them in parallel:

htseq-count \
  --format=bam \
  --order=pos \
  --stranded=reverse \
  --idattr=gene_id \
  --additional-attr=gene_name \
  -n 4 \
  sample1.bam sample2.bam sample3.bam \
  annotation.gtf \
  > gene_counts_htseq.txt

--format=bam matters: htseq-count defaults to expecting SAM text on that positional argument, so pointing it at a BAM without this flag is a very common first error. --order=pos matches BAMs sorted by coordinate, which is what STAR and HISAT2 plus samtools sort normally give you. -n 4 parallelizes across the three input files, and --additional-attr=gene_name tacks on a human-readable gene symbol column alongside the gene_id.

From counts to DESeq2

Both output formats plug directly into differential expression analysis, but each needs a small cleanup step first, since the raw file is not just a counts matrix.

For featureCounts, skip the comment line and drop the annotation columns (Chr, Start, End, Strand, Length), keeping only Geneid and the sample columns:

library(DESeq2)

fc <- read.delim("gene_counts.txt", comment.char = "#")
counts_mat <- as.matrix(fc[, 7:ncol(fc)])
rownames(counts_mat) <- fc$Geneid
colnames(counts_mat) <- c("sample1", "sample2", "sample3")

dds <- DESeqDataSetFromMatrix(
  countData = counts_mat,
  colData = sample_info,
  design = ~ condition
)
dds <- DESeq(dds)

For htseq-count, drop the trailing __no_feature/__ambiguous/etc. rows before building the same DESeqDataSet. Either way, this is exactly the counts-matrix step described in RNA-Seq Analysis: From Raw Reads to Differential Expression: STAR or HISAT2 produces the aligned BAM, featureCounts or htseq-count collapses it to a gene-by-sample matrix, and DESeq2 takes it from there.

Common errors and fixes

Managing the environment

Update both counters in place:

conda update -n bu-counts --override-channels -c conda-forge -c bioconda subread htseq

Export a snapshot so a collaborator (or a future you) can reproduce the exact environment:

conda env export -n bu-counts > bu-counts-env.yml

Remove it cleanly when you no longer need it:

conda remove -n bu-counts --all

Next steps

With a gene counts table in hand from either tool, the natural next stop is differential expression itself — see RNA-Seq Analysis: From Raw Reads to Differential Expression for the full DESeq2 workflow. If your BAMs are not sorted and indexed yet, Install SAMtools and BCFtools covers that step, and if you are still choosing an aligner, Install BWA and Bowtie2 is a good place to start.