STAR quantMode GeneCounts: The Complete Guide to Output, Analysis, and Best Practices

Abstract

This comprehensive guide demystifies the STAR aligner's quantMode GeneCounts output, a critical step in RNA-seq analysis for differential expression. Tailored for researchers and bioinformaticians, we break down the ReadsPerGene.out.tab file structure, explain strandedness selection, and provide a step-by-step workflow for accurate gene counting. The article further tackles common troubleshooting issues, compares STAR's counting to dedicated tools like featureCounts and Salmon, and offers evidence-based recommendations to ensure robust, reproducible results for downstream biomedical and clinical research applications.

Decoding the ReadsPerGene.out.tab File: A Beginner's Guide to STAR's Count Matrix

In the analysis of RNA sequencing (RNA-seq) data, accurate quantification of gene expression is foundational for downstream biological interpretation. The STAR aligner, a widely used tool for mapping RNA-seq reads, offers integrated quantification via the --quantMode GeneCounts parameter. This functionality generates a ReadsPerGene.out.tab file containing four distinct columns of read counts, a feature that is leveraged in large-scale transcriptomic studies and atlases [1] [2]. However, selecting the correct column is crucial, as an incorrect choice can lead to a significant overestimation of background noise or a complete misinterpretation of expression data [3] [4]. This application note provides a detailed guide to the purpose, interpretation, and proper selection of these four columns, framed within robust experimental protocols to ensure analysis accuracy.

Decoding the ReadsPerGene.out.tab File

Column Definitions and Specifications

The ReadsPerGene.out.tab file is the direct output of STAR when run with the --quantMode GeneCounts option. Its structure is standardized, with the first four lines dedicated to summary counts and all subsequent lines dedicated to individual gene counts [4] [5].

Summary Lines (First Four Rows):

• N_unmapped: Total number of unmapped reads.
• N_multimapping: Total number of reads mapped to multiple genomic locations.
• N_noFeature: Total number of reads that did not overlap any gene feature defined in the provided annotation file (GTF/GFF).
• N_ambiguous: Total number of reads that overlapped more than one gene, making their assignment ambiguous [4] [5].

Gene Counts Table:

Table 1: Description of the four columns in the ReadsPerGene.out.tab file.

Column Index	Description	Corresponding htseq-count option	Primary Use Case
1	Gene identifier	-	-
2	Counts for unstranded RNA-seq	-s no	Unstranded (non-strand-specific) libraries
3	Counts for the 1st read strand aligned with RNA	-s yes	Stranded protocols where the first read is sense to the RNA
4	Counts for the 2nd read strand aligned with RNA	-s reverse	Stranded protocols where the second read is sense to the RNA

The Critical Distinction: Mapping vs. Quantification

A fundamental and often misunderstood concept is that STAR's mapping is strand-agnostic. The aligner finds the optimal genomic location for a read regardless of its strand origin. However, the quantification step that follows is entirely strand-aware [3]. The aligner records the strand to which each read was mapped in the BAM file. During quantification, STAR uses the provided gene annotation (GTF file) to determine the canonical strand of a gene. It then cross-references this with the read's mapped strand to assign it to the correct column (column 3 or column 4) based on the library preparation protocol specified by the user [3] [5].

Experimental Protocol: Determining Library Strandedness

Selecting the correct column is impossible without knowing the library preparation kit used. The following protocol provides a data-driven method to confirm the library type.

Protocol 1: Empirical Determination of Strand-Specificity

Purpose: To experimentally verify the strandedness of an RNA-seq library by analyzing the distribution of counts across the four columns.

Principle: In a perfectly stranded library, one of the strand-specific columns (3 or 4) will contain the vast majority of counts aligning to known genes, while the other will contain primarily antisense or background signal [2] [5].

Procedure:

• Run STAR Alignment: Process a representative sample through STAR with the --quantMode GeneCounts option and a comprehensive gene annotation file.
• Generate Count Summary: Execute the following Bash command to sum the total counts for all genes (excluding the summary lines) in each of the three count columns:
  This command works by skipping all lines starting with "N_" and then accumulating the values from columns 2, 3, and 4.
• Interpret Results: Analyze the summed counts [5]:
  • If the Forward Strand and Reverse Strand counts are roughly equal, the library is unstranded (use column 2).
  • If one strand-specific count (e.g., Reverse Strand) vastly exceeds the other, the library is stranded.
  • Compare the high-count strand to your kit's manual (e.g., Illumina TruSeq Stranded Total RNA: use column 4).

Troubleshooting:

• Unexpectedly high antisense counts: This can indicate a library preparation issue or genuine biological antisense transcription, but the dominant strand should still be clear.
• Discrepancies with other tools: While STAR's GeneCounts are designed to match htseq-count defaults, slight variations can occur due to different implementations of overlap rules and handling of ambiguous reads [6]. For critical analyses, consistency in the tool used for final quantification is recommended.

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Key reagents and computational tools for stranded RNA-seq analysis.

Item Name	Function / Description	Example / Specs
Stranded mRNA Library Prep Kit	Isolates mRNA and preserves strand-of-origin information during cDNA synthesis.	Illumina TruSeq Stranded mRNA, NEBNext Ultra II Directional RNA
Reference Genome (FASTA)	The DNA sequence of the organism used as a scaffold for read alignment.	Homosapiens.GRCh38.dna.primaryassembly.fa
Gene Annotation (GTF)	File containing genomic coordinates and strands of all known genes and transcripts.	Gencode v44 comprehensive annotation, Ensembl annotation
STAR Aligner	Spliced Transcripts Alignment to a Reference; performs alignment and quantification.	Version 2.7.10b; requires ~30GB RAM for human genome [7]
High-Performance Computing	Compute instance for resource-intensive alignment steps.	AWS EC2 (e.g., r5ad.2xlarge), HPC cluster [7]

Relationship Between Library Protocol and Count Columns

The following diagram illustrates the logical workflow of how a stranded library preparation protocol determines which of the four columns in the STAR output should be used for final analysis.

Integrated Workflow for Accurate Gene Quantification

This comprehensive workflow chart outlines the key steps from raw sequencing data to a finalized count matrix for differential expression analysis, highlighting the critical decision point for column selection.

The four-column output of STAR's --quantMode GeneCounts is a powerful feature that accommodates both unstranded and stranded RNA-seq libraries. A rigorous understanding of the library preparation chemistry is the most critical factor for accurate analysis. By following the experimental protocols outlined herein—specifically, empirically verifying strandedness and meticulously consulting kit documentation—researchers can confidently select the appropriate count column. This ensures the integrity of gene expression data, forming a reliable foundation for all subsequent differential expression analysis and biological discovery in transcriptomics research.

The STAR (Spliced Transcripts Alignment to a Reference) aligner, when run with the --quantMode GeneCounts parameter, performs read alignment and quantifies gene expression simultaneously. A critical component of its output is the *ReadsPerGene.out.tab file, which begins with several header lines that provide essential summary statistics about the mapping and quantification process [4] [5]. These header lines—N_unmapped, N_multimapping, N_noFeature, and N_ambiguous—offer researchers invaluable metrics for assessing data quality, understanding library preparation properties, and troubleshooting potential issues in RNA-seq experiments [4] [8]. For researchers in drug development and biomedical research, proper interpretation of these metrics is crucial for ensuring the reliability of downstream analyses such as differential expression testing and biomarker identification.

The counting methodology implemented by STAR's GeneCounts follows the same principles as htseq-count with default parameters, using the "union" mode for resolving overlapping features [8] [9]. A read is counted toward a gene if it overlaps (by at least 1 nucleotide) with one and only one gene as defined in the provided annotation file (GTF or GFF) [4] [5]. For paired-end reads, both ends are checked for overlaps with genomic features [5]. The resulting count table includes four columns: gene identifiers followed by three columns corresponding to different strandedness options, enabling researchers to select counts appropriate for their specific library preparation protocol [4] [9].

Comprehensive Explanation of Header Lines

Definition and Interpretation of Each Metric

The header lines in the ReadsPerGene.out.tab file represent different categories of reads that could not be uniquely assigned to genes based on the provided annotation and mapping parameters. Understanding each category is essential for quality assessment and troubleshooting.

N_unmapped: This metric represents reads that could not be aligned to the reference genome under the specified mapping parameters [4] [5]. Unmapped reads can result from various factors including sequencing artifacts, adapter contamination, or the presence of novel sequences not represented in the reference genome. In typical RNA-seq experiments, the unmapped fraction should generally be a small minority of total reads, though this can vary depending on sample quality and reference genome completeness.

N_multimapping: These reads align to multiple locations in the genome with equally good or nearly equally good alignment scores [4] [10]. By default, with standard parameter settings, STAR does not count any alignments from multimapping reads toward gene counts and aggregates them all under this header line [10]. The multimapping fraction tends to be higher in genomes with many paralogous genes or repetitive elements. When using --outFilterMultimapNmax 1 (allowing only uniquely mapped reads), this value should be zero as all multimapping reads would be filtered out [4] [10].

N_noFeature: This category consists of reads that aligned successfully to the genome but do not overlap any annotated genomic feature (e.g., exon) in the provided GTF file [4] [5]. High values typically indicate either incomplete annotation, novel transcripts not present in the annotation, or potential genomic DNA contamination. In stranded protocols, the column with the lowest N_noFeature count usually indicates the correct strand orientation [10].

N_ambiguous: These reads overlap with multiple genes or features and therefore cannot be uniquely assigned to a single gene [4] [5]. This ambiguity arises when genes overlap in genomic coordinates or when reads span exon-exon junctions that are shared between different genes. The "union" counting mode (default for both STAR and htseq-count) resolves these ambiguous cases conservatively by not counting these reads toward any specific gene [4] [8].

Quantitative Relationships Between Metrics

The following table summarizes the key characteristics, potential causes, and optimization strategies for each header metric:

Table 1: Comprehensive Guide to STAR GeneCounts Header Metrics

Metric	Definition	Common Causes	Typical Range	Troubleshooting Approaches
N_unmapped	Reads failing to align to reference	Poor read quality, adapter contamination, novel sequences, incorrect reference	5-15% of total reads	Quality trimming, adapter removal, verify reference compatibility
N_multimapping	Reads mapping to multiple genomic locations	Repetitive regions, paralogous genes, transposable elements	10-30% of total reads [11]	Increase --outFilterMultimapNmax, use transcript-aware aligners like Salmon [11]
N_noFeature	Aligned reads not overlapping annotated features	Incomplete annotation, novel transcripts, genomic DNA contamination	Variable by annotation quality	Use comprehensive annotation (e.g., GENCODE over UCSC [11]), assess gDNA contamination
N_ambiguous	Reads overlapping multiple genes	Overlapping genes, shared exons or junctions	5-15% of aligned reads	Review annotation quality, consider different counting strategies [11]

The relationship between these metrics and total reads follows a logical hierarchy. The sum of all header line values plus the total assigned reads (sum of all gene counts) should approximately equal the total number of input reads, though slight discrepancies may occur due to read filtering during alignment [8]. Specifically:

Total Reads ≈ Nunmapped + Nmultimapping + NnoFeature + Nambiguous + Assigned Reads

It's important to note that different counting tools may categorize reads differently. For instance, featureCounts reports "UnassignedMultiMapping" and "UnassignedAmbiguity" separately, similar to STAR's Nmultimapping and Nambiguous, though the exact classification algorithms may vary slightly between tools [11].

Experimental Protocols for Data Interpretation

Protocol 1: Determining Library Strandedness

Purpose: To determine the strandedness of RNA-seq libraries using STAR GeneCounts output, which is crucial for selecting the appropriate count column for downstream analysis.

Materials and Reagents:

• STAR-aligned BAM files
• ReadsPerGene.out.tab files from STAR with --quantMode GeneCounts
• Bioinformatics workstation with command-line tools
• Reference annotation file (GTF format)

Procedure:

• Execute STAR alignment with the --quantMode GeneCounts parameter to generate the ReadsPerGene.out.tab files [5] [9].

• Extract summary statistics from the output file using command-line tools:

• Calculate total assigned reads for each strandedness column:

• Compute N_genic values for columns 3 and 4 using the formula provided by STAR developer Alexander Dobin [10]: N_genic = TotalReads - N_unmapped - N_multimapping - N_noFeature - N_ambiguous

• Compare Ngenic values between columns 3 and 4. The column with the higher Ngenic value indicates the correct strand orientation for your library [10].

• Verify with known strand-specific genes by checking the distribution of counts for genes with known strand orientation.

Interpretation: The strandedness protocol with the highest Ngenic value and lowest NnoFeature count represents the correct strand orientation for your data [10]. For example, in a study using human A549 cells, researchers observed 153,677 reads in column 3 versus 2,427,536 reads in column 4, clearly indicating a stranded protocol where the reverse column (4) represented the correct strand orientation [5].

Protocol 2: Quality Assessment and Troubleshooting

Purpose: To identify potential issues in RNA-seq experiments by analyzing the distribution of reads across different header metrics.

Materials and Reagents:

• MultiQC tool for aggregating results across multiple samples [8]
• R or Python environment for statistical analysis
• Complete set of ReadsPerGene.out.tab files from all samples

Procedure:

• Generate summary statistics for all samples using MultiQC or custom scripts:

• Create a summary table of all header metrics across samples to identify outliers.

• Calculate percentage distributions for each metric relative to total reads:

• Identify samples with abnormal distributions using pre-defined thresholds:

  • >30% unmapped reads: Potential reference genome or quality issues
  • >40% multimapping: Possible repetitive elements or incomplete genome
  • >50% noFeature: Likely annotation issues or novel transcription
  • >20% ambiguous: Gene model overlaps or complex regions

• Investigate root causes based on the identified issues:

  • For high noFeature: Compare with full annotation sets (GENCODE recommended over UCSC) [11]
  • For high multimapping: Consider allowing multi-mapped reads to contribute to counts with fractional weights [11]
  • For high ambiguous: Evaluate whether current annotation is appropriate for your biological system

Table 2: Research Reagent Solutions for RNA-seq Quality Assessment

Reagent/Resource	Function	Usage Notes
STAR Aligner	Spliced transcript alignment	Use latest version; requires significant RAM for large genomes [5]
GENCODE Annotations	Comprehensive gene annotation	Preferred over UCSC for human/mouse; more complete [11]
MultiQC	Aggregate bioinformatics results	Visualize trends across multiple samples [8]
FastQC	Read quality control	Run before alignment to identify sequencing issues
RSeQC	RNA-seq quality metrics	Provides additional quality metrics like read distribution

Interpretation: A typical high-quality RNA-seq dataset should show relatively consistent proportions across samples for each metric. Significant deviations may indicate technical artifacts that could bias downstream analyses. For example, if most samples have 10-15% multimapping reads but one sample shows 40%, this outlier may require additional investigation or exclusion from downstream analysis.

Visualization and Data Analysis Workflows

To better understand the relationship between read types and the decision process STAR employs during read assignment, the following workflow diagram illustrates the logical sequence:

Diagram 1: STAR Read Assignment Workflow

The strandedness determination protocol can be visualized as a decision flow to guide researchers in selecting the appropriate count column:

Diagram 2: Strandedness Determination Protocol

Advanced Applications in Research Scenarios

Integration with Downstream Analysis Tools

The header metrics from STAR's GeneCounts output provide crucial quality indicators that should inform downstream analytical decisions. In drug development contexts, where reproducibility is paramount, these metrics serve as key quality control checkpoints.

For differential expression analysis, the choice of count column (unstranded, column 3, or column 4) must align with the experimental library preparation method [4] [10]. Selection of the wrong column can introduce substantial noise and potentially mask true biological signals. When analyzing the header metrics:

• For unstranded libraries: Use column 2 counts, recognizing that antisense transcription will be counted toward genes and cannot be distinguished [4] [5].

• For stranded libraries: Use either column 3 or 4 based on the determined strandedness, which provides more accurate quantification by distinguishing sense from antisense transcription [10].

• For complex study designs: Consider that the "antisense" counts (the column not used for primary analysis) may provide valuable information about regulatory antisense transcription relevant to drug mechanisms [4].

In scenarios where novel transcript discovery is a study goal, high N_noFeature rates may not necessarily indicate problems but rather biological novelty. In such cases, researchers can use tools like StringTie or Cufflinks to assemble novel transcripts, then create a merged annotation file and requantify using STAR or other count tools [12] [13] [14].

Troubleshooting Common Issues

Several common issues can be diagnosed through abnormal patterns in the header metrics:

High Nmultimapping in featureCounts: When transitioning from STAR to featureCounts, researchers may observe high "UnassignedMultiMapping" counts [11]. This occurs because featureCounts, by default, does not count multimapping reads. Solution strategies include:

• Allowing reads to map to multiple features in featureCounts (though this requires careful interpretation)
• Using transcript-level quantification tools like Salmon or RSEM that probabilistically assign multimapping reads [11]
• Increasing the --outFilterMultimapNmax parameter in STAR to retain more multimapping reads for downstream analysis

