企业官网建设流程全解析

1. 拿到表达矩阵之后，先想清楚"相关性"到底在问什么

做过几套转录组的人大概都有这么一段经历：差异表达跑完，火山图、热图都出得挺漂亮，导师或者合作方突然问一句——"这些lncRNA和旁边那些mRNA有关系吗？"于是你打开表达矩阵，盯着几千行乘几十列的数字发呆，不知道从哪下手。标题里说的"将lncRNA和mRNA进行相关性分析"，本质上就是回答这个问题：在样本层面，某个lncRNA的表达起伏和某个mRNA的表达起伏是否同步。同步，可能是共调控、可能是同一条通路被激活、也可能是这个lncRNA本身就在顺式或反式地影响这个mRNA。不同步，也不代表没功能，只是这条线索暂时断了。

我先把结论摆在这里：相关性分析在这类项目里属于典型的"下游挖掘"环节，它的价值不在于算出那一个皮尔逊系数，而在于帮你从几千个候选里筛出一个几十到几百对的短名单，供后续做顺式/反式预测、共表达模块划分、以及实验验证的靶点挑选。做得好，它是一条高效的信息压缩管道；做得糙，它就是一堆漂亮的假阳性，而且很麻烦的一点是，这些假阳性往往在生物学上讲得通，会让你误以为发现了新机制。

所以这篇文章想聊的不是"跑通一段代码"，而是这套分析里那些真正影响结论质量的决策点：输入数据该怎么准备才不至于一开始就带偏、样本量小的时候相关系数到底还能不能信、lncRNA这种低表达且噪声大的数据类型在相关性计算里有什么特殊性、以及拿到一堆显著结果之后怎么排优先级才不会挑到一个纯粹由批次效应驱动的"关联"。这些内容在标准流程文档里通常写得很简略，但实际做项目时，它们决定了你的结果能不能被重复、能不能被写进结论。

适合谁看：手上已经有一个表达矩阵（不管是自己跑出来的还是公共数据下载的）、会基本的数据框操作、但还没系统做过lncRNA-mRNA联合分析的人。如果你连TPM和counts的区别都还模糊，建议先补一下上游定量，因为后面的所有讨论都建立在一个前提上——输入的表达量本身是可靠的。这不是客套话，我见过太多相关性分析失败，根子出在定量阶段用了不合适的注释或者样本被污染，跟相关性算法一点关系都没有。

还有一个预期管理的问题。这类分析很容易让人产生"我找到了调控关系"的错觉，但相关性只说明统计上的同步，方向性和因果性都不在其中。我个人的习惯是，在报告里把这类结果统一称为"候选关联对"或者"共表达证据"，绝不写"调控"。这个措辞上的克制，后来帮我省掉了很多解释上的麻烦。

2. 输入矩阵的清洗：三类会悄悄污染相关系数的脏数据

相关性计算本身极度简单，一行函数就能出结果，正因为简单，它几乎不会报错，也就不会提醒你输入有问题。所有的问题都会以"结果看起来不太对劲"的形式暴露出来，而这个时候你往往已经跑完了几十行代码。所以这一步值得花时间，把脏数据在前端拦掉。

2.1 低表达过滤的阈值不能拍脑袋，要跟样本数挂钩

最常被忽略的是低表达基因的处理。一个在30个样本里有28个是0、剩下2个是3和5的lncRNA，它的表达向量几乎是一串零，这种向量跟任何mRNA算出来的相关系数都没有意义，但它极有可能因为那两个非零值的偶然同向而给出一个不小的正相关。更麻烦的是，低表达基因的定量方差极大，技术噪声会主导信号，你算出来的"强相关"很可能只是两个基因的测序噪声恰好同频。

过滤阈值怎么定？我的经验是分两步。第一步用硬阈值做粗筛，要求至少在某一个比例以上的样本中表达量高于某个下限。常用的组合是"至少在50%的样本中counts≥10"，或者在TPM体系下"至少20%样本TPM>1"。第二步是针对lncRNA的特性放宽一点——lncRNA整体表达水平本来就比mRNA低一个量级，如果你用mRNA的标准去卡，最后剩下的lncRNA可能只有几百个，后续分析空间被榨干了。所以我一般会把lncRNA的门槛单独设：比如mRNA要求50%样本CPM>1，lncRNA要求30%样本CPM>0.5，然后把这个差异在方法里写清楚，不要让人觉得你是在凑结果。

