Reading HTSeq Raw Counts

HTSeq is a Python bioinformatics software package, and htseq-count is its command-line program for counting reads by gene. HTSeq raw counts are the unnormalized data table produced by htseq-count.

The table records how many sequence fragments, reads or read pairs, were assigned to each gene in each RNA-seq sample. It is created by aligning reads to a reference genome, then comparing their positions with an annotation containing gene coordinates.

In a count matrix that combines multiple samples, rows are usually genes and columns are samples. Suppose part of the sample_A column looks like this:

gene_id             sample_A
ENSG00000000003     9896

sample_A: a sample name assigned by the researcher
ENSG00000000003: an Ensembl ID that identifies a gene
9896: the number of reads, short sequence fragments produced by the sequencer, judged to have come from this gene

For paired-end data, the two reads are grouped and counted as one pair, so 9896 means that 9,896 read pairs were assigned to this gene. It does not mean that exactly 9,896 RNA molecules were present. The number is called a raw count because it has not yet been corrected for sequencing depth in each sample or for gene length.

The roles of HTSeq and GTF

An aligned read contains only information such as aligned at position 11870 on chr1. The alignment alone does not identify the gene at that position. It must be compared with a separate table of gene locations.

HTSeq is a Python-based bioinformatics software package for processing high-throughput sequencing data. Its htseq-count command counts aligned reads by gene.
GTF (Gene Transfer Format) is an annotation file that records the locations of genes, transcripts, and exons on a genome. Here, annotation does not mean a note that changes the sequence. It means information that assigns biological meaning to genomic coordinates.

For example, suppose an alignment records a read at chr1:1000-1075, and the GTF contains this line:

chr1  source  exon  900  1100  .  +  .  gene_id "GENE1";

The line says that positions 900 through 1100 on chr1 form an exon belonging to GENE1. Because the read position overlaps this interval, htseq-count increments the count for GENE1. In other words, its job is to compare two coordinate tables.

Stage	File or program	Information and role
Input	SAM/BAM/CRAM	Where each read aligned in the genome
Input	GTF	Where each gene and exon lies in the genome
Processing	`htseq-count`	Compare the coordinates in the two files and assign reads to genes
Output	raw counts	Number of reads or fragments assigned to each gene

In RNA-seq, HTSeq usually reads GTF lines marked exon. Multiple exons with the same gene_id are grouped as one gene, and all observations assigned to those exons are summed. GTF: gene annotation files continues with a detailed explanation of the nine GTF columns and gene_id.

The output is a gene × sample matrix

gene_id              sample_1  sample_2
ENSG00000000003           9896      6279
ENSG00000000005             32        49
ENSG00000000419           5129      2406

Combining ordinary bulk RNA-seq results creates a matrix where rows are genes, columns are samples, and values are counts. Input orientation differs across analysis tools, so do not guess the rows and columns. Check the file’s first row and the format required by the tool.

When one count is added

The default union mode checks the set of genes that a read overlaps.

Genes overlapped by the read	Handling
Exactly one	Increment that gene’s count by 1
None	Record in `__no_feature`
Two or more	Record in `__ambiguous` by default

For paired-end data, the two mates are counted not as two separate reads but as one read pair, meaning evidence for one cDNA fragment, so the count increases once. Some options, such as --nonunique=fraction, can produce fractional counts, but ordinary gene-level counts are integers.

Rows beginning with `__` are not genes

HTSeq also totals reads that were not assigned to a gene, grouped by reason. All these special counters begin with __, making them distinguishable from gene rows.

Special counter	Meaning	What to check when the value is high
`__no_feature`	Does not overlap any annotation feature	genome build, chromosome names, strandedness, GTF
`__ambiguous`	Overlaps two or more features, so one cannot be selected	overlapping genes and annotation, overlap mode
`__too_low_aQual`	Excluded because alignment quality is below the minimum	aligner MAPQ and HTSeq `-a` setting
`__not_aligned`	Unaligned record	upstream alignment results
`__alignment_not_unique`	Aligned to multiple genomic locations	repetitive sequences, aligner multi-mapping settings

You may see lines like these at the end of an HTSeq output file:

__no_feature
__ambiguous
__too_low_aQual
__not_aligned
__alignment_not_unique

These rows are not genes, so they are excluded from DESeq2 input. It is still useful to inspect their values by sample before discarding them. If one sample has an unusually high proportion of __no_feature or __alignment_not_unique, first inspect alignment, annotation, and differences in sample quality.

special = genes_by_samples.loc[
    genes_by_samples.index.str.startswith("__")
]
gene_counts = genes_by_samples.loc[
    ~genes_by_samples.index.str.startswith("__")
]

print(special)

A large special counter alone does not mean the sample failed. Expected values differ with the protocol, annotation coverage, and method for handling multi-mapping. Examine differences relative to other samples together with the run options.

Why raw counts should not be compared directly

Sequencing the same tissue more deeply increases the counts of most genes together. Longer genes also provide more opportunities for reads to overlap them. Do not directly compare values for different genes or totals for different samples using raw counts alone.

Value	Correction state	Main use
HTSeq raw count	Uncorrected	Input to count models such as DESeq2 and edgeR
DESeq2 normalized count	Corrected for sample depth	Visualization and checking within-sample patterns
TPM	Corrected for gene length and sample depth	Inspecting relative expression composition within one sample

After receiving raw counts, DESeq2 estimates sequencing depth for each sample and variation for each gene within its model. Because a gene has the same length in every sample when comparing that gene across conditions, counts do not need to be converted to TPM before a DEG test. Converting raw counts into comparable values continues with the different uses of normalized values.

What to check in a received count matrix

Are the values nonnegative integers?
Are the row IDs gene_id, gene symbols, or transcript IDs?
Does each column represent one sample, and do the sample names match the metadata?
Are there __ special counters, and are any values unusually large across samples?
Which genome build and GTF annotation version were used?
Is the data single-end or paired-end, and what were the strandedness and overlap mode?
Are these genuine raw counts rather than normalized values?

A public repository’s processed data classification alone does not provide these answers. Open the file and inspect its rows and values, then check both the processing method on the data page and the paper’s Methods section.

Case: reading the GSE251845 file

GSE251845_htseq_raw_counts.csv.gz from GSE251845 is also an HTSeq count matrix with genes as rows and tissue samples as columns. For example, 24C is the sample name for tumor tissue from patient 24. The number where one gene row meets the 24C column is the number of fragments assigned to that gene in the alignment result for the tissue.

Rows such as __no_feature, __ambiguous, and __alignment_not_unique at the end of the file are not genes. For gene-level differential expression analysis, inspect these rows separately and then exclude them from the count matrix.