Unexpected relationships between columns: In some cases, values in column 3 or 4 may exceed those in column 2, which seems counterintuitive [4]. This occurs because the unstranded counts (column 2) may classify some reads as ambiguous that would be assigned to genes in the stranded columns (3 or 4) [4]. Specifically, reads mapping to the opposite strand of a gene are considered ambiguous in unstranded mode but are assigned to genes in the appropriate stranded column, leading to this apparent discrepancy [4].

Discrepancies between STAR and htseq-count: While STAR's GeneCounts and htseq-count should theoretically produce identical results, small differences can occur due to variations in how filtered reads are handled or different implementations of the counting algorithm [8]. For most consistent results, use the same annotation file that was used during STAR genome indexing.

In RNA sequencing (RNA-seq) analysis, a fundamental step involves quantifying gene expression by counting how many reads map to each gene. The accuracy of downstream differential expression analysis depends heavily on the method used for this assignment. The STAR aligner (Spliced Transcripts Alignment to a Reference) incorporates a built-in counting feature enabled by the --quantMode GeneCounts parameter, which provides a streamlined workflow for gene expression quantification. According to the STAR manual, this feature is designed to produce counts that "coincide with those produced by htseq-count with default parameters" [4]. This application note details the underlying assignment mechanism of STAR, specifically its 1-nucleotide overlap rule, and provides a comprehensive comparison with the popular htseq-count tool, offering practical guidance for researchers and drug development professionals implementing these methods in their transcriptomic analyses.

Core Principles of Read Assignment

STAR's 1-Nucleotide Overlap Rule

STAR employs a specific, well-defined rule for assigning reads to genomic features. When run with the --quantMode GeneCounts option, STAR counts a read toward a gene if it overlaps one and only one gene by at least 1 nucleotide (1nt) [4]. For paired-end sequencing data, this assignment rule considers both ends of the read pair, requiring that both reads meet the overlap criterion with the same single gene.

This implementation mirrors the default behavior of htseq-count, which uses similar logic to resolve feature assignments [4]. The 1nt threshold represents the minimal possible overlap that still indicates a potential biological relationship between the read and the gene.

Output Structure and Strandedness

STAR generates a *ReadsPerGene.out.tab file containing four columns that accommodate different RNA-seq library preparation protocols [4] [5]:

• Column 1: Gene identifiers
• Column 2: Counts for unstranded RNA-seq
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count -s yes)
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count -s reverse)

The appropriate column selection depends entirely on the strandedness of the RNA-seq library preparation protocol used in the experimental design. For non-strand-specific protocols (unstranded), column 2 should be used. For strand-specific protocols, either column 3 or 4 is appropriate based on the specific library construction method [4] [5].

Table 1: Guide to Selecting the Correct Count Column in STAR Output

Library Type	STAR Column	htseq-count Strandedness	Description
Unstranded	2	--stranded=no	Reads counted regardless of strand
Strand-specific	3	--stranded=yes	1st read strand aligned with RNA
Strand-specific	4	--stranded=reverse	2nd read strand aligned with RNA

STAR vs. htseq-count: Implementation Comparison

Algorithmic Similarities and Differences

While both tools aim to solve the same fundamental problem, their implementations differ in several key aspects:

STAR performs read counting during the mapping process, leveraging its internal representation of genomic annotations [4]. This integrated approach offers computational efficiency, as it eliminates the need for a separate counting step after alignment completion.

In contrast, htseq-count operates as a post-alignment tool that processes SAM/BAM files generated by any splice-aware aligner, including STAR [15] [16]. This separation of alignment and counting provides flexibility but requires additional processing time and storage resources.

Both tools default to discarding reads that overlap multiple features (counting them as "ambiguous"), which ensures that expression ratios between samples remain accurate for differential expression analysis, even if total counts are slightly reduced [15].

Handling of Ambiguous Reads

A critical aspect of read counting involves handling reads that ambiguously map to multiple genomic features. The following diagram illustrates the decision logic both tools employ for read assignment:

Diagram 1: Logical workflow for read assignment in STAR and htseq-count showing how reads are categorized based on their overlap with genomic features.

For reads that map to multiple genomic locations (multi-mappers), STAR's default behavior reports the number of reads (not alignments) that mapped uniquely, while htseq-count reports the number of alignments in the SAM file that are of reads that mapped to more than one place [17]. This conceptual difference in counting philosophy can lead to apparent discrepancies in summary statistics, though both approaches maintain biological validity for different analysis goals.

Quantitative Comparison of Output

Under optimal conditions with proper parameter matching, STAR's --quantMode GeneCounts and htseq-count should produce identical or highly similar results [4]. However, several factors can cause discrepancies:

Table 2: Factors Causing Count Discrepancies Between STAR and htseq-count

Factor	Impact on Counts	Solution
BAM file sorting	htseq-count expects name-sorted BAM by default; position-sorted can cause undercounting	Use samtools sort -n for name sorting [18]
Strandedness parameter	Mismatch between tool settings and library type	Verify library protocol and set --stranded accordingly [15]
Annotation version	Different GTF/GFF files used in alignment vs. counting	Use identical annotations for both steps
Alignment filtering	Different MAPQ thresholds or filtering criteria	Ensure consistent alignment quality thresholds

Research indicates that when the same annotations are used and BAM files are properly name-sorted, counts from STAR and htseq-count show strong concordance [18]. The integrated counting approach of STAR offers advantages in processing time and simplicity, while htseq-count provides greater flexibility for testing different counting strategies and parameters.

Experimental Protocol for Read Counting with STAR

Genome Indexing and Alignment

A robust STAR workflow begins with proper genome index generation, which incorporates gene annotations directly into the index structure:

The --sjdbOverhang parameter should be set to the maximum read length minus 1, which for typical 51bp reads is 50 [5]. This parameter defines the length of the genomic sequence around annotated junctions used for alignment and affects splice junction detection accuracy.

Following index generation, alignment with integrated counting proceeds with:

This command generates both a sorted BAM file and the ReadsPerGene.out.tab file containing gene counts, streamlining the workflow from alignment to count matrix generation.

Strandedness Determination

For researchers uncertain about their library's strandedness, the following diagnostic approach using STAR output can help:

Extreme imbalances between columns 3 and 4 indicate strand-specific data, while relatively balanced counts suggest unstranded data [5]. For example, in one analysis, stranded data showed 153,677 forward reads versus 2,427,536 reverse reads mapped to known genes, clearly indicating a stranded protocol where the reverse complement of mRNA was sequenced [5].

Advanced Considerations in Read Assignment

Handling Overlapping Genes

Genomic regions with overlapping genes present particular challenges for read assignment. When a read maps to a location where multiple genes overlap, both STAR and htseq-count (with default settings) will classify it as "ambiguous" and exclude it from counts for any gene [4] [19]. This conservative approach ensures that expression estimates aren't inflated by double-counting, but may lead to undercounting for genes in densely annotated regions.

In eukaryotic genomes, approximately 18% of genes may participate in overlapping transcription units, with average overlapping lengths of 290bp in yeast and exceeding 1kb in mammalian genomes [19]. Reads shorter than these overlapping regions cannot be uniquely assigned, creating quantification challenges, particularly for non-coding RNAs and antisense transcripts.

Specialized tools like IAOseq and MGcount have been developed to address these limitations using probabilistic models that consider reads distribution along transcribed regions [19] [20]. These methods can help resolve expression levels for overlapping genes, though they introduce additional computational complexity.

Special Cases and Edge Conditions

Several special scenarios require researcher attention when interpreting counting results:

• Antisense transcription: When using stranded protocols, the column not selected for analysis (column 3 or 4) represents antisense reads, which may indicate regulatory antisense transcription [4].

• Multi-mapping reads: Reads that align to multiple genomic locations are reported separately from overlapping assignments, and their handling differs between tools [17] [20].

• Incomplete annotations: Reads mapping to unannotated regions are counted as "no_feature," with high values potentially indicating novel transcripts or the need for annotation updates.

The following diagram illustrates the complete categorization of reads in a typical RNA-seq experiment:

Diagram 2: Complete categorization of RNA-seq reads showing all possible assignment outcomes during the counting process.

The Scientist's Toolkit: Essential Research Reagents

Table 3: Key Reagents and Computational Tools for RNA-seq Read Counting

Tool/Resource	Function	Usage Notes
STAR aligner	Spliced alignment of RNA-seq reads	Use --quantMode GeneCounts for integrated counting [4]
HTSeq-count	Post-alignment read counting	Default settings match STAR's counting logic [15]
SAMtools	BAM file manipulation and sorting	Name sorting critical for htseq-count accuracy [18]
Gencode/Ensembl GTF	Gene structure annotations	Use consistent versions across alignment and counting [5]
Reference genome	Genomic sequence for alignment	Must match annotation build (e.g., GRCh38, mm10)
Quality control tools	FastQC, MultiQC	Assess read quality before counting

Troubleshooting Common Issues

Discrepancy Resolution

When facing significant differences between STAR and htseq-count results, researchers should systematically verify:

• File sorting: Ensure BAM files are name-sorted for htseq-count, as position-sorted files can lead to incorrect pairing and undercounting [18].

• Strandedness settings: Confirm that the strandedness parameter matches the experimental library preparation method. The default --stranded=yes in htseq-count will discard half the reads if data is actually unstranded [15] [16].

• Annotation consistency: Verify that identical GTF files are used during STAR genome generation and htseq-count execution.

• Alignment quality filtering: Check that minimum alignment quality thresholds are consistent (STAR's default filtering versus htseq-count's -a parameter).

Interpretation of Special Counters

STAR's ReadsPerGene.out.tab file includes special counters that provide valuable diagnostic information:

• N_unmapped: Reads that could not be aligned to the reference
• N_multimapping: Reads aligned to multiple genomic locations
• N_noFeature: Reads aligned to regions not annotated in the GTF file
• N_ambiguous: Reads overlapping multiple genes

High values in N_noFeature may indicate the presence of novel transcripts, rRNA contamination, or incomplete annotations. Elevated N_ambiguous counts are common in genomes with high gene density or overlapping transcription units.

STAR's implementation of the 1nt overlap rule for read assignment provides a robust, efficient method for gene-level quantification that closely mirrors the established htseq-count approach. The integrated counting workflow reduces computational overhead and potential sources of error from intermediate file handling. For most RNA-seq applications, STAR's --quantMode GeneCounts offers production-ready counting that maintains accuracy while streamlining processing. Understanding the nuances of how both tools handle edge cases—particularly overlapping genes, multi-mapping reads, and strand-specificity—enables researchers to make informed decisions about quantification methods and accurately interpret resulting count data for downstream expression analysis.

The Critical Role of GTF Annotation File in Accurate Read Counting

In RNA sequencing (RNA-seq) analysis, the General Transfer Format (GTF) file serves as the fundamental blueprint that defines the genomic coordinates of all annotated features, including genes, exons, transcripts, and other functional elements [21]. This file plays an indispensable role in the process of read counting, where sequenced RNA fragments are assigned to specific genomic features to quantify gene expression levels. The accuracy and completeness of the GTF file directly determine the reliability of downstream analyses, including differential expression testing and biological interpretation.

The GTF file's structure is meticulously defined, consisting of nine tab-separated columns that provide essential information about each feature [21]:

• seqname: Chromosome or scaffold name
• source: Origin of the annotation
• feature: Type of feature (e.g., gene, exon, CDS)
• start and end: Genomic coordinates
• score: Confidence metric
• strand: Orientation (+ or -)
• frame: Coding frame for CDS features
• attribute: Semicolon-separated metadata providing identifiers and relationships

Within the context of STAR quantMode GeneCounts analysis, the GTF file enables the aligner to classify reads according to their genomic overlaps, generating counts that form the basis for subsequent expression analysis [4] [5].

GTF File Format Specifications and Requirements

Critical Formatting Specifications

For successful compatibility with read counting tools such as STAR, GTF files must adhere to specific formatting standards [22] [21]. The attribute column (column 9) is particularly crucial as it contains the identifiers that establish relationships between features. The gene_id and transcript_id attributes are mandatory for proper gene-level counting, with the format requiring strict use of semicolon separations and quoted values [22].

Table 1: Essential GTF Attribute Requirements for Read Counting

Attribute	Required For	Format Example	Purpose in Read Counting
gene_id	Gene, transcript, exon	gene_id "ENSG00000123456";	Primary key for aggregating counts to genes
transcript_id	Transcript, exon	transcript_id "ENST00000567890";	Links exons to transcripts
gene_name	Gene	gene_name "TP53";	Human-readable gene identifier
gene_biotype	Gene	gene_biotype "protein_coding";	Classifies gene type for filtering

Tools like STAR and featureCounts utilize these attributes to resolve hierarchical relationships—exons to transcripts, transcripts to genes—enabling accurate assignment of reads to the appropriate feature level [4]. The feature column (column 3) must also use standardized terms such as "gene", "exon", "CDS", "fiveprimeUTR", and "threeprimeUTR" to ensure proper interpretation by counting algorithms [22].

Common GTF File Issues and Solutions

Several common issues can compromise read counting accuracy. Seqname mismatches occur when chromosome naming conventions in the GTF file don't match those in the reference genome FASTA file (e.g., "chr1" versus "1"), causing STAR to fail in recognizing features [23]. The solution involves consistent naming across all reference files.

Missing mandatory attributes like gene_id will cause featureCounts to fail, as it relies on these identifiers to aggregate counts [24]. Validation tools can identify and rectify these missing elements before analysis. Incomplete gene models lacking UTR annotations or containing erroneous boundaries can lead to miscalculation of gene spans and misassignment of reads that map to these regions [22].

Integration of GTF Files with STAR quantMode GeneCounts

STAR Alignment and Counting Workflow

The STAR (Spliced Transcripts Alignment to a Reference) aligner incorporates read counting directly into its alignment process through the quantMode GeneCounts parameter [5] [4]. When this option is specified, STAR utilizes the provided GTF annotation during mapping to assign reads to genomic features simultaneously with alignment.

Table 2: STAR quantMode GeneCounts Output Interpretation

Output Column	Strandedness	Description	HTSeq-count Equivalent
Column 1	N/A	Gene identifiers from GTF file	N/A
Column 2	Unstranded	Reads counted regardless of strand	-s no
Column 3	Forward	1st read strand aligns with RNA	-s yes
Column 4	Reverse	2nd read strand aligns with RNA	-s reverse

The selection of the appropriate column from the STAR output is critical and depends entirely on the library preparation protocol used in the experiment [4] [5]. For unstranded protocols, column 2 should be used, while stranded protocols require selection of either column 3 or 4 based on the specific chemistry employed. Misidentification of the strandedness parameter represents one of the most common sources of error in read counting analysis.

Experimental Protocol: Implementing STAR with GTF-Based Counting

Genome Index Generation The initial step requires building a STAR genome index incorporating the GTF annotation file [5]:

The sjdbOverhang parameter should be set to the maximum read length minus 1, which defines the genomic context around splice junctions extracted from the GTF file [5].

Read Alignment and Counting Following index generation, alignment with simultaneous read counting is performed [5] [25]:

This process generates the *ReadsPerGene.out.tab file containing the raw counts for each gene annotated in the GTF file [4]. The counting algorithm follows rules similar to htseq-count: a read is counted if it overlaps (by at least 1 nucleotide) one and only one gene, with both ends of paired-end reads checked for overlaps [4].

Advanced Applications and Methodological Considerations

Specialized Read Counting Scenarios

Beyond standard gene-level counting, GTF files can be engineered to enable specialized counting scenarios. Transposable elements and repetitive regions present particular challenges due to multimapping reads [24]. One solution involves creating aggregated features in the GTF where all elements of a repetitive family are grouped under a single identifier, allowing for family-level expression quantification despite the impossibility of distinguishing individual genomic loci.

For intronic and non-coding RNA quantification, the GTF file must contain appropriate feature types. While standard annotations focus on protein-coding exons, comprehensive analysis requires inclusion of features such as "ncRNA", "snoRNA", "miRNA", and explicit intronic regions [22]. These can be added to standard GTF files through computational prediction or integration with specialized databases.

The SAF (Site Annotation Format) provides an alternative to GTF for featureCounts when analyzing custom genomic intervals not present in standard annotations [24]. This format requires only five columns (GeneID, Chr, Start, End, Strand) without the hierarchical structure of GTF, making it suitable for counting reads in user-defined regions such as chromatin accessibility peaks or custom genomic bins.

Quality Control and Troubleshooting

Rigorous quality assessment of GTF-based read counting involves examining the count distribution across feature types and the proportion of reads assigned to different categories. STAR's ReadsPerGene.out.tab file includes summary statistics at the beginning that report key metrics [4] [5]:

• N_unmapped: Reads failing to align to the genome
• N_multimapping: Reads aligning to multiple genomic locations
• N_noFeature: Reads aligning to genomic regions not annotated in GTF
• N_ambiguous: Reads overlapping multiple features in the GTF

High proportions of N_noFeature can indicate poor annotation quality, sample contamination, or the use of an inappropriate GTF file for the experimental context. Elevated N_ambiguous counts suggest overlapping feature definitions in the GTF file that require resolution.