注意：不要用"过滤后保留了多少个基因"来反向调阈值。这在项目里很常见——发现结果太少就把阈值调低，结果多了就把阈值调高。这种操作会让整个分析失去统计意义，写进方法学部分也会被审稿人抓住。

这里还有个容易混淆的点：过滤应该在标准化之前还是之后？我的做法是先用counts做过滤（因为counts的方差结构更符合过滤所需的判断），再对过滤后的矩阵做标准化（TPM/CPM/VST）。如果你先标准化再过滤，标准化的缩放因子会被那些本该被剔除的低表达基因拉偏，尤其是CPM这种基于总量归一的方法，一堆零值基因会稀释分母的合理性。

2.2 样本层面的离群点比基因层面的问题更致命

基因过滤解决的是"每个变量够不够格参与计算"，样本离群点解决的是"每个观测值是否可信"。一个被污染的样本、一个降解严重的样本、或者一个样本分组标错的样本，它对所有基因的影响是一致的，所以它的存在会系统性地抬高或压低大量基因对的相关性，制造出一批假关联。而且这种污染在整体相关性上体现得很隐蔽，因为它的影响是全局的、同向的。

我一般会做三层检查。第一层是样本间相关性热图，看有没有某个样本跟其他所有样本的相关性都明显偏低（比如均值低出0.1以上）。第二层是PCA或者MDS，看离群方向和离群距离，同时观察它是否偏离了自己的分组中心——有时候样本在PC1上偏离，但恰好和另一个分组靠得近，这往往意味着标签有误。第三层是对可疑样本单独检查表达分布，看它的基因检出数、中位数表达水平、以及有没有大量基因处于极端值。

要不要删？这取决于偏离程度和生物学解释。如果偏离是可以被技术原因解释的（比如某个样本的RIN值明显低），而且删掉之后不影响分组平衡，那删。如果偏离可能反映真实的生物学异质性（比如某个个体本身状态特殊），那更稳妥的做法是保留但在敏感性分析里做去留对比，看核心结论是否稳健。这个对比做起来不费劲，但能显著提升结论的可信度。

2.3 批次效应必须在相关性计算之前处理，不能之后

这一条我想单独说，因为它是最容易被低估的。假设你的样本分两批测序，而这两个批次恰好和你的分组高度重叠（这在真实项目里非常常见，因为没有谁会特意打乱批次和分组的对应关系），那么所有基因在两个批次间的差异都会表现为"组间差异"，进而所有基因的表达向量都会因为分组结构而呈现出共同的变化模式。这时候你算lncRNA和mRNA的相关性，算的其实是"它们是否都受批次影响"，而不是"它们之间是否有调控关系"。结果就是你会得到一大堆跨全基因组的、毫无拓扑结构的强相关，看起来信号极强，实际上全是系统偏倚。

判断方法很直接：把批次作为协变量，看PCA上前几个主成分和批次的相关性，如果批次解释了很大的方差比例，就必须处理。处理方法上，如果批次和分组完全混杂（即每个批次只包含一个分组），那说实话没有干净的解法，只能承认这个混杂无法在统计上消除。如果批次和分组有部分交叉，可以用ComBat这类方法做校正，或者更简单粗暴但有时更稳妥的做法是把批次放进下游模型作为协变量（这对相关性分析不适用，但对后续验证有用，因为相关性计算是两两之间的，不太好直接塞协变量）。

如果确实需要校正后再算相关性，我一般会对比校正前后的相关性分布。校正之后，如果全基因组相关性中位数从0.6降到0.2左右，而少数目标基因对的相关性保持稳定，那说明校正做对了——去掉了共同的系统偏倚，保留了特异信号。反之，如果所有相关性都均匀下降到接近0，那你可能把真实信号也一起削掉了，得回头检查校正模型。

3. 相关系数怎么选：皮尔逊、斯皮尔曼和这中间的那些坑

工具链上这块其实没有什么神秘的东西，R里的cor()、Hmisc::rcorr()、或者Python里的scipy.stats.pearsonr，一行就能出结果。真正需要动脑的是选哪个系数，以及为什么。这个选择会实质性地改变你的结果集，不是风格问题。

3.1 皮尔逊假设线性，而表达量数据经常不听话

皮尔逊相关系数衡量的是线性关系，它默认两个变量服从近似正态分布、没有极端离群值、方差大致齐性。转录组数据在这三条上经常全部踩雷：counts是右偏的，标准化之后即使接近正态也常有长尾，而且某些基因（尤其是低表达基因）会在个别样本上出现极高的值，这个值会主导整个相关系数。

