企业官网建设流程全解析

开头部分，直接切入。

做RNA-seq数据分析的人，十有八九都绕不开一个环节：转录本定量。从测序仪下机到拿到count矩阵，中间经过质控、比对、定量三步，而featureCounts几乎是目前处理"从比对结果到表达量计数"这一步最顺手的工具之一。它不花哨，不追求花活，但胜在快、稳、内存占用低，而且精度在绝大多数场景下都不输给其他方案。

这篇文章不打算照搬官方文档。我想以一次完整的实操为例，从比对产物的检查开始，一直聊到featureCounts的参数选择、结果解读、常见坑，再到下游差异分析的数据衔接。如果你手头正好有BAM文件需要变成表达矩阵，或者刚入门RNA-seq分析，想在"定量"这一步踩少一点坑，这篇内容应该正好对得上你的需求。

1. 转录本定量：从比对到计数的完整链路

1.1 为什么在"定量"这一环选featureCounts

RNA-seq定量的大方向其实就两类：一类是alignment-free的，比如Salmon、kallisto，它们直接基于转录组序列做伪比对，速度快但依赖转录组索引和注释完整性；另一类是基于alignment的相对传统路线，先让reads比对到基因组，再统计reads落在哪些基因/转录本上。featureCounts属于后者，而且用起来非常干脆，输入是BAM文件加注释GTF/GFF，输出是计数矩阵和统计摘要。

我之所以在大多数项目里坚持用featureCounts，原因很朴素。第一是速度快，它对reads的比对位置做了高效哈希索引，实测单线程情况下也能在几分钟内处理几千万对reads，加多线程之后基本不会成为流程瓶颈。第二是内存占用小，几十G的BAM文件，默认配置下跑完也就占几个G内存，这对没有高配服务器的团队特别友好。第三是输出干净，主输出文件是一个标准的tab分隔计数矩阵，summary文件把成功比对、未比对、多映射、重复等情况单独统计，便于做质控。这三条优点叠加下来，featureCounts在教学、科研和工业流程里都成了默认选项。

不过选择它也有代价。featureCounts本质上是基于基因注释的计数，它没法直接处理注释文件里没收录的新转录本结构。如果你的课题高度依赖新异构体发现，那可能还需要结合转录本拼接的结果做补充；但大多数差异表达分析场景下，基因层面的定量完全够用。这也是为什么我在这篇文章开头先把适用边界说清楚：明确你想要的是"基因表达量"，而不是"异构体级别精确量化"。

1.2 比对环节决定了定量的天花板

很多人把注意力全放在featureCounts参数上，忽略了上游比对质量的拖累。实际上，featureCounts只能统计"已经正确比对到参考基因组上的reads"，比对错位、多映射、嵌合比对，这些都会直接影响最终定量结果。定量这一步再怎么优化，也只能在比对质量给定的前提下做文章。所以，我建议每一次跑featureCounts之前，都先花十几分钟把比对结果检查一遍。

比对工具的选择也很有讲究。目前RNA-seq主流是STAR和HISAT2，两者都倾向于转录组感知的拼接比对，可以跨越内含子。STAR速度极快但内存占用高；HISAT2在表型上更温和一些。如果你的数据是来自核糖体RNA去除后的文库，比对率一般能在80%到90%以上；如果低于70%，别急着跑定量，先回头看看参考基因组版本、注释文件版本、以及样本是否存在污染或降解。引用参考基因组版本混乱是很多项目里最隐蔽的坑——样品处理、比对、定量、差异分析全流程都要统一在同一个版本上，否则后期的基因ID比对很容易对不上。

1.3 核心工作流总览

为了方便后面具体展开，我这里先把完整流程列出一个清单式的路线，后面每个环节都会单独细讲。

