做甲基化分析这几年,我见过太多人踩坑。
特别是拿到 GEO 数据后,面对那一堆密密麻麻的 ID,头都大了。
很多人第一反应就是去搜“怎么注释”,然后点进一堆广告链接。
结果下了一堆软件,装了一堆依赖包,最后跑出来全是 NA。
真的想骂人。
今天我就把压箱底的经验掏出来,不整那些虚的。
咱们直接说怎么用最笨但最稳的方法,搞定 geo 2r差异甲基化区域注释。
先说个真实案例。
上个月有个粉丝找我,说他跑了三个月的脚本,注释结果跟文献对不上。
我一看他的代码,好家伙,直接拿 probe ID 去对基因名。
这能对上才怪。
甲基化探针和基因之间,哪有那么多一一对应。
有的探针在启动子区,有的在增强子,有的在基因间区。
你直接当基因名用,不报错才怪。
所以,第一步,千万别急着注释。
先看清楚你的数据格式。
是 beta 值?还是 M 值?
如果是 Illumina 450k 或者 EPIC 芯片,那 probe ID 肯定是以 ILMN_开头的。
这时候,你手里得有个“地图”。
这个地图,就是探针注释文件。
很多人不知道去哪下,其实最简单的方法,就是去 GEO 页面找 Supplementary Data。
通常作者都会上传一个 annotation file。
如果作者没上传,别慌。
去 Bioconductor 找对应的包。
比如 HumanMethylation450kanno.ilmn12.hg19 这种。
下载下来,里面全是探针对应的基因信息。
这里有个坑,很多人直接用 merge 函数合并数据。
结果发现数据量变大了,因为一个基因对应多个探针。
这时候你要决定,是取平均?还是取最大值?
我一般建议,如果是做差异甲基化区域(DMR),最好先做 DMR 预测,再注释。
而不是先注释再找差异。
因为先注释,你会丢失很多非编码区的调控信息。
这才是 geo 2r差异甲基化区域注释 的核心逻辑。
别被那些标题党误导,说只要注释就能发现 biomarker。
那是扯淡。
第二步,清洗数据。
把那些不在染色体上的探针,剔除掉。
把那些跨染色体映射的探针,也剔除。
这一步很关键,不然你后面画图,点都点不到正确的位置。
我见过有人把线粒体 DNA 的探针也放进来,结果热图全是红色,吓死人。
第三步,关联基因。
这里有个小技巧。
不要只关联 TSS 区域。
要把 promoter, 5'UTR, 3'UTR, gene body 都考虑进去。
特别是 gene body 区域的甲基化,有时候跟基因表达是正相关的。
这点很多人会忽略。
我在处理一个肺癌数据集时,就发现好几个抑癌基因,启动子没甲基化,但基因体甲基化很高,表达反而低。
这说明啥?
说明调控机制很复杂。
这时候,如果你只用简单的 geo 2r差异甲基化区域注释 方法,就会漏掉关键信息。
第四步,可视化。
用 ggplot2 画个 Manhattan plot。
看看你的差异位点,是不是富集在某个染色体区域。
如果散得像撒胡椒面,那大概率是你注释错了,或者批次效应没去除干净。
最后,我想说,做生物信息,心态要稳。
别指望一键生成完美结果。
每一步都要检查,每一步都要思考。
你付出的每一分细心,都会在最后的图表里体现出来。
别偷懒,别复制粘贴别人的代码而不看含义。
这才是做科研该有的样子。
希望这篇干货,能帮你省下几个熬夜的夜晚。
毕竟,头发比代码值钱。
加油吧,同行们。
本文关键词:geo 2r差异甲基化区域注释