Figure 1: GTF-driven RNA-seq analysis workflow with STAR quantMode, showing critical decision points for accurate read counting.

Research Reagent Solutions

Table 3: Essential Tools for GTF-Based Read Counting Analysis

Tool/Resource	Function	Application Context
STAR	Spliced alignment with integrated counting	Primary alignment and quantification [5]
featureCounts	Read summarization for genomic features	Alternative counting engine [24]
ENSEMBL GTF	Curated annotation resource	Standardized gene models [21]
GENCODE GTF	Comprehensive annotation	Human and mouse studies [5]
AGAT Toolkit	GTF/GFF manipulation and conversion	Format troubleshooting and conversion [24]
RStudio	Downstream statistical analysis	Differential expression [25]
SAMtools	BAM file processing	Alignment file manipulation [5]
FASTQC	Read quality control	Pre-alignment quality assessment [25]

From Mapping to Count Table: Your Step-by-Step Workflow for Differential Expression

Essential STAR Command Parameters for Activating quantMode GeneCounts

This application note provides a comprehensive technical overview of the STAR RNA-seq aligner's quantMode GeneCounts feature, which enables simultaneous read alignment and gene-level quantification. We detail the essential command parameters, output interpretation, and experimental protocols required for successful implementation in pharmaceutical and clinical research settings. The GeneCounts function provides a computationally efficient alternative to separate alignment and quantification steps, generating count data directly compatible with downstream differential expression analysis tools like DESeq2. This guide covers critical considerations for stranded versus unstranded library protocols, parameter optimization, and quality control metrics to ensure data integrity for drug discovery pipelines.

The quantMode GeneCounts option, introduced in STAR release 2.4.2a, integrates read counting directly into the alignment process, providing a streamlined workflow for RNA-seq expression analysis [26]. This feature counts reads per gene during mapping according to the same "union" mode used by htseq-count with default parameters, where a read is counted if it overlaps one and only one gene by at least one nucleotide [26]. For paired-end experiments, both ends are checked for overlaps with gene features [26]. This integrated approach eliminates the need for separate counting steps while maintaining compatibility with established bioinformatics pipelines for drug target identification and biomarker discovery.

A key advantage of GeneCounts is its production of counts that coincide with those generated by htseq-count, providing researchers with a validated counting methodology [26]. The counting is performed only on uniquely mapped reads, excluding multi-mappers from the final counts unless specifically included through alternative parameter configurations [26]. The implementation specifically counts reads that overlap exonic regions of genes, meaning reads aligning wholly within intronic space are not counted toward gene expression [26]. This exonic-focused counting provides a more accurate representation of mature transcript abundance, which is particularly valuable in pharmaceutical research investigating expression changes in response to therapeutic interventions.

Essential Command Parameters

Core Parameter Specification

Implementing quantMode GeneCounts requires specific command-line parameters that control both alignment behavior and counting functionality. The most fundamental parameters for activating and configuring gene counting are detailed in the table below.

Table 1: Essential STAR Parameters for GeneCounts Implementation

Parameter	Value	Function	Required
--quantMode	GeneCounts	Activates read counting per gene during alignment	Yes
--sjdbGTFfile	/path/to/annotations.gtf	Provides gene annotations for counting	Yes
--runThreadN	<number>	Specifies processor cores for parallelization	Recommended
--readFilesCommand	zcat or gunzip -c	Enables reading of compressed input files	Conditional
--outFileNamePrefix	<output_prefix>	Defines naming convention for output files	Recommended
--outSAMtype	BAM SortedByCoordinate	Specifies sorted BAM output for downstream analysis	Optional

A typical STAR command incorporating these essential parameters appears as follows:

Advanced Parameter Configurations

For enhanced junction discovery and counting accuracy, particularly in complex transcriptomes relevant to disease states, additional parameters can be incorporated:

The --sjdbOverhang parameter should be set to the read length minus 1, which is critical for accurate splice junction detection [5]. For example, with 150bp reads, this value should be 149. The --twopassMode Basic enables novel junction discovery in the first pass that is then incorporated into the second alignment pass, improving mapping accuracy in transcriptomic studies investigating novel splice variants in disease pathways [9].

Output Interpretation and Data Structure

ReadsPerGene.out.tab Format

STAR generates gene count data in the ReadsPerGene.out.tab file, which contains four columns corresponding to different strandedness options [4] [5] [9]. The file structure is consistent across experiments and facilitates proper selection of count data based on library preparation methodology.

Table 2: ReadsPerGene.out.tab File Structure

Column	Content	HTSeq-count Equivalent	Usage
1	Gene identifier	N/A	Gene IDs from the provided GTF file
2	Unstranded counts	--stranded no	For non-strand-specific protocols
3	Counts for 1st read strand	--stranded yes	For strand-specific protocols where the first read aligns with RNA
4	Counts for 2nd read strand	--stranded reverse	For strand-specific protocols where the second read aligns with RNA

The file begins with several summary lines prefixed with "N_" that provide mapping statistics [4] [5]:

• N_unmapped: Total number of unmapped reads
• N_multimapping: Reads that aligned to multiple genomic locations
• N_noFeature: Reads that did not overlap any gene feature
• N_ambiguous: Reads that overlapped multiple genes and could not be uniquely assigned

Strandedness Selection Protocol

Proper selection of the appropriate column from the ReadsPerGene.out.tab file is critical for accurate gene expression quantification. The following decision pathway and experimental validation method ensures correct column selection:

Diagram 1: Strandedness selection workflow

To experimentally determine the correct column for analysis, researchers should calculate the total counts across all genes for columns 3 and 4 (excluding the N_ summary lines) [26]. For example:

If one column contains significantly higher counts (typically >80-90% of the total), this indicates a stranded protocol, and the larger count column should be selected [5] [26]. For Illumina TruSeq Stranded protocols, column 4 (corresponding to htseq-count's --stranded reverse) is typically appropriate [26]. When both columns have similar totals, the protocol is likely unstranded, and column 2 should be used [26].

Experimental Protocols

Complete RNA-seq Analysis Workflow

The following protocol describes an end-to-end workflow for gene expression analysis using STAR's quantMode GeneCounts, optimized for pharmaceutical research applications:

Pre-mapping Steps:

• Quality Control: Perform QC on raw FASTQ files using FastQC
• Adapter Trimming: Remove adapter sequences using Trimmomatic or Cutadapt
• Genome Index Preparation: Generate STAR genome indices with gene annotations

STAR Alignment with GeneCounts:

Post-mapping Analysis:

• Count Matrix Compilation: Extract the appropriate column from ReadsPerGene.out.tab across all samples
• Downstream Analysis: Import count matrix into DESeq2 or edgeR for differential expression analysis
• Quality Assessment: Verify mapping statistics from Log.final.out files

Genome Index Generation Protocol

Proper genome index generation is a prerequisite for successful alignment and counting. The following protocol generates a comprehensive STAR index:

The --sjdbOverhang parameter should be set to the maximum read length minus 1, which for modern sequencing platforms is typically 100-150 [5] [9]. This parameter determines the length of the genomic sequence around annotated junctions that is used for constructing the splice junction database, directly impacting mapping accuracy across splice boundaries.

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Resources

Resource	Specification	Function in Workflow
STAR Genome Index	Species-specific with GTF annotations	Reference for alignment and gene feature identification
Reference Annotations	GTF format from GENCODE, Ensembl, or RefSeq	Defines gene models for read counting
High-Performance Computing	16+ CPU cores, 64+ GB RAM	Enables parallel processing of large RNA-seq datasets
Quality Control Tools	FastQC, MultiQC	Assesses read quality and mapping statistics
Differential Expression Tools	DESeq2, edgeR, limma-voom	Identifies statistically significant expression changes
Stranded RNA-seq Library Kit	Illumina TruSeq Stranded, NEBNext Ultra II	Preserves strand information for transcript origin

Troubleshooting and Quality Assurance

Common Issues and Solutions

Discrepancies between strandedness columns: When the sum of columns 3 and 4 does not equal column 2, this typically results from reads mapping to genes that overlap on opposite strands [26]. In such cases, the unstranded count (column 2) considers these reads ambiguous and excludes them, while stranded counting assigns them to separate genes based on strand orientation [26]. This is expected biological behavior rather than a technical artifact.

Low feature assignment rates: If the N_noFeature counts are excessively high (>50% of total reads), verify that the annotation file used for counting matches the reference genome and contains comprehensive gene models. For non-model organisms, consider using --sjdbGTFfeatureExon and --sjdbGTFtagExonParentGene parameters to customize feature types and parent gene identifiers [27].

Count discrepancies with htseq-count: While STAR GeneCounts is designed to replicate htseq-count output, minor differences may occur due to variations in how alignment information is parsed [6] [26]. For absolute consistency with htseq-count, researchers can use the Aligned.sortedByCoord.out.bam file as direct input to htseq-count, though this eliminates the efficiency advantage of integrated counting.

Quality Control Metrics

Successful GeneCounts implementation should yield the following quality metrics:

• Uniquely mapped reads: >70% of total reads [9]
• Gene feature assignment: >60% of uniquely mapped reads
• Strand specificity: >90% of reads in one column for stranded protocols [5]
• Multimapping rate: <10% of total reads (unless using specific multi-mapper counting)

The Log.final.out file provides comprehensive mapping statistics for assessing these metrics, including unique mapping rates, splice junction detection, and insertion/deletion patterns that inform overall data quality for downstream analysis in drug discovery pipelines.

How to Determine Your Library's Strandedness from the Count Data

In RNA sequencing (RNA-Seq) analysis, accurately interpreting gene expression data requires knowing whether your library preparation was stranded or unstranded. Stranded RNA-Seq preserves the original orientation of transcripts, enabling precise detection of antisense transcription and enhancing transcript annotation [28]. The STAR aligner, through its --quantMode GeneCounts function, provides a built-in solution to determine this crucial parameter from your data itself [10] [4]. This application note details a robust protocol for analyzing the STAR output file ReadsPerGene.out.tab to definitively determine your library's strandedness, a critical step for ensuring the accuracy of downstream differential expression analysis.

When run with the --quantMode GeneCounts option, STAR quantifies reads per gene and generates a *ReadsPerGene.out.tab file. This file contains four columns corresponding to different strandedness assumptions [4]:

• Column 1: Gene identifier
• Column 2: Counts for unstranded RNA-seq
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count -s yes)
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count -s reverse)

The core principle for determining strandedness is that the correct strand-specific column (3 or 4) will yield the highest number of reads confidently assigned to genes and the fewest reads assigned to no feature [10].

Materials and Equipment

Research Reagent Solutions and Essential Materials

Table 1: Key reagents, tools, and software required for the analysis.

Item Name	Function/Description
STAR Aligner	Performs alignment and read counting with --quantMode GeneCounts [4].
ReadsPerGene.out.tab file	STAR output containing counts for different strandedness modes.
Computing Environment	Linux/server environment with adequate resources to run STAR.
Text processing tools	awk, grep, or scripting language (e.g., Python, R) for data extraction.

Step-by-Step Computational Protocol

Step 1: Generate the Count Data Using STAR

Align your RNA-Seq reads with STAR, ensuring you include the --quantMode GeneCounts option and provide a reference annotation file via --sjdbGTFfile during genome generation or mapping [4]. A typical command structure is:

Step 2: Inspect the Output File

Examine the header of your *ReadsPerGene.out.tab file. The first lines report overarching statistics before the gene-specific counts begin [10] [4]:

These lines are followed by rows of gene IDs and their corresponding counts in columns 2, 3, and 4.

To objectively determine strandedness, calculate the total number of reads assigned to genomic features (N_genic) for each strandedness column (2, 3, and 4). N_genic is calculated as [10]: N_genic = Total_Reads - (N_unmapped + N_multimapping + N_noFeature + N_ambiguous)

Extract the values for N_unmapped, N_multimapping, N_noFeature, and N_ambiguous from the file header. The Total_Reads can be derived from the alignment report in the Log.final.out file or approximated as the sum of all assigned and unassigned reads.

Step 4: Interpret the Results and Determine Strandedness

Compare the calculated N_genic values for columns 2, 3, and 4, and the N_noFeature values.

Table 2: Interpretation guide for strandedness determination.

Observation	Interpretation	Recommended Column
Column 3 has the highest N_genic and lowest N_noFeature.	Library is stranded, and the "1st read strand" corresponds to the RNA strand.	Column 3
Column 4 has the highest N_genic and lowest N_noFeature.	Library is stranded, and the "2nd read strand" corresponds to the RNA strand (reverse).	Column 4
Column 2 has the highest N_genic, with Columns 3 & 4 being roughly equal and having high N_noFeature.	Library is unstranded.	Column 2

The diagnostic workflow for this decision process is illustrated below.

Troubleshooting and Data Validation

• Unexpectedly High N_noFeature: If the N_noFeature counts are high across all columns, verify that the chromosome names in your GTF annotation file match those in the reference genome used for alignment [10].
• Ambiguous Results: In rare cases, the differences between columns might be subtle. As the STAR developer notes, the column with the highest N_genic value points to the correct strand, and the ratio of N_genic between the top two columns provides an upper bound for the strand error in the library [10].
• Understanding Count Discrepancies: Do not expect the sum of columns 3 and 4 to equal column 2. This is because a read that overlaps a gene on one strand and also overlaps another gene on the opposite strand is considered "ambiguous" for unstranded counting (column 2) but can be assigned unambiguously to one gene when using the correct stranded column (3 or 4) [4].

Determining the strandedness of an RNA-Seq library is a foundational step in data analysis. The protocol outlined here, leveraging the built-in counting features of the STAR aligner, provides a reliable and computationally efficient method to deduce this information directly from the data. Correctly selecting the count column based on this diagnostic analysis ensures that downstream applications, such as differential expression testing, accurately reflect the true biology of the transcriptome.

The --quantMode GeneCounts feature in the STAR aligner is a powerful tool that integrates read alignment with transcript quantification, streamlining the RNA-seq analysis workflow. When this option is activated, STAR generates a *ReadsPerGene.out.tab file containing read counts for each gene based on the provided genome annotation [4] [26]. This file contains four columns that correspond to different strandedness options, allowing researchers to select the appropriate counts based on their library preparation protocol [9] [5]. The selection of the correct column is not merely a procedural step but is critical for generating accurate gene expression data, as an incorrect choice can lead to misquantification, particularly for genes that overlap with others on the opposite strand [29].

The fundamental principle behind this counting method is that a read is assigned to a gene if it overlaps (by at least 1 nucleotide) with one and only one gene [26]. For paired-end reads, both ends are checked for overlaps [26]. This counting strategy coincides with the "union" mode of the popular htseq-count tool, providing consistency with established bioinformatics workflows [26].

Understanding the Columns: Structure and Interpretation

The ReadsPerGene.out.tab file is structured with a consistent format that begins with summary lines followed by gene-specific counts. The first four lines provide crucial summary statistics [4]:

• N_unmapped: Total number of unmapped reads.
• N_multimapping: Number of reads mapped to multiple genomic locations.
• N_noFeature: Reads that do not overlap any annotated gene.
• N_ambiguous: Reads that overlap multiple genes, making unambiguous assignment impossible.

Following these summary lines, the file lists counts for individual genes in four columns with specific strandedness interpretations [9] [5]:

Table 1: Description of Columns in ReadsPerGene.out.tab

Column	Strandedness	htseq-count Equivalent	Recommended Protocol
1	Gene ID	Not applicable	Not applicable
2	Unstranded	--stranded no	Non-strand-specific protocols
3	1st read strand aligned with RNA	--stranded yes	Stranded protocols where the first read is sense
4	2nd read strand aligned with RNA	--stranded reverse	Stranded protocols where the first read is antisense

It is important to note that only uniquely mapped reads are counted in these columns [26]. Multi-mapping reads are excluded from gene counts, though they are reported in the N_multimapping summary line [26].

Protocol-Specific Column Selection Guide

Decision Workflow for Column Selection

The following diagram illustrates the systematic process for determining the correct column for your RNA-seq data:

Selection Criteria Based on Library Type

Unstranded RNA-seq (Column 2)

For non-strand-specific protocols, use the counts in column 2 [5]. In these protocols, cDNA is synthesized without preserving strand information, resulting in reads deriving from either strand without distinction [29]. The unstranded count assigns reads to genes regardless of their original transcriptional strand [4]. This approach is sufficient for many applications but can lead to ambiguous assignments for genes that overlap on opposite strands [29].

Stranded RNA-seq (Columns 3 or 4)

Stranded protocols preserve the information about which DNA strand was originally transcribed, allowing for more accurate assignment of reads to their correct transcriptional unit [29]. The choice between column 3 and 4 depends on the specific stranded protocol used:

• Column 3 (--stranded yes): Use when the first read of the pair is aligned to the same strand as the RNA transcript (sense orientation) [9] [5].
• Column 4 (--stranded reverse): Use when the first read of the pair is aligned to the opposite strand of the RNA transcript (antisense orientation) [9] [5]. This is the case for the Illumina TruSeq Stranded protocol [26].

