生信分析里10个静默算错的坑，跑完才发现就晚了

生信分析里10个静默算错的坑，跑完才发现就晚了

不报错的错误，往往最贵

生信分析最痛苦的，不是报错，而是不报错。流程跑到了终点，图画完了，汇报的时候被问了一句“你这个分组是不是反了”，前几天的工作全部白费。这类问题有一个共同的名称：静默的错误——正常的退出和没有warning的日志，但结果是错误的。

今年很多避坑类的文章都在说同一个道理——不是不会写代码拖慢了进度，而是不知道哪一步悄悄改变了结论。下面10个地方按流程顺序排，每条讲清楚哪里错了，会导致什么问题，如何预防。它们大多是参数默认值、样本顺序、ID 映射这类“看起来不重要”的细节，恰恰是新手跑完才发现的坑。

坑一：样本表与表达矩阵各自为政。由GEO下载的矩阵的列名经常是GSM号，而在临床信息表中是另一种顺序。用行号拼起来，用cbind粘上去，就算样本的顺序变了，也不会报错。后果很直接：分组标签错乱，差异基因方向整体颠倒，甚至“富集”出几条看似合理的通路。防范的办法很简单-所有表都使用唯一的样本ID进行merge，join前后打印维度和分组频数以及前几行，在跑差异之前随便抽取3个样本回原始数据核对一下表达值。这步花五分钟，可省去三天返工。

坑二：参考基因组与注释文件不匹配。使用hg38，而定量使用的是hg19的GTF，又或者genome是来自UCSC，而GTF来自于Ensembl，一个是chr开头的，另一个不是。运气好的话，工具会直接报错；运气不好的话，就是一大批reads落在“没有注释”的地方，定量值普遍偏低，你以为数据质量不高。芯片数据中为平台注释版本错误，探针被映射至旧注释，基因名张冠李戴。将genome、GTF的版本、来源和下载日期写入一个README中，比对完成后查看assigned的比例是否与同类数据相当，基本可以排除。

坑三：链特性和双端方向弄反。dUTP建库链特异性定量按非链特异计算，会导致反义链上的转录本被计算在内，lncRNA与反义转录本定量系统性偏差；链方向设置错误，正反义颠倒。双端数据R1、R2互换也不报错，比对率和定量都很诡异。这类问题需要询问建库的公司，或查看方法部分。不确定可以运行一次链特异性推断以检查它是否与预期一致。

坑四：多比对reads、重复序列默认处理。同样的基因，不同的定量工具对于multi-mapping reads的处理策略也是不同的，有些直接忽略，有些则保留并平分权重。换个工具，换个默认设置，表达量可能差一截，而你只会觉得“这批数据怎么这么不稳定”。PCR重复也是一样，去重在一些流程中是常规操作，在某些流程中会降低高表达基因的强度。的工具版本和重要参数需要记录，更换工具时最少用一份数据对几个管家基因进行交叉验证。

坑五：用FPKM/TPM做差异分析输入。这是最经典、也最容易“反正都能跑通”的错误。DESeq2、edgeR等工具是基于原始counts的均值-方差关系的统计模型，输入TPM相当于破坏了模型的前提，跑是能跑完，但得到的p值是不可信的。FPKM和TPM也不能混着比：FPKM是先除长度，再除总的reads，每个样本加起来的总和并不相等，因此跨样本无法比较；TPM是先归一，再除长度，这样每个样本加起来的总和是相同的，样本之间可以比较，但是它依旧无法代替counts来进行检验。其次，FPKM依赖于你所使用的转录本集合，注释的变化会导致其值的变化，所以也不适合用于保存比较。

坑六：基因ID转换这里丢东西了。EnsemblID后带有版本号，去掉版本号则无法匹配；一个symbol对应多个EnsemblID，使用求和聚合会将不同转录本的表达量叠加为一个不存在的较高表达量，使用求平均数则会抹平真实差异；一个Ensembl对应多个symbol，使用merge会产生重复行，后续的统计结果将重复计数。还有更阴险的，为了“干净”，把表达量=0的行整行删除，恰好删除的基因是某一通路集中的，后面的富集结果就不对了。转换前后的各看一次行和重复数、NA数，这是最廉价的防线。

坑七：设计矩阵忽略了批次和配对信息。最坏的例子就是批次和分组是混杂的，第一批全是对照，第二批全是对照，那这个时候，所有的算法都是无能为力的，差异表达的基因其实就是测序批次的差异。轻一点的批次没有写在design中，把批次效应计算到残差中，p值偏保守或者激进，自己也不知道。配对样本也一样，同一个体处理前后的、同一只鼠左右的都要采用配对设计或把个体当作协变量纳入模型。拿到样本先画一个分组*批次的表，看看有没有空白的格子，有批次的话写进design，不要只在讨论中提及。

坑八：将技术重复误认为是生物学重复。同一个样本分两管建库，高度一致只能说明建库测序稳定，并不代表个体差异。但人们往往把它们当作两个重复进行分析，导致方差低估和假阳性倍增。另一种极端是每组只有一例也照样跑出p值-那不是统计，是算术，工具只是没有拒绝你而已。在样本表中需要专门设置一列区分生物重复和技术重复，做差异分析之前打印每组样本数量，少于3个不要强行解释通路。

坑九：只看p值不看校正后的padj或事后调阈值。测试20,000基因一次，p < 0.05，仅凭运气就应该有1,000个“显著”。DESeq2的结果表中padj才是经过多重检验校正的，但是很多人按照pvalue筛选，筛选后发现基因过少，就把阈值从0.05放宽到0.1，log2FC从1降到0.58，直到结果“好看”为止。事后调整并非绝对禁止，但若要写在方法部分，最好要在独立的数据中加以验证。更稳妥的方法是先确定前把的阈值，运行完毕不再更改，在汇报时给出padj和效应量的结果，而不是仅仅汇报一个基因数量。

坑十：富集分析的背景集和ID类型不符。ORA类的富集默认拿整个基因组做背景，但是你的实验只是检测到了部分基因，所以会造成系统的高估。正确做法是把“这次分析中实际有资格被检出的基因”作为背景集。其他几种常见的错误：人源基因列表用了小鼠注释库、Ensembl ID直接扔给只接受symbol的tool、上下调基因混合做不分方向的ORA。富集前确定三要素：物种正确性、ID类型正确性、背景是否为实际检测基因，再开始进行通路富集。

静默错误没法靠“小心一点”避免，只能靠固定动作。下面这些检查点，全部做一遍也就半小时，但能拦住上面绝大部分问题。建议把它们写成脚本，每次分析结束自动打印一遍；如果用的是 PaperOK 的生信模式这类会记录每一步输入、参数和产物的工具，回查时对着记录逐条核对会轻松很多，公开数据也能直接检索下载，省掉一部分手工搬运。但记录归记录，判断仍然得你自己做。

这些坑之所以叫静默，就是因为它不会给你任何提示。工具只是工具，它不会帮你判断样本是否配型错误、注释是否统一、输入是否使用counts。将上述动作固化，远比使用一个更“智能”的工具有效。

流程跑通只是及格，算对了才算数