1. 原始测序数据（FASTQ）
2. 质控与过滤（FastQC + trimmomatic/fastp）
3. 比对到参考基因组（STAR/HISAT2），输出SAM/BAM
4. 对BAM进行排序、压缩、建立索引（samtools sort/index）
5. 检查比对质量（比对率、插入片段分布、reads覆盖均匀性）
6. 准备注释文件（GTF/GFF，务必与比对所用参考一致）
7. 运行featureCounts进行定量，得到count矩阵
8. 质控summary、评估统计指标
9. 转换成TPM/FPKM，进入DESeq2/edgeR或limma流程做差异分析
10. 可视化（PCA、热图、火山图）

这个流程的核心环节就两个：比对和定量。比对决定你能看到什么，定量决定你能数清楚多少。两者互相制约，缺一不可。

2. 比对结果的质量核查与格式预处理

2.1 SAM/BAM文件到底要怎么看

拿到比对软件生成的BAM文件，别直接丢给featureCounts。我见过不少新手上来就跑featureCounts，结果counts全是一堆0，回头一查，是BAM没排序，或者注释和比对版本不一致。先花几分钟检查BAM文件的基本信息，能省去后面一大截排查时间。

常用命令是samtools和samtools stats。比如我想快速看一眼BAM头信息和比对统计，可以这样：

# 查看BAM头中参考序列信息 samtools view -H sample.bam | head -50 # 生成详细的比对统计报告，包括总reads、比对率、重复率、插入片段分布等 samtools stats sample.bam > sample.bam.stats

在samtools stats的输出里，我最关注几个关键字段：

• raw total sequences：BAM中总reads数，双端数据会每个mate各计一次，所以总数是read pair的两倍；
• reads mapped：成功比对的reads数量；
• reads mapped and paired：双端中配对且正确比对的reads数，这个比例越高越好；
• insert size average、insert size standard deviation：插入片段长度均值和标准差，能帮助判断文库质量；
• reads duplicate：重复reads的比例，高重复率通常意味着文库复杂度低，定量的时候要警惕。

如果双端reads的insert size均值明显偏离文库目标（例如目标300bp，实际却达到800bp），要考虑是不是样本降解严重，或者比对时参数没配对好方向。这一阶段发现问题，比跑到定量完才发现结果没法解释要划算得多。

2.2 侧记：mummer的序列比对结果怎么看

这里不得不提一个经常和RNA-seq比对混淆的话题。有时候项目里并不只是做转录组定量，还会顺带做基因组层面的共线性分析或者变异检测，这时候会用MUMmer这类全基因组比对工具。MUMmer的输出结果格式和STAR/HISAT2的SAM/BAM完全不同，很多第一次接触的人会对着delta文件发懵。

MUMmer输出的delta文件，每一行记录的是两条序列之间的一个局部比对的比对位置、长度和错配信息。简单说，它是把两个基因组之间所有相似片段以"块"的形式列出来。理解delta文件的关键不在于逐行读完，而是用delta-filter、show-coords、show-aligns这些工具去提取有用信息。比如show-coords能把delta转成类似表格格式，给出两个基因组的比对起点、终点、覆盖率和相似度。如果你关心两个基因组之间某个基因区域是否保守，就在show-coords的输出里按位置区间去过滤，再结合dot plot图看整体的共线性趋势。

我在这里提这个，是因为在实际项目中，我遇到过有人把BAM比对结果和MUMmer比对结果混为一谈。它们是两种完全不同的东西：前者是短测序reads和单个参考基因组的比对，用来定位每一条read；后者是长序列或者整个基因组与另一个基因组的比对，用来研究结构变异和进化关系。弄清楚它们的差异，能帮你更快定位问题出在哪一环。

2.3 排序、压缩、索引的必要性

featureCounts本身不要求BAM一定按坐标排序，它甚至可以直接处理未排序的SAM文件，但实际项目中我强烈建议统一用samtools sort处理一遍。原因有三个。一是排序后的BAM在后续其他分析里几乎都通用，比如IGV可视化、call variant、提取bamCoverage信号图，全部要求坐标排序。二是featureCounts在多线程模式下处理排序后的BAM会更稳定，磁盘IO顺序读也行。三是如果你需要多次运行不同参数的定量，就没必要反复重跑排序步骤。