Empirical Verification Methods

When protocol information is unavailable, empirical methods can determine the correct column:

• Compare Total Counts: Sum the counts across all genes in columns 3 and 4 (excluding the summary lines). The column with the significantly higher total count corresponds to the correct strand orientation [26].
• Inspect Known Genes: Examine specific genes with known strand orientation. For a correctly chosen stranded column, counts for genes on the opposite strand should be minimal or represent genuine antisense transcription [5].
• Check N_noFeature: The N_noFeature values differ between columns 3 and 4. The correct stranded column typically has a lower N_noFeature count [26].

Important Considerations and Common Issues

Mathematical Relationships Between Columns

A common point of confusion arises when column 2 is not equal to the sum of columns 3 and 4 [4] [26]. This occurs due to ambiguous reads that overlap genes on opposite strands [26]. In unstranded counting (column 2), these reads are marked as ambiguous and excluded from individual gene counts. In stranded counting (columns 3 and 4), the same reads can be assigned to different genes based on their strand information [26]. Therefore, the relationship is typically: Column 2 ≤ Sum of Columns 3 and 4 [26].

Advantages of Stranded RNA-seq

Stranded protocols provide significant advantages for transcriptome analysis [29]:

• Resolution of Overlapping Genes: Approximately 11,000 genes (19%) in GENCODE annotations overlap with genes on the opposite strand [29]. Stranded RNA-seq unambiguously assigns reads in these regions.
• Antisense Transcription Detection: Enables discovery and quantification of natural antisense transcripts [29].
• More Accurate Quantification: Reduces misassignment of reads to wrong genes, improving differential expression analysis accuracy [29].

Table 2: Impact of Strandedness on Read Assignment

Feature	Unstranded RNA-seq	Stranded RNA-seq
Ambiguous Read Rate	~6.1%	~2.9%
Opposite Strand Overlap	Cannot resolve	Accurately resolves
Antisense Detection	Limited	Enabled
Differential Expression	Potentially inaccurate for overlapping genes	More accurate

Differences with htseq-count

While STAR's GeneCounts is designed to match htseq-count output, slight differences may occur due to [6]:

• Handling of alignment boundaries and overlap definitions
• Variations in multi-mapping read filtering
• Different implementations of the "union" mode

For maximum consistency with htseq-count, some protocols recommend running htseq-count separately on the BAM output, though STAR's integrated counting is generally sufficient for most applications [4].

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Reagents and Kits for RNA-seq Library Preparation

Reagent/Kit	Function	Strandedness	Key Features
TruSeq Stranded mRNA Kit	PolyA-selected mRNA library prep	Stranded (use Column 4)	dUTP second-strand marking; high compatibility with FFPE samples
dUTP Second-Strand Marking Method	cDNA synthesis with strand marking	Stranded	Incorporates dUTP in second strand; enzymatic degradation before PCR
Oligo(dT) Primers	mRNA enrichment	Varies	Selects polyadenylated transcripts; used in most mRNA-seq protocols
Ribosomal RNA Depletion Kits	rRNA removal	Varies	Preserves non-polyadenylated transcripts; broader transcriptome coverage
RNase H	RNA degradation	Protocol-dependent	Degrades RNA strand in RNA-DNA hybrids; used in cDNA synthesis

Selecting the correct column from STAR's ReadsPerGene.out.tab output is fundamental for accurate RNA-seq quantification. Researchers must understand their library preparation protocol and apply the systematic selection process outlined in this guide. Stranded RNA-seq protocols generally provide superior accuracy, particularly for complex genomes with abundant overlapping transcripts, and represent the recommended approach for future RNA-seq studies [29]. By following these guidelines and understanding the underlying principles, researchers can ensure the generation of reliable gene expression data for downstream analysis.

Integrating STAR GeneCounts into a DESeq2 or edgeR Differential Expression Pipeline

This application note provides a comprehensive protocol for integrating STAR RNA-seq aligner output, specifically the quantMode GeneCounts results, into established differential expression (DE) analysis pipelines using DESeq2 or edgeR. The STAR aligner generates read counts per gene that serve as direct input for count-based statistical models in DE analysis. We detail the complete workflow from experimental design to result interpretation, focusing on the critical steps of data import, quality control, statistical modeling, and visualization. This protocol emphasizes proper handling of STAR's multi-column count output based on library strandedness and integrates robust analysis practices to ensure biologically meaningful results for researchers and drug development professionals.

RNA sequencing (RNA-seq) has become the primary method for transcriptome analysis, enabling detailed investigation of gene expression patterns across biological conditions [30]. A fundamental goal in many RNA-seq experiments is identifying differentially expressed genes between experimental groups, which provides crucial insights into molecular mechanisms and potential therapeutic targets [31].

The STAR (Spliced Transcripts Alignment to a Reference) aligner provides integrated read counting through its --quantMode GeneCounts option, generating counts that coincide with those produced by htseq-count with default parameters [4] [5]. These counts serve as ideal input for count-based statistical frameworks like DESeq2 and edgeR, which rely on raw, unnormalized count data to accurately model biological variation and test for differential expression [32] [33].

This protocol bridges the gap between alignment and differential expression analysis, providing a standardized approach for leveraging STAR's counting capabilities within robust statistical frameworks.

Background

STAR quantMode GeneCounts Output

When run with --quantMode GeneCounts, STAR outputs read counts per gene into a ReadsPerGene.out.tab file with four columns corresponding to different strandedness options [4] [5]:

• Column 1: Gene identifier
• Column 2: Counts for unstranded RNA-seq
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count option -s yes)
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count option -s reverse)

The appropriate column selection is critical and depends on the library preparation protocol used. For unstranded protocols, column 2 should be used. For stranded protocols, either column 3 or 4 is selected based on the specific strand orientation of the library preparation method [4].

Statistical Foundations of Differential Expression Analysis

DESeq2 and edgeR employ similar statistical approaches for differential expression analysis, both using the negative binomial distribution to model count data while accounting for biological variability [33]. These methods internally account for library size differences, making it essential to provide raw count data rather than pre-normalized values [32]. The analysis workflow typically includes normalization, dispersion estimation, statistical testing, and multiple testing correction.

Materials and Methods

Research Reagent Solutions

Table 1: Essential reagents and computational tools for RNA-seq differential expression analysis

Category	Item	Function
Alignment & Counting	STAR aligner	Spliced alignment of RNA-seq reads to reference genome with integrated counting [5]
Reference Files	Genome sequence (FASTA)	Reference genome for read alignment [5]
	Gene annotation (GTF)	Gene models for read assignment and counting [5]
Differential Expression	DESeq2	Statistical analysis of differential expression using negative binomial models [34] [33]
	edgeR	Statistical analysis of differential expression with empirical Bayes methods [35]
Quality Control	FastQC	Quality assessment of raw sequencing reads [31]
	Trimmomatic/fastp	Read trimming and adapter removal [30]
Programming Environment	R	Statistical computing environment for analysis [34] [35]
	Bioconductor	Repository for bioinformatics packages including DESeq2 and edgeR [32]

The following diagram illustrates the complete workflow from raw sequencing data to differential expression results:

Strandedness Determination and Column Selection

Proper selection of the count column from STAR's output requires understanding the library preparation protocol. The following decision diagram guides this critical step:

Table 2: Guide to selecting the appropriate column from STAR ReadsPerGene.out.tab

Library Type	Protocol Examples	STAR Column	htseq-count Equivalent
Unstranded	Standard Illumina, non-stranded	Column 2	Default parameters
Stranded (forward)	dUTP, NSR, NNSR	Column 4	-s reverse
Stranded (reverse)	ScriptSeq, Ligation	Column 3	-s yes

To empirically verify strandedness, researchers can examine counts from known strand-specific genes or use the following approach to compare column totals [5]:

Protocol: DESeq2 Analysis with STAR Counts

Data Import and Preprocessing

The first step involves importing STAR count data into R and creating a DESeq2 dataset object:

DESeq2 Differential Expression Analysis

With the count data properly formatted, proceed with the DESeq2 analysis pipeline:

Quality Control and Visualization

Comprehensive quality assessment is critical for validating analysis results:

Protocol: edgeR Analysis with STAR Counts

Data Import and Normalization

The edgeR pipeline shares similarities with DESeq2 but employs different normalization strategies:

Differential Expression Testing

edgeR provides multiple testing approaches; here we demonstrate the quasi-likelihood method:

Results Interpretation and Downstream Analysis

Both DESeq2 and edgeR generate comprehensive result tables containing:

Table 3: Key columns in differential expression results

Column	Description	Interpretation
baseMean	Average normalized count across all samples	Overall expression level
log2FoldChange	Log2 fold change between conditions	Effect size (biological significance)
pvalue	Nominal p-value for statistical test	Unadjusted significance
padj	Multiple testing-adjusted p-value (FDR)	Statistical significance after correction

Visualization Strategies

Effective visualization enhances interpretation of differential expression results:

Critical Analysis Parameters

Table 4: Key parameters for differential expression analysis

Parameter	DESeq2 Function	edgeR Function	Considerations
Fold Change Threshold	lfcThreshold in results()	logFC in topTags()	Biological relevance vs. statistical power
FDR Cutoff	alpha in results()	Default 0.05 in topTags()	Balance false discoveries vs. missed findings
Independent Filtering	Automatic in results()	filterByExpr()	Reduces multiple testing burden

Troubleshooting and Quality Assessment

Common Issues and Solutions

• Low alignment rates: Verify genome index compatibility with read length and quality-trim reads [30]
• Sample clustering by batch rather than condition: Include batch effects in design formula (~ batch + condition)
• Inflation of low p-values: Check dispersion estimates and consider additional covariates
• Unexpected strandedness results: Validate with positive control genes known to be strand-specific

Quality Metrics

Table 5: Expected ranges for key quality metrics

Metric	Target Range	Potential Issues
STAR Alignment Rate	>70% of reads	Poor RNA quality, adapter contamination
Exonic Mapping Rate	>60% of aligned reads	DNA contamination, ribosomal RNA
DESeq2 Mean-SD Trend	Smooth decreasing curve	Outliers, inadequate replication
PCA Sample Clustering	Grouping by experimental condition	Batch effects, sample mislabeling

This protocol provides a standardized approach for integrating STAR's quantMode GeneCounts output into robust differential expression analysis pipelines using DESeq2 or edgeR. The key to success lies in proper experimental design, appropriate strandedness determination, careful quality control, and informed interpretation of statistical outputs. By following this comprehensive workflow, researchers can confidently identify biologically meaningful expression changes that advance understanding of molecular mechanisms and support drug development efforts.

The integration of STAR alignment and counting with established statistical methods represents an efficient and reproducible workflow for RNA-seq differential expression analysis, providing a solid foundation for transcriptomic investigations across diverse biological contexts.

Solving Common quantMode GeneCounts Issues: Discrepancies, Low Counts, and Customization

A nuanced understanding of the STAR --quantMode GeneCounts output is crucial for accurate gene expression quantification. A common point of confusion arises when the sum of counts in columns 3 and 4 does not equal the value in column 2. This discrepancy is not an error but a direct consequence of how reads overlapping genes on opposite strands are classified under different strandedness modes. This application note demystifies the structure of the ReadsPerGene.out.tab file, provides a definitive protocol for selecting the correct count column, and outlines a robust framework for troubleshooting and quality control to ensure the integrity of downstream differential expression analysis.

Decoding the ReadsPerGene.out.tab File

The --quantMode GeneCounts option in STAR directs the aligner to output a *ReadsPerGene.out.tab file alongside the standard alignment files. This file contains raw read counts per gene, which serve as the primary input for differential expression tools like DESeq2 and edgeR [36].

The file is structured into four columns, each representing counts based on a different assumption about the library's strandedness [4] [5]:

• Column 1: Gene identifier.
• Column 2: Counts for unstranded RNA-seq. A read is counted toward a gene regardless of its mapping strand.
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count -s yes).
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count -s reverse).

In addition to gene counts, the file includes several summary lines that provide a high-level overview of the mapping outcomes [4]:

• N_unmapped: Reads that could not be aligned to the genome.
• N_multimapping: Reads that aligned to multiple genomic locations. It is critical to note that multi-mapped reads are not counted toward any gene and are entirely relegated to this category [26] [10].
• N_noFeature: Reads that aligned to genomic regions not annotated as a gene in the provided GTF file.
• N_ambiguous: Reads that overlap multiple genes or overlap genes on opposite strands. The treatment of these reads is the primary source of the discrepancy between columns.

The Core Principle: How Strandedness and Read Ambiguity Affect Counting

The central question of why Column 2 ≠ Column 3 + Column 4 is explained by the concept of read ambiguity [26] [37].

The Ambiguity in Unstranded Counting (Column 2)

In unstranded counting, a read that overlaps features on opposite strands cannot be assigned with confidence to a single gene. Consider a read that maps to a genomic location where two genes, GeneA on the + strand and GeneB on the - strand, overlap. From the perspective of the unstranded protocol, this read could have originated from either gene. To avoid double-counting, STAR follows the same logic as htseq-count and classifies this read as "ambiguous" (N_ambiguous). Consequently, it is not counted for either GeneA or GeneB in Column 2 [4] [37].

The Resolution in Stranded Counting (Columns 3 & 4)

Stranded counting leverages the library preparation protocol information. The same read that was ambiguous for Column 2 can now be assigned based on its strand of origin. If the read aligns to the + strand, it will be counted for GeneA in the appropriate stranded column (Column 3 or 4). If it aligns to the - strand, it will be counted for GeneB [4] [37]. Therefore, in stranded columns, reads that were excluded from Column 2 due to strand ambiguity are now incorporated into specific gene counts.

Table 1: How Read Fate Depends on Counting Method and Strandedness

Read Mapping Scenario	Fate in Column 2 (Unstranded)	Fate in Columns 3 & 4 (Stranded)
Read maps uniquely to a gene on one strand.	Counted for the gene.	Counted for the gene in the corresponding strand-specific column.
Read maps to a location where genes on the + and - strands overlap.	Ambiguous: Added to N_ambiguous, not counted for any gene.	Resolved: Counted for the gene on the strand to which the read aligns.
Read does not map to any annotated gene.	Added to N_noFeature.	Added to N_noFeature.

This principle is visualized in the following diagram, which illustrates how a single read mapping to overlapping genes on opposite strands is classified differently.

Figure 1: Resolution of read assignment in unstranded versus stranded counting. A read that is ambiguous and excluded from gene counts in unstranded mode can be uniquely assigned to a gene in the correct stranded column.

Experimental Protocol: Column Selection and Data Preparation for Differential Expression

Selecting the incorrect column for analysis is a critical error that will compromise all downstream results. The following step-by-step protocol ensures the correct column is chosen and the data is properly formatted for differential expression analysis.

Step 1: Determine Your Library's Strandedness

If your library preparation method is known (e.g., Illumina TruSeq Stranded mRNA kit), this information is definitive. For such reverse-stranded libraries, you should use Column 4 [36] [5]. If the protocol is unknown, you can determine it empirically from the ReadsPerGene.out.tab file itself.

Procedure:

• For a given sample, calculate the total number of reads assigned to genes for each stranded column. This can be done by summing all the counts in Column 3 (excluding the N_* lines) and separately summing all the counts in Column 4.
• Alternatively, and more simply, compare the N_noFeature values for Column 3 and Column 4 [10]. The correct strand will have a significantly lower N_noFeature count because more reads will fall within the exonic regions of genes annotated on the correct strand.
• The column with the lower N_noFeature count (and higher gene-assigned count) is the one to use for your analysis [5] [10].

Table 2: Empirical and A Priori Methods for Determining the Correct Count Column

Method	Procedure	Interpretation
A Priori (Known Protocol)	Check experimental records for the library kit used.	- Unstranded: Use Column 2.- Stranded (e.g., Illumina TruSeq Stranded): Use Column 4 [36].
Empirical (N_noFeature Check)	Compare the N_noFeature values in Columns 3 and 4 of your ReadsPerGene.out.tab file.	The column with the lower N_noFeature value is the correct one for your library [10].
Empirical (Total Genic Reads)	Sum all gene counts (excluding N_* lines) in Columns 3 and 4.	The column with the higher total genic count is the correct one [10].

Step 2: Generate the Count Matrix for Differential Expression

Differential expression tools like DESeq2 require a unified count matrix across all samples. The following command-line protocol, adapted from the BiocoreCRG RNA-seq course, demonstrates how to create this matrix from STAR output [36].

Application Note: This example assumes a reverse-stranded library (using Column 4).

The final raw_counts_A549_matrix.txt file will be a tab-delimited matrix where the first column is the gene identifier and each subsequent column contains the raw counts for that gene in each sample, ready for input into DESeq2 or a similar tool.

The overall workflow, from sequencing reads to a count matrix, is summarized below.

Figure 2: A recommended workflow for processing RNA-seq data with STAR and generating a count matrix for differential expression analysis.

The Scientist's Toolkit: Essential Research Reagents and Materials

