step0-install_packagesstep1-download_data1.背景知识2.关于GEO中的几个文件说明A.2种familyfile(s)B.SeriesMatrixFile(s)C.GSE64634_RAW.tar3.关于下载镜像问题4.关于探针id转换idmap1包——针对bioconductor包附加小功能——对基因名字进行注释(`annoGene`)idmap2包——万能芯片探针ID注释平台包(提取soft文件)idmap3包——idmap1和idmap2都无法注释的平台AnnoProbe包5.整体代码这里抄step2-checkstep3-DEGstep4-anno-go-kegg真正代码detailedplot综合显示图(更推荐)step5-anno-GSEAstep6-anno-GSVAstep7-visualization致谢
简单安装R包和设置镜像代码:
而且我发现,最近NCBI更新后,GEO数据库也更新了!
点击后:
可以看到官网的代码和我们目前用的代码基本一致。
下面一个个的来说吧!
也就是SOFTformattedfamilyfile(s)和MINiMLformattedfamilyfile(s),通过字面,我们可以理解这是2种不同格式的探针说明文件。为什么我会这样说呢?因为我下载了soft文件来查看:
这时首行的内容:
就是告诉你这个数据是GPL570这个平台测序得到的:
中间有很多行关于GSM的,可能记录的是用到这个平台测序的sample:
接着就是一系列探针信息了:
我们可以通过在R种提取信息对我们得到的矩阵种探针做注释。
这里面包含了3部分资料:数据属性+患者信息+表达矩阵
数据属性在最前面的几行,和患者信息之间有一个空行,但是他们都是以“!”开头:
接下来是患者信息:
紧接着患者信息下一行会有个提示:!series_matrix_table_begin,然后下面就是表达矩阵信息了:
这个就是我们需要的表达矩阵信息了,矩阵中每一行是一个基因(也就是一个探针),每一列是一个样本(GSM)
根据常识猜了猜想:
X和Y应该代表着芯片上的位置,这个可能和探针有对应关系。
MEAN是信号的平均值,STDV是信号的标准差。
NPIXELS是像素点吧。
既然知道了这是原始数据,jimmy又给发了代码,感觉不学习下有点对不起自己,那就跟着代码过一遍吧:
可以看到,整个过程非常的简单,就是利用了oligo这个R包而已。读入所有的cel文件后,利用rma()函数将数据进行了normalization,得到的是一个ExpressionSet对象!!
后面会比较我们利用rawdata得到的结果和我们直接下载的SeriesMatrixFile(s)文件之间的区别!
通过学习,发现实际上jimmy是开发了一个R包,或者说包装了一个函数吧,如果不想了解具体原理,那么用法如下:
具体用法也非常简单:
所以可以猜到,应该是jimmy事先用循环的方式帮我们下载好了很多GEO数据,并做成了Rdata格式的文件!非常的良苦用心了!!
因为GEO数据矩阵中,横坐标都是探针的代号,如下图:
我们只看这些代号并不能理解具体是什么基因,于是这就需要我们做id转换了:将探针的代号转为基因symbol。
这里又要提到jimmy“开发”的一个R包了:
老规矩,还是先学习下:
(具体的R包名称是bioc_package后加“.db”)
因为很多时候,用户是找不到这些GPL平台对应的R包,或者下载安装困难,其实仅仅是需要探针与基因对应关系,没有必要去下载安装这几十个M的包。于是jimmy就开发了idmap1这个R包来帮我们:
安装方法(这里好像只能用方法1,因为在载入包的同时附带有一个变量p2s_df加载,如果用方法2和3没办法得到这个变量):
关于我们得到的ExpressionSet对象,可以通过gset@annotation得到我们需要的注释平台信息。
前面说到了这个变量p2s_df,像我这么喜欢资源的人,当然也要保存一份在本地了。大家有兴趣的可以自己write到本地保存。哈哈哈哈哈哈哈哈哈哈(jimmy应该不会打我吧)!
同样的,还是idmap1这个R包里的函数,如果安装过这个R包就无须再安装了,如果没有安装又想用这个功能,还真的没有办法,因为这些数据存在于这个包的自带变量中(humanGTF、mouseGTF、ratGTF):
R语言parse函数与eval函数的字符串转命令行及执行操作:
parse()函数能将字符串转换为表达式expression;eval()函数能对表达式求解
这次是根据[soft文件](#####A.2种familyfile(s))进行提取信息得到的!
如果是我自己来处理这样的文件,我应该会分2步:
结果证明这样有错误,但是,具体原因有空再去找吧。下面看正确的读入方法——借助GEOqueryR包工具:
一般来说,大家关心的其实就是探针的ID,以及基因的symbol列。有了这个变量后,就可以按照R语言基本操作来提取我们需要的信息了。
注意:我检查了得到的结果,里面存在有的探针ID对应2个基因,如下图:
虽然不知道这些代表着什么意思,但是,我将这个数据和bioconductor包里的hgu133plus2.db数据做了比较,结果是这样的:自己提取的结果中如果是一个ID对应2个基因,那么这个探针在bioconductor包里基本上找不到数据。而其他一个ID对应2个基因的结果均和bioconductor包一致。
当然了,我这只是单独来看一个平台的探针,而在idmap2jimmy已经帮我们整理好了,直接用就行了!!
安装方法:
(比较慢)
查看支持的平台:
如果我们拿到的soft注释文件中是序列信息,那么我们该怎么做呢?
应该是先将序列比对到参考基因组上,然后通过提取基因注释文件中的数据得到基因symbol!
而在idmap3包中,jimmy已经帮我们做好了!!说他宠粉也是真的了!!我都懒得做,他居然还写了个循环来完成了这个事。
使用方法:
heatmap_top1000_sd.png
cor_top500_mad.png
MA.png
heatmap_top200_DEG.png
GO系列结果过于冗余:
npc_VS_normal_dotplot_gene_diff_BP.png
npc_VS_normal_dotplot_gene_diff_CC.png
npc_VS_normal_dotplot_gene_diff_MF.png
npc_VS_normal_dotplot_gene_down_BP.png
npc_VS_normal_dotplot_gene_down_CC
npc_VS_normal_dotplot_gene_down_MF.png
npc_VS_normal_dotplot_gene_up_BP.png
npc_VS_normal_dotplot_gene_up_CC.png
npc_VS_normal_dotplot_gene_up_MF.png
gene_down_GO_all_barplot.png
gene_up_GO_all_barplot.png
因为用到的这个样本用GSVA没有得到显著性结果,所以没有图出来,具体也没有深究,有需要日后再仔细研究吧
气泡图:
Gene-ConceptNetwork图:每一个小蓝圈表示一个基因,其颜色表示FC值,每个KEGGterm圈的大小由里面包含基因的数目决定。
成环:更炫酷了,但是感觉图形展示不方便
不成环:信息展示更有力吧
EnrichmentMap图:这里和上面的图类似,只不过不再显示具体的基因,而是直接画出每个term和term之间的关系,每个圆圈代表着一个term,圆圈大小代表着有多少个基因,颜色表示p值。
如果term和term之间有共同的基因,那么就会连接起来,聚在一起。
Heatmap-likefunctionalclassification:
和我们常规的热图不太像,这里纵轴是每个KEGG通路,横轴是涉及到的基因。颜色表示FC值。
上面所有的代码都来自生信技能树曾老板jimmy的帮助,同时我在测试运行的过程中又进行了部分改进和增补。
就用曾老板亲自编辑的感谢词来感谢吧:
FatYangthankDr.JianmingZeng(UniversityofMacau),andallthemembersofhisbioinformaticsteam,biotrainee,forgenerouslysharingtheirexperienceandcodes!