一个标准的预处理命令链大致是这样：

# SAM转BAM并排序 samtools view -bS sample.sam | samtools sort -@ 8 -m 4G -o sample.sorted.bam # 建立索引（部分工具需要，featureCounts不强制但建议） samtools index -@ 8 sample.sorted.bam

需要注意，samtools sort的-m参数控制的是每个线程的最大内存用量，不是总内存。如果机器只有32G内存，-@ 8 -m 4G意味着最多可能吃掉32G内存，容易导致OOM。实际生产里我更习惯给每个线程2G，比如-@ 12 -m 2G，刚好压住内存上限，速度也不慢。

有另一个小技巧，就是给BAM添加RG标签，也就是read group信息。如果后续要用GATK流程处理变异，或者做一些需要合并多个样本的定量比较，RG标签是必须的。即便只是跑featureCounts，建议在比对时顺手加上RG标签，免得下游要重新处理一遍。

3. featureCounts实操：参数详解与命令模板

3.1 安装与版本选择

featureCounts是Subread软件包的一部分，安装方式比较灵活。最省事的办法是用conda：

conda install -c bioconda subread

或者去SourceForge下载源码自行编译。我建议优先使用conda，原因有二：一是依赖关系处理得干净，不会出现本地库冲突；二是版本管理方便，做项目复现的时候能锁定版本。featureCounts版本号看起来影响不大，但在不同版本之间，默认参数和行为有细微变化，比如某些版本对链特异性参数的解释做了调整。如果你长期维护同一套流程，最好把subread版本固定在某个已知稳定版本，写进环境配置文件里。

用conda安装完成后，可以直接验证版本：

featureCounts -v

看到类似featureCounts v2.0.6的输出，就说明环境没问题了。如果你用的是服务器且没有root权限，conda基本是最顺滑的选择，因为可以安装到个人目录下。

3.2 GTF/GFF注释文件的准备

注释文件的质量直接决定定量的准确性。featureCounts支持GTF、GFF以及SAF三种格式，其中GTF和GFF是直接从Ensembl、UCSC、GENCODE下载的标准格式，SAF格式则是一种极简的四列格式：GeneID、Chr、Start、End。日常使用中，我绝大多数时候直接传GTF文件给-a参数，不额外转SAF，除非我需要自定义定量区间。

这里有个关键点：注释文件必须和比对时的参考基因组版本配套。比如你用Ensembl的GRCh38参考基因组跑的STAR，那么注释也建议用Ensembl对应版本的GTF，不要混用UCSC的注释。染色体命名格式不一致是常见的坑，比如Ensembl用的是chr1、chr2，而某些版本可能直接用1、2，特征计数阶段一个对不上，全部counts为0也不奇怪。

下载GTF时，我还习惯做一步"过滤性压缩"：只用exon行来做定量，因为featureCounts默认-t exon -g gene_id，这样统计的是每个基因的所有外显子区域。GTF里有很多其他feature类型，比如CDS、utr、gene、transcript，如果不在-t里指定，软件不会用到。这也是featureCounts默认行为很对人胃口的地方，你不用自己手动提取外显子区间，它自己会处理多外显子基因的合并问题。

3.3 核心命令模板与实际参数选择

下面这是我在双端RNA-seq项目里最常用的模板之一：

featureCounts \ -T 8 \ -a annotation.gtf \ -o counts.txt \ -t exon \ -g gene_id \ -s 2 \ -p --countReadPairs \ -Q 10 \ -B \ -C \ sample1.sorted.bam sample2.sorted.bam sample3.sorted.bam

逐项说明一下这些参数的作用。