Successful execution of the RNA-seq analysis pipeline depends on the following key components.

Table 3: Essential Research Reagent Solutions for RNA-seq Analysis with STAR

Item	Function / Role in the Workflow
STAR Aligner	Spliced Transcripts Alignment to a Reference. Performs fast, accurate alignment of RNA-seq reads and can simultaneously generate read counts per gene [4] [5].
Reference Genome (FASTA)	The genomic sequence for the target organism. Provides the reference coordinates for aligning the sequenced reads.
Annotation File (GTF/GFF)	File containing the coordinates and strandedness of all known genes, transcripts, and exons. Critical for both the alignment and the gene counting steps. Must be from the same source and version as the reference genome [26] [5].
High-Performance Computing (HPC) Cluster	STAR is computationally intensive for both genome indexing and read alignment. An HPC environment with sufficient RAM and multiple CPU cores is typically required.
Differential Expression Tool (e.g., DESeq2, edgeR)	R/Bioconductor package that takes the raw count matrix as input and uses statistical models to identify differentially expressed genes between experimental conditions [36].

Troubleshooting and Quality Control

This section addresses common issues and questions that arise when working with STAR gene counts.

Extreme Cases: Column 2 is Zero While Columns 3 and 4 are High

In rare cases, you might observe a gene where the value in Column 2 is 0, but the values in Columns 3 and 4 are very high (e.g., in the thousands) [37]. This is an extreme manifestation of the ambiguity principle. It occurs when two genes are perfectly or near-perfectly overlapping on opposite strands. In the unstranded view (Column 2), all reads mapping to this locus are considered ambiguous. However, in the stranded view (Columns 3 and 4), the reads are neatly partitioned between the two genes based on their strand of origin [37]. The solution is to verify the gene annotation in a genome browser and confidently use the correct stranded column.

Consistency with HTSeq-Count

While STAR's --quantMode GeneCounts is designed to produce counts that "coincide with those produced by htseq-count with default parameters" [26], minor differences can sometimes occur. If perfect parity with an independent htseq-count run is absolutely required, it is advisable to use the htseq-count results. However, for most applications, the counts generated directly by STAR are perfectly valid and more efficient to produce.

The non-equivalence of Column 2 and the sum of Column 3 + Column 4 in STAR's output is a deliberate and logical feature of the counting algorithm, not a software bug or data error. This behavior stems from the differential handling of reads that map to overlapping genes on opposite strands. By understanding the principles of strandedness and read ambiguity, researchers can correctly select the appropriate count column based on their experimental design. The protocols and troubleshooting guide provided herein empower scientists to confidently generate accurate gene count matrices, ensuring a solid foundation for all subsequent differential expression and biomarker discovery efforts.

In the analysis of RNA-sequencing data, accurate read quantification is a critical step for downstream differential expression analysis. A common approach involves using aligners like STAR (Spliced Transcripts Alignment to a Reference), which offers integrated counting via the --quantMode GeneCounts option, ostensibly producing counts that "coincide with those produced by htseq-count with default parameters" [26]. However, researchers frequently encounter substantial discrepancies between the gene counts generated by STAR and those from standalone counting tools such as featureCounts or HTSeq-count [38] [6]. These inconsistencies can arise from subtle differences in parameters, annotation file handling, and library preparation settings. This Application Note details the primary causes of these discrepancies and provides validated protocols to harmonize results, ensuring reliable and reproducible gene expression data for the research and drug development community.

Understanding the Tools and Their Outputs

STAR's Integrated Quantification

The STAR aligner incorporates a two-step process for read alignment: seed searching and clustering, stitching, and scoring [39] [40]. When run with --quantMode GeneCounts, it counts reads per gene concurrently with alignment. The output file ReadsPerGene.out.tab contains four columns [4] [26]:

• Column 1: Gene identifier.
• Column 2: Counts for unstranded RNA-seq.
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count -s yes).
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count -s reverse).

A critical point is that STAR, by default, counts only uniquely mapping reads and considers a read to belong to a gene if it overlaps (by at least 1 nucleotide) one and only one gene, following the "union" mode of HTSeq-count [26] [15]. Reads that are multi-mapping are excluded from gene counts and are reported in the N_multimapping tally [10].

Standalone Counting Tools

• HTSeq-count: This tool provides several overlap resolution modes, with union being the most commonly recommended [41] [15]. It is considered more conservative in its counting, as it may classify a higher proportion of reads as ambiguous if they overlap more than one feature [41]. A crucial setting for paired-end data is the -r parameter, which specifies the sort order of the input BAM file (name or pos); an incorrect setting here is a documented source of count discrepancies [42].
• featureCounts: Part of the Subread package, featureCounts is noted for its high speed (~20x faster than HTSeq-count) and slightly more liberal counting strategy [41]. For paired-end data, it is essential to use the --countReadPairs parameter to count fragments (i.e., read pairs) instead of individual reads; omitting this will result in featureCounts reporting approximately double the counts of STAR and HTSeq-count [38].

Table 1: Key Characteristics of Read Counting Tools

Tool	Primary Function	Counting Strategy	Speed	Key Parameter for PE Data
STAR (--quantMode)	Alignment & Quantification	Counts unique mappers, "union" mode	Fast (integrated)	N/A (inherently counts fragments)
HTSeq-count	Quantification only	Conservative, higher ambiguous rates [41]	Slow	-r pos (for position-sorted BAM) [42]
featureCounts	Quantification only	Liberal, lower ambiguous rates [41]	Very Fast [41]	--countReadPairs [38]

Strandedness Settings

A major source of count discrepancies is the incorrect specification of library strandedness.

• Problem: Selecting the wrong column in STAR's output or using an incorrect strandedness parameter (-s in featureCounts/HTSeq-count) will lead to drastically low counts, as reads are assigned to the wrong strand or excluded [4].
• Diagnosis: Researchers can determine the correct strandedness by comparing the total counts assigned to genes (the N_genic value) for columns 3 and 4 in the STAR output. The column with the higher N_genic value indicates the correct strand [10]. Tools like infer_experiment.py from the RSeQC package can also be used for validation [38].
• Solution: Ensure consistency across your workflow. The table below maps common stranded protocol types to the correct settings for each tool.

Table 2: Correcting for Library Strandedness

Library Type	STAR Column	HTSeq-count -s	featureCounts -s
Unstranded	2	no	0
Stranded (e.g., dUTP, Illumina TruSeq)	4	reverse	2
Reversely Stranded	3	yes	1

Annotation File Issues

The structure and content of the Gene Transfer Format (GTF) file can directly impact counting.

• Problem: A documented issue involves the GTF file's attribute used to group features into genes. The default STAR parameter --sjdbGTFtagExonParentTranscript is set to transcript_id. If a GTF file has empty transcript_id attributes, STAR cannot correctly assign reads to genes, leading to misquantification [38].
• Solution: If the GTF file uses gene_id as the primary feature identifier, the STAR command should be modified with --sjdbGTFtagExonParentTranscript gene_id during the genome generation or mapping step [38]. This ensures reads are grouped by the correct gene identifier.

Paired-End Read Counting

A fundamental difference in how paired-end data is handled can cause a near two-fold difference in counts.

• Problem: By default, STAR and HTSeq-count for paired-end data count a successfully aligned read pair as a single fragment (one count) [38] [42]. However, the default behavior of featureCounts is to count reads, not pairs. This results in featureCounts producing counts that are approximately twice as high as STAR's if not corrected [38].
• Solution: When running featureCounts on paired-end data, always include the --countReadPairs flag. This ensures it counts fragments, making its output directly comparable to STAR and HTSeq-count [38].

Read Sorting and Multi-Mapping

• Problem (HTSeq-count): For paired-end data, HTSeq-count requires the input BAM file to be sorted in a specific order. The default expected order is by name, but if the BAM is sorted by position (a common output from aligners like STAR) and the -r pos parameter is not specified, HTSeq-count will produce incorrect, often significantly lower, counts [42].
• Solution: Always specify the sort order of your input BAM file when using HTSeq-count with the -r parameter (-r pos for most STAR-generated BAMs) [42].
• Note on Multi-mappers: All tools discussed here, by default, exclude reads that map to multiple locations in the genome (multi-mappers) from gene-level counts. This is a deliberate design choice for differential expression analysis to maintain the accuracy of expression ratios between samples [15].

Recommended Experimental Protocol

Verified Workflow for Consistent Counts

The following step-by-step protocol ensures consistent results between STAR and featureCounts, based on a resolved case study [38].

• STAR Genome Indexing:

• STAR Alignment and Quantification:

  • After running, use the ReadsPerGene.out.tab file and the strandedness diagnosis method (Section 3.1) to select the correct count column (e.g., column 4 for stranded protocols).

• featureCounts Quantification (from STAR BAM):

Validation and Troubleshooting

• Cross-Tool Validation: After running both pipelines, compare counts for a few specific genes. The counts should be highly correlated, if not identical.
• Visual Inspection in IGV: For genes with persistent discrepancies, manually inspect the read alignments in a genome browser like IGV to understand how the different tools might be assigning or excluding reads [42].
• Check Special Counters: Always examine the special counters at the bottom of the output files (__no_feature, __ambiguous, N_multimapping, etc.) in STAR and HTSeq-count. Significant differences in these categories can pinpoint the source of the problem.

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools

Item	Function/Description
STAR Aligner	Ultrafast splice-aware aligner for RNA-seq data; includes integrated read counting [39] [40].
featureCounts	High-speed read quantification program; part of the Subread package, ideal for high-throughput workflows [41].
HTSeq-count	Python-based script for precise read counting; widely used and cited in differential expression studies [41] [15].
Reference GTF File	Gene annotation file (e.g., from GENCODE or Ensembl); must be consistent between alignment and quantification steps.
RSeQC Package	Contains infer_experiment.py, a tool to empirically determine the strandedness of your RNA-seq library [38].
Integrative Genomics Viewer (IGV)	High-performance visualization tool for interactive exploration of large genomic datasets; essential for manual validation of ambiguous reads [42].

Workflow and Decision Diagrams

The following diagram illustrates the logical workflow for diagnosing and resolving counting discrepancies.

Diagram 1: A systematic diagnostic workflow for resolving gene count discrepancies between STAR and other quantification tools.

Achieving consistent gene counts between STAR's --quantMode GeneCounts and standalone tools like featureCounts and HTSeq-count is readily attainable by systematically addressing common pitfalls. The key parameters requiring rigorous validation are library strandedness, annotation file attributes, and for paired-end data, the counting unit (fragments vs. reads). By adhering to the detailed protocols and troubleshooting guide provided in this Application Note, researchers can ensure the accuracy and reproducibility of their RNA-seq quantification, forming a solid foundation for all subsequent differential expression and biomarker discovery efforts in drug development.

The Spliced Transcripts Alignment to a Reference (STAR) is a widely used RNA-seq aligner that can perform mapping and quantification in a single step. A key feature for expression analysis is the --quantMode GeneCounts option, which directs STAR to count the number of reads per gene during the alignment process. These counts form the foundational data for downstream differential expression analysis [4].

When using --quantMode GeneCounts, STAR outputs a file named ReadsPerGene.out.tab with four columns [4]:

• Column 1: Gene ID
• Column 2: Counts for unstranded RNA-seq
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count option -s yes)
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count option -s reverse)

The correct column must be selected based on the strandedness of the RNA-seq library preparation protocol to ensure accurate gene-level quantification [4].

The Challenge of Non-Standard GTF/GFF3 Files

STAR requires gene annotations in GTF or GFF3 format to perform gene counting. The accuracy of --quantMode GeneCounts depends entirely on the correct interpretation of parent-child relationships between exons, transcripts, and genes within these annotation files. Non-standard files, particularly many GFF3 files, often present two major challenges:

• Inconsistent Exon-to-Gene Linking: Many GFF3 files specify the parent transcript for each exon but omit the direct parent gene. Since gene-level counting requires knowing which gene an exon belongs to, this missing link can prevent proper gene assignment [43].
• Non-Standard Feature Types: The exon feature, which STAR uses by default to define counting intervals, might be labeled differently (e.g., CDS). Furthermore, the identifiers used to link features (e.g., Parent, ID, gene_id) can vary significantly between annotation sources [44].

Failure to correctly configure STAR for these non-standard files often results in a warning like: WARNING: ... no gene_id for line: ... and can lead to a significant number of reads being categorized as N_noFeature in the output, indicating a failure to associate reads with genes [43].

Key Parameters for Handling Non-Standard Annotations

To address these challenges, STAR provides specific parameters that allow users to define how to extract the necessary information from the annotation file.

Parameter Specifications

Table 1: Key STAR Parameters for Non-Standard GTF/GFF3 Files

Parameter	Default Value	Function	Use Case
--sjdbGTFfeatureExon	exon	Specifies the feature type in the GTF/GFF3 file that corresponds to an exon.	Use when exons are labeled as CDS or another term instead of "exon" [43].
--sjdbGTFtagExonParentTranscript	parent	Specifies the attribute tag that links an exon to its parent transcript.	Use when the link uses a tag other than "Parent", such as transcript_id [44].
--sjdbGTFtagExonParentGene	gene_id	Specifies the attribute tag that links an exon to its parent gene.	Critical for GFF3 files where exons only link to a transcript; used with --quantMode GeneCounts [43] [44].

Parameter Interaction Logic

The relationship between these parameters and the annotation file's structure is logical but must be precisely defined. The following diagram illustrates the decision-making workflow STAR uses to assign exons to genes based on user-provided parameters.

Experimental Protocol: Resolving Parent-Gene Links in GFF3

This protocol provides a step-by-step methodology for generating a STAR genome index and performing alignment with gene counting when working with a non-standard GFF3 file where exons list only a transcript as a parent.

Problem Diagnosis and File Inspection

• Initial Run and Log Check: Run the STAR genome generation step with a basic command. Upon completion, carefully inspect the Log.out file for warnings such as WARNING: ... no gene_id for line: ... [43].

• Inspect GFF3 Structure: Examine the structure of your GFF3 file using command-line tools.

  Example output of a problematic GFF3 structure:

  In this example, the exon has a Parent=mRNA_1 but no direct link to gene_1.

Solution A: Using STAR Parameters (Direct Method)

If inspection reveals that the transcript features (e.g., mRNA) themselves contain the gene ID in their attributes, you can resolve the hierarchy using STAR's parameters directly.

• Identify Key Tags: From the GFF3 inspection, determine:

  • The feature type for exons (e.g., exon or CDS).
  • The tag used for the exon's parent (e.g., Parent).
  • The tag in the transcript's attributes that contains the gene ID (e.g., Parent pointing to gene_1 in the mRNA line above).

• Build Genome Index with Explicit Parameters:

  In this command, --sjdbGTFtagExonParentGene Parent instructs STAR to find the gene ID by looking at the Parent attribute of the transcript that the exon belongs to.

Solution B: Converting GFF3 to GTF (Recommended)

The developer of STAR recommends converting GFF3 files to GTF format using gffread from the Cufflinks package, as this often resolves parent-gene linkage issues more robustly [43].

• Software Requirement: Ensure gffread is installed.

• Execute GFF3 to GTF Conversion:

• Build Genome Index with Converted GTF:

  The GTF format typically has a more straightforward structure with gene_id and transcript_id attributes directly available for each exon, simplifying the process for STAR.

Alignment and Quantification

After generating the genome index with either method, proceed with the alignment and quantification step.

The Scientist's Toolkit: Essential Research Reagents and Materials

Table 2: Key Resources for RNA-seq Analysis with STAR

Category	Item / Software	Function / Description	Key Parameters / Attributes
Alignment & Quantification	STAR Aligner	Maps RNA-seq reads to a reference genome and can perform inline gene expression counting.	--quantMode GeneCounts, --sjdbGTFfile [4]
Genome Reference	Reference Genome (FASTA)	The reference sequence to which RNA-seq reads are aligned.	Must be consistent with the annotation file.
Genome Annotation	Gene Annotation (GTF/GFF3)	File containing coordinates and hierarchies of genomic features (genes, transcripts, exons).	Source (e.g., Ensembl, RefSeq), version.
File Conversion	gffread (Cufflinks)	Utility for converting GFF3 files to GTF format, resolving parent-gene linkage issues.	-T (output GTF) [43]
Downstream Analysis	HTSeq-count	Alternative read-counting tool; used for validation or if quantification is performed post-alignment.	-s (strandedness) [4]
Critical STAR Parameters	--sjdbGTFtagExonParentGene	Solver for GFF3 issues: Specifies the GFF attribute for the exon's parent gene.	Value is often ID or Parent [43] [44]
	--sjdbGTFfeatureExon	Defines which feature type in the annotation file represents an exon.	Default: exon; Alternative: CDS [43]

The STAR aligner, when run with the --quantMode GeneCounts option, performs read alignment and quantifies gene expression simultaneously. A key output file, ReadsPerGene.out.tab, contains the raw count data used in downstream differential expression analysis. The file includes four columns: column 1 contains gene identifiers, column 2 provides counts for unstranded RNA-seq, column 3 for the 1st read strand aligned with RNA (corresponding to htseq-count -s yes), and column 4 for the 2nd read strand aligned with RNA (corresponding to htseq-count -s reverse) [4] [5].

