# From BAM to Counts: A Practical Guide to Gene-Level Quantification with featureCounts


## Key Takeaways

- **Annotation-Alignment Compatibility is Critical:** The GTF/GFF annotation file must precisely match the reference genome build and chromosome naming convention used during BAM alignment; mismatches lead to zero or severely underestimated gene counts.
- **Paired-End and Strandness Settings are Paramount:** Utilizing the `-p` flag for paired-end data ensures fragments, not individual reads, are counted, preventing expression inflation; correctly specifying library strandness (unstranded, stranded, or reverse-stranded) is essential to avoid a large fraction of unassigned reads.
- **Multi-mapping Reads Require Deliberate Handling:** By default, featureCounts discards reads mapping to multiple genomic locations; the `-M` option can include these, with fractional counting (`--fraction`) offering a compromise to avoid overinflating counts in repetitive regions.
- **Summary Statistics are Diagnostic Tools:** The `.summary` file provides crucial metrics on read assignment, including total reads, successfully assigned reads, and reasons for unassigned reads (e.g., multi-mapping, no overlap), enabling rapid identification of annotation or parameter issues.
- **Gene-Level Counts are Standard for Differential Expression:** featureCounts excels at generating gene-level counts, serving as the foundational matrix for differential expression analysis, but it inherently loses isoform-specific information.

---

RNA sequencing experiments generate alignment files in BAM format that contain the genomic positions of sequenced reads. These alignments must be converted into a count matrix that records how many reads map to each gene before differential expression analysis can proceed. featureCounts is a widely used program that performs this gene-level quantification step efficiently. This article explains how to run featureCounts correctly, how to choose annotation files and parameters, how to handle paired-end reads and strandness, and how to produce a count matrix suitable for downstream analysis.

The intended reader is a researcher or laboratory professional who has already aligned RNA-seq reads to a reference genome and now needs to generate gene-level counts. The guidance assumes basic familiarity with the Unix command line and with BAM file formats. The practical outcome is a reproducible command sequence that produces a count matrix with gene identifiers in rows and sample names in columns, ready for import into differential expression tools.

## What featureCounts Does in the RNA-Seq Workflow

The RNA-seq analysis workflow proceeds through several stages. Raw sequencing reads are first assessed for quality, then aligned to a reference genome, and then quantified at the level of genes or transcripts. The quantification step is the bridge between alignment and biological interpretation. Without a count matrix, differential expression analysis cannot begin.

featureCounts belongs to the Subread package and is designed to assign aligned reads to genomic features such as genes, exons, or other annotated regions. The program takes BAM files as input and produces a table of counts. Each row in the output table corresponds to a feature from the annotation file, and each column corresponds to one input sample. The count value represents the number of reads that were successfully assigned to that feature.

The program is widely used because it is computationally efficient and can process multiple samples in a single run. It also provides detailed summary statistics that help researchers assess how well the quantification step performed. These statistics include the total number of reads, the number of successfully assigned reads, and the number of reads that could not be assigned due to various reasons such as ambiguous mapping or lack of overlap with annotated features.

The choice of quantification tool depends on the research question. Gene-level counts are appropriate for standard differential expression analysis where the unit of interest is the gene. Transcript-level quantification is more complex and requires different tools that estimate isoform abundances. For most experiments that ask which genes change expression between conditions, gene-level counts from featureCounts are sufficient and provide a solid foundation for downstream analysis.

## At a Glance

The table below summarizes the key decisions and parameters for running featureCounts in a typical RNA-seq analysis.

| Decision Point | Default Setting | Recommended Practice | Consequence of Getting It Wrong |
|---|---|---|---|
| Annotation file | None, must be supplied | Use GTF or GFF file matching the exact genome build and chromosome naming convention used for alignment | Reads will not be assigned to features, producing zero or very low counts |
| Paired-end handling | Single-end mode | Use `-p` for paired-end data to count fragments instead of individual reads | Fragment counts will be doubled, inflating expression values |
| Strandness | Unstranded (0) | Set to stranded (1) or reverse-stranded (2) according to library preparation protocol | Large fraction of reads will be unassigned, reducing the usable count matrix |
| Multi-mapping reads | Discarded | Keep default for standard analysis, use `-M` for repetitive element studies | Repetitive region expression will be underestimated or overestimated depending on choice |
| Counting mode | Meta-feature level | Use default gene-level counting for differential expression analysis | Counts will not correspond to genes, complicating downstream interpretation |