所以我的常规做法是：要么对表达量做log2转换（TPM/CPM加个小常数后取log），要么直接用斯皮尔曼。log转换是必须的，这一点上我不太接受"TPM已经是标准化的所以不用转"这种说法——标准化解决的是量纲和测序深度，不解决分布形态。你可以在log转换后再算皮尔逊，这时它的行为会稳定得多。

至于斯皮尔曼，它本质上是对秩次做皮尔逊，对单调非线性关系敏感，对离群值鲁棒。对于lncRNA这类噪声较大的数据，斯皮尔曼往往更稳。但它有个副作用：当样本数很少、且存在大量并列值（比如低表达基因里一堆零）时，秩次的处理会引入偏差，相关系数的分布会变得很奇怪。所以我一般会两种都算，然后看两套结果的交集。取交集是个很朴素的稳健性策略，虽然简单，但在报告里很好解释——"在两种相关性度量下均显著"这个表述比单纯说"相关性显著"要扎实得多。

3.2 样本量小于15的时候，别急着看p值

这个坑我踩过，而且踩得挺深。相关性的显著性检验对样本量极度敏感。在n=6的情况下，|r|>0.81大概就能到p<0.05；在n=50的情况下，|r|>0.28就够了。这意味着在小样本项目里，你很容易得到一批"显著"的关联，但它们对应的相关系数实际上很高却不稳定——换一个样本重测，可能就掉到0.3去了。

更麻烦的是多重检验。假设你有2000个lncRNA和15000个mRNA，全部两两计算就是三千万对。哪怕用BH校正，在小样本下你依然可能得到成百上千的"显著对"，而这些对里有多少是真的，没人知道。我的处理策略是双阈值：既要求校正后的p值达标，也要求|r|超过一个经验下限（大样本时0.6，小样本时0.8）。这个|r|阈值不是从统计表里查的，是从实践里攒的——我做过几次同一批样本的重测，发现真实共表达的基因对，跨批次的重现相关系数大多在0.6以上；低于这个水平的对，重现率非常低。

另外提一句，样本量小的时候，与其纠结单一基因对的相关性，不如把注意力放到模块层面。共表达模块（几个基因聚成一群，整体表达模式一致）在样本量小的情况下鲁棒性远好于单对基因，因为聚类用了全部基因的信息，等效于做了降噪。WGCNA那套思路在这类场景里比逐对计算更靠谱。

3.3 偏相关能救命，但也可能把真信号一起剔掉

当两个基因的相关性是由第三个因素驱动的（比如一个公共的上游转录因子，或者细胞组成比例），偏相关可以把这个第三变量的影响剔除掉。听起来很美好，实际操作上有两个陷阱。

第一，你得知道应该控制什么。如果控制变量选错了，你会得到一个无意义的残差相关。常见的合法控制变量包括：已知的批次、样本的细胞组成估计值、以及一些显著的技术协变量（比如RIN、测序深度）。但如果你把某个中间通路的基因也塞进去当协变量，那就等于在问"去掉中介之后还有没有直接效应"，这已经是另一个问题了。

第二，偏相关的计算在小样本下极其不稳定。它本质上是在做回归取残差，样本量小时残差的估计误差会放大。我一般只在样本量足够（n≥20）且有明确混杂因素的情况下才用偏相关，并且一定同时报告原始相关，让读者知道偏相关改变了多少。

实际报告里，我会用这样一张表来呈现关键结果，比单纯堆p值清楚得多：

4. 从几千万对里筛出能看的那几百对：策略与取舍

算完之后你会面对一个巨大的矩阵，直接看是看不动的。这一步的核心是分层筛选和排序，把注意力集中在最可能成立的那部分结果上。

4.1 先按基因组位置分层：顺式候选和反式候选分开处理

lncRNA和mRNA的关联在生物学上有两类非常不同的成因。一类是位置上的邻近——lncRNA在基因组上靠近某个mRNA（比如在它的启动子上游几kb内，或者在其内含子区域），这种情况下两者的表达相关很可能是顺式的局部转录调控或者共享调控元件导致的。另一类是远距离的共表达，可能是通过某个共享通路或反式调控因子。

把这两类混在一起排优先级是不合理的。我的做法是先按坐标把lncRNA和mRNA做位置配对（一般取转录起始位点距离在100kb以内作为顺式候选窗口，这个窗口大小可以按物种调整），然后：

