# Reading HTSeq Raw Counts

> How htseq-count uses BAM and GTF files to assign reads to genes, and how to read a raw count matrix and special counters.

**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:

```text
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.

:::note[Raw count does not mean FASTQ]
FASTQ contains the raw reads produced by a sequencer. An HTSeq raw count is processed data made after FASTQ reads have been aligned and assigned to genes. Here, `raw` means a **count before normalization**.
:::

## 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:

```text
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](/en/reference/sam-bam/) | Where each read aligned in the genome |
| Input | [GTF](/en/reference/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](/en/reference/gtf/) continues with a detailed explanation of the nine GTF columns and `gene_id`.

## The output is a gene × sample matrix

```text
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:

```text
__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.

```python
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](/en/lessons/rna-seq-quantification/) 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.

### Official documentation

- [HTSeq documentation: counting reads by feature with htseq-count](https://htseq.readthedocs.io/en/latest/htseqcount.html)
- [HTSeq tutorial: the internal process for assigning reads to features](https://htseq.readthedocs.io/en/latest/tutorials/exon_example.html)