-T 8是设置线程数，这里8个线程适合普通的12核服务器，如果你的机器核数更多，可以往上调，但不要超过物理核心数，否则反而会拖慢速度。

-a annotation.gtf是注释文件路径，注意要用绝对路径，尤其跑大规模流程时，不要依赖相对路径去猜。

-t exon表示统计exon这个feature类型，如果注释文件里还想统计其他类型，可以在这里另写，但一般分析都不用动。

-g gene_id是告诉featureCounts用GTF的哪个attribute字段作为基因标识符。Ensembl GTF里这一列通常是gene_id "ENSG00000000001"，所以默认就是按基因ID计数。如果你想要转录本水平定量，可以把-g改成transcript_id，同时-t保持exon，这样计数单位就是转录本，需要注意reads在多转录本共享外显子时的哈希分配逻辑。

-s 2是链特异性参数，这个值的选择要依据文库制备方式而定。很多商业化的链特异性文库（比如dUTP）要用-s 2，即反向链。普通的非链特异性文库则用-s 0。搞错链特异性会带来严重问题：reads被计入反义链或错误链，最终差异分析完全无法解释。我经常见到新手把-s 2套用在非链特异性数据上，结果一半以上的genes出现奇怪的表达模式。

-p --countReadPairs是告诉featureCounts输入是paired-end数据，并且按read pair（即片段）计数，而不是按单条read计数。如果这里是单端数据，保留-p会直接报错。双端数据单位是fragment，单端是read，这个区别在后期的库大小归一化时影响不大，但计数数值本身差了一倍。

-Q 10是设置最低比对质量阈值，低于这个质量值的reads会被丢弃。默认值是0，但实际建议至少给到10，可以把低质量比对过滤掉。对于质量极低的数据，我会给到20，不过要注意可能误伤部分真实但质量偏低的reads。

-B表示只统计成对且都正确比对的reads，-C表示丢弃那些比对位置冲突的reads（比如两个mate比对到不同染色体上的情况）。这两个参数在双端数据分析时加上，能让定量结果更干净，但也意味着比对率报告会略低一些，这个需要在看summary时心里有数。

命令执行后，主要输出counts.txt和counts.txt.summary。counts.txt第一列是Geneid，之后每个样本占一列，数值是该基因上比对上的fragment数。summary文件则统计了每个样本的比对总数、成功计数reads、没有特征的reads、多映射reads、无法比对reads等分类，是所有后续质控的第一步。

3.4 多线程与内存策略

featureCounts的多线程效率在常见RNA-seq工具里属于中等偏上。实测下来，8线程和16线程对运行时间的改善不是线性的，因为线程增多后，内存带宽和IO会成为瓶颈。处理几千万对reads的数据，8线程通常就能在几分钟内完成，如果数据量特别大，比如全转录组（long RNA）+多重复样本，可以开到16线程，但内存可能从3G涨到6G左右，注意服务器内存余量。

一个比较务实的方法是先在单个样本上跑通整条命令，记录耗时和内存峰值，再决定是否批量并行多个样本。比如服务器有32核，与其跑featureCounts -T 32处理一个样本，不如拆成4个管道同时各跑-T 8，这样每个样本都能更快产出，整体效率更高。不过这要求你的磁盘IO能抗住并发读写，否则反而会互相拖慢。

还有一个小提醒：featureCounts默认会把临时文件写到临时目录，如果服务器/tmp空间不足，可能导致运行失败。可以在跑之前先export TMPDIR=/path/to/tmp，或者手动指定一个空间充足的工作目录。这种问题在大型样本上很少发生，但在海量小文件或者磁盘配额紧张的环境里很常见。

4. 定量结果的解读与归一化

4.1 输出文件有哪些东西

featureCounts跑完后，工作目录下一般会出现两个文件：counts.txt和counts.txt.summary。看名字很简单，但实际上counts.txt里包含好几段信息。打开文件，前面21行是注释信息行，每行以#开头，记录的是命令行、版本、输入文件等。这些行在导入R时要注意跳过，否则会报错。

