featureCounts实战指南：从BAM到基因表达矩阵的全流程参数详解

做RNA-seq数据分析的人大概率绕不过featureCounts这个工具。它快、稳、内存占用低是当前从BAM文件到基因表达计数矩阵最常用的方案之一。但我在实际带项目、审流程时发现很多人对featureCounts的使用还停留在跑一条命令出来一个文件的层面比对完直接丢进去默认参数跑完拿到counts就开始跑差异分析。等到结果不对劲才回头排查发现是链特异性设错了、双端read被当成单端数了、多比对read的处理完全没考虑——而这些坑featureCounts本身不报错只会悄悄给你一份偏掉的结果。这篇文章我想从一次完整的实战项目出发把从比对到下游分析这条链路拆开讲透为什么比对阶段就要为定量做准备、featureCounts内部到底怎么数read、参数怎么调才不是照抄命令、输出结果怎么解读和换算、怎么无缝衔接到DESeq2这类差异分析工具。内容会尽量贴近实际操作把我踩过的坑和验证过的经验一并写出来。1. 从比对到定量为什么featureCounts只认BAM以及比对阶段就埋下的坑featureCounts本身不比对序列它的输入是已经比对好的BAM/SAM文件。这看起来简单但问题恰恰出在这里很多人以为只要能比对上的read就是好read实际上比对的参数、输出方式、排序方式直接影响featureCounts能不能正确识别和计数。1.1 比对器选型与排序要求featureCounts不挑剔比对器STAR、HISAT2、Subread aligner甚至bowtie2的产物它都能处理。但有一点很关键——输入BAM必须按基因组坐标排序不能是比对完成后的原始顺序。STAR默认输出就是按坐标排序的但HISAT2默认输出是reads的顺序需要额外用name sort再转coordinate sort或者直接用samtools sort处理。我在实际项目里的做法是统一走STAR因为它的两步比对模式对剪接位点的处理更成熟而且默认输出就是coordinate-sorted省一步工序。如果你用了HISAT2建议比对后加一条命令samtools sort - 8 -o sample.sorted.bam sample.bam samtools index sample.sorted.bam排序不只是为了让后续IGV可视化顺手更关键的是featureCounts在计数时需要顺序扫描BAM文件中的比对记录只有coordinate-sorted的输入才能保证正确性。1.2 比对质量过滤在比对端做还是定量端做这是个很容易被忽略的问题。STAR比对会输出很多多比对reads和低质量比对如果不在比对阶段标记就会把选择压力全部推给featureCounts的-Q参数。我的建议是分两层控制第一层在比对阶段STAR用--outFilterMismatchNmax 10和--outFilterMultimapNmax 20这类参数先把明显不合理的比对过滤掉第二层在定量阶段featureCounts的-Q 10过滤掉比对质量低于10的read。为什么是10而不是0因为MAPQ0在STAR里通常代表多比对read或比对到重复区域的read这些read在基因组上有多处位置可以匹配直接纳入计数会把重复序列区域如rRNA簇、转座子的counts虚高。但要小心许多非模式生物的参考基因组注释不完整过度过滤反而损失真实信号。常规转录组分析中-Q 10是一个稳妥的折中值。1.3 比对软剪辑soft-clip与剪接reads的处理RNA-seq数据里有大量跨越剪切位点的reads比对软件会把read的一部分比对到外显子A跨过内含子再把剩余部分比对到外显子B中间这段就变成了剪辑状态。这类reads是转录本定量的核心证据featureCounts天然支持这类跨越内含子的read计数——它依据GTF注释中的exon坐标只要read比对到的区间和外显子区间有重叠就计入该基因。这里有个细节featureCounts对soft-clip的处理默认是忽略被clip掉的部分只看实际比对到基因组上的部分。也就是说如果一个read有10bp被clip掉了剩下的140bp比对到了外显子上它依然会计数只是那10bp不影响重叠判断。这个设计在大多数场景下是合理的我很少去改它。2. featureCounts的核心计数逻辑读懂-s、-p、-M这些参数背后的原理featureCounts的计数逻辑本质上是read-overlap counting——把read的基因组坐标和GTF注释里每个基因的exon坐标做重叠判断如果一个read与某个基因的任意外显子有重叠默认要求至少1bp就把它归属于这个基因。听起来很简单但各种参数会在不同环节改变归属判断的结果。2.1 链特异性-s参数为什么如此重要链条特异性是新手犯错误最多的地方。它有三个取值0非链特异性、1链特异性read与注释链相同、2链特异性read与注释链相反。绝大多数建库试剂盒做的是dUTP法也就是保留第二条链的信息。这种数据下双端read中第一个readread1实际比对到了转录本的相反链上所以正确参数是-s 2。而老式的链特异性建库如Illumina TruSeq Stranded mRNA早期版本是read1与转录本同链对应-s 1。如果你用的是非链特异性试剂盒才选-s 0。我见过一个真实案例某课题组用链特异性试剂盒建库却用-s 0跑了整个流程结果差异基因数量翻了一倍多——因为反义链的reads也会被计入正义链基因大量本来是噪声的信号混了进来。更隐蔽的是很多基因本身存在天然反义转录本natural antisense transcripts这些reads在-s 0下会被错误地计入对应基因造成表达量的系统性偏差。判断自己数据是哪种链方向最直接的办法是比对完后看featureCounts的summary文件如果你设了-s 2正常情况下Assigned比例应该超过70%如果不到50%且Unassigned_Unmapped不高大概率是链方向设反了赶紧换成-s 1跑一遍对比一下——这个动作只需几分钟能避免整个下游分析在错误的基础上白跑。2.2 双端read计数-p和--countReadPairs的配合双端测序中一个片段fragment会测出两条readsread1和read2。如果不用-p参数featureCounts会把这两条reads当作两个独立的计数单位一个片段被数两次——对于表达量高的基因这个double counting效应会直接翻倍而且在基因长度不同时引入非线性误差。正确做法是加-p让featureCounts把来自同一fragment的两条reads合并成一个计数单位。加-p之后如果两条reads都比对到了同一个基因就计1如果只有一条比对上了也计1因为该fragment仍然有转录证据只有当两条reads比对到不同基因时才会进入Ambiguity或其他category。featureCounts新版本里还有一个参数--countReadPairs效果等同于-p的显式声明。跑的时候建议用featureCounts --help确认你用的版本支持哪些写法老版本可能只认-p。2.3 多映射reads与多重基因-M和-O的真实影响默认情况下featureCounts会忽略multi-mapping reads也就是比对到基因组多个位置的read。它的逻辑是我无法可靠判断这个read来自哪个拷贝干脆不计入任何基因。对大多数基因而言这是对的否则重复基因家族看家基因、组蛋白基因的表达量会被严重高估。但在某些场景下忽略多映射read会丢失信息。比如想评估转座子相关基因、核糖体蛋白基因的表达这些基因家族成员序列高度相似几乎不可能被唯一比对。这时可以加-MfeatureCounts会把多映射read分配给所有与之重叠的基因之一实际上是每个基因都尝试分配会累加到每一个可能匹配的基因上。-O参数则允许一个read同时分配给多个重叠的基因overlapping genes。人类基因组里不少基因存在gene overlap比如A基因的外显子嵌在B基因的内含子里。默认情况下这样的read会标记为Ambiguity而不分配给任何基因加-O后会同时计数到所有重叠基因。这两个参数都要谨慎使用因为它们都会让不同基因之间的counts产生关联性影响后续差异分析的独立性假设。我的默认策略是常规差异表达分析不加-M和-O只做特定基因家族如MHC、嗅觉受体表达评估时才加。要手动加上并明确说明否则审稿人会问。2.4 计数单位与meta-feature结构为什么-G和-t要用exonfeatureCounts的-t默认为exon-g默认为gene_id它的含义是把所有exon当作meta-feature的一部分汇总到gene_id这一层。也就是说featureCounts默认输出的是基因水平的表达量。如果你想看转录本异构体transcript isoform水平的表达量需要把-g改成transcript_id。但这里有个重要提示-g transcript_id并不能真正解决异构体定量问题。因为一个多外显子基因的多个异构体共享大量外显子基于overlap的计数法无法区分reads具体来自哪个异构体。真正的异构体定量应该用RSEM、Salmon、Kallisto这类基于转录组序列的算法。featureCounts做转录本水平的quantification结果在异构体层面不可靠。所以我的建议是用featureCounts时老老实实做基因水平定量这是它的绝对优势区域异构体层面的分析交给专门的工具。3. 实操阶段一条命令跑完定量的完整过程与参数调优这一节直接给出我在真实项目中跑featureCounts的完整操作流程以及每一步为什么要这样设参数。3.1 输入文件准备featureCounts需要两个输入BAM文件和GTF注释文件。BAM文件的准备细节在第一章讲过补充一点同一批样本用同一个参考基因组版本和同一种比对参数否则跨样本比较的表达量没有意义。如果一个项目分批次测序建议把所有人的数据合到一起重新比对避免批次间的比对差异干扰下游分析。GTF注释文件的版本选择也很关键。很多人从Ensembl下载GTF只用最新版本其实应该根据参考基因组的版本配套选择。Ensembl的GTF和UCSC的参考基因组序列版本号往往不对应比如Ensembl release 110对应的GRCh38和UCSC的hg38虽然坐标体系基本相同但基因注释的细节有所不同。强烈建议下载Ensembl自己的全基因组fasta文件和对应的同名GTF避免版本错位。3.2 完整的featureCounts命令与参数解析我比较典型的一条命令长这样featureCounts \ -a Homo_sapiens.GRCh38.110.gtf \ -o counts.txt \ -T 12 \ -p \ -s 2 \ -Q 10 \ -g gene_id \ -t exon \ -C \ sample1.sorted.bam sample2.sorted.bam sample3.sorted.bam逐个说下参数参数值含义与我的判断依据-aGTF文件注释来源必须与参考基因组配套-ocounts.txt输出文件前缀会生成counts.txt和counts.txt.summary-T 12线程数多线程能带来近似线性加速但注意I/O可能成为瓶颈-p双端计数按fragment计数避免double counting-s 2链特异性dUTP法建库的标准设置具体看建库试剂盒说明-Q 10MAPQ阈值过滤低质量/多比对reads-g gene_idmeta-feature基因水平汇总-t exonfeature类型只统计落在exon上的reads-C排除嵌合比对避免同一fragment的两条reads比对到不同染色体时被错误处理3.3 多批次样本的批量运行写法实际项目里样本少则十几个多则上百个不可能手动一个个写文件名。我习惯用shell循环或者xargs批量处理for bam in $(cat bamlist.txt); do sample$(basename $bam .sorted.bam) featureCounts \ -a Homo_sapiens.GRCh38.110.gtf \ -o ${sample}.counts.txt \ -T 8 -p -s 2 -Q 10 -g gene_id -t exon \ $bam done如果样本量特别大更推荐把所有bam文件一次性传给featureCounts它内部会对多文件做合并计数输出一个合并的counts矩阵比循环单独跑再合并更省时省力featureCounts -a annotation.gtf -o all_counts.txt -T 16 -p -s 2 *.sorted.bam3.4 运行后的第一件事检查summary文件featureCounts运行完和counts.txt同目录下会生成一份.summary文件这是我最先看的东西。它的格式长这样Status sample1.sorted.bam Assigned 21456789 Unassigned_Unmapped 456789 Unassigned_MultiMapping 1234567 Unassigned_NoFeatures 234567 Unassigned_Ambiguity 345678 Unassigned_MappingQuality 12345 Unassigned_Secondary 0 Unassigned_Nonjunction 0 Unassigned_Duplicate 0 Unassigned_Overlapping_Length 0 Unassigned_Chimera 0不同版本的featureCounts列出的Unassigned类别略有差异但核心的Assigned比例是所有版本都有的。正常的人类转录组数据Assigned比例应该在60%-85%之间。如果低于50%逐个排查Unassigned_Unmapped很高比对环节出问题了或者样品与参考基因组不匹配物种搞错了很常见Unassigned_MultiMapping很高reads过多比对到多个位置可能建库时rRNA去除不彻底或参考基因组里重复序列过多Unassigned_NoFeatures很高reads能比对到基因组但落不到注释的exon区域。可能是注释太旧或者样品里有参考基因组未注释的新转录本也可能是比对到了内含子区域未剪接的pre-mRNA污染Unassigned_Ambiguity很高reads落在多个基因的重叠区域。链特异性设错也会造成这个比例上升我通常要求项目里的Assigned比例稳定在70%以上才继续下游分析如果个别样本低于这个线会优先复查链特异性和比对质量而不是硬着头皮往下走。4. 输出结果解读从count矩阵到TPM/FPKM的换算逻辑featureCounts直接输出的是raw count——每个基因在所有样本里观测到的fragment数量。这个数字本身不适合用来做基因间比较因为不同基因的长度不同长得越长的基因测到read的概率天然越大。4.1 counts.txt文件结构标准的输出文件是tab分隔的文本前6行是注释信息以#开头记录运行的命令和参数从第7行开始是计数矩阵Geneid Chr Start End Strand Length sample1.bam sample2.bam ENSG00000223972 1 11869 14409 14409 123 456 ENSG00000227232 1 14404 29570 - 29570 234 567注意这个Length列它不是基因的全长而是该基因所有外显子合并后的总长度不考虑异构体时的unique exonic length。很多人在换算RPKM时把Length列理解成转录本长度其实featureCounts这里的Length是计算时用的外显子合并长度和实际的转录本全长会有差异。4.2 为什么推荐TPM而不是FPKM/RPKMRNA-seq领域目前达成的共识是跨样本比较用TPMTranscripts Per Million跨基因比较、差异分析用raw count。RPKM/FPKM的问题是它先把read数除以基因长度得到RPK再除以所有基因的RPK总和得到FPKM这个两步归一化会让不同样本的归一化结果受整套转录组组成影响同一样本里高表达基因会压低其他基因的数值。TPM的算法顺序刚好相反先计算每个基因的read per kilobaseRPK再用所有基因的RPK总和作为归一化因子scaling factor这样每个样本的TPM总和恒定约一百万样本间可直接比较。这也是为什么TCGA等大型项目建议用TPM做表达量比较。从featureCounts的counts矩阵换算TPM我一般在R里完成逻辑很简单library(tidyverse) counts - read.delim(counts.txt, row.names 1, comment.char #) gene_length - counts$Length # 去掉长度信息列保留计数矩阵 count_matrix - counts %% select(starts_with(sample)) # 转换为TPM rpk - count_matrix / (gene_length / 1000) tpm - t(t(rpk) / (colSums(rpk) / 1e6)) tpm_df - as.data.frame(tpm) write.csv(tpm_df, tpm_matrix.csv)这段代码的核心逻辑就是把每个基因的count数除以它的外显子长度kb得到RPK再除以每个样本的RPK总和乘以一百万得到TPM。TMM归一化trimmed mean of M-values是edgeR/limma里常用的另一种策略它适合差异分析前做但它不是表达量的绝对度量不建议作为表达量比较的数值来用。4.3 哪些场景用count、哪些用TPM这个分界线必须清晰差异表达分析DESeq2、edgeR、limma使用raw count作为输入这些工具内部有自己的归一化逻辑DESeq2的median-of-ratiosedgeR的TMM不需要你事先转换成TPM主成分分析、聚类热图、相关性分析这类看整体表达谱关系的可视化推荐用TPM比较单个基因的表达高低也用TPM。我最常看到的一个错误就是有人把count矩阵归一化成TPM之后丢进DESeq2导致所有样本的总测序深度信息被抹平差异分析的统计检验失去意义。这个坑在审稿时经常遇到也是新手最容易混淆的地方。5. 下游差异分析的衔接DESeq2需要什么样的干净数据featureCounts的输出可以直接对接各种差异分析工具但这里面还是有不少需要注意的地方。5.1 从counts.txt构建DESeq2输入DESeq2需要的输入是一个raw count矩阵 样本分组信息矩阵行为基因列为样本数值必须是整数raw count不是归一化后的浮点数。直接把featureCounts的counts.txt读进来就能用但必须注意去掉Chr, Start, End, Strand, Length这些注释列基因名作为行名不能有重复样本名要规范建议统一加前缀sampleID_condition这种格式构建DESeq2对象的示范代码library(DESeq2) # 读入count矩阵 count_data - read.delim(counts.txt, row.names 1, comment.char #) %% select(starts_with(sample)) # 构建样本信息表 col_data - data.frame( row.names colnames(count_data), condition factor(c(rep(control, 3), rep(treatment, 3))) ) # 构建DESeq2对象 dds - DESeqDataSetFromMatrix( countData round(count_data), colData col_data, design ~ condition ) # 运行差异分析 dds - DESeq(dds) # 提取差异结果 res - results(dds, contrast c(condition, treatment, control)) res_df - as.data.frame(res) %% rownames_to_column(gene_id) %% filter(padj 0.05, abs(log2FoldChange) 1)这里round()是防患未然因为某些上游工具会输出非整数的count值DESeq2要求输入整数矩阵。5.2 过滤低表达基因在DESeq2之前还是之后DESeq2内部做independent filtering时会自动过滤极低表达基因所以理论上你可以不做额外的过滤。但为了统计学稳健性我通常会在构建DESeqDataSetFromMatrix之前先过滤掉在超过半数样本里count 10的基因keep - rowSums(count_data 10) (ncol(count_data) / 2) count_data_filtered - count_data[keep, ]这个过滤可以帮助去掉那些sample量极低、无法可靠估计离散度的基因减少多重检验校正的负担。但注意不要设太激进的阈值否则会丢掉一些表达量低但生物学上有意义的基因比如转录因子。5.3 从DESeq2回到表达量可视化分析做完之后做热图或者折线图展示差异基因的表达模式时又要回到TPM。一个标准流程是用DESeq2做差异分析拿到差异基因列表用featureCounts的counts矩阵算TPM把差异基因的TPM值做z-score行标准化然后画热图这种做法既利用了DESeq2的统计严谨性也保留了TPM的直观可比性。注意在报告里要写清楚归一化方法。5.4 常见的结果不显著排查思路很多人在这一步会碰到这种情况差异分析跑完一个显著基因都没有或者显著基因列表奇怪到没法解释。我一般从这几个方向排查链特异性反了没有反了的话正义链基因信号和反义链转录本信号纠缠在一起差异会被稀释测序深度够不够一般认为每个样本的total reads至少在2000万以上取决于物种和预算低于这个线低表达基因的统计功效极差biological replicates是不是太少了3 vs 3是底线2 vs 2基本做不出可靠结果分组设计里是不是混入了批次效应如果对照组和处理组分别在不同的上机批次测的批次效应会完全掩盖处理效应即使DESeq2加了~ batch condition也只是亡羊补牢6. 我整理过的一批高频报错与排查经验最后这部分是我在实际项目里收集到的、反复出现的一些问题和解决思路分享出来希望对大家有用。6.1 GTF文件格式问题导致的莫名其妙结果featureCounts对GTF的要求很严格如果GTF文件有问题它要么报错要么给出奇怪的结果。最常见的几种feature类型不匹配GTF第九列里gene_id或transcript_id字段缺失。有的厂商自定义GTF用的是geneId而不是gene_id用-g gene_id会报Attribute gene_id is not found。GTF里exon的ID不唯一同一个gene_id下多个exon有相同的exon_id不算错但如果gene_id重复出现而坐标不连续featureCounts会把它们合并成一个大meta-feature。GTF和BAM的染色体命名不匹配GTF用chr1而BAM用1或者反过来结果就是所有read都匹配不上exonUnassigned_NoFeatures高达100%。这是新手最容易犯的错误之一。解决方案拿到GTF后先less看一眼前几行用grep chr1确认染色体命名方式拿一个单基因的BAM和GTF做小范围验证确认计数结果符合预期再全量跑。6.2 双端read计数异常总量看着对但结果对不上如果跑完发现总counts数接近所有样本read数的一倍多半是没加-p被double counting了。featureCounts加-p之后计数单位从read变成fragment总counts自然会小于总read数。这是最容易自查的一个指标。另外有些情况下--countReadPairs和-p同时用featureCounts会提示参数冲突或直接报错。新版featureCounts2.0以上更建议使用--countReadPairs替代-p但旧版1.x只认-p。建议跑之前先featureCounts -v看版本再决定参数写法。6.3 Unassigned_NoFeatures比例过高从注释和比对两个方向找问题现象是summary文件里Unassigned_NoFeatures占到30%甚至更多Assigned掉到60%以下。优先确认参考基因组和注释是不是同一来源同一版本。我遇到过某个公开数据集的比对是用GRCh38做的但下载的GTF是hg19的导致大量reads落不进任何exon。这种情况不会报错因为部分基因在两个版本间坐标是保守的但整体比例会明显异常。如果注释没问题检查一下BAM里reads的insert size分布。如果insert size特别长超过基因的平均内含子长度说明建库本身质量可能有问题DNA污染或者RNA降解都会导致reads在基因组上的分布偏离预期。6.4 多线程跑的时候系统崩溃资源限制问题-T 16看起来很美但如果你的机器内存只有16GBAM文件又大多线程可能直接OOM。featureCounts的每个线程需要独立的内存缓冲区处理reads样本多、注释复杂时内存占用会指数上涨。我的建议是线程数不要超过CPU核心数的70%单次处理的样本数控制在10个以内大项目分批跑完再合并。跑大矩阵前先用一个小样本测试内存占用再决定加多少线程。6.5 和StringTie/RSEM结果差异大不一定是谁错了有人在项目里同时用了featureCounts和StringTie做定量发现同一样本同一基因的表达量差异可能达到数倍于是怀疑某个工具有bug。实际上这很可能不是bug而是两者的统计口径不同featureCounts基于比对到基因外显子区域的read计数不做转录本区分StringTie是基于组装转录本的模型来估计表达量对低表达基因可能产生更平滑的估计值但过度依赖组装质量RSEM/Salmon这类基于转录组序列定量的工具对多异构体基因的分配方式和featureCounts完全不同不要指望不同工具给出完全一致的数值关键的是看差异分析结果的趋势是否一致。如果同一处理差异的方向和倍数在两个工具下都一致说明结论可靠如果差异方向都对不上那就要检查是不是bias太严重。6.6 链特异性复查的快速方法如果你不确定该用-s 1还是-s 2但手头有比对好的BAM和GTF可以用一个非常快的检验取一条已知高表达基因的reads看它们的flag字段里的strand信息和GTF注释的strand是否一致。更简单的做法是分别用-s 1和-s 2各跑一次同一个样本看哪个Assigned比例高。真实数据里正确方向的Assigned比例通常会明显高于错误方向一般差10个百分点以上。这个方法不严谨但非常高效适合快速确认。经验和心得小结featureCounts本身是个足够傻瓜的工具一条命令就能跑完定量但真正决定结果质量的是转录组流程里那些看似不起眼的选择链特异性、双端处理、比对过滤器、GTF版本、多映射过滤。我个人在跑新项目时已经形成了一个固定习惯每次定量完先花两分钟看summary文件的Assigned比例和各Unassigned分布确认都在合理范围才继续。这个习惯至少帮我避免过两三次整个项目推倒重来级别的灾难。建议新手也把summary文件的检查当成流程里不可省略的一步。最后分享一个小技巧如果你用featureCounts跑完定量想在论文方法部分写得规范可以看它生成的counts.txt前几行注释里面完整保存了运行命令直接在material and methods里引用即可省去手动记录参数的过程。