• 对位置邻近的对，相关性只作为辅助证据，重点看方向和距离。即使相关性中等（0.4~0.5），只要位置合理，也值得进候选名单。
• 对远距离的对，要求更高的相关性阈值和更严格的显著性，因为随机共表达的概率更高。

这样做的好处是优先级有了生物学依据，而不是单纯按p值排序。按纯统计量排序的结果，排在最前面的往往是一批功能上毫无关联的高表达管家基因对，这种结果汇报出去很尴尬。

4.2 用表达模式聚类降低逐对比较的噪声

逐对相关性有个固有缺陷：它只看两个基因，信息量小，对噪声敏感。换个角度，先把所有基因的样本表达谱做聚类，把表达模式相近的基因聚成若干模块，然后看哪些模块里同时富集了lncRNA和mRNA。这种模块层面的一致性比单对相关性稳健得多，尤其是在样本量不大的时候。

具体做法是：把过滤后的表达矩阵（log转换、标准化）做样本间或基因间的层次聚类，切出几十个模块，然后用模块内的平均表达谱（或者第一个主成分，即模块特征基因）去和其他模块算相关性。模块内如果lncRNA和mRNA都占一定比例，那这个模块内部的基因对就是高优先级的候选，因为它们不仅彼此相关，而且和一堆其他基因同步变化，说明背后有个共同的调控程序在驱动。

这一步的额外价值在于，它天然引入了"多方证据"的概念。一个lncRNA-mRNA对如果既在两两相关里显著，又在同一个共表达模块里，那它的可信度要比单纯的逐对显著高不少。我在筛选时会给这类对加分，排序时往前放。

4.3 结果排序别只看p值，加几个"人工加权项"

这是我最想强调的一节，因为纯按统计量排序几乎必然导致结果没有生物学可读性。我的排序表里通常包含这样几列：

• 校正后的p值（基础门槛）
• 相关系数绝对值（强度门槛）
• 是否顺式（位置优先）
• 是否同模块（一致性加分）
• 两个基因的表达水平是否都足够（避开低表达噪声对）
• 是否存在已知的功能注释线索（如lncRNA邻近已知免疫基因等）

这几项加权之后排出来的名单，可解释性会好很多。加权的方法不需要很复杂，简单排序打分就行，甚至手工挑都行——反正最后进实验验证的也就十几对，前面的排序只要保证这几对确实排得靠前就达到目的了。

提示：排序之前一定要把结果表导出成带基因名、坐标、相关方向、强度、显著性的宽表存下来。后面写报告、做图、跟合作方沟通，全指望这张表。临时再算一遍很容易因为环境或者随机种子不一致导致数字对不上，这种低级错误在项目后期返工时非常耗时间。

5. 那些让我返工过的具体问题

这一节讲几个实际踩过的坑，都是标准流程里不会提的，但每个都让我多花了一两天。

5.1 注释版本不匹配，配对全部错位

有一次我拿到一个lncRNA注释文件和一个mRNA注释文件，分别来自不同的数据库版本，染色体命名一个是chr1，一个是1，坐标也是基于不同的参考版本。我直接按基因名做配对，结果算出来一批"位置邻近"的候选，位置关系完全是错的。这个问题肉眼极难发现，因为基因名是对的、坐标看起来也合理，只是参考基因组版本不对。

从那以后我加了一条硬性检查：在配对之前，先核对两个注释文件的染色体命名规则、坐标体系版本、以及基因ID类型（Ensembl ID还是Symbol，带不带版本号）。如果是从不同来源拼起来的注释，最好统一用同一套GTF重新导出。这个检查只花几分钟，但能避免后续全部返工。

5.2 相关系数的方向在报告里被写反了

正相关和负相关在报告里的表述很容易出错，尤其是当你在代码里对某些基因做了取反操作（比如log转换时加了负号，或者某些方法内部做了秩次反转）之后。我在一次内部汇报里把一批负相关写成了正相关，因为中间做了一次矩阵转置，导致方向和基因的对应关系整体翻转。

这个问题的根源在于，相关性分析里基因对的"方向"是由行名和列名的对应关系决定的，一旦中间有任何一步重排（转置、排序、子集提取），很容易错位。我的做法是：算完相关性后立刻做一次抽样验证——随机抽5对，手动从原始矩阵里把两个基因的表达值拉出来，肉眼看一下走势，确认正负号没错。这个检查非常土，但非常有效。