从第22行开始是表格正文，前六列分别是Geneid、Chr、Start、End、Strand、Length。最后一列Length是基因的外显子合并总长度，这个值在计算TPM/FPKM时非常关键。之后的每列代表每个样本的count数。这里有个容易踩坑的点：同一GTF下，一个基因可能在多个位置有重叠区间，featureCounts在输出中会留一行记录它的合并区间，Chr、Start、End这列是合并后的坐标，不是单一外显子坐标。

counts.txt.summary则是一张容易被人忽略的质控表。它把每个样本的reads分成几个类别：Total、Assigned、Unassigned_Unmapped、Unassigned_Secondary、Unassigned_MappingQuality、Unassigned_NoFeatures、Unassigned_Overlapping_Length、Unassigned_Ambiguity等。我最关注的是Assigned占比，合格样本一般要大于60%，质量好的可以到80%以上。如果Assigned比例很低，比如只有40%，那就必须排查是注释不匹配还是比对质量问题。

4.2 定量质量怎么看才靠谱

经验不足的人往往只看Total和Assigned两个数字，但这不够。我总结了一套快速判断标准：

• Unassigned_NoFeatures比例高：说明大部分reads比对到了注释文件没有记录的区域，可能是注释版本太旧，或者比对到线粒体、rRNA区域的比例偏高；
• Unassigned_Ambiguity比例高：说明reads落在多个基因重叠区域无法唯一判定，常见于基因密集区域或重复区域，这个比例一般不会太高，如果高到20%以上，要考虑是否是链特异性参数设错导致reads同时落在正负链基因上；
• Unassigned_MappingQuality比例高：说明大量reads比对质量低于-Q阈值，这种情况常见于参考基因组污染或者样品来源物种不匹配。

拿到BAM文件后，还可以用featureCounts自带的-v或者RSeQC工具来做更细的可视化质控。但不管用哪个工具，核心就一句话：Assigned比例并非越高越好，你需要看的是被丢掉的reads各自去了哪一类，结合文库类型和物种来判断是否合理。

4.3 从count矩阵到TPM/FPKM转换

featureCounts输出的是原始count数。这个数值受基因长度、测序深度、文库大小影响，直接拿去比较样本间表达量会失真。最常见做法是先转成TPM，再进入各差异分析工具。TPM的转换公式比较直白：TPM = (reads_per_kb / sum(reads_per_kb)) * 1e6，其中reads_per_kb = count / (gene_length_in_bp / 1000)。这个公式的核心含义是把基因长度归一化后，再按所有基因的总转录本量做比例归一化。

我自己写过一个简单的R函数来做转换：

counts_to_tpm <- function(counts, feature_length) { rate <- counts / (feature_length / 1000) tpm <- t( t(rate) / colSums(rate) ) * 1e6 return(tpm) } counts <- read.table("counts.txt", header=TRUE, row.names=1, skip=1) tpm <- counts_to_tpm(counts[, 7:ncol(counts)], counts$Length)

这里有个新手的常见误区：FPKM和TPM的语义差别。FPKM是片段每千碱基每百万reads，它的分母用的是总reads，而TPM分母用的是归一化后的转录本总量。TPM的核心优点是在不同样本间，所有基因的TPM总和相同，更适合样本间比较。现在的差异分析工具如DESeq2、edgeR本身并不需要TPM，它们吃原始count，做内部的库大小归一化；TPM更多用于表达量展示和富集分析前的输入。所以流程上千万别混淆：差异分析用count，展示/比较绝对表达量用TPM。

5. 常见问题与排查技巧

5.1 为什么Assigned比例那么低

这是featureCounts使用中最高频的问题。Assigned比例低，首先要看summary表格里的Unassigned分类，再对应排查。

