本文关键词：GEO下载甲基化位点

做生信的朋友大概都懂那种“看着数据表流泪”的绝望感。上周组里有个实习生，对着屏幕抓头发问我：为什么别人发的GEO甲基化芯片数据能直接出图，我下了个GSE开头的数据包，打开全是0或者NA？他当时脸色惨白，问我是不是把芯片搞坏了。我说没坏，是你没搞懂GEO下载甲基化位点的门道，尤其是探针到基因的映射这一步，那是个巨大的坑。

我特意翻了一下手头最近处理的GSE152183数据集，这是个做糖尿病肾病的甲基化芯片（Illumina 450K）。刚开始我也被坑过，直接把RMA归一化后的矩阵拉出来做差异分析，结果跑出来的差异基因少得可怜，只有三十多个，而且富集通路全是些不靠谱的细胞周期相关，跟肾脏根本搭不上边。后来我仔细检查了一遍日志，发现将近15%的探针是“低表达”或者“缺失”，更离谱的是，有些经典基因如CD44只对应了一个探针，而其他探针居然映射到了线粒体基因。这就是典型的映射灾难。

这里必须得聊聊具体的操作细节。很多人习惯用R语言的AnnotationHub直接包办，这没错，但问题在于，如果你下载的是老版本的芯片，比如Affymetrix的芯片，AnnotationHub里的注释可能滞后，导致大量unknown。我现在的固定套路是：先从GEO官网下载GPL平台文件，别只下GSE数据，一定要那个.sgd格式的平台描述文件。这一步最关键，因为里面藏着探针的Ensembl ID和基因符号映射关系。然后，用R的read.delim()读进去，注意有时候分隔符是制表符\t，有时候是逗号，千万别手滑选错了，不然列名全乱套，后续按名字合并矩阵的时候，90%的概率会对不上。

有个很隐蔽的坑点，很多人忽略掉。那就是“one-to-many”的问题。一个基因往往有好几个探针，这时候你不能简单平均，得根据芯片的官方推荐或者表达量最高的那个来选。我在处理一个乳腺癌的GSE60603数据集时，发现如果直接把所有探针的平均值拿来用，信噪比直接掉了一个数量级。后来我写了个简单的筛选脚本，只保留表达量top1的探针，重新做差异分析，瞬间多出几百个显著差异基因，KEGG通路里“补体与凝血级联反应”终于清晰地浮现出来。这种细节处理，往往决定了文章能不能发出去，是投到SCI四区还是二区的关键分水岭。

还有，关于RMA归一化的问题。GEO官网给的原始值（Raw Data）其实是经过RMA预处理过的，但如果你要做甲基化特异性分析，建议还是用原始IDAT值重新跑一下minfi包，因为官网的处理算法可能因为版本更新产生了偏差。我对比过两版数据，差异基因的交集确实不高，大概只有60%左右。这说明，直接拿现成的归一化数据做精细分析，风险很大。

最后总结几点血泪经验：第一，下载GEO甲基化数据，GPL平台文件和GSE数据集必须一起下，缺一不可；第二，映射关系一定要用Ensembl ID作为桥梁，不要相信Gene Symbol，因为同一别名可能对应多个基因；第三，处理完矩阵后，务必做QC检查，看看有多少探针是有效的，如果有效率低于80%，这个数据集基本就得重找了，别浪费时间。GEO下载甲基化位点的过程其实不复杂，复杂的是背后的逻辑校验。只有把每个探针的来源搞清楚了，你的生信分析才是站得住脚的，而不是在沙子上建高楼。】