上次更新了一个GEO数据框的帖子,包含了GEO分析的主要组成部分,其中有些区域为了代码复用,不是很容易懂。现在对其进行注释,并增加KEGG分析的内容,力求能让他变成最有诚意的GEO数据分析教程,文末还录制了一个友好的导学视频,希望物尽其用。 以下是正文:GEO数据库包罗万象,对于每个领域的科研工作者很有帮助。 遇到比较困难的地方就使用别人的小工具来代替。比如,ID转换,GO分析等。直到去年我才能够完全使用R语言实现。 这两年,网上的教程帖子井喷一样涌现,大有教师超过学员的趋势。 但是,GEO的分析跟TCGA不一样,里面要考虑的细节比较多。
看文献得到GSE号,这是数据的名片,GSE32575, 在这个网址输入这个号我们可以知道这个芯片的信息 点击进去后,会有这个芯片的介绍 认真读一读就知道作者为什么要做这个芯片。 如果还不是很清楚,可以往下看,可以看到实验设计,甚至还有这个芯片当时发表的文章。下载下来读一读,可以看看当时这个芯片的哪些信息被使用了,还有哪些没被使用,我们可以提出新的科研假设,用别人的信息去支持。而这也是芯片数据挖掘的目的:
在往下看还有平台信息GPL6102,这个信息可以帮助我们来注释文件。 图中的48表示这个芯片做了48个样本,每个样本都写明是什么样本。 假如点击 Analyze with GEO2R 按钮,就可以按照以下帖子介绍的方法,无代码分析这个芯片。无代码GEO芯片分析图文教程 黑色区域是代码块,可以一整块一整块复制到Rstudio的脚本中,然后一行行运行,代码框中的#号表示注释,也就是这一行的内容不会被当成代码执行。 关于如何配置自己的R语言环境,可以看这篇帖子。 捣鼓好了之后,我们就可以往下操作了,而当你完成这个帖子的那一刻,你也会一脚踏进数据挖掘的大门,至少说,见了世面,以后不会被别人那些云山雾罩的话给忽悠。 1.加载R包获取数据。## 加载R包 以后只要更改GSE号,就可以直接下载别的GEO数据,很方便。 2.通过pData函数获取分组信息## 获取分组信息,点开查阅信息 3.通过exprs函数获取表达矩阵并校正exprSet=exprs(gset) 可以先简单看一下整体样本的表达情况,其中 boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2) 需要人工校正一下,用的方法类似于Quntile Normalization 这里我们用limma包内置的一个函数 library(limma) 4.判断是否需要进行数据转换我们通过分组信息已经知道,前12个样本本次不需要,所以先去掉 exprSet = as.data.frame(exprSet)[,-seq(1,12)] 打开现在的表达矩阵,是这个样子的 肉眼看到数字很大,这时候需要log转换(选log2)。 此外还有一个脚本可以自动判断是否需要log转换,这个脚本来自于GEO2R,我改写了一点点。 ex <> 这个语句会自动判断,如果已经log话就跳过,并提示
如果没有log话,他自动log2,并且提示:
这时候再来看这个数据就规则多了 但是我们发现一个问题,这个表达数据行名是探针,不是我们熟悉的基因名称如果TP53,BRCA1这样的,所以我们需要注释他。 5.获取注释信息这步会很顺利,也会很难,是技术上的限速环节。 platformMap <>'platformMap.txt') 我们根据上述帖子里面提到的platformMap,找到对应的R包是:
这里面的GPL6102是平台名称,我们在一开始已经在网页中知道了,如果忘记了,不要紧,这行代码可以自动获取 index = gset@annotation 我们也发现,表格中对应的 'illuminaHumanv2',跟我们需要的R包 'illuminaHumanv2.db',中间相差个“.db”,如果是单个,就自己加上,我们也可以使用代码的方法,自动获取: index = gset@annotation 安装R包并加载 ## 第一句是为了增加镜像 获取探针对应的symbol信息 ## 获取表达矩阵的行名,就是探针名称 制作probe2symbol转换文件 probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F) ![]() 有了这个文件,我们就可以进行探针转换了 有些平台是没有对应的R包的,可以自己下载平台信息来做 skr!GEO芯片数据的探针ID转换 有些GEO平台的探针转换比较麻烦 正则表达式是我们认识世界的哲学 6.探针转换与基因去重一个基因会对应对个探针,有些基因名称会是重复的,这些都需要处理。 library(dplyr) 现在数据变成这个样子 ![]() 而且,行数从原来的48701变成18962行。 7.差异分析这里是limma包的核心环节,也是理解以后其他差异分析工具如Deseq2的基础。 我们选取格式比较简单的。如果没有配对信息,差异分析这样做: design=model.matrix(~ group_list) 其中allDiff得到6799行。 pairinfo = factor(rep(1:18,2)) 得到的差异分析结果allDiff_pair有7650个,比直接做差异分析要多接近1000个。 ![]() 这里配对和非配对的差异分析,只相差一步,也就是是否有配对信息 ![]() 8.作图验证(非必要)差异分析的结果,配对和非配对理解起来比较抽象,我们用图片直观地感受一下。 首先基因名称需要在列,所以需要用t函数,行列转置。 data_plot = as.data.frame(t(exprSet)) ![]() 作图展示: library(ggplot2) ![]() 看起来杂乱一片,根本得不到有用信息,我们加入他的配对信息,再来作图: library(ggplot2) ![]() 突然间飞到枝头变凤凰。而这些基因就是普通差异分析会遗漏的部分,配对差异分析的时候他们又变得有意义。 再来看看真正有差异的基因吧,比如: library(ggplot2) ![]() 我们还可以批量地作图,这段不必要,可以自己一张张画,我只是展示一下如何批量画出差异最大的8个基因: library(dplyr) ![]() 而配对样本的展示,我们用TCGA的数据演示过很多次 9.后续分析差异分析后还有很多操作, library(pheatmap) ![]() 画热图的意义在哪? ![]() 第二看差异基因:差异基因提取出来的热图,就应当呈现横竖两条线,把表格分成四个象限,也就是差异基因有高有低,这才符合常识。 我们还可以画一个火山图 library(ggplot2) ![]() 火山图画起来简单,难的是如何定制化展示。比如在图上有不同的颜色,不同的点来表示数据。有什么意义呢?我其实理解的也不是很透彻,但是考虑到他的横坐标是变化倍数,纵坐标是p值的负数,那么p值越小,倍数越大的基因就会出现在左右上方。 数据是死的,而图形化展示是活的,火山图可能是大家觉得比较好的展示差异基因的方法。从这个图中我们还可以看出,两边是否对称,比如有的药物就是抑制基因转录的,那么加了药物的两组,可能会出现,上调和下调基因不对称的情形,这个从数据中可以分析出来,但是用图就一目了然了。 此外还有耳熟能详的GO分析,KEGG分析,我觉得都是名气很大用的很少的聚类方法。 Y叔的神包clusterprofiler可以做出特别好看的图。 比如GO分析,有三个组成部分,分别是CC,cellular compartment 细胞组分,BP,biological process,生物进程,MF,molecular function,分子功能。 ## 加载R包 ![]() 如何看这个图,如何解释?看图容易解释难。简单说来,一看颜色,二看大小。颜色越红p值越小越可信,点越打表示富集的基因比例越多,越有可能被验证。 举个例子,BP那一项中,有个叫neutrophil mediate immunity的GO通路被富集了,看字面意思应该是中性粒细胞介导的免疫反应,联系到这到实验设计,就是一开始的summary和文章部分,我觉得这个是合理的,因为他研究的就是长期慢性炎症对肥胖的影响,这里我们发现,炎症导致的免疫反应可能可以解释这个事情。还有很多其他的GO通路,需要自己根据背景去解释。去发现线索。 还有一个是广为人知的KEGG通路富集分析,Y叔的神包依然可以胜任。 EGG - enrichKEGG(gene= gene$ENTREZID, ![]() 这一看,全是大咖类的通路,每一个都是身价百万,说明啥问题,这篇文章的实验设计能够从多个方面解释。细胞周期,免疫反应,凋亡,等等,假如是我们自己的结果,我们就需要用自己的背景去思考,这么多通路是并联的还是串联的,是以此发生的级联反应,还是单兵作战呢?我觉得大部分情况下,我们会尝试用一个通路去说明。这是需要背景知识的。 假如我们对凋亡感兴趣,我们还可以看看我们有多少基因被富集在了凋亡通路,Y叔的神包clusterprofiler也是可以做的。而且做的很云淡风轻。 test <> ![]() 我们发现想看的凋亡通路牌号是,hsa04210,那我们就用下面的代码浏览一下,会跳出一个网页。 browseKEGG(EGG, 'hsa04210') ![]() 其中标红的就是富集到的通路。 那么这时候我们就会想,差异基因是有高有低的,而通路富集的时候缺失了这一个信息。一种解决方案是,分别用高表达和低表达基因去做这个分析。我觉得信息是冗余和偏误的,因为任何一个通路的激活,都是伴随着基因的高表达和低表达,本身这就是一个整体,没必要分开分析。 另外一种方法是,先总体分析,然后再标注出上调和下调的基因。那么Y叔的神包也可以做的。 ![]() 好像上了一点档次,留给大家自己去实现。 那么我们这个GEO的分析就讲完了。来看看我们获得了哪些图片呢? ![]() ![]() 只要是参加过线下课的同学,或者只要R语言基本能力过关的同学,轻而易举就能重复这个帖子的所有内容。 即使是零基础,安装这个帖子准备好工作环境,也能完成。 学习R语言,从这一课开始 如果继续做GESA,基因富集分析也是可以的,我自己觉得这个分析不需要人为设定差异基因的阈值,比如我们认为大于2倍就叫差异基因,这是不科学的,GSEA分析更加靠谱,可以代替GO和KEGG分析。 ![]() ![]() ![]() 假如我们通过以上的方法最终锁定了一两个基因,那么我们还可以对这个基因进行批量的相关性分析。 那么现在的问题是,技术上实现了,理论上该如何整合。R语言充其量只是个工具,神器不神,主要看人。假如你现在已经读过很多文献,做过很多实验,R语言一定有用,因为你首先想的会是用R语言解决自己的科学问题。 假如你现在还没有怎么读过文献,先不要急着要学习R语言。经常有学生加我微信,上来就说我要学习生信,我一听心里都很慌张,因为我知道,就我掌握的东西离生信还离得远,我就是从科研人员的角度功利地去运用生信知识解决问题。所以学习生信之前,先冷静几天,不要拿自己的爱好跟别人的专业去比较。 不分类别地去多看10分左右的文献可能是更好的策略,再高分的就高屋建瓴谈哲学对初学者不友好,只能学其行不能学其神,低分的容易把初学者带入科研皮条客的大坑,10分左右的文献,起手不高,论证严谨,大部分情况下比许多导师还可靠(这里用词很谨慎)。 有时候我会偏执地认为,读文献是性价比最高的提升科学思维的方法,一个科研人员的水平高低可以简单地用文献阅读量来比较。然后在回过头来看看这个帖子凝练自己的科学问题。有了科学问题我们就有了灵魂,才知道如何在海量资源中筛选有用的数据。 下周,我会更新一个长长的TCGA的帖子,这样我的2018年就这么过去了。
不需要任何赞赏,如果觉得这个帖子好看,就在右下角点击好看,这是微信的新功能。 |
|