如果是Unassigned_NoFeatures偏高，十有八九是注释不匹配。比如参考基因组和GTF版本不配套。有一个快速排查法：随机取几条Assigned为0的高表达基因（通常是核糖体蛋白基因、组蛋白基因），在IGV里打开BAM文件和GTF，看reads到底落在哪里。如果reads密集覆盖在基因外显子区但GTF完全没显示，那很可能就是GTF版本错误。如果reads覆盖在基因附近但featureCounts没有统计，那可能是链特异性参数反了。

我还遇到过一种情况是BAM文件里包含了多个染色体上的reads，但GTF只覆盖了主要染色体，比如线粒体、decoy contig没有注释条目。此时Unassigned_NoFeatures比例会虚高，解决办法是在比对之前就把参考序列过滤成主要染色体，比如chr1到chr22、chrX、chrY、chrM。

5.2 链特异性参数到底选几

链特异性是新手最容易犯错的参数。RNA-seq文库制备现在几乎都带链信息，但不同试剂盒的标记方式不统一。Illumina TruSeq Stranded（dUTP法）产生的reads，其cDNA第一条链来自RNA的反转录，测序得到的read1方向与RNA相反，因此应该用-s 2。而有些老式试剂盒，比如以前的SciClone或者某些接头方法，可能需要-s 1。

实际上最稳妥的办法是拿一个已知单链高表达基因来做验证。比如在人类样本中选取一个已知只从正链转录的基因（如GAPDH），看看-s 0、-s 1、-s 2三种设置下双端的counts哪一组能正确反映其表达方向。也可以用RSeQC的infer_experiment.py工具，输入BAM和GTF，它能估计链特异性系数并给出推荐参数：

infer_experiment.py -r annotation.bed -i sample.sorted.bam

输出会告诉你“Fraction of reads explained by”两条链的比例，如果两条链比例接近，是没链特异性（-s 0）；如果反义链比例大于0.8，就是你想要的-s 2。这个工具虽然老，但判断方向这件事上从来没有过时。

5.3 多映射reads该怎么处理

重复区域或多拷贝基因（比如rRNA基因簇、组蛋白基因）会让一条read能比对到基因组多个位置，这叫多映射reads。featureCounts默认会丢弃这些reads，不计入任何基因，以避免多重计数造成混淆。但在某些特定分析中，比如研究rRNA相关基因的表达，把这些reads完全丢掉会导致表达量被严重低估。

featureCounts给了一个-M参数，开启后多映射reads会被计数，但会均匀分配到所有可能的位置。这个参数在转录组结构分析中有争议，建议在大多数差异表达分析里保持默认关闭，除非你有明确理由。还可以结合--fraction参数，让多映射reads的计数按比例拆分，而不是每个位点都计一次，这样更温和一些。

在基因组富集分析中，曾有研究指出多映射reads的过度排除会导致重复区域相关基因的信号丢失，这是不是bug，是人可能没注意到的统计偏差。所以在处理多映射reads之前，先问问自己：我关心的基因是否落在重复区域？如果是，那就要在比对时使用-M --fraction；如果不是，默认丢弃对你没影响。

5.4 样本与样本之间的基因ID对不齐怎么办

这个坑我印象很深。早期我处理一批公共数据集，每个样本的GTF来源各不相同，导致count矩阵汇总时，有的基因叫ENSG00000…，有的叫NM_001…，有的甚至混用了基因symbol。这种情况在批量下载公共数据时特别常见。最直接的解决办法，是生成一个统一的gene_id到symbol的映射表。

建映射表的过程又会遇到一个经典问题："怎样用函数比对两列打乱的数据并找出不重复的数据"。比如一款注释文件里有gene_id和gene_symbol两列，但顺序被故意打乱，或者有大量重复的symbol对应不同gene_id。这时候用Excel的VLOOKUP或R里的merge都可以，但VLOOKUP有个很尴尬的点：重复值只会匹配到第一个结果，无法看到所有映射关系。我建议用R的dplyr包来处理，清晰且可回溯：