In this file, several summary lines begin with "N_" which provide crucial diagnostic information about reads that were not assigned to genes [4]:

• N_unmapped: Reads that could not be aligned to the reference genome.
• N_multimapping: Reads that aligned to multiple genomic locations.
• N_noFeature: Reads that aligned to genomic regions not corresponding to any annotated gene feature in the provided GTF/GFF file.
• N_ambiguous: Reads that overlapped multiple genes and could not be uniquely assigned [4].

A high proportion of N_noFeature reads, especially when accompanied by low gene counts, indicates that a significant portion of your data is not contributing to gene expression measurements, potentially reducing statistical power and reliability of results. This guide provides a systematic approach to diagnosing and resolving the underlying causes.

Diagnostic Workflow for High N_noFeature and Low Counts

The following diagram outlines a logical, step-by-step diagnostic process for investigating high N_noFeature and low gene counts in your STAR output. This workflow guides you from initial data checks to specific investigations of common issues.

Diagnostic Checks and Protocols

Verifying Library Strandedness Configuration

Objective: To ensure the correct column from the ReadsPerGene.out.tab file is selected based on the experimental library preparation protocol. Using the wrong column is a primary cause of low counts and a high N_noFeature rate.

Background: In a standard stranded RNA-seq protocol (e.g., dUTP method), the sequence of Read 2 corresponds to the original RNA strand. Therefore, Read 2 maps to the sense strand of the gene. During STAR alignment, reads are checked for overlaps with gene annotations on both genomic strands. The N_noFeature count will be artificially inflated if reads are assessed for overlap against the incorrect strand [5].

Protocol:

• Determine Your Library Type: Consult your lab notes or sequencing facility documentation to confirm whether the library preparation was unstranded, stranded-forward (e.g., Standard Illumina), or stranded-reverse (e.g., dUTP method).
• Select the Correct Count Column:
  • For unstranded libraries, use column 2.
  • For stranded-forward libraries (where Read 1 is sense), use column 3.
  • For stranded-reverse libraries (where Read 2 is sense, common with dUTP), use column 4 [4] [5].
• Perform a Strandedness Check: You can computationally validate strandedness using the ReadsPerGene.out.tab file itself. As demonstrated in an example analysis [5], the following command compares the total counts in columns 3 and 4 for all genes:
  • Interpretation: In a correctly identified stranded-reverse library, you expect the count in column 4 to be vastly greater than in column 3 [5]. If they are similar, you may be using the wrong strandedness parameter.

Inspecting Annotation (GTF/GFF) Compatibility

Objective: To verify that the gene annotation file used during STAR genome generation or mapping is compatible with your genome assembly and sequencing reads.

Background: The N_noFeature category consists of reads that map to the genome but do not overlap any feature (e.g., exon) defined in the annotation file [4]. This occurs if there is a mismatch between the genome assembly and the annotations, if the annotations are incomplete, or if novel transcripts are expressed.

Protocol:

• Check Genome-Annotation Version Match: Ensure the versions of your reference genome FASTA file and annotation GTF file match (e.g., both are GRCh38.p13 from Ensembl). Using GENCODE annotations with an Ensembl genome of the same version is usually acceptable, but mixing major versions (e.g., GRCh37 with GRCh38 annotations) will cause widespread assignment failures.
• Validate GTF File Content: Inspect the first few lines of your GTF file to confirm it contains gene and exon features.
  Look for lines where the third column (feature type) is gene or exon.
• Investigate Read Mapping with IGV: For a deeper dive, visually inspect where unassigned reads are mapping.
  • Extract a subset of unmapped reads using samtools:

  • Alternatively, randomly sample reads flagged as "noFeature" by their SAM tags (this requires a custom script).
  • Load your sorted BAM file, the annotation GTF file, and the genome into the Integrative Genomics Viewer (IGV). Navigate to genomic coordinates with high read coverage but no annotated features. This can reveal whether reads are mapping to unannotated regions, such as novel isoforms, non-coding RNAs, or genomic regions absent from your GTF file.

Checking Read Quality and Adapter Content

Objective: To identify issues with raw read quality or adapter contamination that may prevent reads from mapping to gene features, even if they are of transcriptional origin.

Background: Sequencing artifacts, such as low-quality bases or undetected adapter sequences, can prevent STAR from aligning reads to the reference genome or can cause them to align to incorrect, intergenic locations, inflating the N_noFeature statistic [45].

Protocol:

• Run FastQC: Perform quality control on your raw FASTQ files using FastQC.

• Interpret FastQC Report: Pay close attention to:
  • "Per base sequence quality": Ensures base qualities are generally high (Q>30) across all cycles.
  • "Adapter Content": Determines if adapter sequences are present in your data. A significant presence (>1%) requires trimming.
  • "Sequence Length Distribution": Confirms read lengths are as expected.
• Trim Adapters and Low-Quality Bases: If the FastQC report indicates issues, use a tool like Trimmomatic or Cutadapt [45].
  • Example with Trimmomatic (PE data):

• Re-run STAR and FastQC: Map the trimmed reads with STAR and run FastQC again on the trimmed files to confirm the issues have been resolved.

Reviewing STAR Mapping Parameters

Objective: To ensure that key parameters controlling mapping sensitivity and accuracy are appropriately set, and not inadvertently causing valid reads to be discarded or misclassified.

Background: STAR's mapping behavior is controlled by numerous parameters. Overly stringent parameters can cause reads derived from transcripts to fail to align or to align poorly, pushing them into the unmapped or noFeature categories [46].

Protocol:

• Check Critical Alignment Parameters: Review your STAR command for the following parameters:
  • --outFilterScoreMinOverLread and --outFilterMatchNminOverLread: These control the minimum alignment score and matched bases relative to read length. The defaults are often sufficient, but increasing them too much can discard valid, shorter alignments.
  • --alignIntronMin and --alignIntronMax: Define the minimum and maximum intron sizes. The defaults (21 and 0, which means 1000000) are suitable for most eukaryotes. Setting --alignIntronMin too high can prevent the detection of small introns.
  • --alignSJoverhangMin and --alignSJDBoverhangMin: Control the minimum overhang for unannotated and annotated splice junctions. The default 5 is reasonable, but increasing it may miss some junctions.
• Consult the STAR Log File: Examine the Log.final.out file from your STAR run. Key metrics to check:
  • Unmapped reads: Should typically be low (<10-20%).
  • Mismatch rate per base: Should be low (~0.5-1.5%).
  • Deletion and Insertion rates: Should be very low. A high percentage of unmapped reads could indicate issues with read quality, adapter contamination, or overly strict parameters.
• Consider Two-Pass Mapping: For improved mapping of novel splice junctions, which can reduce N_noFeature by assigning reads to previously unannotated exons, use the --twopassMode Basic option [4]. This mode performs two alignment passes, using the splice junctions discovered in the first pass as annotations for the second pass.

The Scientist's Toolkit: Research Reagent Solutions

Table 1: Essential reagents and software tools for RNA-seq analysis and troubleshooting.

Item Name	Function / Purpose	Example / Note
STAR Aligner	Spliced alignment of RNA-seq reads to a reference genome.	Requires a genome index built with --runMode genomeGenerate [46].
High-Quality Annotation File	Provides genomic coordinates of genes and exons for read quantification.	Use from GENCODE or Ensembl; ensure version matches genome assembly [5].
FastQC	Quality control tool for high-throughput sequence data.	Identifies adapter contamination, low-quality bases, and other technical issues [45].
Trimmomatic / Cutadapt	Read trimming tool to remove adapters and low-quality bases.	Critical for cleaning data before alignment to improve mapping rates [45].
SAMtools	Utilities for manipulating and viewing alignments in SAM/BAM format.	Used for sorting, indexing, and extracting subsets of BAM files for analysis [5].
Integrative Genomics Viewer (IGV)	High-performance visualization tool for genomic data.	Essential for visually diagnosing mapping issues and inspecting noFeature reads [45].
R/Bioconductor Packages	Statistical analysis and visualization of count data (e.g., DESeq2, edgeR).	Used for downstream differential expression analysis after count quantification [45].

STAR quantMode vs. Other Quantification Tools: A Strategic Comparison for Your RNA-seq Study

Table of Contents

• Executive Summary
• Tool Overview and Core Characteristics
• Detailed Methodologies and Protocols
• Performance and Benchmarking Insights
• The Scientist's Toolkit
• Conclusions and Strategic Recommendations

In RNA-seq analysis, the step of read quantification is critical for accurate gene expression measurement. Among the plethora of available tools, alignment-based counters like STAR GeneCounts, RSEM, and Subread featureCounts are widely used, each with distinct philosophies and performance characteristics. STAR GeneCounts offers integrated quantification during alignment, providing HTSeq-compatible counts. Subread featureCounts is a fast and efficient standalone counter, while RSEM employs a sophisticated probabilistic model to estimate expression at both gene and transcript levels, particularly handling multi-mapped reads. Benchmarking studies reveal that while these tools show high correlation in overall results, their accuracy varies, especially for genes with low or high expression levels and for complex isoforms. Tools like featureCounts, which discard multi-mapping reads, can underpredict expression for repetitive genes, whereas RSEM's assignment strategy can improve accuracy in such contexts [47] [48] [49]. The choice of tool thus depends on the specific research goals, such as the need for speed, transcript-level resolution, or handling of ambiguous reads. This application note provides a detailed comparison, experimental protocols, and data-driven recommendations to guide researchers and drug development professionals in selecting and implementing the optimal quantification strategy.

The following table summarizes the core attributes, advantages, and limitations of the three featured quantification tools.

Table 1: Core Characteristics of STAR GeneCounts, RSEM, and featureCounts

Feature	STAR --quantMode GeneCounts	RSEM (RNA-Seq by Expectation Maximization)	Subread featureCounts
Primary Function	Integrated read alignment and gene-level quantification [4] [50]	Transcript- and gene-level abundance estimation from BAM files [47] [49]	Fast, standalone gene-level read summarization [51]
Core Algorithm	Direct counting based on genomic overlap (mimics HTSeq) [4]	Expectation-Maximization (EM) algorithm to resolve multi-mapping reads [47] [49]	Efficient chromosome hashing and feature overlap algorithm [51]
Handling of Multi-mapped Reads	Does not count them; reports total number separately as N_multimapping [10]	Probabilistically assigns them to all compatible transcripts, influencing gene-level estimates [49]	Can be excluded (default) or included with a fractional count [51]
Isoform Resolution	No (Gene-level only) [50]	Yes (Transcript-level is its primary function) [49]	No (Gene-level only) [51]
Key Advantage	Convenience and speed; generates counts simultaneously with alignment [50]	High accuracy for transcript-level quantification; sophisticated handling of ambiguities [50] [49]	Extremely fast and memory-efficient; simple, transparent counting logic [51] [48]
Main Limitation	Does not attempt to resolve multi-mapping or transcript ambiguities [4]	Computationally intensive; requires more resources and steps [50]	Simple counting may underperform for genes with complex paralogs or isoforms [49]

The workflow relationships and fundamental differences in how these tools assign reads are illustrated below.

Detailed Methodologies and Protocols

Protocol for STAR--quantMode GeneCounts

STAR's integrated quantification is activated via a single command-line flag during the alignment step. Its output requires careful selection of the column corresponding to the library strandedness.

• Typical Command:

  The --quantMode GeneCounts option instructs STAR to count reads per gene during alignment [4].

• Output Interpretation: STAR generates a *ReadsPerGene.out.tab file. The first four lines are summary statistics, followed by gene counts in four columns [4]:

  • Column 1: Gene identifier.
  • Column 2: Unstranded counts.
  • Column 3: Counts for forward strand (for -s yes in HTSeq).
  • Column 4: Counts for reverse strand (for -s reverse in HTSeq).

• Determining Strandedness: To select the correct column, calculate the total assigned reads (N_genic = TotalReads - N_unmapped - N_multimapping - N_noFeature - N_ambiguous) for columns 3 and 4. The column with the higher N_genic value corresponds to the correct strand specificity of your library [10].

Protocol for SubreadfeatureCounts

featureCounts is a standalone tool known for its speed and is often run on a BAM file generated by STAR or other aligners.

• Typical Command:

  • -s 2: Indicates a reverse-stranded library (use -s 0 for unstranded, -s 1 for forward-stranded).
  • -p: Specifies that fragments (paired-end reads), not individual reads, will be counted.
  • -t exon & -g gene_id: Define the feature and attribute in the GTF file used for counting [51].

• Key Considerations:

  • By default, featureCounts does not count multi-mapping reads, which can lead to lower counts than expected if such reads are present [51].
  • It includes an option (--countReadPairs) to count read pairs instead of reads for paired-end data, which is crucial for accurate fragment counting [51].

Protocol for RSEM

RSEM is a more complex tool that requires a preparation step to build a transcriptome reference and then calculates expression. It can use its own alignment method or work with a STAR-generated BAM file.

• Workflow with STAR Alignment:

  • Prepare Reference: First, build a custom RSEM reference from a FASTA and GTF file.

  • Calculate Expression: Use STAR as the aligner within RSEM.
    This workflow leverages STAR's alignment accuracy while utilizing RSEM's powerful quantification model [50] [49].

• Output: RSEM outputs several files, with the *.genes.results file containing gene-level expected counts, TPM (Transcripts Per Million), and FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values. These expected counts are suitable for input to differential expression tools like DESeq2 after rounding [50].

Performance and Benchmarking Insights

Independent benchmarking studies provide critical insights into the relative performance of quantification tools, evaluating them against gold standards like qRT-PCR and simulated data with known "ground truth."

Table 2: Performance Benchmarks from Comparative Studies

Benchmark Metric	STAR GeneCounts / HTSeq	RSEM	featureCounts / HTSeq	Context and Notes
Correlation with qRT-PCR (R²)	Highest correlation (0.85-0.89) reported for HTSeq [47]	Correlates well, but slightly lower than HTSeq in some studies [47]	High correlation, similar to other count-based methods [48]	Measures how well RNA-seq expression estimates match validation data [47].
Accuracy vs. Simulation (MARDS)	N/A	More accurate for transcript-level quantification [49]	Simple counting can be less accurate for complex genes [49]	MARDS (Mean Absolute Relative Difference) measures deviation from simulated truth; lower is better [49].
Sensitivity to Expression Level	Reliable for medium to high abundance genes [48]	Can quantify lowly expressed isoforms, but accuracy is affected [49]	May be less sensitive for very low or high expression genes [48]	Performance varies with expression level; complex isoforms are challenging for all [48].
Impact of Multi-mapped Reads	Reads are discarded, potentially lowering counts for some genes [10]	Reads are probabilistically assigned, improving completeness [49]	Reads are typically discarded (default), similar to STAR [51]	A key differentiator; assignment strategies recover otherwise lost information [49].

The following diagram synthesizes the decision logic for choosing the appropriate tool based on common research scenarios and priorities.

The Scientist's Toolkit

Successful implementation of RNA-seq quantification pipelines relies on a suite of reliable software and reference materials. The table below lists key reagents and computational tools essential for this field.

Table 3: Essential Research Reagents and Tools for RNA-seq Quantification

Item	Function / Description	Example / Note
STAR Aligner	Spliced Transcripts Alignment to a Reference; performs alignment and can generate integrated gene counts [4] [48]	https://github.com/alexdobin/STAR
Subread Package	Contains the featureCounts program for fast read summarization [51] [48]	http://subread.sourceforge.net/
RSEM	Estimates gene and transcript abundances using an EM algorithm for accurate quantification [47] [49]	https://deweylab.github.io/RSEM/
Reference Genome	The DNA sequence of the organism under study; required for alignment [48]	e.g., GRCh38 for human; GRCm38 for mouse
Annotation File (GTF/GFF)	File containing genomic coordinates of genes, transcripts, and exons; essential for read assignment [4] [48]	e.g., From GENCODE or Ensembl
Salmon	A fast, alignment-free tool for transcript-level quantification; often used as a modern alternative or complement [50] [49]	https://combine-lab.github.io/salmon/
DESeq2 / edgeR	R/Bioconductor packages for differential expression analysis following count quantification [52] [48]	Accepts count matrices from all tools listed

Based on the comprehensive comparison and benchmarking data, the following strategic recommendations are proposed:

