软件知识分享 | featureCounts:高效的基因组特征读段计数工具

一、featureCounts 的开发背景

      featureCounts 是由澳大利亚 墨尔本沃尔特与伊丽莎·霍尔医学研究所（Walter and Eliza Hall Institute of Medical Research, WEHI） 的 Wei Shi（施伟） 及其生物信息学团队开发的。

      它是 Subread 软件包中的一个核心组件，主要用于高通量测序数据（如 RNA-Seq、ChIP-Seq 等）的定量分析，能够高效、准确地将测序 reads 分配到基因组特征（如基因、外显子）上。

发表：2014 年发表于 Oxford 旗下 Bioinformatics，是 RNA-seq 基因表达定量领域引用最广泛的工具之一。

开发初衷：RNA-seq 定量分析的核心步骤是将比对后的 Reads 分配到基因组特征中生成表达计数。其他工具在大样本、高深度测序场景下，普遍存在运行速度慢、内存占用高、对多线程支持不佳等问题。featureCounts 通过引入染色体分块哈希（chromosome-block hashing）和特征树搜索（feature tree search）算法，将运算效率提升一到两个数量级。

支持：单端（SE）/ 双端（PE）RNA-seq 测序数据

输入文件：已完成基因组比对的 BAM / SAM 文件

二、核心模块

      featureCounts 虽为一款命令行工具，但其功能模块设计层次分明，涵盖输入处理、特征注释过滤、计数算法、多文件联合与输出控制。

1. 输入模块（Input Processing）

BAM/SAM 文件读取：支持单个或多个 BAM 文件同时处理。

格式兼容：原生支持 .bam或.sam 格式，可以同时读取多个文件统一输出。

2. 注释模块（Annotation Engine）

GTF/GFF 特征注释：以常见的 GTF/GFF 格式注释文件定义基因组分区特征。

特征类型与属性 ID：两个参数直接决定计数对象。

#参数：-t / --fType  #指定计数特征类型（feature type），默认 exon-g / --gType  #用于分组聚合的特征属性 ID（attribute type），默认 gene_id

3. 计数算法模块（Counting Algorithm）

读段分配规则：链特异性建库支持，通过 -s 参数设定链特异性策略

#参数：    -s 0 #（默认）非链特异性-s 1 #链特异性（有义链）-s 2 #反向链特异性

单端 vs 双端策略：

单端（SE）：每条 Read 视为独立计数单位双端（PE）：默认每个配对 Fragment 计为 1 次（而非 2 次），通过 -p 启用，结合 -B（两端均比对）、-C（排除不一致比对）提高准确性

多比对读段处理：-M 参数启用多比对读段计数。默认仅计唯一比对（MAPQ ≥ 255 或比对次数 = 1）。启用后，每条多比对读段分配 1/n 的权重（分数计数）。

最小 MAPQ 过滤：-Q 参数设定最小比对质量阈值，默认 0，推荐 RNA-seq 设 -Q 15。

重叠读段计数方法：

-O #（允许重叠特征）默认一个 Read 只分配给一个特征，开启后一个 Read 可分配给多个重叠特征--minOverlap #最小重叠碱基数，默认 1

分数计数（Fractional counting）：启用 -M 和 -O 后，featureCounts 会自动使用分数分配策略，避免多重计数冲突。

4. 多文件联合定量模块（Multi-file Summarization）

featureCounts 可同时读入多个 BAM 文件，生成联合计数矩阵。

#参数：featureCounts -a annotation.gtf -o all_counts.txt \sample1.bam sample2.bam sample3.bam sample4.bam#跨文件特征过滤：--minMQS #为每个文件设定独立的 MAPQ 阈值。

5. 输出控制与报告模块（Output & Reporting）

核心输出：文本计数矩阵文件，第一列到第二列依次为 Geneid、Chr、Start、End、Strand、Length 注释列，后续每一列代表一个样本的计数。

统计摘要输出：featureCounts 运行结束后会在 STDERR 输出详细的比对分配统计摘要，包含总读段数、成功分配读段数及比例、因不同原因未分配的读段数及比例（无特征、模糊匹配、多匹配、短片段等）。

通过 -R 参数可将每条读段的分配状态输出到单独文件，供深度排查。

三、使用示例与实战

使用案例：RNA-seq 基因表达定量

# 单端数据：按基因计算表达矩阵featureCounts -T 8 -t exon -g gene_id \  -a Homo_sapiens.GRCh38.109.gtf \  -o gene_counts.txt \  sample1.bam sample2.bam sample3.bam# 双端链特异性数据（Illumina TruSeq Stranded Total RNA）featureCounts -T 8 -t exon -g gene_id \  -p -s 2 -B -C \  -a Homo_sapiens.GRCh38.109.gtf \  -o stranded_gene_counts.txt \  sample1.bam sample2.bam

1. 计数矩阵（主输出文件）

一个标准的 .txt 表格文件，以 Tab 分隔，布局参考如下：

Geneid  Chr     Start   End     Strand  Length  sample1.bam  sample2.bamENSG00000141510  chr17  7661779 7687538 -       27921   1285          1102ENSG00000142168  chr19  45430116 45432199 +       573     432           398...

第一行包含所有样本文件名作为列名。该格式无需任何转换即可直接读入 R 的 DESeq2、edgeR、limma 包。

2. 运行统计摘要（STDERR 输出）

这些统计摘要可以直接接入 MultiQC 自动化汇总，在批量样本中快速识别偏离阈值的异常样本。

3. 读段分配日志文件（可选，-R C 参数启用）

每条读段的分配状态记录，用于排查个体读段分配失败原因。

使用经验分享：不知道文库链特异性可以先用 RSeQC 推断链类型再给 featureCounts 设置对应参数。

项目地址：

https://subread.sourceforge.net/

参考文献：