library(dplyr) df1 <- read.table("gencode_genes.tsv", header=TRUE) df2 <- read.table("counts_genes.tsv", header=TRUE) merged <- left_join(df2, df1, by = c("gene_id" = "gene_id")) # 找出在df2中但不在df1中的gene_id，即无法映射的ID unmatched <- df2 %>% anti_join(df1, by = "gene_id")

这样处理完，你就知道哪些基因ID在两个文件里对不上，哪些在df1里有重复记录。之后再用distinct()过滤重复symbol，或者按基因ID去重后，再重建count矩阵。这一步看起来琐碎，但做不好直接导致下游差异分析结果不可复现。

5.5 批量处理时的脚本技巧

如果你有几十个样本，当然不想一个个手敲featureCounts命令。可以用一个简单循环把位置参数拼出来：

counts_bam=$(ls /path/to/bams/*.sorted.bam | tr '\n' ' ') featureCounts -T 12 -a annotation.gtf -o counts.txt -t exon -g gene_id \ -s 2 -p --countReadPairs -Q 10 -B -C ${counts_bam}

这一行把目录下所有sorted.bam一次性传给featureCounts，它会自动按样本名区分列。注意通配符展开时不要混入其他bam，否则会导致某列是无效文件。写完脚本后再加一层判空保险，比如检查${counts_bam}里至少有一个文件名，避免跑了个寂寞。

6. 下游分析衔接：差异表达与可视化

6.1 让count矩阵直接对接DESeq2

定量做完，下一步最常见就是差异表达分析。以DESeq2为例，它的输入其实非常简洁：一个count矩阵加一个共有一列样本信息的colData。featureCounts的输出可以直接整理成下面这种格式：

library(DESeq2) # 读取featureCounts输出，跳过注释行 countdata <- read.table("counts.txt", header=TRUE, row.names=1, skip=1) countdata <- countdata[, 7:ncol(countdata)] # 准备样本信息 condition <- factor(c("control", "control", "treat", "treat")) coldata <- data.frame(row.names = colnames(countdata), condition) # 构建DESeq2对象并运行 dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ condition) dds <- DESeq(dds) res <- results(dds, alpha = 0.05)

这个流程需要特别注意的地方是，DESeq2会自己对count做归一化，不要再手动转TPM或CPM。如果先转成TPM再塞给DESeq2，它的内部离散度估计和方差稳定化步骤都会出问题，结果会失真。我一直觉得，把"归一化"这件事交给响应的统计模型去做，不要自己做太多额外加工，是最安全的。

6.2 一个完整示意的运行输出

跑完DESeq2后，你会拿到一个results对象，包含log2FoldChange、p值、padj等。按照padj < 0.05且log2FoldChange绝对值大于1的阈值，就能筛出候选差异基因。之后可以用ggplot2画火山图、PCA图，或者用pheatmap画热图。

我个人习惯在进入富集分析之前，先做两个检查：第一，确认总样本数不要少于3 vs 3，否则统计功效不够，差异基因数量会非常少；第二，检查PCA图上样本聚类是否跟实验分组匹配。如果对照组和实验组没有明显分开，可能是批次效应太强，这种情况要先做批次校正或者用limma的removeBatchEffect，再继续下游。这些操作虽然不在featureCounts的直接范畴内，但属于"从定量到下游分析"这个完整流程里绕不开的一环。

6.3 跨工具联合使用的经验

在有些场景里，我还会用featureCounts的结果辅助验证alignment-free定量工具的结果。比如用Salmon跑完转录本定量，再把结果汇总到基因层面，然后和featureCounts的基因计数做相关性分析。如果两者在大多数样本上的相关系数很低，基本可以断定某个环节出了问题。这类交叉验证在大型多组学项目里很有价值，因为单靠一种工具的定量结果有盲区。

有一回，一个项目中featureCounts和Salmon结果差异很大。后来发现是Salmon使用的转录本数据库版本比STAR的参考基因组版本新了一整版，两边的基因注释ID混着用，导致比对阶段和定量阶段的基准根本不同，后续所有联合分析都失真。从那以后，我在任何流程里都坚持记录所有工具的版本号、索引版本、注释版本，最好能生成一份yaml或json格式的配置文件存到项目目录。这不算麻烦，但在你复现结果或者处理审稿意见时，能省下大量时间和沟通成本。

6.4 进一步拓展：从基因水平到转录本水平

featureCounts默认在基因水平计数，但它也支持转录本水平。只需要把-g参数从gene_id改成transcript_id，其他逻辑不变。但要注意，转录本水平定量的复杂度和解释难度都会上升。同一基因的多个转录本共享外显子，新外显子组合信息可以造成一个reads同时对应多个异构体。featureCounts对这种共享外显子的reads有局部处理策略，但没法做到精确的异构体分辨。如果你想做异构体水平的表达分析，更稳妥的是Salmon或者RSEM这类基于转录本的定量方法，它们会结合序列组成来分配reads归属。

因此我的建议是：绝大多数差异表达项目，用featureCounts做基因水平定量，简单透明且稳定性好；只有当课题核心是异构体切换时，才考虑转录本水平，并通过特定方法交叉验证。

7. 实战场景中的小技巧与心得

7.1 一个常见参数组合速查表

说了这么多，我把常用参数的组合场景做成一个小表，方便你们在实际跑数据时快速翻查。

这些组合并不是死的，不同分析目标会微调，但如果拿不准，就先用这个表里的配置跑一遍，绝对不会犯大方向错误。

7.2 再谈一次"函数比对两列打乱数据"的实际案例

前面第5.4里用R处理的，是一个"两列数据打乱后找出不重复项"的场景。这个坑不只在featureCounts的下游会遇到，在构建注释映射表的时候更常见。比如我从Ensembl下载的GTF，与另一个数据库导出的gene symbol表，顺序往往是不一致的，同一个ENSEMBL ID可能出现多次，也可能完全找不到对应symbol。对应的处理思路是：right_join会保留右侧独有的，outer_join会保留两边的所有记录，而你要灵活组合这些操作，才能看清全貌。

另一种常见情况是Excel环境中用VLOOKUP去匹配乱序两列，最后匹配出一堆#N/A。我建议直接切换到R或者pandas处理，原因是表格数据在做键值匹配时，特别是多对一、一对多关系里，可视化地逐个检查处理结果，能避免VLOOKUP匹配到错误候选值的无声失败。我当年写的最多的一句话是"下次不要用VLOOKUP匹配基因ID"，现在还是想多说一句：匹配之前，先看清楚键的唯一性和重复度。

7.3 featureCounts之外的补充选择

最后聊两句横向选择的问题。featureCounts不是目前市面上唯一的定量工具，但它是应用面最广的。你可以用HTSeq-count，功能类似但速度非常慢；也可以用RSEM，它做转录本分配更精细，但耗时高很多；还可以走Salmon、kallisto这条alignment-free路线。如果你处理的是超大规模样本，比如几百个样本的群体转录组，featureCounts配合STAR的two-pass模式通常能扛住；如果样本量特别巨大且不需要解析新异构体，Salmon的速度优势就会体现出来。

选型这件事，没有绝对最优，只有最适合当前数据规模和科学问题。特色实验如单细胞RNA-seq通常不走featureCounts，而是用STARsolo或者Cell Ranger内部的CR计数，因为它们需要处理UMI和barcode。但在bulk RNA-seq和大多数常规转录组课题里，featureCounts永远是那个可靠的下限保障。

所以我的习惯是：别神化任何工具，先跑一遍默认配置，看一眼summary，检查每条归类，再根据实际问题做精细调整。定量分析这件事，耐心和细致比炫酷参数重要得多。