• For Standard Gene-level Differential Expression: Pipelines using STAR with --quantMode GeneCounts or featureCounts are highly robust and widely accepted. They are simple, fast, and produce reliable count matrices for tools like DESeq2 and edgeR [48]. The choice between them can be based on convenience (STAR's integration) versus the need for a standalone counter (featureCounts).

• For Transcript-level Analysis or Complex Genes: RSEM is the superior choice when the research question requires quantification of transcript isoforms or when working with gene families with high sequence similarity. Its probabilistic model for resolving multi-mapping reads provides more accurate and complete expression estimates in these challenging scenarios [49].

• For Maximum Accuracy and Modern Workflows: Consider using a hybrid approach, such as STAR alignment followed by RSEM quantification, or adopting alignment-free tools like Salmon or Kallisto. These alignment-free tools have been shown to be both fast and accurate, often outperforming traditional alignment-based counting in benchmarks [50] [49].

Ultimately, the quantification tool should not be selected in isolation but as a component of a holistic RNA-seq analysis strategy, considering the alignment software, reference annotations, and the biological questions driving the research.

RNA sequencing (RNA-seq) data analysis fundamentally relies on determining the origin of sequenced fragments. The choice of method for this task significantly impacts the speed, resource requirements, and accuracy of downstream results such as transcript abundance estimation and differential expression analysis. Exact alignment tools, like STAR, perform base-by-base alignment of reads to a reference genome, providing precise genomic coordinates. In contrast, pseudoalignment tools, such as Salmon and Kallisto, use lightweight algorithms to rapidly determine transcript compatibility of reads without exact base-level alignment [53] [54]. This article explores the technical distinctions between these approaches, focusing on their implications for quantifying gene expression, with specific attention to the context of STAR's quantMode GeneCounts output.

Core Computational Concepts and Differences

Foundational Algorithms and Workflows

The fundamental difference lies in the core algorithm and the primary reference used for mapping.

• STAR (Exact Alignment): STAR uses a traditional alignment-based approach. It employs a maximal mappable seed search to map reads to a reference genome, identifying all possible positions for each read [55] [56]. It is splice-aware, meaning it can identify exon-intron boundaries, making it suitable for detecting novel splice junctions and genetic variations [55]. Its workflow typically involves first aligning reads to the genome and then counting reads per gene as a secondary step [53].

• Salmon & Kallisto (Pseudoalignment): These tools use alignment-free, lightweight algorithms that work directly from a reference transcriptome (in FASTA format). They avoid the computationally intensive base-to-base alignment [54].

  • Kallisto introduces the concept of pseudoalignment, which rapidly determines the set of transcripts compatible with a given sequenced fragment using a pre-built index of k-mers (default k=31) without determining its exact position [57] [58].
  • Salmon utilizes a related but distinct quasi-mapping approach and incorporates a sophisticated dual-phase parallel inference algorithm. A key enhancement in recent Salmon versions is selective alignment, which uses a more sensitive scheme for finding potential mapping loci and scores them, bridging a gap between pure pseudoalignment and traditional alignment [59] [60].

The diagram below illustrates the fundamental workflow differences between these approaches.

Key Feature and Performance Comparison

The following table summarizes the critical differences between the tools, influencing researchers' choices based on their experimental goals and computational resources.

Table 1: Comparative analysis of STAR, Kallisto, and Salmon.

Feature	STAR	Kallisto	Salmon
Core Method	Exact alignment to the genome [4] [55]	Pseudoalignment to the transcriptome [57] [58]	Quasi-mapping / Selective alignment to the transcriptome [59] [60]
Primary Output	Genomic alignments (BAM), gene-level read counts [4] [55]	Transcript-level abundances (TPM, estimated counts) [55] [58]	Transcript-level abundances (TPM, estimated counts) [59] [60]
Speed	Slower; can take hours per sample [53] [55]	Very fast; minutes per sample [53] [57]	Very fast; comparable to Kallisto [60]
Memory Usage	High (e.g., ~30 GB for human genome) [53] [56]	Low (can run on a laptop) [53] [57]	Low to Moderate [60]
Strengths	Detects novel junctions, variants; provides genomic context [53] [55]	Extremely fast and resource-efficient for quantification [53] [57]	Fast and feature-rich with advanced bias correction [60]
Limitations	Computationally intensive; quantification is not its primary strength [53]	Limited to annotated transcriptome; no genomic discovery [53] [54]	Limited to annotated transcriptome; no genomic discovery [54]

STAR quantMode GeneCounts: An Exact Alignment Quantification Model

The STAR --quantMode GeneCounts option provides a direct method for obtaining gene-level counts from exact alignments, which serves as a useful reference point for understanding quantification output.

Mechanism and Output Interpretation

When using --quantMode GeneCounts, STAR counts reads that overlap (by at least 1 nucleotide) one and only one gene, based on the provided GTF/GFF annotation file [4]. This method produces a ReadsPerGene.out.tab file. A key feature of this output is its handling of strandedness, which is crucial for accurate interpretation in strand-specific protocols [4].

Table 2: Interpretation of STAR GeneCounts output columns.

Column	Content	Recommended Use
Column 1	Gene ID	Gene identifier from the annotation file.
Column 2	Counts for unstranded RNA-seq	Use for non-strand-specific protocols.
Column 3	Counts for the 1st read strand aligned with RNA	Use for strand-specific protocols (e.g., -s yes in htseq-count).
Column 4	Counts for the 2nd read strand aligned with RNA	Use for strand-specific protocols (e.g., -s reverse in htseq-count).

The sums of columns 3 and 4 may not equal column 2 because unstranded counts (column 2) consider reads overlapping genes on opposite strands as "ambiguous" and discard them, whereas the stranded columns (3 and 4) assign these reads to a specific gene based on strand [4].

Technical Considerations

• Input Requirements: This mode requires gene annotations (--sjdbGTFfile) provided during genome generation or at the mapping stage [4].
• Read Assignment: It uses an "intersection-nonempty" rule similar to htseq-count with default parameters. A read is counted if it overlaps one and only one gene [4].
• Limitations: The counting model is straightforward and does not employ statistical models to resolve multi-mapping reads across isoforms or paralogous genes in a principled manner. Reads that map to multiple genes are discarded [53] [4].

Protocols for Quantification

STAR Alignment and Gene Counting

This protocol generates gene-level counts via exact genomic alignment.

• Genome Generation: First, generate a STAR genome index.

• Alignment & Quantification: Map reads and count gene overlaps in a single step.
  The gene counts will be in sample_ReadsPerGene.out.tab [4].

Kallisto Pseudoalignment and Quantification

This protocol uses pseudoalignment for rapid transcript-level quantification.

• Indexing: Build a kallisto index from a transcriptome file.

• Quantification: Run the pseudoalignment and quantification.
  The output abundance.tsv file in the output_dir contains transcript abundance estimates in TPM and estimated counts [61] [58].

Salmon Quasi-mapping and Quantification

This protocol uses quasi-mapping with selective alignment for accurate, bias-aware quantification.

• Indexing: Build a Salmon index. Using a decoy-aware transcriptome is recommended.

• Quantification: Quantify reads with bias correction and validation mappings enabled.
  The output quant.sf file in the salmon_quant directory contains the abundance estimates [59] [54].

Successful execution of these computational protocols requires specific, well-curated input files and resources.

Table 3: Essential materials and resources for RNA-seq alignment and quantification.

Item	Function	Example Source & Notes
Reference Genome	Exact alignment with STAR requires a reference genome sequence as the alignment target.	ENSEMBL, UCSC Genome Browser; provided as a FASTA file (e.g., GRCh38.primary_assembly.genome.fa).
Reference Transcriptome	Pseudoalignment with Salmon and Kallisto requires a set of known transcript sequences.	ENSEMBL (e.g., Homo_sapiens.GRCh38.cdna.all.fa.gz); must match the organism and assembly version of the genome [54].
Genome Annotations	File linking genomic coordinates to gene and transcript models. Essential for STAR's gene counting and for assigning transcriptomic data to genomic coordinates.	GTF or GFF3 format from ENSEMFL or GENCODE. The choice between "comprehensive" and "filtered" annotations can impact results, especially in scRNA-seq [56].
Sequencing Reads	The raw data from the RNA-seq experiment to be analyzed.	Typically gzipped FASTQ files. Quality control (e.g., with FastQC) is recommended before alignment/quantification.
Salmon/Kallisto Index	A pre-processed, searchable representation of the reference transcriptome, enabling rapid pseudoalignment.	Generated by the salmon index or kallisto index commands. Pre-built indices for common organisms are available via refgenie [59].
STAR Genome Index	A pre-processed representation of the reference genome, enabling fast and accurate splice-aware alignment.	Generated by the STAR --runMode genomeGenerate command. Requires significant memory and time to build.

The choice between exact alignment with STAR and pseudoalignment with Salmon or Kallisto is fundamental and depends on the research objectives. STAR's exact alignment is indispensable for discovery-based genomics, including the identification of novel splice junctions, fusion genes, and genetic variants. Its integrated GeneCounts provides a straightforward, alignment-based quantification metric. In contrast, Salmon and Kallisto offer a paradigm of speed and efficiency for transcript quantification, excelling in resource-limited environments and large-scale studies where the goal is accurate abundance estimation rather than genomic discovery. Furthermore, Salmon's integrated bias models can improve the accuracy of abundance estimates and the reliability of subsequent differential expression analysis [60]. Understanding these core differences allows researchers to strategically select the tool that best aligns with their experimental design and biological questions.

In RNA-seq data analysis, determining the abundance of gene expression is a fundamental step. The process of assigning sequencing reads to genes, known as "counting," can be approached through different methodologies, primarily categorized as integrated counting and advanced ambiguity resolution. Integrated counting, as exemplified by STAR's quantMode GeneCounts, offers a streamlined workflow by performing alignment and quantification simultaneously [4]. In contrast, advanced ambiguity resolution methods, such as those employed by RSEM and Kallisto, utilize more sophisticated statistical models to handle the complex issue of reads that map to multiple genes or isoforms [50]. The choice between these approaches represents a critical trade-off between analytical convenience and resolution power, directly impacting the accuracy and interpretability of downstream differential expression analyses. This application note delineates the operational parameters, advantages, and limitations of each method to guide researchers in selecting the optimal quantification strategy for their specific experimental context within drug development and basic research.

Key Concepts and Definitions

Integrated Counting

Integrated counting refers to quantification methods built directly into alignment tools. The STAR aligner's --quantMode GeneCounts is a prime example, which counts the number of reads per gene while mapping [4]. Its output, the *ReadsPerGene.out.tab file, provides counts for unstranded and strand-specific protocols based on the assumption that a read should be counted if it overlaps one and only one gene [4]. This method is conceptually equivalent to the popular htseq-count tool and generates a straightforward table of counts suitable for input into differential expression packages like DESeq2 [50].

Advanced Ambiguity Resolution

Advanced ambiguity resolution encompasses tools that employ specialized algorithms to handle the inherent uncertainty in assigning reads to features. This is particularly valuable for resolving:

• Multi-mapping reads: Reads that align equally well to multiple genomic locations.
• Isoform-level quantification: Determining the specific transcript isoforms from which reads originated [50]. Tools such as RSEM (which often uses STAR alignments) and Kallisto (a pseudoaligner) fall into this category. They use probabilistic models to fractionally assign reads to genes or isoforms, thereby providing more nuanced expression measures than simple integer counts [50].

Comparative Analysis: Integrated Counting vs. Advanced Ambiguity Resolution

The table below summarizes the core characteristics, advantages, and disadvantages of integrated counting and advanced ambiguity resolution methods.

Table 1: Comparison of Integrated Counting and Advanced Ambiguity Resolution Methods

Feature	Integrated Counting (e.g., STAR GeneCounts)	Advanced Ambiguity Resolution (e.g., RSEM, Kallisto)
Core Principle	Counts reads overlapping a single gene during alignment [4].	Uses probabilistic models to resolve multi-mapping reads and assign reads to isoforms [50].
Primary Output	Integer read counts per gene.	Estimated counts or abundances, often at the transcript level.
Speed & Convenience	High; quantification occurs during alignment, streamlining the workflow [50].	Varies; RSEM can be slower, while Kallisto is very fast due to pseudoalignment [50].
Resolution	Gene-level only [50].	Transcript-/isoform-level, which can be aggregated to gene-level [50].
Handling of Ambiguity	Simple; a read is counted only if it maps uniquely to one gene. Ambiguous reads are categorized and excluded [4].	Sophisticated; probabilistically distribreads among compatible genes or transcripts, which is considered superior for accuracy [50].
Ease of Interpretation	High; simple count data is straightforward to understand and explain [50].	Lower; expression values are more abstract and require understanding of the underlying statistical model [50].
Best Suited For	Rapid gene-level differential expression analysis where isoform information is not critical.	Studies focusing on alternative splicing, isoform usage, or when the highest quantification accuracy is required [50].

Experimental Protocols

Protocol 1: Gene-Level Quantification with STAR quantMode

This protocol describes a standard workflow for generating gene-level count data using the STAR aligner's integrated counting function.

Research Reagent Solutions:

• STAR Aligner: Spliced Transcripts Alignment to a Reference, used for fast and accurate alignment of RNA-seq reads [4].
• Reference Genome Index: A species-specific genome index built by STAR using a reference genome sequence and annotation file.
• Annotation File (GTF/GFF): A file containing genomic coordinates of genes and transcripts, required for the --sjdbGTFfile parameter [4].

Methodology:

• Alignment and Quantification: Execute the STAR command with the --quantMode GeneCounts flag. This instructs STAR to count reads per gene during the alignment process.

• Output Interpretation: STAR generates a *ReadsPerGene.out.tab file. Select the appropriate column based on your library preparation protocol:
  • Column 2: Use for unstranded RNA-seq data.
  • Column 3: Use for stranded data where the first read strand is the same as the RNA strand (e.g., -s yes in htseq-count).
  • Column 4: Use for stranded data where the first read strand is the reverse of the RNA strand (e.g., -s reverse in htseq-count) [4].
• Downstream Analysis: The selected column of counts can be used directly as input for differential expression analysis with tools like DESeq2 or edgeR.

Protocol 2: Transcript-Level Quantification with RSEM

This protocol leverages STAR for alignment but uses RSEM for quantification, providing more accurate resolution of ambiguous reads and transcript-level abundances.

Research Reagent Solutions:

• RSEM (RNA-Seq by Expectation Maximization): A software package for estimating gene and isoform expression levels from RNA-seq data.
• STAR Aligner: Used here in a specific mode to generate alignments suitable for RSEM.
• Reference Transcriptome: A FASTA file of all transcript sequences for the organism of interest.

Methodology:

• STAR Alignment for RSEM: Run STAR using the -quantMode flag with the TranscriptomeSAM option. This produces a BAM file where reads are aligned to the transcriptome, which is the required input for RSEM.

• RSEM Calculation: Use the rsem-calculate-expression function to process the transcriptome BAM file.

• Output Analysis: RSEM generates *.genes.results and *.isoforms.results files containing expected counts and TPM (Transcripts Per Million) values. The expected counts in the *.genes.results file, after rounding, can be used for gene-level differential expression analysis in DESeq2 [50].

Workflow Visualization

RNA-seq Quantification Workflow Comparison

The choice between integrated counting and advanced ambiguity resolution is not a matter of one being universally superior, but rather depends on the specific research question and experimental goals.

For rapid, gene-level differential expression analysis where the primary aim is to identify which genes are up- or down-regulated, the integrated counting provided by STAR quantMode GeneCounts offers an excellent balance of speed, simplicity, and reliability. Its straightforward count data is easy to interpret and integrates seamlessly with standard analysis pipelines [50].

When the research objective requires transcript-level resolution, involves genes with multiple isoforms, or demands the highest possible accuracy in quantification by resolving multi-mapping reads, advanced ambiguity resolution with tools like RSEM or Kallisto is the recommended approach. Although it adds a step to the workflow and produces more complex output, the ability to dissect isoform-specific expression can be critical in many drug development contexts, such as understanding the mechanisms of action or the effects of splicing-modifying compounds [50].

Researchers should align their choice of quantification method with the biological complexity they seek to uncover. For many standard gene-level analyses, the convenience of integrated counting is sufficient. For deeper investigation into transcriptomic diversity, the power of advanced ambiguity resolution is indispensable.

RNA sequencing (RNA-Seq) has become a foundational technology for transcriptomic research, enabling large-scale inspection of mRNA levels in living cells and leading to insights on gene expression, biomarker discovery, and therapeutic target identification [62]. A critical step in the RNA-Seq workflow is the alignment of sequenced reads to a reference genome and the subsequent quantification of gene-level abundances. The STAR aligner (Spliced Transcripts Alignment to a Reference) is a widely used software that performs this function with high accuracy, capable of processing large volumes of RNA-seq data [7]. A significant feature of STAR is its --quantMode GeneCounts function, which generates a ReadsPerGene.out.tab file containing counts of reads overlapping each gene. These counts are calculated during the mapping process and are designed to coincide with those produced by htseq-count with default parameters, providing a direct input for differential expression analysis tools like DESeq2 [4] [36]. This application note provides a structured framework for selecting the appropriate quantification tool based on a researcher's specific biological questions and available computational resources.

Understanding STAR quantMode GeneCounts Output

Structure and Interpretation of the Output File

When STAR is run with the --quantMode GeneCounts option, it produces a *ReadsPerGene.out.tab file. This file contains four columns that correspond to different strandedness options, which must be selected based on the library preparation protocol used for the RNA-seq experiment [4].

• Column 1: Gene identifier.
• Column 2: Counts for unstranded RNA-seq.
• Column 3: Counts for the 1st read strand aligned with RNA (equivalent to htseq-count option -s yes).
• Column 4: Counts for the 2nd read strand aligned with RNA (equivalent to htseq-count option -s reverse).

The appropriate column must be selected for downstream analysis based on the strandedness of the sequencing library. For unstranded protocols, column 2 is used. For stranded protocols, either column 3 or 4 is selected depending on the specific kit chemistry. Incorrect selection will result in a substantial overestimation of "no feature" counts and an inaccurate representation of gene expression [4].

Key Statistical Considerations for Output Validation

It is important to understand the relationship between the columns in the output file. The sum of columns 3 and 4 does not necessarily equal column 2. This occurs because reads considered "ambiguous" for unstranded analysis (column 2) might be assigned to a specific gene when a strandedness option is applied (columns 3 or 4). Specifically, reads that overlap a feature on the opposite strand are considered ambiguous for unstranded counting but are assigned to genes when using the stranded columns. This explains why values in column 3 or 4 can sometimes be higher than the value in column 2 for the same gene [4].

Figure 1: STAR GeneCounts output interpretation workflow. Researchers must select the appropriate column based on their library preparation protocol.

Decision Framework for Quantification Tool Selection

Comparative Analysis of Quantification Approaches

The choice of quantification method involves trade-offs between analytical accuracy, computational requirements, and suitability for specific research goals. A recent large-scale multi-center benchmarking study involving 45 laboratories and 140 bioinformatics pipelines highlighted significant variations in performance, particularly for detecting subtle differential expression [1].

Table 1: Tool Selection Framework Based on Research Question and Resources

Research Question	Recommended Tool	Computational Demand	Accuracy Profile	Ideal Use Cases
Gene-level DE with reference genome	STAR + GeneCounts	High (RAM-intensive) [7]	High accuracy for known genes [4] [36]	Clinical diagnostics, validation studies [1]
Transcript-level or metagenomic analysis	Pseudoaligners (Salmon, Kallisto)	Low to moderate [7]	Fast with accuracy trade-offs [7]	Exploratory analysis, isoform quantification [7]
Rapid screening or resource-constrained environments	HISAT2 + featureCounts [62]	Moderate	Established, reliable [62]	Beginner-friendly protocols, standard differential expression [62]
Large-scale clinical biomarker discovery	STAR (optimized cloud workflow) [7]	Very High (cloud-scale)	High consistency for subtle expression [1]	Multi-site studies, pharmacogenomics [7] [1]

Experimental and Resource Considerations

When designing an RNA-seq study, researchers must consider multiple experimental and resource factors that influence tool selection:

• Sample Size and Throughput Requirements: For processing tens or hundreds of terabytes of RNA-sequencing data, a scalable, cloud-native architecture is recommended. Cloud-based optimization of STAR can provide significant execution time and cost reduction through techniques like early stopping (23% reduction in total alignment time) and appropriate instance selection [7].

• Subtle Differential Expression Detection: When identifying clinically relevant subtle differential expressions (e.g., between disease subtypes or stages), the accuracy and reproducibility of the tool become critical. A recent study showed greater inter-laboratory variations in detecting subtle differential expressions, with experimental factors including mRNA enrichment and strandedness emerging as primary sources of variation [1].

• Computational Resources: STAR typically requires a large amount of RAM (tens of GiBs, depending on the reference genome size) and high-throughput disks to scale efficiently with increasing thread numbers [7]. For environments with limited computational resources, pseudoaligners or HISAT2 may be more practical alternatives [7] [62].

Figure 2: Decision tree for selecting RNA-seq quantification tools based on research goals and constraints.

Experimental Protocol: End-to-End Workflow with STAR GeneCounts

Sample Preparation and Sequencing Considerations

The accuracy of STAR GeneCounts quantification is highly dependent on proper experimental design and sample preparation. Based on multi-center benchmarking studies, the following recommendations ensure reliable results:

• RNA Quality and Integrity: Use high-quality RNA samples (RIN > 8) to minimize degradation artifacts. The Quartet project reference materials provide quality standards for clinical transcriptomic studies [1].

• Library Preparation: Select stranded mRNA-seq protocols to accurately assign reads to their transcript of origin. Unstranded protocols may result in ambiguous assignments for genes overlapping on opposite strands. The benchmarking study identified mRNA enrichment and library strandedness as primary sources of inter-laboratory variation [1].

• Replication and Sequencing Depth: Include a minimum of three biological replicates per condition for robust statistical power in differential expression analysis. Sequence to a depth of 20-30 million reads per sample for standard gene-level differential expression analysis.

• Spike-in Controls: Incorporate ERCC RNA spike-in controls to monitor technical performance and normalize for sample-specific biases, particularly in studies detecting subtle expression differences [1].

Detailed Computational Protocol for STAR Alignment and Quantification

This protocol provides a step-by-step procedure for generating gene counts from raw FASTQ files using STAR aligner with quantMode GeneCounts.

Step 1: Software Installation and Genome Index Preparation

Step 2: Read Alignment and Gene Counting with STAR

Step 3: Output Processing for Downstream Analysis

After running STAR, the sample_ReadsPerGene.out.tab file will contain the gene counts. For differential expression analysis with DESeq2, create a count matrix by compiling the appropriate column from all samples:

Step 4: Differential Expression Analysis with DESeq2

Research Reagent Solutions and Essential Materials

Table 2: Essential Research Reagents and Computational Resources for RNA-seq Quantification

Item	Function	Example Specifications	Considerations
STAR Aligner [7]	Spliced alignment of RNA-seq reads to reference genome	Version 2.7.10b or newer	Requires ~30GB RAM for human genome; benefits from high-throughput storage
Reference Genome [36]	Reference sequence for read alignment	GENCODE v29 with comprehensive gene annotation	Ensure consistency between FASTA and GTF files used for indexing and alignment
SRA Toolkit [7]	Access and conversion of public sequencing data	prefetch, fasterq-dump utilities	Direct access to NCBI SRA database; files hosted on AWS us-east-1 [7]
DESeq2 [36]	Differential expression analysis	R/Bioconductor package	Uses negative binomial generalized linear models; most genes not DE [36]
Computational Instance [7]	Cloud computing resource for alignment	AWS EC2 optimized for STAR	Cost-efficient instance types identified; spot instances verified applicable [7]
ERCC Spike-in Controls [1]	Technical standards for performance assessment	92 synthetic RNAs with known concentrations	Enables quality assessment at subtle differential expression levels [1]
Quality Control Tools [62]	Assessment of raw read quality	FastQC, Trimmomatic	Identifies adapter contamination and low-quality sequences requiring trimming

Validation and Quality Control Procedures

Performance Metrics and Benchmarking

To ensure the reliability of gene quantification results, implement the following quality control procedures:

• Signal-to-Noise Assessment: Calculate PCA-based signal-to-noise ratio (SNR) values to discriminate the quality of gene expression data. In multi-center studies, SNR values for samples with small biological differences averaged 19.8 (range 0.3-37.6) compared to 33.0 (range 11.2-45.2) for samples with large biological differences [1].

• Ground Truth Validation: Validate absolute gene expression measurements using reference datasets when available. In benchmarking studies, correlations with TaqMan datasets were higher for the Quartet samples (average r=0.876) than for MAQC samples (average r=0.825), highlighting the importance of appropriate reference materials [1].

• Cross-Laboratory Consistency: For multi-site studies, implement standardized protocols for experimental processes including mRNA enrichment, library strandedness, and bioinformatics pipelines, as these have been identified as primary sources of variation [1].

Troubleshooting Common Issues

• Low Correlation Between Replicates: Check RNA quality, library preparation consistency, and potential sample swaps. Examine PCA plots for batch effects and outliers.

• High Multimapping Rates: Consider increasing the --outFilterMultimapNmax parameter or using more stringent filtering. Alternatively, explore tools that probabilistically assign multimapping reads.

• Unexpected Strandedness Patterns: Verify library preparation protocol and select the appropriate column from the ReadsPerGene.out.tab file. For stranded protocols, one of the stranded columns (3 or 4) should contain the majority of reads.

• Memory Limitations: For resource-constrained environments, consider using a subset of the genome annotations or switching to less memory-intensive aligners like HISAT2 for initial exploratory analysis [62].

The selection of an appropriate quantification tool for RNA-seq analysis requires careful consideration of research questions, experimental design, and computational resources. STAR with quantMode GeneCounts provides a robust, accurate solution for gene-level quantification, particularly suited for studies requiring high sensitivity for detecting subtle differential expression, such as clinical biomarker discovery and pharmacogenomics. The decision framework presented here enables researchers to systematically evaluate trade-offs between accuracy, computational demand, and analytical needs, while the detailed experimental protocol ensures reproducible, high-quality results. As RNA-seq continues to transition into clinical diagnostics, adherence to best practices in both experimental execution and computational analysis becomes increasingly critical for generating reliable, actionable insights from transcriptomic data.

Conclusion

STAR's quantMode GeneCounts provides a highly convenient and accurate method for generating gene-level counts directly during alignment, producing results that align closely with htseq-count. Success hinges on correctly selecting the strandedness column based on your library preparation protocol and understanding the inherent trade-offs compared to more complex quantification tools. For standard differential gene expression analysis, it is a robust and efficient choice. Future directions involve tighter integration with transcript-level quantifiers for isoform-resolution studies and continued optimization for cloud-based, high-throughput processing of large-scale clinical RNA-seq datasets in translational research.

References

• 1. A real-world multi-center RNA-seq benchmarking study ... [https://www.nature.com/articles/s41467-024-50420-y]
• 2. a repository of uniformly processed RNA sequencing data - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC6446219/]
• 3. How to perform strand-aware (RNA-Seq) mapping using ... [https://www.biostars.org/p/407788/]
• 4. RNA-seq: Explain STAR quantMode geneCounts values [https://www.biostars.org/p/218995/]
• 5. Mapping with STAR [https://biocorecrg.github.io/RNAseq_course_2019/alnpractical.html]
• 6. quantMode GeneCounts gives different counts with HTSeq ... [https://github.com/alexdobin/STAR/issues/974]
• 7. Performance Analysis and Optimization of the STAR ... [https://arxiv.org/html/2506.12611v1]
• 8. Implementing STAR 's GeneCounts in Generals Statistics... Output [https://github.com/MultiQC/MultiQC/issues/992]
• 9. Alignment-based method | Guide to RNA-seq Analysis [https://ycl6.gitbook.io/guide-to-rna-seq-analysis/raw-read-processing/mapping/alignment-based-method]
• 10. --quantmode GeneCounts multi-mapped and strand ... [https://groups.google.com/g/rna-star/c/mkooNLzyJYc]
• 11. Unassigned Multimapping in featurecounts - reference-annotation... [https://help.galaxyproject.org/t/unassigned-multimapping-in-featurecounts/2399]
• 12. STAR直接就可以输出readsCount，为什么还需要 ... [https://blog.csdn.net/qazplm12_3/article/details/120359482]
• 13. STAR直接就可以输出readsCount，为什么还需要 ... [https://cloud.tencent.com/developer/article/1936443]
• 14. STAR直接就可以输出readsCount，为什么还需要featurecounts？ [https://jishu.proginn.com/doc/57976476928ece6c9]
• 15. Counting reads in features with htseq-count [https://htseq.readthedocs.io/en/release_0.11.1/count.html]
• 16. counting reads within features — HTSeq 2.0.3 documentation [https://htseq.readthedocs.io/en/master/htseqcount.html]
• 17. STAR and HTSeq: unique alignment inconsistency [https://groups.google.com/g/rna-star/c/A2BdNlcn9Us]
• 18. STAR --quantMode GeneCounts vs HTSeq-count [https://www.biostars.org/p/484626/]
• 19. inferring abundance of overlapping genes using RNA-seq data [https://pmc.ncbi.nlm.nih.gov/articles/PMC4331702/]
• 20. MGcount: a total RNA-seq quantification tool to address multi ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04544-3]
• 21. GFF/GTF File Format [http://www.ensembl.org/info/website/upload/gff.html]
• 22. Annotating Genomes with GFF3 or GTF files - NCBI - NIH [https://www.ncbi.nlm.nih.gov/genbank/genomes_gff]
• 23. Ensembl gene annotation gtf for rat problem with RNA STAR [https://help.galaxyproject.org/t/ensembl-gene-annotation-gtf-for-rat-problem-with-rna-star/14818]
• 24. Counting RNA-seq reads for more features than GTF files ... [https://www.biostars.org/p/9570587/]
• 25. Transcriptomics / RNA-seq Alignment with STAR / Hands-on [https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/rna-seq-bash-star-align/tutorial.html]
• 26. new release: counting reads per gene [https://groups.google.com/g/rna-star/c/gZRJx3ElRNo]
• 27. --quantMode GeneCounts options? [https://groups.google.com/g/rna-star/c/C7FpTyCJKLI]
• 28. Stranded vs Non-Stranded RNA-Seq [https://blog.genewiz.com/stranded-vs-non-stranded-rna-seq]
• 29. Comparison of stranded and non-stranded RNA-seq ... [https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-015-1876-7]
• 30. A comprehensive workflow for optimizing RNA-seq data analysis [https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-024-10414-y]
• 31. Data analysis pipeline for RNA-seq experiments [https://pmc.ncbi.nlm.nih.gov/articles/PMC6373869/]
• 32. RNA-seq workflow - gene-level exploratory analysis and ... [https://www.bioconductor.org/help/course-materials/2016/CSAMA/lab-3-rnaseq/rnaseq_gene_CSAMA2016.html]
• 33. rnaseq: Usage [https://nf-co.re/rnaseq/latest/docs/usage/differential_expression_analysis/theory]
• 34. From Gene Counts to Differential Expression - DESeq2 ... [https://ashleyschwartz.com/posts/2023/05/deseq2-tutorial]
• 35. Using DESeq2 and edgeR in R [https://gtpb.github.io/ADER18F/pages/tutorial1.html]
• 36. Differential expression analysis [https://biocorecrg.github.io/RNAseq_course_2019/differential_expression.html]
• 37. -quantMode GeneCounts output big differences between ... [https://groups.google.com/g/rna-star/c/qFLxINhF7cI]
• 38. Strandedness issues in STAR quantification [https://www.biostars.org/p/9613715/]
• 39. Alignment with STAR | Introduction to RNA-Seq using high ... [https://hbctraining.github.io/Intro-to-rnaseq-hpc-O2/lessons/03_alignment.html]
• 40. STAR: ultrafast universal RNA-seq aligner - PubMed Central [https://pmc.ncbi.nlm.nih.gov/articles/PMC3530905/]
• 41. featureCounts or htseq-count? - CVR Bioinformatics [https://bioinformatics.cvr.ac.uk/featurecounts-or-htseq-count/]
• 42. HTSeq-count vs. FeatureCounts: Significant Difference in ... [https://www.biostars.org/p/9604786/]
• 43. --sjdbGTFtagExonParentGene [https://groups.google.com/g/rna-star/c/M0q8M5FscA4]
• 44. Option for exon Parent feature type? · Issue #387 [https://github.com/alexdobin/STAR/issues/387]
• 45. a practical guide to RNA sequencing data analysis [https://pmc.ncbi.nlm.nih.gov/articles/PMC12599993/]
• 46. Mapping RNA-seq Reads with STAR - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC4631051/]
• 47. Benchmarking RNA-Seq quantification tools - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC5003039/]
• 48. A comparison of transcriptome analysis methods with ... [https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-022-08465-0]
• 49. Evaluation and comparison of computational tools for RNA ... [https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-017-4002-1]
• 50. STAR quantMode vs RSEM vs Kallisto [https://bioinformatics.stackexchange.com/questions/11470/star-quantmode-vs-rsem-vs-kallisto]
• 51. Difference between read counts STAR vs featureCounts [https://www.biostars.org/p/401442/]
• 52. Systematic comparison and assessment of RNA-seq ... [https://www.nature.com/articles/s41598-020-76881-x]
• 53. Could you explain the difference between STAR, ... [https://www.biostars.org/p/400121/]
• 54. Introduction to RNA-Seq using high-performance computing [https://hbctraining.github.io/Intro-to-rnaseq-hpc-salmon/lessons/04_quasi_alignment_salmon.html]
• 55. Kallisto vs. STAR: Alignment and Quantification of Bulk ... [https://www.elucidata.io/blog/kallisto-vs-star-alignment-and-quantification-of-bulk-rna-seq-data]
• 56. Comparative analysis of common alignment tools for single ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8848315/]
• 57. About - Pachter Lab [https://pachterlab.github.io/kallisto/about]
• 58. Manual - Pachter Lab [https://pachterlab.github.io/kallisto/manual]
• 59. The mapping-based mode of Salmon runs in two phases [https://salmon.readthedocs.io/en/latest/salmon.html]
• 60. Salmon: fast and bias-aware quantification of transcript ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5600148/]
• 61. Pseudoalignment with Kallisto - Protocols [https://protocols.hostmicrobe.org/compute/kallisto]
• 62. A Guide to Basic RNA Sequencing Data Processing and ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12067304/]