夜雨聆风一、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
参数 | 说明 |
-a | 参考注释 GTF/GFF 文件路径 |
输出计数矩阵文件名 | |
线程数 | |
计数特征类型 | |
聚合用的属性 ID | |
链特异性模式（0/1/2） | |
双端模式 | |
双端均需比对成功 | |
排除不一致配对 | |
计入多比对读段 | |
允许跨特征重叠分配 | |
最小 MAPQ 阈值 | |
最小重叠碱基数 | |
输出每条读段的分配状态 |
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/
参考文献：
