研究方向为食管癌的同学老师们一定知道GSE53625这个数据集,来源于《LncRNA profile study reveals a three-lncRNA signature associated with the survival of esophageal squamous cell carcinoma patients》这篇文章。据我所知这应该是目前食管鳞癌最大的数据集,并且分类简单明确,就是癌和正常组织,临床数据完整,个人感觉非常适合作为数据挖掘的开始部分。因此学完了生信技能树的生信入门课程后我决定从这个数据集入手。
一、下载数据
就是小洁老师给的代码,刚开始想直接从R里面下载,奈何学校的网实在太差,这个数据集104MB,硬是下载不下来,只能从GEO里找到表达矩阵用迅雷下的,结果还是下了1个多小时才下下来(学校的网真是欲哭无泪)。完事之后仔细检查了一下自主下载和代码下载的文件名和大小完全一样。接着提取表达矩阵、临床数据,saveRdata。
二、基因注释
按照小洁老师给的代码接下来就是分组和基因注释了,如前所述,这个数据集分组超级简单,第二列就给出了包含normal和cancer的样本名。当我以为一切顺利的时候问题就来了,当我用http://www.bio-info-trainee.com/1399.html这个网站找注释文件的时候,发现GPL18109搜不到,只能用第二种方法,在GEO里面下载soft文件,如下
可以看到这里面只有ID列并没有我所熟知的symbol列,只有对应的基因序列,难道要走自主注释?当时就崩溃了。搜了一下好像有软件可以用基因序列对应到gene symbol,难道我做个基因转换还得再学一个陌生的软件?于是我就在生信技能树公众号使用了关键词“基因序列”居然搜索到了这个数据集相关的推文,并且也给出了前辈们做好了的基因注释文件https://mp.weixin.qq.com/s/pMVNYn-kMljavQvnhvkS_w。
三、差异基因
该找的文件都已经找到,接下来就是关键的差异基因分析了。按照小洁老师给的代码认真整理完数据之后,激动的运行了limma,LogFC选的1,p_value选的0.05,之后table了一下change列,发现居然没有差异基因!!!!
LogFC调整到0.5之后才勉强有几百个差异基因。这里已经彻底懵逼,甚至用一秒钟时间怀疑了一下小洁老师的代码:)。为了验证我的结果,我搜了一下用这个数据集已经发表的数据挖掘文章,找到了一篇挖掘到3000多个差异基因的文章,结果相差这么多肯定不对劲,这时候我想到了GEO2R这个网页工具,并且下载了差异表达矩阵,再用代码标出上调下调基因,结果得到了跟上面那篇文章一模一样的结果,之后一切顺利,火山图,热图,GO,KEGG。
因为还想用这个原始表达矩阵做一些其他的分析,我决定找一下问题所在。这时我注意到第一步的数据处理有一个log表达量的步骤,难道是log之后表达量变小了差异就没有了?(当时还不知道limma包必须要log哈)之后我跳过了这一步重新运行所有代码,奇迹般地得到了正确的结果,当回过头来看了一下表达矩阵
突然意识到这表达量为什么这么小,结合刚才跳过了log步骤得到了正确的结果。。。。
最后真诚的希望小洁老师以后上课的时候能讲一下这个数据集,对于新手来说真的很不友好!