5.3 把小样本的高相关当成了强证据

前面提过，这里再展开说一下那个具体的教训。有个项目只有8个样本（4对4），我算出一批|r|>0.9的基因对，当时很兴奋，觉得信号极强。后来补做了几个样本，这些对里有超过一半的相关系数掉到了0.5以下。原因很简单，8个样本里只要有两三个样本的表达值碰巧同向，相关系数就会被拉得很高，而这种"碰巧"的概率并不低。

从那以后我给自己定了一个规矩：样本数低于10的时候，相关性结果只做探索性描述，不写进结论，也不拿去挑验证靶点。如果样本实在少，我会退而求其次，用表达模式的定性描述（比如"在处理组中一致上调"）替代相关系数，这种表述虽然不够精确，但至少不会被误读。

5.4 忽略了lncRNA自身的位置和链方向

做顺式配对时，lncRNA的链方向很重要。反义链上的lncRNA和它的重叠mRNA，表达相关性的生物学含义与同向的情况不同——反义转录本有时会抑制正义链的转录，这时你会看到负相关；有时又会因为共享调控元件而正相关。如果不区分链方向，把它们混在一起统计，结论会自相矛盾。

我在配对表里加了链方向这一列之后，很多原本看着矛盾的结果就解释通了。这个细节在读文献时经常被一句"邻近的lncRNA"带过，但实际做分析时它是必填项。

6. 结果交付前的自查清单

跑完分析、排完序、挑出候选之后，别急着写结论。下面这几项检查是我每次都会做一遍的，能拦掉大部分会在后面被质疑的问题。

第一项，重跑一次随机对照。把所有样本的标签随机打乱，重新算一遍相关性，看能出多少"显著对"。这个数字就是你结果集里假阳性的基准水位。如果你的真实分析里显著对的数量和随机打乱后差不多，那说明你找到的东西很大程度上是噪声。这个对照做起来很快，但说服力极强，写进报告里也显得严谨。

第二项，换一种标准化方法验证核心结果。如果主分析用的是TPM，就用VST或者log-CPM再算一遍，看排名前二十的候选对是否依然在列。稳健的关联对不会因为标准化方法变化而消失。这一项能同时验证你对低表达过滤和log转换的处理是否合理。

第三项，检查候选对在单个样本层面的表现。具体来说，把几个重点候选对的散点图画出来，看看是不是有某个样本在驱动整个相关性。如果去掉一个点之后相关系数大幅下降，那这个对就不够稳。这个方法在统计上叫影响点诊断，做起来就是画个图的事，但对判断结果的可靠性非常有用。

第四项，核对候选mRNA的功能注释是否和你的研究场景自洽。如果一个关于免疫应答的项目，筛出来的高优先级候选mRNA全是嗅觉受体或者角蛋白家族，那大概率是哪里出了问题，得回头查。这不是过度解读，而是基本的合理性检查。

第五项，把完整的参数记录下来：过滤阈值、标准化方法、相关性度量、多重检验校正方式、阈值设定、样本剔除情况。这些信息在半年后你自己都未必记得，更别说合作者或者后续接手的人。我的习惯是在项目目录里放一个参数说明文件，每改一次参数就更新一行，附上修改理由。这个小习惯在项目复盘时价值极高。

还有一个经验性的判断：如果一个lncRNA和几十个mRNA都呈现强相关，先别高兴，这更可能说明这个lncRNA的表达谱本身受某个全局因素影响（比如它在某个样本里因为技术原因异常高），而不是它真的调控了几十个靶标。真正特异的调控关系通常表现为"少数强关联加一些弱关联"，而不是"一大片一样强的关联"。看到过度广泛的关联，第一反应应该是去查这个lncRNA的表达分布，看它是不是在个别样本上异常。

最后分享一个我在实际操作中的体会：相关性分析在这类项目里的正确定位，是"生成可验证的假设"，而不是"证明调控关系"。把它当成一个高效的候选筛选器，用完就往下走——做顺式/反式预测、做靶基因富集、设计敲低实验——不要在相关系数本身上反复纠结。我见过不少人花大量时间调整阈值、换算法、调参数，试图让结果"更显著"，但如果上游的样本质量和定量本身有问题，下游怎么调都是白费。把时间花在样本检查、注释核对、以及后续的功能验证设计上，回报率要高得多。