## Preparing the Input Files

### BAM Files from Alignment

The input to featureCounts is one or more BAM files produced by an aligner. Common aligners include STAR, HISAT2, and Bowtie2. The BAM files must be sorted by genomic position for featureCounts to process them efficiently. Most aligners produce sorted BAM files by default, but if the output is unsorted, the files must be sorted using a tool such as samtools before running featureCounts.

The BAM files should contain only reads that passed the alignment quality filter. Reads that failed to align or that aligned to multiple locations with low quality can be filtered before quantification. However, the decision to filter multi-mapping reads is not always straightforward and depends on the research question. This topic is discussed in more detail in the section on multi-mapping reads.

### Annotation Files in GTF or GFF Format

featureCounts requires an annotation file that defines the genomic coordinates of features. The most common formats are GTF (Gene Transfer Format) and GFF (General Feature Format). These files contain records for genes, transcripts, exons, and other genomic elements. The annotation file must match the reference genome version that was used for alignment. Using an annotation file from a different genome build will produce incorrect counts because the coordinates will not correspond to the alignment positions.

Annotation files are available from several public sources. The [NCBI](https://www.ncbi.nlm.nih.gov/) provides reference genome assemblies and associated annotation files for many organisms. The Ensembl project, hosted by the [European Bioinformatics Institute](https://www.ebi.ac.uk/training), also provides comprehensive annotation files for a wide range of species. These files can be downloaded directly from the respective websites.

The choice of annotation file affects the results. A gene-level annotation file that defines genes as the union of all their exons will produce counts that reflect total gene expression. An annotation file that defines each transcript separately will produce transcript-level counts, which are more complex to interpret. For standard gene-level analysis, the gene-level annotation is appropriate.

### Checking Annotation and Alignment Compatibility

A common source of error is a mismatch between the annotation file and the reference genome used for alignment. Chromosome names must match exactly. For example, if the alignment used a genome where chromosomes are named "chr1", "chr2", and so on, the annotation file must use the same naming convention. If the annotation file uses "1", "2", and so on without the "chr" prefix, featureCounts will not assign any reads to features.

The genome build version must also match. If the alignment was performed against GRCh38 but the annotation file is from GRCh37, the coordinates will be shifted and counts will be incorrect. Checking the genome build and chromosome naming convention before running featureCounts can save considerable time and prevent the need to rerun the analysis.

## Core Principles of Read Assignment

### How featureCounts Assigns Reads to Features

featureCounts works by comparing the genomic coordinates of each aligned read with the coordinates of features in the annotation file. A read is assigned to a feature if it overlaps the feature's genomic interval. The program uses a set of rules to determine which feature receives the read when overlaps are ambiguous.

For paired-end reads, featureCounts can use information from both mates. The default behavior is to count a fragment if either mate overlaps a feature. The program also offers options to require both mates to overlap the same feature or to use the fragment as a single unit. The choice of paired-end handling affects the count values and should be made based on the experimental design.

The program distinguishes between reads that overlap a single feature and reads that overlap multiple features. Reads that overlap multiple features are considered ambiguous and are not counted unless the researcher specifies a different behavior. The summary statistics report the number of ambiguous reads, which provides useful information about the complexity of the annotation and the mapping quality.

### The Role of the Annotation Hierarchy

GTF files contain a hierarchy of features. Genes contain transcripts, and transcripts contain exons. featureCounts can use different levels of this hierarchy for counting. The default is to count at the meta-feature level, which means that exons belonging to the same gene are combined into a single feature. This approach produces gene-level counts.

The program identifies meta-features using the "gene_id" attribute in the GTF file. All exons that share the same gene_id are combined into one meta-feature. The counts reported for each gene_id represent the total number of reads that overlap any exon of that gene.

This approach has an important consequence. Reads that map to regions where two genes overlap are assigned to both genes if the overlap is not resolved by the program's rules. In practice, such overlaps are rare in most genomes, but they can occur in regions with complex gene structures.

### Strandness and Its Effect on Counting

RNA-seq libraries can be prepared in different ways that affect the strand of the sequenced reads. In a stranded library, the reads retain information about which DNA strand was transcribed. In an unstranded library, this information is lost. The strandness of the library must be specified to featureCounts so that it can correctly assign reads to features on the correct strand.

The program offers three settings for strandness: unstranded, stranded, and reverse-stranded. The correct setting depends on the library preparation protocol. Many commercial kits produce reverse-stranded libraries, where the first read in the pair corresponds to the reverse complement of the original RNA. The protocol documentation should specify the strandness of the library.

Using the wrong strandness setting will result in a large fraction of reads being unassigned. The summary statistics will show a high proportion of reads that do not overlap features on the specified strand. If this occurs, the strandness setting should be checked and corrected.

## Running featureCounts

### Basic Command Structure

The basic command to run featureCounts is straightforward. The program takes the annotation file, the BAM files, and a set of options that control the counting behavior. A minimal command looks like this:

```
featureCounts -a annotation.gtf -o counts.txt sample1.bam sample2.bam
```

This command uses the annotation file "annotation.gtf", processes the two BAM files, and writes the count matrix to "counts.txt". The output file contains the count matrix, and a companion file with the ".summary" suffix contains the summary statistics.

The program writes additional files depending on the options used. The count matrix file has a header line that describes the columns, followed by one row per feature. The first six columns contain annotation information such as gene ID, chromosome, start, end, strand, and length. The remaining columns contain the counts for each sample.

### Specifying Paired-End Reads

Paired-end reads require the `-p` option to be specified. Without this option, featureCounts treats each read in a pair as an independent fragment, which will double-count fragments and produce incorrect results. The command for paired-end data is:

```
featureCounts -p -a annotation.gtf -o counts.txt sample1.bam sample2.bam
```

The program also offers options to control how paired-end fragments are counted. The default behavior counts a fragment if either mate overlaps a feature. The `-B` option requires both mates to overlap the same feature. The `-P` option requires both mates to overlap the feature and to be properly paired according to the aligner's definition.

The choice of paired-end options depends on the analysis goals. Requiring both mates to overlap a feature is more stringent and reduces the number of counted fragments. This approach may be appropriate when the alignment quality is uncertain or when the researcher wants to minimize false assignments.

### Selecting the Counting Mode

featureCounts offers several counting modes that control how reads are assigned to features. The default mode counts reads that overlap any exon of a gene. The `-M` option allows multi-mapping reads to be counted. The `-O` option allows reads that overlap multiple features to be assigned to all of them.

The default behavior is to discard multi-mapping reads and reads that overlap multiple features. This conservative approach is appropriate for most analyses because it avoids counting reads that cannot be uniquely assigned. However, some research questions require the inclusion of multi-mapping reads, particularly when studying repetitive elements or gene families with high sequence similarity.

The choice of counting mode should be documented in the methods section of any publication. Different modes will produce different count values, and the results of differential expression analysis can be affected by the choice.

### Running Multiple Samples in One Command

featureCounts can process multiple BAM files in a single command. This approach is efficient because the annotation file is loaded only once. The command lists all BAM files after the options:

```
featureCounts -p -a annotation.gtf -o counts.txt sample1.bam sample2.bam sample3.bam sample4.bam
```

The output file will have one column per sample. The sample names are taken from the BAM file names, with the directory path and the ".bam" extension removed. The researcher can rename the columns in the output file to match the experimental design.

Processing all samples in a single run ensures that the same parameters are used for all samples. This consistency is important for downstream analysis because it avoids introducing batch effects from different parameter settings.

## Handling Paired-End Reads Correctly

### Fragment Counting versus Read Counting

Paired-end sequencing produces two reads from each fragment. The two reads are sequenced from opposite ends of the same fragment. When counting gene expression, the fragment is the biological unit of interest, not the individual reads. Counting each read separately would double-count each fragment and inflate the expression values.

featureCounts handles this by default when the `-p` option is specified. The program treats the two mates as a single fragment and counts the fragment once. The count value for each gene represents the number of fragments that overlap that gene.

The distinction between read counts and fragment counts is important when comparing results across experiments. Some tools report read counts, while others report fragment counts. The methods section of a publication should specify which type of count was used.

### Handling Discordant Pairs

Discordant pairs are read pairs where the two mates align to different chromosomes or in an unexpected orientation. These pairs can arise from structural rearrangements, alignment errors, or chimeric transcripts. featureCounts treats discordant pairs according to the options specified.

By default, the program counts a fragment if either mate overlaps a feature, even if the pair is discordant. The `-P` option requires the pair to be properly paired according to the aligner's definition. Using this option will discard discordant pairs.

The decision to include or exclude discordant pairs depends on the research question. For standard gene expression analysis, discordant pairs are usually excluded because they may represent artifacts. For studies of structural variation or fusion genes, discordant pairs are informative and should be retained.

### The Effect of Insert Size

The insert size is the length of the fragment between the two reads. featureCounts does not require the insert size to be specified, but the program uses the alignment information to determine which features are overlapped by the fragment.

For fragments that span multiple exons, the two mates may align to different exons of the same gene. featureCounts counts the fragment if either mate overlaps an exon of the gene. This approach correctly assigns the fragment to the gene even when the mates align to different exons.

The insert size distribution can be checked after alignment to verify that the library preparation worked correctly. An unusually narrow or wide insert size distribution may indicate a problem with the library preparation or the alignment.

## Understanding Strandness

### Why Strandness Matters

The strandness of an RNA-seq library determines which strand of the DNA the reads correspond to. In a stranded library, the reads can be assigned to the correct strand of the gene. In an unstranded library, a read that maps to a gene could come from either the sense or the antisense strand.

The strandness setting in featureCounts controls which strand is used for counting. The program offers three options: unstranded (0), stranded (1), and reverse-stranded (2). The correct option depends on the library preparation protocol.

Using the wrong strandness setting will produce incorrect counts. In a stranded library, reads from the sense strand will not be counted if the reverse-stranded option is used, and vice versa. The summary statistics will show a high proportion of unassigned reads, which is a clear sign that the strandness setting is wrong.

### Determining the Strandness of Your Library

The strandness of a library is determined by the library preparation protocol. Most commercial kits specify the strandness in their documentation. The researcher should consult the protocol documentation to determine the correct setting.

If the strandness is unknown, it can be inferred from the data. One approach is to run featureCounts with each of the three strandness settings and compare the proportion of assigned reads. The setting that produces the highest proportion of assigned reads is likely the correct one.

Another approach is to examine the alignment of reads to a set of known genes. In a stranded library, reads should align predominantly to the sense strand of genes. In an unstranded library, reads should align equally to both strands.

### Common Strandness Settings for Popular Kits

Many popular library preparation kits produce reverse-stranded libraries. In these libraries, the first read in the pair corresponds to the reverse complement of the RNA. The correct strandness setting for these libraries is reverse-stranded (2).

Some kits produce stranded libraries where the first read corresponds to the RNA sequence itself. The correct setting for these libraries is stranded (1).

Unstranded libraries are less common but are still used in some protocols. The correct setting for these libraries is unstranded (0).

The protocol documentation should be checked to confirm the strandness. If the documentation is unavailable, the strandness can be inferred from the data as described above.

## Managing Multi-Mapping Reads

### What Multi-Mapping Reads Are

Multi-mapping reads are reads that align to multiple locations in the genome with the same or similar quality. These reads arise from repetitive elements, gene families with high sequence similarity, or recent duplication events. The aligner assigns a mapping quality score that reflects the confidence in the alignment.

By default, featureCounts discards multi-mapping reads. The summary statistics report the number of reads that were discarded because they mapped to multiple locations. This conservative approach avoids counting reads that cannot be uniquely assigned to a single gene.

The proportion of multi-mapping reads varies by organism and by genomic region. Genomes with large amounts of repetitive sequence will produce more multi-mapping reads than genomes with less repetitive sequence.

### When to Include Multi-Mapping Reads

Some research questions require the inclusion of multi-mapping reads. Studies of repetitive elements, transposons, or gene families with high sequence similarity cannot rely on uniquely mapped reads alone because most reads from these regions will be multi-mapping. The challenge of quantifying multi-mapping reads is a known limitation in sequence-based analysis, particularly when studying regions with high sequence homology.

The `-M` option in featureCounts allows multi-mapping reads to be counted. When this option is used, a multi-mapping read is counted for every feature that it overlaps. This approach can inflate the counts for genes in repetitive regions, but it provides a more complete picture of expression from these regions.

The decision to include multi-mapping reads should be made based on the research question and should be documented in the methods. Including multi-mapping reads will change the count values and may affect the results of differential expression analysis.

### The Fractional Counting Option

featureCounts offers a fractional counting option that assigns a fraction of a multi-mapping read to each feature it overlaps. The `-M` option combined with the `--fraction` option assigns a weight of 1 divided by the number of features overlapped by the read.

This approach avoids inflating the counts for genes in repetitive regions. Each multi-mapping read contributes a total of one count, distributed across the features it overlaps. The fractional counting option is a compromise between discarding multi-mapping reads entirely and counting them fully for every feature.

The fractional counting option is appropriate when the researcher wants to include multi-mapping reads but does not want to inflate the expression values for genes in repetitive regions. The choice between full counting and fractional counting should be documented in the methods.

## Quality Control and Summary Statistics

### Interpreting the Summary File

featureCounts produces a summary file that reports the number of reads in each category. The summary file has the same base name as the output file with ".summary" appended. The file contains rows for each category and columns for each sample.

The categories include total reads, successfully assigned reads, reads that failed to align, reads that aligned to multiple locations, reads that aligned to no feature, reads that aligned to multiple features, and reads with ambiguous mapping. The proportions of reads in each category provide a quick assessment of the quantification quality.

A high proportion of successfully assigned reads indicates that the annotation file and the parameters are appropriate. A low proportion of assigned reads suggests a problem with the annotation file, the strandness setting, or the alignment quality.

### Expected Proportions of Assigned Reads

The proportion of reads that are successfully assigned to features depends on several factors. These include the quality of the annotation, the strandness setting, the read length, and the proportion of reads from repetitive regions.

For a well-annotated genome and a properly configured run, the proportion of assigned reads is typically high. Values above 70 percent are common for mammalian genomes. Lower values may indicate problems with the annotation or the parameters.

The summary statistics should be examined for every sample. Samples with unusually low assignment rates should be investigated. The cause may be a problem with the library preparation, the alignment, or the quantification parameters.

### Using Summary Statistics for Troubleshooting

The summary statistics provide diagnostic information that can be used to troubleshoot problems. If the proportion of reads that failed to align is high, the alignment step should be reviewed. If the proportion of reads that aligned to no feature is high, the annotation file may be incorrect or the strandness setting may be wrong.

If the proportion of reads that aligned to multiple features is high, the annotation may contain overlapping features. This situation can occur in regions with complex gene structures or when the annotation contains redundant entries.

The summary statistics should be recorded for each analysis run. These records are useful for comparing results across runs and for documenting the quality of the quantification step in publications.

## Output Files and Their Structure

### The Count Matrix File

The count matrix file produced by featureCounts has a specific structure. The first line starts with a "#" character and contains the command that was used to run the program. This line is useful for documentation and reproducibility.

The second line is the header line. It starts with "Geneid" followed by the annotation columns and then the sample names. The annotation columns include "Chr", "Start", "End", "Strand", and "Length". The sample names are derived from the BAM file names.

Each subsequent line corresponds to one feature. The first column contains the gene ID from the annotation file. The next columns contain the chromosome, start position, end position, strand, and length of the feature. The remaining columns contain the counts for each sample.

### The Summary File

The summary file contains the read counts for each category. The first column contains the category names, and the remaining columns contain the counts for each sample. The categories are described in the section on summary statistics.

The summary file is useful for quality control and for documenting the quantification step. The proportions of reads in each category can be calculated from the counts in the summary file.

### Preparing the Count Matrix for Downstream Analysis

The count matrix produced by featureCounts can be imported into R for downstream analysis. The standard approach is to read the file into a data frame and then extract the count columns. The gene IDs in the first column become the row names of the count matrix.

The count matrix should be checked for quality before proceeding with differential expression analysis. Genes with very low counts across all samples are usually filtered out. The filtering threshold depends on the analysis method and the experimental design.

The count matrix should also be checked for sample-level issues. Samples with very low total counts may have failed library preparation or sequencing. These samples should be examined and possibly excluded from the analysis.

## Parameter Selection and Tradeoffs

### Choosing Between Gene-Level and Transcript-Level Counts

The choice between gene-level and transcript-level counts depends on the research question. Gene-level counts are appropriate for standard differential expression analysis where the unit of interest is the gene. Transcript-level counts are needed when the research question involves isoform switching or alternative splicing.

featureCounts produces gene-level counts by default. The program can also produce counts at other levels by using the `-f` option, which counts at the feature level instead of the meta-feature level. Feature-level counts correspond to individual exons or other annotated features.

Transcript-level quantification requires different tools that estimate isoform abundances from the read data. These tools are more complex and require more computational resources than gene-level counting.

### The Effect of Read Length

The read length affects the proportion of reads that can be assigned to features. Longer reads are more likely to span exon boundaries and can be assigned to genes with higher confidence. Shorter reads are more likely to map to repetitive regions or to overlap multiple features.

The read length is determined by the sequencing platform and the library preparation protocol. The researcher cannot change the read length after sequencing, but the read length should be considered when interpreting the results.

### The Effect of Annotation Complexity

The complexity of the annotation affects the proportion of reads that can be assigned to features. Annotations with many overlapping features will produce more ambiguous reads than annotations with well-separated features.

Some annotations contain multiple entries for the same gene, such as when different transcript models are included. These redundant entries can cause reads to be assigned to multiple features. The `-O` option allows reads to be assigned to all overlapping features, which can be useful when the annotation contains redundant entries.

The choice of annotation file is an important decision that affects the results. The annotation should be appropriate for the organism and the research question. Using a comprehensive annotation that includes all known transcripts will produce different results than using a conservative annotation that includes only well-supported transcripts.

## Reproducibility and Documentation

### Recording Parameters for Each Run

The parameters used for each featureCounts run should be recorded for reproducibility. The command line that was used is recorded in the first line of the count matrix file. This record is useful for reproducing the analysis and for documenting the methods in publications.

The parameters that should be recorded include the annotation file, the strandness setting, the paired-end options, and the counting mode. The version of featureCounts and the version of the annotation file should also be recorded.

### Using Workflow Management Tools

Workflow management tools can help ensure that the quantification step is reproducible. These tools automate the execution of analysis steps and record the parameters and inputs for each step.

The [nf-core](https://nf-co.re/docs) project provides community-developed pipelines that include RNA-seq quantification steps. These pipelines follow community standards and are designed to be reproducible. The documentation for nf-core pipelines describes the parameters and the expected inputs and outputs.

[Galaxy](https://training.galaxyproject.org/) provides a web-based platform for running bioinformatics analyses. The Galaxy Training Network offers tutorials that cover RNA-seq analysis, including quantification with featureCounts. These tutorials provide step-by-step instructions that can be followed by researchers who prefer a graphical interface.

### Version Control for Analysis Code

The code used for the analysis should be under version control. This practice ensures that changes to the analysis code are tracked and that the exact code used for a particular analysis can be recovered.

The [Carpentries](https://carpentries.org/lessons) offers lessons on version control with Git. These lessons cover the basics of tracking changes to files and collaborating with others. Version control is an essential component of reproducible research.

## Common Failure Patterns and Troubleshooting

### Low Assignment Rates

A low proportion of assigned reads is the most common problem in gene-level quantification. The cause is often a mismatch between the annotation file and the alignment. The chromosome naming convention and the genome build should be checked first.

The strandness setting is another common cause of low assignment rates. If the library is stranded and the wrong strandness setting is used, most reads will not be assigned to features on the correct strand. The summary statistics will show a high proportion of reads that aligned to no feature.

The annotation file may also be incorrect or incomplete. The annotation should be checked to ensure that it covers the regions where the reads align. If the annotation is missing genes or exons, reads that align to those regions will not be assigned.

### Unexpected Count Distributions

The distribution of counts across genes can reveal problems with the quantification. Genes with extremely high counts may represent artifacts such as reads from ribosomal RNA or other abundant contaminants. Genes with zero counts across all samples may be absent from the annotation or may not be expressed in the samples.

The count distribution should be examined for each sample. Samples with unusual distributions may have failed library preparation or sequencing. These samples should be investigated before proceeding with downstream analysis.

### Discrepancies Between Replicates

Biological replicates should produce similar count distributions. Large discrepancies between replicates may indicate a problem with one of the samples. The summary statistics for each sample should be compared to identify samples with unusual assignment rates.

The count matrix can be examined for sample-level outliers. Principal component analysis or hierarchical clustering can reveal samples that are very different from the others. These samples should be investigated before proceeding with differential expression analysis.

## Limitations of Gene-Level Quantification

### Loss of Isoform Information

Gene-level quantification does not distinguish between different isoforms of the same gene. A gene with multiple isoforms will produce a single count that represents the total expression of all isoforms. This limitation is acceptable for many research questions but is important to recognize.

If the research question involves isoform switching or alternative splicing, gene-level counts are not sufficient. Transcript-level quantification is needed to estimate the abundance of each isoform separately. Resources such as [Cortexa](https://doi.org/10.1186/s12859-024-05919-y) demonstrate the value of integrating expression and alternative splicing data for understanding complex regulatory processes, but such analyses require specialized processing beyond standard gene-level counting.

### Ambiguous Read Assignment

Some reads cannot be uniquely assigned to a single gene. These reads may overlap multiple genes or may map to regions that are not annotated. The proportion of ambiguous reads depends on the complexity of the annotation and the read length.

The default behavior of featureCounts is to discard ambiguous reads. This conservative approach avoids inflating the counts for genes in complex regions. However, the discarded reads may contain information that is relevant to the research question.

### The Effect of Annotation Quality

The quality of the annotation directly affects the quality of the counts. An annotation that is incomplete or contains errors will produce counts that do not accurately reflect gene expression.

The annotation should be appropriate for the organism and the research question. For well-studied organisms, the annotation is generally of high quality. For less-studied organisms, the annotation may be incomplete, and the counts should be interpreted with caution.

## Safety and Regulatory Context

### Data Management and Privacy

RNA-seq data may contain sensitive information, particularly if the samples come from human subjects. The data should be managed according to the relevant regulations and institutional policies. The [NCBI](https://www.ncbi.nlm.nih.gov/) provides resources for data submission and access that follow established standards.

Researchers should be aware of the data sharing requirements of their funding agencies and journals. Many journals require that sequencing data be deposited in public databases. The data should be prepared for submission according to the database requirements.

### Computational Resource Considerations

featureCounts is computationally efficient, but the analysis of large datasets requires adequate computational resources. The memory and disk space requirements depend on the number of samples and the size of the genome.

The analysis should be run on a machine with sufficient resources to avoid crashes or excessive runtime. The resource requirements should be considered when planning the analysis.

### Professional Escalation Criteria

Some problems require professional escalation. If the proportion of assigned reads is very low and the cause cannot be identified, the analysis should be reviewed by a bioinformatics specialist. If the annotation file appears to be incorrect, the annotation provider should be contacted.

If the results of the quantification are inconsistent with the experimental design, the analysis should be reviewed before proceeding with downstream analysis. The review should include an examination of the alignment quality, the annotation file, and the quantification parameters.

## Frequently Asked Questions

### What is the difference between gene-level and transcript-level counts?

Gene-level counts represent the total expression of a gene, combining all of its isoforms. Transcript-level counts represent the expression of individual isoforms. Gene-level counts are produced by featureCounts by default and are appropriate for standard differential expression analysis. Transcript-level counts require specialized tools that estimate isoform abundances and are needed when the research question involves isoform switching or alternative splicing.

### How do I know which strandness setting to use?

The strandness setting depends on the library preparation protocol. Most commercial kits specify the strandness in their documentation. If the strandness is unknown, run featureCounts with each of the three settings and compare the proportion of assigned reads. The setting that produces the highest proportion of assigned reads is likely the correct one.

### Should I include multi-mapping reads in my analysis?

The decision to include multi-mapping reads depends on the research question. The default behavior is to discard multi-mapping reads, which is appropriate for most analyses. Studies of repetitive elements or gene families with high sequence similarity may require the inclusion of multi-mapping reads. The `-M` option enables this behavior, and the `--fraction` option assigns fractional counts to avoid inflating expression values.

### Why is the proportion of assigned reads low?

A low proportion of assigned reads can have several causes. The most common causes are a mismatch between the annotation file and the alignment, an incorrect strandness setting, and an incomplete annotation. Check the chromosome naming convention and the genome build first, then verify the strandness setting, and finally examine the annotation file for completeness.

### Can I run featureCounts on multiple samples at once?

Yes, featureCounts can process multiple BAM files in a single command. List all BAM files after the options. The output file will have one column per sample. Processing all samples in a single run ensures that the same parameters are used for all samples, which is important for consistency.

### What should I do if the counts for my replicates are very different?

Large discrepancies between replicates may indicate a problem with one of the samples. Compare the summary statistics for each sample to identify samples with unusual assignment rates. Examine the count distributions and consider whether the sample should be excluded from the analysis.

### How do I prepare the count matrix for differential expression analysis?

The count matrix produced by featureCounts can be imported into R for downstream analysis. Read the file into a data frame and extract the count columns. Filter out genes with very low counts across all samples. Check the count matrix for sample-level issues before proceeding with differential expression analysis.

### Where can I find annotation files for my organism?

Annotation files are available from several public sources. The [NCBI](https://www.ncbi.nlm.nih.gov/) provides reference genome assemblies and associated annotation files for many organisms. The Ensembl project, hosted by the [European Bioinformatics Institute](https://www.ebi.ac.uk/training), also provides comprehensive annotation files. The annotation file must match the reference genome version used for alignment.

## Related Bioinformatics Guides

- [RNA-Seq vs ChIP-Seq: Complementary Approaches for Gene Regulation](/knowledge/bioinformatics/rna-seq-vs-chip-seq-complementary-approaches-for-gene-regulation)
- [Gene Set Enrichment Analysis in R: A Practical Tutorial for Interpreting Omics Data](/knowledge/bioinformatics/gene-set-enrichment-analysis-in-r-a-practical-tutorial-for-interpreting-omics-data)
- [RNA-Seq vs Microarray: Choosing the Right Gene Expression Profiling Platform](/knowledge/bioinformatics/rna-seq-vs-microarray-choosing-the-right-gene-expression-profiling-platform)
- [Single-Cell RNA Sequencing Quality Control: A Practical Guide to Filtering and Metrics](/knowledge/bioinformatics/single-cell-rna-sequencing-quality-control-a-practical-guide-to-filtering-and-metrics)
- [RNA-Seq vs qPCR: Validation and Comparison](/knowledge/bioinformatics/rna-seq-vs-qpcr-validation-and-comparison)

## Related Clinical & Scientific Guides

* [A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data](/knowledge/bioinformatics/a-practical-guide-to-detecting-antimicrobial-resistance-genes-in-shotgun-metagenomic-data)
* [Computational Immunology: Modeling the Immune System](/knowledge/bioinformatics/computational-immunology-modeling-the-immune-system)
* [How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices](/knowledge/bioinformatics/how-to-set-hard-filters-for-germline-variant-calling-a-practical-guide-to-gatk-best-practices)


## References and Further Reading

- [NCBI Data Resources](https://www.ncbi.nlm.nih.gov/). National Center for Biotechnology Information.
- [EMBL-EBI Training](https://www.ebi.ac.uk/training). European Bioinformatics Institute.
- [Bioconductor](https://bioconductor.org/). Bioconductor Project.
- [Galaxy Training Network](https://training.galaxyproject.org/). Galaxy Project.
- [nf-core Documentation](https://nf-co.re/docs). nf-core.
- [The Carpentries Lessons](https://carpentries.org/lessons). The Carpentries.
- [HERVOminer: a sequence similarity-based approach for recognizing endogenous retrovirus origin of the peptidome.](https://doi.org/10.1038/s41698-026-01370-9). 2026.
- [Cortexa: a comprehensive resource for studying gene expression and alternative splicing in the murine brain.](https://doi.org/10.1186/s12859-024-05919-y). 2024.
- [Integrated omics reveals disease-associated radial glia-like cells with epigenetically dysregulated interferon response in multiple sclerosis.](https://doi.org/10.1016/j.neuron.2025.09.022). 2025.
- [Kif15 orchestrates neuronal-microglial communication via CX3CL1 to impede nerve regeneration.](https://doi.org/10.1016/j.jbc.2026.113090). 2026.
- [Reproducibility of PD patient-specific midbrain organoid data for &lt,i&gt,in vitro&lt,/i&gt, disease modeling.](https://doi.org/10.1016/j.isci.2025.113541). 2025.
- [A cell-intrinsic glucocorticoid biosynthesis and sensing circuit maintains a homeostatic Th17 cell state.](https://doi.org/10.1016/j.immuni.2026.05.015). 2026.
- [Tutorial: RNA-seq differential expression & pathway analysis with Sailfish, DESeq2, GAGE, and Pathview](https://doi.org/10.6084/M9.FIGSHARE.1619655.V1). 2015.
- [Current best practices in single-cell RNA-seq analysis: a tutorial](https://doi.org/10.15252/msb.20188746). Molecular Systems Biology, 2019.
- [Galaxy Training Material for single-cell RNA-seq tutorial with Plant Datasets](https://doi.org/10.5281/ZENODO.4597857). 2021.
- [Become Competent in Generating RNA-Seq Heat Maps in One Day for Novices Without Prior R Experience.](https://doi.org/10.1007/978-1-0716-1084-8_17). Methods in molecular biology, 2021.
- [Variant Calling from RNA-seq Data Using the GATK Joint Genotyping Workflow.](https://doi.org/10.1007/978-1-0716-2293-3_13). Methods in molecular biology, 2022.

> This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.