featureCounts: An efficient general-purpose program for assigning sequence reads to genomic features

Yang Liao; Gordon K Smyth; Wei Shi

Source-linked AI summary

featureCounts: An efficient general-purpose program for assigning sequence reads to genomic features

Yang Liao, Gordon K Smyth, Wei Shi

arXiv:1305.3347v2May 15, 2013q-bio.GNq-bio.QM

TL;DR

Counting reads assigned to genomic features is central to many sequencing analyses, yet general-purpose summarization tools have been limited. The paper introduces featureCounts, which supports diverse RNA and DNA sequencing applications and is reported to be an order of magnitude faster and far more memory-efficient than other methods while maintaining high accuracy concordance.

Problem
Read counting is required for many genomic analyses, but only a handful of general-purpose tools were available for summarizing reads across genomic features.
Method
featureCounts assigns RNA- or DNA-sequencing reads to genomic features using optimized feature-search structures, with support for multithreading and configurable summarization options.
Results
FeatureCounts showed high concordance in summarization accuracy while running an order of magnitude faster on average and using substantially less memory than other methods.
Takeaways & Limitations
FeatureCounts provides a general-purpose, computationally efficient option for summarizing reads in large sequencing studies.

Abstract

from arXiv · show

Next-generation sequencing technologies generate millions of short sequence reads, which are usually aligned to a reference genome. In many applications, the key information required for downstream analysis is the number of reads mapping to each genomic feature, for example to each exon or each gene. The process of counting reads is called read summarization. Read summarization is required for a great variety of genomic analyses but has so far received relatively little attention in the literature. We present featureCounts, a read summarization program suitable for counting reads generated from either RNA or genomic DNA sequencing experiments. featureCounts implements highly efficient chromosome hashing and feature blocking techniques. It is considerably faster than existing methods (by an order of magnitude for gene-level summarization) and requires far less computer memory. It works with either single or paired-end reads and provides a wide range of options appropriate for different sequencing applications. featureCounts is available under GNU General Public License as part of the Subread (http://subread.sourceforge.net) or Rsubread (http://www.bioconductor.org) software packages.

1 Introduction

Read summarization counts reads overlapping genomic features, enabling biological interpretation and count-based analyses across sequencing applications. The paper introduces featureCounts as an efficient, versatile program designed to handle diverse sequencing data and feature types.

Motivation: Read counting summarizes mapped reads over exons, genes, promoter regions, gene bodies, and other genomic intervals for downstream statistical analyses.These counts are required by many count-based methods for differential expression.
Challenges: Read-counting tools must support DNA and RNA sequencing, single- and paired-end reads, and reads containing insertions, deletions, or fusions.These requirements make practical read counting more subtle than it may initially appear.
Interpretation: Gene-level RNA-seq counts summarize gene expression but do not distinguish isoforms when multiple transcripts arise from the same gene.Estimating individual isoform expression is therefore a separate challenge from assigning reads to genes.
Related work: Existing read-counting software was limited despite the importance of read counting, motivating a general-purpose solution.The introduction notes that only a handful of general-purpose read-counting tools were available.
Contribution: featureCounts quantifies RNA- or DNA-sequencing reads for any genomic feature using chromosome hashing, feature blocking, other efficiency strategies, and multithreading.It is available as a Unix command or a function.

2 Data formats and inputs

featureCounts accepts aligned reads in SAM or BAM format together with genomic features in GFF or SAF format, with automatic read-format detection. It supports single or paired reads, strand-specific counting, mapping-quality filtering, and summarization at feature or meta-feature levels.

2.1 Input data: Inputs comprise one or more aligned-read files in SAM or BAM format and a genomic-feature list in GFF or SAF format.SAM or BAM format is detected automatically, and the alignment and annotation should correspond to the same reference.
2.1 Input data: SAF provides a simpler annotation format with five required columns: feature identifier, chromosome name, start position, end position, and strand.These columns provide the minimal sufficient information for read quantification.
2.1 Input data: Feature identifiers are assumed to be unique in both GFF and SAF annotations.This follows the convention used in the commonly used GTF refinement of GFF.
2.1 Input data: Reference-sequence counts vary with the application, from chromosome-scale references to much larger sets of contigs or transcriptome sequences.Incomplete or low-quality assemblies can make each contig a reference sequence, while RNA-seq may use transcriptome alignments.
2.1 Input data: Strand-specific counting is supported when strand information is provided, and users can require a minimum mapping-quality score for assigned reads.Mapping results usually include mapping-quality scores for mapped reads.
2.2 Single and paired-end reads: Paired reads define a fragment, so featureCounts counts fragments rather than reads when paired-end data are used.Pairs are assumed to occupy consecutive SAM or BAM positions by default; name sorting can enforce this arrangement.
2.3 Features and meta-features: FeatureCounts summarizes reads at feature or meta-feature levels, where features are intervals and meta-features group biologically related features such as exons within genes.Features sharing a feature identifier in GFF or SAF are assigned to the same meta-feature.

3 Algorithm

featureCounts assigns reads or fragments to genomic features by comparing their mapped bases, including gaps, with feature regions and calling any overlap of at least 1 bp a hit. It accelerates assignment by hashing reference sequences and narrowing searches through a bins–blocks–features hierarchy.

Read assignment: featureCounts calls a hit when any read or fragment overlaps a feature by at least 1 bp, accounting for insertions, deletions, exon junctions, and fusions.Assignment compares the mapping location of every base with the genomic region spanned by each feature.
Read assignment: At the meta-feature level, a hit is assigned when a read or fragment overlaps any component feature, and overlaps within the same meta-feature count once.Thus, an exon-spanning read is counted once for its gene even when it overlaps multiple exons.
Read assignment: Users can exclude multi-overlap reads or count them for every overlapped feature, with the choice depending on the experiment type.The paper recommends excluding reads or fragments overlapping more than one gene in RNA-seq experiments.
Data structure: A hash table rapidly matches reference sequence names in SAM files and GFF annotations, after which reads and features are analyzed separately for each reference sequence.This is particularly useful when the reference contains many sequences.
Data structure: The algorithm sorts features by start position, divides each reference sequence into nonoverlapping 128kb bins, and groups consecutive features within each bin into blocks.The number of blocks per bin is the square-root of that bin’s feature count, rounded up, making block and bin counts nearly equal.
Data structure: The hierarchical bins–blocks–features structure speeds assignment by narrowing candidate regions before comparing the query read with individual features.Queries proceed from genomic bins to blocks in overlapping bins and then to features in overlapping blocks.

4 Implementation

featureCounts is implemented in C with memory-efficient processing and configurable multithreading. The same compiled code is exposed through an R wrapper that integrates directly with R-based statistical analysis workflows.

Implementation: The Unix featureCounts command is written entirely in C and minimizes memory use by retaining only stage-specific feature annotation data.Its C implementation also supports user-specified multithreading, with one thread as the default.
Implementation: The Rsubread featureCounts function wraps the same compiled C code, preserving its efficiency while providing an R programming interface.It returns an R data object that can be passed directly to analysis software such as edgeR or limma.

5 Performance on RNA-seq data

On SEQC RNA-seq data, featureCounts matched summarizeOverlaps for single-end gene counts and counted more paired-end fragments than competing methods. It was also considerably faster and more memory efficient, while recovering additional fragments that were likely genuinely expressed.

Dataset and setup: The evaluation used 6.8 million pairs of 101bp reads from the SEQC project, aligned to GRCh37 and annotated with 25,702 genes and 225,071 exons.Counts were summarized at the gene level using RefSeq annotation build 37.2.
Single-end reads: featureCounts and summarizeOverlaps yielded identical counts for every gene when evaluated on single-end reads.The comparison used only the first read from each paired-end read pair.
Paired-end fragments: summarizeOverlaps counted far fewer paired-end fragments because it required both ends to map successfully, unlike featureCounts and htseq-count.Fragments with only one mapped end can still have relatively high mapping confidence at 101bp read length.
Paired-end fragments: featureCounts assigned some fragments that htseq-count treated as ambiguous by considering how many reads from each fragment overlapped each gene.More than 86% of these additional fragments belonged to genes already supported by at least 100 unambiguous fragments assigned by both methods.
Computational performance: featureCounts was >10folds faster than the other programs and used less memory; summarizeOverlaps still used 20 times as much memory when run chromosome by chromosome.All programs were run using a single CPU without multithreading.

6 Performance on ChIP-seq data

On a mouse H3K27me3 ChIP-seq dataset, featureCounts summarized paired-end fragments over broad gene regions while allowing multi-overlap assignments. It matched countOverlaps exactly and substantially outperformed competing tools in speed, memory use, and fragment counts.

Dataset and setup: The evaluation used 15 million pairs of 35bp reads from mouse mammary stem cells, with fragments summarized over each gene body plus the immediately upstream 3kb region.Reads were mapped to mm9, and only properly paired fragments between 50bp and 500bp were evaluated.
Counting results: featureCounts and countOverlaps produced identical counts for every gene, whereas coverageBED and htseq-count assigned fewer fragments.coverageBED used only the first read of each fragment, while htseq-count presumably omitted multi-overlap fragments.
Counting results: htseq-count counted 7–8% fewer fragments than featureCounts, and changing from union to intersection-nonempty mode did not substantially close the gap.featureCounts was configured to count multi-overlap fragments, consistent with the evaluation’s assignment rule.
Performance comparison: About 5× faster and about 10× less memory than the next most efficient tool, featureCounts achieved the best overall efficiency.Running countOverlaps chromosome by chromosome reduced its peak memory, but it remained more than 100× higher than featureCounts.

7 Performance when the number of reference sequences is large

featureCounts was evaluated on RNA-seq reads simulated from the budgerigar genome, whose annotation comprised 16,204 genes across 2,850 scaffolds. It was fastest and most memory-efficient, while matching summarizeOverlaps’ counts and retaining a speed advantage.

Experimental dataset: 16,204 genes were annotated across 2,850 scaffolds in the budgerigar genome assembly used for the large-reference-sequence evaluation.The assembly was generated in the Assemblathon 2 project.
Performance comparison: featureCounts was fastest and used the least memory among the compared gene-level summarization methods.The comparison included featureCounts, summarizeOverlaps, and htseq-count, with htseq-count run in union mode.
Counting accuracy: featureCounts and summarizeOverlaps produced the same counts, whereas htseq-count produced slightly fewer counts.Counts were compared after summarizing the simulated reads at the gene level.
Performance comparison: featureCounts increased its speed advantage over summarizeOverlaps in this large-reference-sequence evaluation.The program maintained its efficiency advantage over the other methods.

8 Theoretical analysis of algorithmic complexity

The theoretical analysis compares the time and memory growth of featureCounts with other read-summarization algorithms as the number of reads and genomic features increase. featureCounts has the lowest time complexity, while htseq-count has the best space complexity and featureCounts remains close behind.

Complexity framework: The analysis derives theoretical time and space-growth expressions for featureCounts and competing algorithms as reads, features, bins, and overlapping within-bin features increase.Theoretical complexity abstracts away hardware, operating-system, and other implementation factors.
Time complexity: featureCounts has the lowest theoretical time complexity among the compared algorithms.Its hierarchical within-bin search is more efficient than coverageBED’s sequential search, while htsed-count’s red-black tree search has higher complexity because log f is typically larger than the square-root of features per bin.
Time and space complexity: coverageBED uses more bins and typically has more within-bin features than featureCounts, increasing its theoretical complexity.The compared quantities are b2 > b1 and k2 > k1.
Space complexity: htseq-count has the best theoretical space complexity, while featureCounts is close behind; BEDTools, countOverlaps, and summarizeOverlaps require more space under the stated dependencies.featureCounts remains close to htseq-count because its bin count b1 is usually small compared with f.
Implementation effects: Practical efficiency gains for featureCounts additionally arise from high-performance C programming and direct memory manipulation.Observed running time and memory usage depend on implementation efficiency as well as inherent algorithmic complexity.

9 Discussion

The discussion presents featureCounts as an accurate read summarization method with substantially lower computational cost than alternatives. It also emphasizes multithreading, configurable overlap handling, and implementation in Subread and Rsubread.

9 Discussion: FeatureCounts showed high concordance with alternative methods in summarization accuracy while being an order of magnitude faster on average.The comparison found a large difference in computational cost despite similar summarization accuracy.
9 Discussion: FeatureCounts supports multithreaded processing, making it useful for summarizing data from large sequencing studies.The reported results used a single thread, while the program also supports multithreaded processing.
9 Discussion: Users can control whether reads overlapping multiple features or meta-features are counted, depending on the experiment type.The discussion distinguishes overlap-counting choices between RNA-seq and genomic DNA sequencing experiments.
9 Discussion: FeatureCounts is implemented in Subread and Rsubread, with an R interface that calls the compiled C program.The Rsubread implementation has the same speed and memory usage as the entirely C-written Subread package.

11 FUNDING

The work was supported by Australian NHMRC funding, Victorian State Government infrastructure support, and Australian Government NHMRC IRIIS.

11 FUNDING: Funding included an NHMRC Project Grant (1023454), an NHMRC Fellowship to GKS, Victorian State Government Operational Infrastructure Support, and Australian Government NHMRC IRIIS.The passage identifies these sources as supporting or enabling the work.