其他
GEO数据挖掘流程——代码版（方便抄袭）
GEO数据挖掘流程——代码版（方便抄袭）
微信公众号：生信小知识
关注可了解更多的教程及单细胞知识。问题或建议，请公众号留言;
内容目录
step0-install_packagesstep1-download_data1.背景知识2.关于GEO中的几个文件说明A. 2种family file(s)B. Series Matrix File(s)C. GSE64634_RAW.tar3.关于下载镜像问题4.关于探针id转换idmap1包——针对bioconductor包附加小功能——对基因名字进行注释（`annoGene`）idmap2包——万能芯片探针ID注释平台包（提取soft文件）idmap3包——idmap1和idmap2都无法注释的平台AnnoProbe包5. 整体代码这里抄step2-checkstep3-DEGstep4-anno-go-kegg真正代码detailed plot综合显示图（更推荐）step5-anno-GSEAstep6-anno-GSVAstep7-visualization致谢
step0-install_packages
简单安装R包和设置镜像代码：
1# 镜像设置
2if (T) {
3 rm(list = ls())
4 options()$repos
5 options()$BioC_mirror
6 options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
7 options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
8 options()$repos
9 options()$BioC_mirror
10}
11
12# https://bioconductor.org/packages/release/bioc/html/GEOquery.html
13# BiocManager安装的R包
14if (T) {
15 if (!requireNamespace("BiocManager", quietly = TRUE))
16 install.packages("BiocManager")
17
18 pkgs=c("KEGG.db","GSEABase","GSVA","clusterProfiler","GEOquery","limma","impute","genefu","org.Hs.eg.db","hgu133plus2.db")
19 for (pkg in pkgs) {
20 if (!requireNamespace(pkg, quietly = TRUE))
21 BiocManager::install(pkg,ask = F,update = F)
22 }
23}
24
25# 直接通过install.package函数安装的R包
26if (T) {
27 options()$repos
28 pkgs=c("FactoMineR", "factoextra",'WGCNA',"ggplot2", "pheatmap","ggpubr")
29 for (pkg in pkgs) {
30 if (!requireNamespace(pkg, quietly = TRUE))
31 install.packages(pkg)
32 }
33}
34
35
36# 载入R包
37if (T) {
38 library("FactoMineR")
39 library("factoextra")
40 library(GSEABase)
41 library(GSVA)
42 library(clusterProfiler)
43 library(genefu)
44 library(ggplot2)
45 library(ggpubr)
46 library(hgu133plus2.db)
47 library(limma)
48 library(org.Hs.eg.db)
49 library(pheatmap)
50}
step1-download_data
1.背景知识
这里通过使用GEOqueryR包来帮助我们下载数据，比手动登录GEO数据库一个个点击要方便很多！
而且我发现，最近NCBI更新后，GEO数据库也更新了！
点击后：
点击后：
可以看到官网的代码和我们目前用的代码基本一致。
2.关于GEO中的几个文件说明
下面一个个的来说吧！
A. 2种family file(s)
也就是SOFT formatted family file(s)和MINiML formatted family file(s)，通过字面，我们可以理解这是2种不同格式的探针说明文件。为什么我会这样说呢？因为我下载了soft文件来查看：
这时首行的内容：
就是告诉你这个数据是GPL570这个平台测序得到的：
中间有很多行关于GSM的，可能记录的是用到这个平台测序的sample：
接着就是一系列探针信息了：
我们可以通过在R种提取信息对我们得到的矩阵种探针做注释。
B. Series Matrix File(s)
这里面包含了3部分资料：数据属性+患者信息+表达矩阵
数据属性在最前面的几行，和患者信息之间有一个空行，但是他们都是以“!”开头：
接下来是患者信息：
紧接着患者信息下一行会有个提示：!series_matrix_table_begin，然后下面就是表达矩阵信息了：
这个就是我们需要的表达矩阵信息了，矩阵中每一行是一个基因（也就是一个探针），每一列是一个样本（GSM）
C. GSE64634_RAW.tar
这个我一开始不知道是什么东西，后来问了问jimmy，然后他给我发了个资料你要挖的公共数据集作者上传了错误的表达矩阵肿么办。大致瞅了瞅，知道这个应该是芯片的原始数据，然后也把这个文件给下载下来看了看：
根据常识猜了猜想：
X和Y应该代表着芯片上的位置，这个可能和探针有对应关系。
MEAN是信号的平均值，STDV是信号的标准差。
NPIXELS是像素点吧。
既然知道了这是原始数据，jimmy又给发了代码，感觉不学习下有点对不起自己，那就跟着代码过一遍吧：
1# BiocManager::install(c( 'oligo' ),ask = F,update = F)
2library(oligo)
3# BiocManager::install(c( 'pd.hg.u133.plus.2' ),ask = F,update = F)
4library(pd.hg.u133.plus.2)
5
6# 读入cel文件数据
7celFiles <- list.celfiles('./',listGzipped = T)
8celFiles
9affyRaw <- read.celfiles( celFiles )
10# 提取矩阵并做normalization
11eset <- rma(affyRaw)
12eset
13# 检查数据
14dat=exprs(eset)
15dim(dat)
16boxplot(dat,las=2)
可以看到，整个过程非常的简单，就是利用了oligo这个R包而已。读入所有的cel文件后，利用rma()函数将数据进行了normalization，得到的是一个ExpressionSet对象！！
后面会比较我们利用rawdata得到的结果和我们直接下载的Series Matrix File(s)文件之间的区别！
3.关于下载镜像问题
jimmy曾经发表过GEO数据库中国区镜像横空出世推文，因为我每次下载都是用vpn，但是怕以后万一没有vpn了，所以还是学习下，以防万一嘛！再说了，可以学习下新知识，何乐而不为呢？
通过学习，发现实际上jimmy是开发了一个R包，或者说包装了一个函数吧，如果不想了解具体原理，那么用法如下：
1# 安装方法1：
2library(devtools)
3install_github("jmzeng1314/GEOmirror")
4library(GEOmirror)
5# 安装方法2：
6source('http://raw.githubusercontent.com/jmzeng1314/GEOmirror/master/R/geoChina.R')
7# 安装方法3（源码）：
8geoChina <- function(gse='GSE2546',mirror='tercent'){
9 suppressPackageStartupMessages(library(GEOquery))
10 options(download.file.method="libcurl")
11 # eSet=getGEO('GSE2546', destdir=".", AnnotGPL = F, getGPL = F)
12 # http://49.235.27.111/GEOmirror/GSE2nnn/GSE2546_eSet.Rdata
13 # gse='GSE2546';mirror='tercent'
14 gse=toupper(gse)
15 down=ifelse(as.numeric(gsub('GSE','',gse))<1000,
16 paste0('/GEOmirror/GSEnnn/',gse,
17 '_eSet.Rdata'),
18 paste0('/GEOmirror/',
19 gsub('[0-9][0-9][0-9]$','nnn',gse),'/',gse,
20 '_eSet.Rdata'))
21
22 if(mirror=='tercent'){
23 up='http://49.235.27.111'
24 }
25 tpf=paste0(gse, '_eSet.Rdata')
26 download.file(paste0(up,down),tpf,mode = "wb")
27 suppressWarnings(load(tpf))
28 return(gset)
29}
具体用法也非常简单：
1eSet=geoChina('GSE2345')
通过安装方法3（源码）我们可以看出这个包的原理：
通过访问下载 http://49.235.27.111/GEOmirror/GSE2nnn/GSE2546_eSet.Rdata
所以可以猜到，应该是jimmy事先用循环的方式帮我们下载好了很多GEO数据，并做成了Rdata格式的文件！非常的良苦用心了！！
4.关于探针id转换
因为GEO数据矩阵中，横坐标都是探针的代号，如下图：
我们只看这些代号并不能理解具体是什么基因，于是这就需要我们做id转换了：将探针的代号转为基因symbol。
idmap1包——针对bioconductor包
这里又要提到jimmy“开发”的一个R包了：
第一个万能芯片探针ID注释平台R包——37个芯片平台
老规矩，还是先学习下：
一般重要的芯片在R的bioconductor里面都是有包的。但是需要我们单独下载对应的包。具体关系如下：
（具体的R包名称是bioc_package后加“.db”）
1gpl bioc_package title
2GPL201 hgfocus [HG-Focus] Affymetrix Human HG-Focus Target Array
3GPL96 hgu133a [HG-U133A] Affymetrix Human Genome U133A Array
4GPL571 hgu133a2 [HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array
5GPL97 hgu133b [HG-U133B] Affymetrix Human Genome U133B Array
6GPL570 hgu133plus2 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
7GPL13667 hgu219 [HG-U219] Affymetrix Human Genome U219 Array
8GPL8300 hgu95av2 [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array
9GPL91 hgu95av2 [HG_U95A] Affymetrix Human Genome U95A Array
10GPL92 hgu95b [HG_U95B] Affymetrix Human Genome U95B Array
11GPL93 hgu95c [HG_U95C] Affymetrix Human Genome U95C Array
12GPL94 hgu95d [HG_U95D] Affymetrix Human Genome U95D Array
13GPL95 hgu95e [HG_U95E] Affymetrix Human Genome U95E Array
14GPL887 hgug4110b Agilent-012097 Human 1A Microarray (V2) G4110B (Feature Number version)
15GPL886 hgug4111a Agilent-011871 Human 1B Microarray G4111A (Feature Number version)
16GPL1708 hgug4112a Agilent-012391 Whole Human Genome Oligo Microarray G4112A (Feature Number version)
17GPL13497 HsAgilentDesign026652 Agilent-026652 Whole Human Genome Microarray 4x44K v2 (Probe Name version)
18GPL6244 hugene10sttranscriptcluster [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]
19GPL11532 hugene11sttranscriptcluster [HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array [transcript (gene) version]
20GPL6097 illuminaHumanv1 Illumina human-6 v1.0 expression beadchip
21GPL6102 illuminaHumanv2 Illumina human-6 v2.0 expression beadchip
22GPL6947 illuminaHumanv3 Illumina HumanHT-12 V3.0 expression beadchip
23GPL10558 illuminaHumanv4 Illumina HumanHT-12 V4.0 expression beadchip
24GPL6885 illuminaMousev2 Illumina MouseRef-8 v2.0 expression beadchip
25GPL81 mgu74av2 [MG_U74Av2] Affymetrix Murine Genome U74A Version 2 Array
26GPL82 mgu74bv2 [MG_U74Bv2] Affymetrix Murine Genome U74B Version 2 Array
27GPL83 mgu74cv2 [MG_U74Cv2] Affymetrix Murine Genome U74 Version 2 Array
28GPL339 moe430a [MOE430A] Affymetrix Mouse Expression 430A Array
29GPL6246 mogene10sttranscriptcluster [MoGene-1_0-st] Affymetrix Mouse Gene 1.0 ST Array [transcript (gene) version]
30GPL340 mouse4302 [MOE430B] Affymetrix Mouse Expression 430B Array
31GPL1261 mouse430a2 [Mouse430_2] Affymetrix Mouse Genome 430 2.0 Array
32GPL8321 mouse430a2 [Mouse430A_2] Affymetrix Mouse Genome 430A 2.0 Array
因为很多时候，用户是找不到这些GPL平台对应的R包，或者下载安装困难，其实仅仅是需要探针与基因对应关系，没有必要去下载安装这几十个M的包。于是jimmy就开发了idmap1这个R包来帮我们：
安装方法（这里好像只能用方法1，因为在载入包的同时附带有一个变量p2s_df加载，如果用方法2和3没办法得到这个变量）：
也可以直接下载：https://codeload.github.com/jmzeng1314/idmap1/legacy.tar.gz/master
1# 安装方法1：
2library(devtools)
3install_github("jmzeng1314/idmap1")
4library(idmap1)
5# 安装方法2：
6source('https://raw.githubusercontent.com/jmzeng1314/idmap1/master/R/getIDs.R')
7# 安装方法3（源码）：
8getIDs <- function(gpl){
9 # gpl='gpl570'
10 gpl=toupper(gpl)
11 if(!gpl %in% unique(p2s_df$gpl)){
12 stop('your gpl is not in our annotation list.')
13 }
14 return(p2s_df[p2s_df$gpl==gpl,])
15}
具体用法也非常简单：
1ids=getIDs('gpl570')
2head(ids)
关于我们得到的ExpressionSet对象，可以通过gset@annotation得到我们需要的注释平台信息。
前面说到了这个变量p2s_df，像我这么喜欢资源的人，当然也要保存一份在本地了。大家有兴趣的可以自己write到本地保存。哈哈哈哈哈哈哈哈哈哈（jimmy应该不会打我吧）！
附加小功能——对基因名字进行注释（`annoGene`）
同样的，还是idmap1这个R包里的函数，如果安装过这个R包就无须再安装了，如果没有安装又想用这个功能，还真的没有办法，因为这些数据存在于这个包的自带变量中（humanGTF、mouseGTF、ratGTF）:
1# 安装方法2(无效)：
2source('https://raw.githubusercontent.com/jmzeng1314/idmap1/master/R/annoGene.R')
3# 安装方法3（源码）(无效)：
4annoGene <- function(IDs,ID_type,species='human',out_file){
5 if(length(unique(IDs))<1){
6 stop("You should give me some genes to be annotated!!!")
7 }
8 if(!ID_type %in% c("ENSEMBL" ,"SYMBOL")){
9 stop("We only accept ENSEMBL or SYMBOL !!!")
10 }
11 if(species=='human'){
12 GTF <- humanGTF
13 }else if(species=='mouse'){
14 GTF <- mouseGTF
15 }else if(species=='rat'){
16 GTF <- ratGTF
17 }else{
18 stop("We only accept human or mouse, or rat, ")
19 }
20 res <- GTF[eval(parse(text=paste0("GTF$",ID_type))) %in% IDs, ]
21
22 missIds <- IDs[!(IDs %in% eval(parse(text=paste0("res$",ID_type))))]
23 missIdsPercentage = round((length(missIds)/length(IDs))*100,2)
24 if(length(missIds)!=0){
25 warning(
26 paste0(missIdsPercentage ,"% of input IDs are fail to annotate... ")
27 # 5.29% of input gene IDs are fail to map...
28 )
29 }
30 if (hasArg(out_file)) {
31 results=res
32 if(grepl('.html$',out_file)){
33 Ensembl_prefix <- "https://asia.ensembl.org/Homo_sapiens/Gene/Summary?g="
34 href = paste0(Ensembl_prefix, results$ENSEMBL)
35 results$ENSEMBL = paste0("<b><a target=\"_black\" href=", shQuote(href), ">", results$ENSEMBL, "</a></b>")
36
37 symbol_prefix <- "http://www.ncbi.nlm.nih.gov/gene?term="
38 href = paste0(symbol_prefix, results$SYMBOL)
39 results$SYMBOL = paste0("<b><a target=\"_black\" href=", shQuote(href), ">", results$SYMBOL, "</a></b>")
40
41 y <- DT::datatable(results, escape = F, rownames = F)
42 DT::saveWidget(y,file = out_file)
43 }else if(grepl('.csv$',out_file)){
44 write.csv(results,file =out_file )
45 }else{
46 stop("We only accept csv or html format !!!")
47 }
48
49 }
50 return(res)
51}
这里补充学习了2个函数：
R语言parse函数与eval函数的字符串转命令行及执行操作：
parse()函数能将字符串转换为表达式expression；eval()函数能对表达式求解
idmap2包——万能芯片探针ID注释平台包（提取soft文件）
第二个万能芯片探针ID注释平台R包——有122个GPL
这次是根据[soft文件](#####A. 2种family file(s))进行提取信息得到的！
如果是我自己来处理这样的文件，我应该会分2步：
用linux中的grep -v命令将表头以“#”和“！”开头的行去掉
1grep -v ^! GPL570.soft |grep -v ^# |grep -v ^\^ >GPL570.txt
用R读入数据
1tmp=read.table("data/GPL570.txt",header = T,sep = '\t',fill = T)
结果证明这样有错误，但是，具体原因有空再去找吧。下面看正确的读入方法——借助GEOqueryR包工具：
1library(GEOquery)
2gpl <- getGEO('GPL570', destdir="data/")
3GPL=Table(gpl)
一般来说，大家关心的其实就是探针的ID，以及基因的symbol列。有了这个变量后，就可以按照R语言基本操作来提取我们需要的信息了。
注意：我检查了得到的结果，里面存在有的探针ID对应2个基因，如下图：
虽然不知道这些代表着什么意思，但是，我将这个数据和bioconductor包里的hgu133plus2.db数据做了比较，结果是这样的：自己提取的结果中如果是一个ID对应2个基因，那么这个探针在bioconductor包里基本上找不到数据。而其他一个ID对应2个基因的结果均和bioconductor包一致。
当然了，我这只是单独来看一个平台的探针，而在idmap2jimmy已经帮我们整理好了，直接用就行了！！
安装方法：
（比较慢）
也可以直接下载：https://codeload.github.com/jmzeng1314/idmap2/legacy.tar.gz/master
1library(devtools)
2install_github("jmzeng1314/idmap2")
3library(idmap2)
使用方法：
1library(idmap2)
2ids=get_soft_IDs('gpl570')
查看支持的平台：
1gpl_list[,1:4] # gpl_list是R包自带的变量
idmap3包——idmap1和idmap2都无法注释的平台
第三个万能芯片探针ID注释平台R包——idmap1和idmap2都无法注释的平台
如果我们拿到的soft注释文件中是序列信息，那么我们该怎么做呢？
应该是先将序列比对到参考基因组上，然后通过提取基因注释文件中的数据得到基因symbol！
而在idmap3包中，jimmy已经帮我们做好了！！说他宠粉也是真的了！！我都懒得做，他居然还写了个循环来完成了这个事。
安装方法：
（比较慢）
也可以直接下载：https://codeload.github.com/jmzeng1314/idmap3/legacy.tar.gz/master
1library(devtools)
2install_github("jmzeng1314/idmap3")
3library(idmap3)
使用方法：
1library(idmap3)
2ids=idmap3::get_pipe_IDs('GPL21827')
查看支持的平台：
1gpl_list[,1:4] # gpl_list是R包自带的变量
AnnoProbe包
芯片探针序列的基因注释已经无需你自己亲自做了——一个汇总
安装方法：
https://codeload.github.com/jmzeng1314/AnnoProbe/legacy.tar.gz/master
1library(devtools)
2install_github("jmzeng1314/AnnoProbe")
3library(AnnoProbe)
使用方法：
1gpl='GPL21827'
2probe2gene=idmap(gpl,type = 'pipe')
5. 整体代码这里抄
1rm(list = ls()) ## 魔幻操作，一键清空~
2options(stringsAsFactors = F)
3# 下面代码只需要修改file的名字
4# 下载的文件会自动保存到./data/下
5# 得到的gset是一个ExpressionSet对象
6
7
8# 只需要修改file的名字，即可下载得到相应的geo数据
9# 获取患者信息，需要自行修改
10file='GSE64634'
11# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64634
12
13if (T) {
14 # 数据下载
15 if (T) {
16 library(GEOquery)
17 # 这个包需要注意两个配置，一般来说自动化的配置是足够的。
18 #Setting options('download.file.method.GEOquery'='auto')
19 #Setting options('GEOquery.inmemory.gpl'=FALSE)
20 fdata=paste0(file,"_eSet.Rdata")
21 fpath=paste0("data/",fdata)
22 if(!file.exists(fpath)){
23 gset <- getGEO(file, destdir="data/",
24 AnnotGPL = F, ## 注释文件
25 getGPL = F) ## 平台文件
26 gset=gset[[1]]
27 save(gset,file=fpath) ## 保存到本地
28 }
29 load(fpath)
30 }
31 gset
32
33 # 获取患者信息，这里需要自行修改
34 if (T) {
35 pd=pData(gset)# 根据disease state:ch1一列得知分组信息
36 group_list=c(rep('normal',4),rep('npc',12))# nasopharyngeal carcinoma NPC 鼻咽癌
37 table(group_list)
38 }
39
40 # 对数据进行normalization
41 if (T) {
42 dat=exprs(gset)
43 dim(dat)
44
45 dat[1:4,1:4]
46 boxplot(dat,las=2)
47 dat=dat[apply(dat,1,sd)>0,]# 去除都是0的探针
48 dat[dat<0]=1
49 boxplot(dat,las=2)
50
51 #dat=log2(dat+1)
52 #boxplot(dat,las=2)
53 library(limma)
54 dat=normalizeBetweenArrays(dat)
55 boxplot(dat,las=2)
56 }
57
58
59 # 探针注释2种方法，推荐方法2
60 # 方法1:比较麻烦，而且不方便，一般不用这种方法
61 if(F){
62 library(GEOquery)
63 #Download GPL file, put it in the current directory, and load it:
64 gpl <- getGEO('GPL570', destdir="data/")
65 GPL=Table(gpl)
66 colnames(GPL)
67 head(GPL[,c(1,11)]) ## you need to check this , which column do you need
68 probe2gene=GPL[,c(1,11)]
69 head(probe2gene)
70 save(probe2gene,file='probe2gene.Rdata')
71 }
72
73 # 方法2:用hgu133plus2.db这个R包比较方便
74 if (T) {
75 library(hgu133plus2.db)
76 ids=toTable(hgu133plus2SYMBOL)
77 head(ids)
78 ids=ids[ids$symbol != '',]
79 ids=ids[ids$probe_id %in% rownames(dat),]# 过滤没法注释的探针
80
81 dat[1:4,1:4]
82 dat=dat[ids$probe_id,]# 调整顺序，让dat的顺序和ids中的一致
83
84 ids$median=apply(dat,1,median)
85 ids=ids[order(ids$symbol,ids$median,decreasing = T),]# 按照基因名、中位数大小排序
86 ids=ids[!duplicated(ids$symbol),]# 只保留相同symbol中中位数最大的探针
87 dat=dat[ids$probe_id,]# 调整顺序，让dat的顺序和ids中的一致
88 rownames(dat)=ids$symbol# id转换
89 dat[1:4,1:4]
90 }
91}
92
93save(dat,group_list,file = 'data/step1-output.Rdata')
step2-check
画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转置。格式要求data.frame。
关于scale：在基因表达量的数据中，有些基因表达量极低，有些基因表达量极高，因此把每个基因在不同处理和重复中的数据转换为平均值为0，方差为1的数据，所以在这里也需要先转置。
1# PCA检查
2# PCA，图会保存在pic/all_samples_PCA.png
3if (T) {
4 # 载入数据
5 if (T) {
6 rm(list = ls()) ## 魔幻操作，一键清空~
7 options(stringsAsFactors = F)
8
9 load(file = 'data/step1-output.Rdata')
10 print(table(group_list))
11 # 每次都要检测数据
12 dat[1:4,1:4]
13 }
14
15 # PCA，图会保存在pic/all_samples_PCA.png
16 if (T) {
17 exprSet=dat
18 # 过滤标准
19 if (T) {
20 dim(exprSet)
21 # 过滤标准需要修改,目前是保留每一个基因在>5个人中表达量>1
22 exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
23 boxplot(apply(exprSet,1 , mean))
24 dim(exprSet)
25 # 过滤标准需要修改,目前是去除每一个基因在>5个人中表达量>12
26 exprSet=exprSet[!apply(exprSet,1, function(x) sum(x>12) > 5),]
27 dim(exprSet)
28 }
29 # 去除文库大小差异normalization，同时利用log做scale
30 exprSet=log(edgeR::cpm(exprSet)+1)
31 dat=exprSet
32 dat=as.data.frame(t(dat)) # 画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换。格式要求data.frame
33 library("FactoMineR")# 计算PCA
34 library("factoextra")# 画图展示
35
36 dat.pca <- PCA(dat, graph = F)
37 # fviz_pca_ind按样本 fviz_pca_var按基因
38 fviz_pca_ind(dat.pca,
39 geom.ind = "point", # c("point", "text)2选1
40 col.ind = group_list, # color by groups
41 # palette = c("#00AFBB", "#E7B800"),# 自定义颜色
42 addEllipses = T, # 加圆圈
43 legend.title = "Groups"# 图例名称
44 )
45 ggsave('pic/all_samples_PCA.png')
46 }
47}
48
49
50# 挑选1000个SD最大的基因画表达量热图
51# 结果图存放在pic/heatmap_top1000_sd.png中
52if (T) {
53 # 载入数据
54 if (T) {
55 rm(list = ls())
56 load(file = 'data/step1-output.Rdata')
57 dat[1:4,1:4]
58 }
59
60 # 挑选1000个SD最大的基因画表达量热图
61 # 结果图存放在pic/heatmap_top1000_sd.png中
62 if (T) {
63 cg=names(tail(sort(apply(dat,1,sd)),1000))# 找到SD最大的1000个基因
64 library(pheatmap)
65 headmap.dat=dat[cg,]
66 # 把每个基因在不同处理和重复中的数据转换为平均值为0，
67 # 方差为1的数据，所以在这里也需要先转置(针对基因)。
68 headmap.dat=t(scale(t(headmap.dat)))
69 headmap.dat[headmap.dat>2]=2
70 headmap.dat[headmap.dat< -2]= -2
71
72 # 分组信息设置
73 ac=data.frame(group=group_list)
74 rownames(ac)=colnames(headmap.dat)
75
76 ## 可以看到TNBC具有一定的异质性，拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
77 pheatmap(headmap.dat,show_colnames =T,show_rownames = F,
78 annotation_col=ac,
79 filename = 'pic/heatmap_top1000_sd.png')
80 }
81}
82
83
84# 过滤标准需要修改
85# 利用top500_mad基因画相关性热图热图
86# 结果图存放在pic/cor_top500_mad.png中
87if (T) {
88 # 导入数据
89 if (T) {
90 rm(list = ls())
91 load(file = 'data/step1-output.Rdata')
92 dat[1:4,1:4]
93 exprSet=dat
94 }
95
96 # 利用所有基因画相关性热图
97 if (T) {
98 ac=data.frame(group=group_list)
99 rownames(ac)=colnames(exprSet)
100 pheatmap::pheatmap(cor(exprSet),
101 annotation_col = ac,
102 show_rownames = F,
103 filename = 'pic/cor_all_genes.png')
104 }
105
106 # 过滤标准
107 if (T) {
108 dim(exprSet)
109 # 过滤标准需要修改,目前是保留每一个基因在>5个人中表达量>1
110 exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
111 boxplot(apply(exprSet,1 , mean))
112 dim(exprSet)
113 # 过滤标准需要修改,目前是去除每一个基因在>5个人中表达量>12
114 exprSet=exprSet[!apply(exprSet,1, function(x) sum(x>12) > 5),]
115 dim(exprSet)
116 }
117
118 # 数据normalization处理
119 if (T) {
120 # 去除文库大小差异normalization，同时利用log做scale
121 exprSet=log(edgeR::cpm(exprSet)+1)
122 # 取top500 mad 的基因画图
123 exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
124 }
125
126 # 利用top500_mad基因画相关性热图热图
127 if (T) {
128 M=cor(exprSet)
129 pheatmap::pheatmap(M,
130 show_rownames = F,
131 annotation_col = ac,
132 filename = 'pic/cor_top500_mad.png')
133 }
134}
all_samples_PCA.png
heatmap_top1000_sd.png
cor_top500_mad.png
step3-DEG
原始版火山图：
plot(logFC,-log10(P.Value))原始版MA图：
plot(AveExpr,logFC)
1# 载入数据,检查数据
2if (T) {
3 rm(list = ls())
4 options(stringsAsFactors = F)
5 library(ggpubr)
6 load(file = 'data/step1-output.Rdata')
7 # 每次都要检测数据
8 dat[1:4,1:4]
9 table(group_list)
10 boxplot(dat[1,]~group_list) #按照group_list分组画箱线图
11
12 # boxplot的美化版
13 bplot=function(g){
14 df=data.frame(gene=g,stage=group_list)
15 p <- ggboxplot(df, x = "stage", y = "gene",
16 color = "stage", palette = "jco",
17 add = "jitter")
18 # Add p-value
19 p + stat_compare_means()
20 }
21}
22# 利用定义好的函数检查数据
23bplot(dat[1,])
24bplot(dat[2,])
25bplot(dat[3,])
26bplot(dat[4,])
27
28
29# limma
30library(limma)
31# 方法1：不制作比较矩阵，简单
32# 但是做不到随心所欲的指定任意两组进行比较
33if (T) {
34 design=model.matrix(~factor( group_list ))
35 fit=lmFit(dat,design)
36 fit2=eBayes(fit)
37 ## 上面是limma包用法的一种方式
38 options(digits = 4) #设置全局的数字有效位数为4
39 topTable(fit2,coef=2,adjust='BH')
40}
41
42# 方法2：制作比较矩阵
43# 可以随心所欲的指定任意两组进行比较
44if (T) {
45 design <- model.matrix(~0+factor(group_list))
46 colnames(design)=levels(factor(group_list))
47 head(design)
48 exprSet=dat
49 rownames(design)=colnames(exprSet)
50 design
51
52 # 比较矩阵
53 # 这个矩阵声明，我们要把 npc 组跟 Normal 进行差异分析比较
54 contrast.matrix<-makeContrasts("npc-normal",
55 levels = design)
56 contrast.matrix
57
58 deg = function(exprSet,design,contrast.matrix){
59 ##step1
60 fit <- lmFit(exprSet,design)
61 ##step2
62 fit2 <- contrasts.fit(fit, contrast.matrix)
63
64 fit2 <- eBayes(fit2) ## default no trend !!!
65 ##eBayes() with trend=TRUE
66
67 ##step3
68 # 有了比较矩阵后，coef=1，而number=Inf是把所有结果都打印出来
69 tempOutput = topTable(fit2, coef=1, number =Inf)
70 nrDEG = na.omit(tempOutput)
71 #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
72 head(nrDEG)
73 return(nrDEG)
74 }
75
76 deg = deg(exprSet,design,contrast.matrix)
77
78 head(deg)
79}
80
81
82save(deg,file = 'data/deg.Rdata')
83
84load(file = 'data/deg.Rdata')
85head(deg)
86bplot(dat[rownames(deg)[1],])
87## for volcano and MA plot
88# 结果存放在pic/volcano.png和pic/MA.png
89if(T){
90 nrDEG=deg
91 head(nrDEG)
92 attach(nrDEG)
93 # 原始版火山图
94 plot(logFC,-log10(P.Value))
95 library(ggpubr)
96 df=nrDEG
97 df$y= -log10(P.Value)
98 ggscatter(df, x = "logFC", y = "y",size=0.5)
99 # 定义logFC=2为阈值
100 df$state=ifelse(df$P.Value>0.01,'stable',
101 ifelse( df$logFC >2,'up',
102 ifelse( df$logFC < -2,'down','stable') )
103 )
104 table(df$state)
105 df$name=rownames(df)
106 head(df)
107 ggscatter(df, x = "logFC", y = "y",size=0.5,color = 'state')
108 ggscatter(df, x = "logFC", y = "y", color = "state",size = 0.5,
109 label = "name", repel = T,
110 #label.select = rownames(df)[df$state != 'stable'] ,
111 label.select = c('TTC9', 'AQP3', 'CXCL11','PTGS2'), #挑选一些基因在图中显示出来
112 palette = c("#00AFBB", "#E7B800", "#FC4E07"))
113 ggsave('pic/volcano.png')
114
115 # MA图
116 ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
117 df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
118 ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
119 table(df$p_c )
120 ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
121 palette = c("green", "red", "black") )
122 ggsave('pic/MA.png')
123}
124
125## for heatmap
126if(T){
127 load(file = 'data/step1-output.Rdata')
128 # 每次都要检测数据
129 dat[1:4,1:4]
130 table(group_list)
131 x=deg$logFC
132 names(x)=rownames(deg)
133
134 # cg中存放着变化上升和下降的前100个基因名
135 cg=c(names(head(sort(x),100)),
136 names(tail(sort(x),100)))
137 library(pheatmap)
138 n=t(scale(t(dat[cg,])))
139
140 n[n>2]=2
141 n[n< -2]= -2
142 n[1:4,1:4]
143 pheatmap(n,show_colnames =F,show_rownames = F)
144 ac=data.frame(group=group_list)
145 rownames(ac)=colnames(n) #将ac的行名也就分组信息（是‘no TNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
146 pheatmap(n,show_colnames =F,
147 show_rownames = F,
148 cluster_cols = F,
149 annotation_col=ac,filename = 'pic/heatmap_top200_DEG.png') #列名注释信息为ac即分组信息
150
151
152}
153
154write.csv(deg,file = 'data/deg.csv')
volcano.png
MA.png
heatmap_top200_DEG.png
step4-anno-go-kegg
关于基因名和ID之间的各种转换：
bitr(geneID, fromType, toType, OrgDb, drop = TRUE)例如：
1head(unique(deg$symbol))
2# [1] "LOC388780" "SLC6A6" "RGS22" "AKR1D1" "AOC1" "FOXJ1"
3bitr(unique(deg$symbol), fromType = "SYMBOL",
4 toType = c( "ENTREZID"),
5 OrgDb = org.Hs.eg.db)
这里将KEGG和GO富集分析函数自定义为3个函数——
run_kegg、run_go、kegg_plot每次运行前先运行下面代码：
1# KEGG pathway analysis
2# 具体结果数据在data/下，图在pic/下
3run_kegg <- function(gene_up,gene_down,geneList=F,pro='test'){
4 gene_up=unique(gene_up)
5 gene_down=unique(gene_down)
6 gene_diff=unique(c(gene_up,gene_down))
7 ### over-representation test
8 # 下面把3个基因集分开做超几何分布检验
9 # 首先是上调基因集。
10 kk.up <- enrichKEGG(gene = gene_up,
11 organism = 'hsa',
12 #universe = gene_all,
13 pvalueCutoff = 0.9,
14 qvalueCutoff =0.9)
15 head(kk.up)[,1:6]
16 kk=kk.up
17 dotplot(kk)
18 barplot(kk)
19 kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
20 head(kk)
21 write.csv(kk@result,paste0('data/',pro,'_kk.up.csv'))
22
23 # 首先是下调基因集。
24 kk.down <- enrichKEGG(gene = gene_down,
25 organism = 'hsa',
26 #universe = gene_all,
27 pvalueCutoff = 0.9,
28 qvalueCutoff =0.9)
29 head(kk.down)[,1:6]
30 kk=kk.down
31 dotplot(kk)
32 barplot(kk)
33 kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
34 write.csv(kk@result,paste0('data/',pro,'_kk.down.csv'))
35
36 # 最后是上下调合并后的基因集。
37 kk.diff <- enrichKEGG(gene = gene_diff,
38 organism = 'hsa',
39 pvalueCutoff = 0.05)
40 head(kk.diff)[,1:6]
41 kk=kk.diff
42 dotplot(kk)
43 barplot(kk)
44 kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
45 write.csv(kk@result,paste0('data/',pro,'_kk.diff.csv'))
46
47
48 kegg_down_dt <- as.data.frame(kk.down)
49 kegg_up_dt <- as.data.frame(kk.up)
50 down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.01,];down_kegg$group=-1
51 up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.01,];up_kegg$group=1
52 #画图设置, 这个图很丑，大家可以自行修改。
53 kegg_plot(up_kegg,down_kegg)
54
55 ggsave(g_kegg,filename = paste0('pic/',pro,'_kegg_up_down.png') )
56
57 ### GSEA
58 if(geneList){
59
60 ## GSEA算法跟上面的使用差异基因集做超几何分布检验不一样。
61 kk_gse <- gseKEGG(geneList = geneList,
62 organism = 'hsa',
63 nPerm = 1000,
64 minGSSize = 20,
65 pvalueCutoff = 0.9,
66 verbose = FALSE)
67 head(kk_gse)[,1:6]
68 gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
69 gseaplot(kk_gse, 'hsa04110',title = 'Cell cycle')
70 kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
71 tmp=kk@result
72 write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
73
74
75 # 这里找不到显著下调的通路，可以选择调整阈值，或者其它。
76 down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
77 up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
78
79 g_kegg=kegg_plot(up_kegg,down_kegg)
80 print(g_kegg)
81 ggsave(g_kegg,filename = paste0(pro,'_kegg_gsea.png'))
82 }
83}
84
85# GO database analysis
86# 具体结果数据在data/下，图在pic/下
87run_go <- function(gene_up,gene_down,pro='test'){
88 gene_up=unique(gene_up)
89 gene_down=unique(gene_down)
90 gene_diff=unique(c(gene_up,gene_down))
91 g_list=list(gene_up=gene_up,
92 gene_down=gene_down,
93 gene_diff=gene_diff)
94 # 因为go数据库非常多基因集，所以运行速度会很慢。
95 if(T){
96 go_enrich_results <- lapply(g_list, function(gene) {
97 lapply( c('BP','MF','CC') , function(ont) {
98 cat(paste('Now process ',names(gene),ont ))
99 ego <- enrichGO(gene = gene,
100 #universe = gene_all,
101 OrgDb = org.Hs.eg.db,
102 ont = ont ,
103 pAdjustMethod = "BH",
104 pvalueCutoff = 0.99,
105 qvalueCutoff = 0.99,
106 readable = TRUE)
107
108 print( head(ego) )
109 return(ego)
110 })
111 })
112 save(go_enrich_results,file =paste0('data/',pro, '_go_enrich_results.Rdata'))
113
114 }
115 # 只有第一次需要运行，就保存结果哈，下次需要探索结果，就载入即可。
116 # load(file=paste0('data/',pro, '_go_enrich_results.Rdata'))
117
118 n1= c('gene_up','gene_down','gene_diff')
119 n2= c('BP','MF','CC')
120 for (i in 1:3){
121 for (j in 1:3){
122 fn=paste0("pic/",pro, '_dotplot_',n1[i],'_',n2[j],'.png')
123 cat(paste0(fn,'\n'))
124 png(fn,res=150,width = 1080)
125 print( dotplot(go_enrich_results[[i]][[j]] ))
126 dev.off()
127 }
128 }
129}
130
131
132
133# 合并up和down的基因KEGG结果，放在一幅图里展示
134kegg_plot <- function(up_kegg,down_kegg){
135 dat=rbind(up_kegg,down_kegg)
136 colnames(dat)
137 dat$pvalue = -log10(dat$pvalue)
138 dat$pvalue=dat$pvalue*dat$group
139
140 dat=dat[order(dat$pvalue,decreasing = F),]
141
142 g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
143 geom_bar(stat="identity") +
144 scale_fill_gradient(low="blue",high="red",guide = FALSE) +
145 scale_x_discrete(name ="Pathway names") +
146 scale_y_continuous(name ="log10P-value") +
147 coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
148 ggtitle("Pathway Enrichment")
149 print(g_kegg)
150}
真正代码
1# 载入数据
2if (T) {
3 rm(list = ls())
4 options(stringsAsFactors = F)
5
6 load(file = 'data/deg.Rdata')
7 head(deg)
8}
9
10# 数据预处理
11## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
12logFC_t=1.5
13# 预处理1
14if (T) {
15 deg$g=ifelse(deg$P.Value>0.05,'stable',
16 ifelse( deg$logFC > logFC_t,'UP',
17 ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
18 )
19 table(deg$g)
20 head(deg)
21 deg$symbol=rownames(deg)
22 library(ggplot2)
23 library(clusterProfiler)
24 library(org.Hs.eg.db)
25 df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
26 toType = c( "ENTREZID"),
27 OrgDb = org.Hs.eg.db)
28 head(df)
29 DEG=deg
30 head(DEG)
31
32 DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')
33 head(DEG)
34 save(DEG,file = 'data/anno_DEG.Rdata')
35}
36# 预处理2
37if (T) {
38 gene_up= DEG[DEG$g == 'UP','ENTREZID']
39 gene_down=DEG[DEG$g == 'DOWN','ENTREZID']
40 gene_diff=c(gene_up,gene_down)
41 gene_all=as.character(DEG[ ,'ENTREZID'] )
42 data(geneList, package="DOSE")
43 head(geneList)
44 boxplot(geneList)
45 boxplot(DEG$logFC)
46
47 geneList=DEG$logFC
48 names(geneList)=DEG$ENTREZID
49 geneList=sort(geneList,decreasing = T)
50}
51
52# detailed plot
53if (T) {
54 source('kegg_and_go_up_and_down.R')
55 run_kegg(gene_up,gene_down,pro='npc_VS_normal')
56 # 需要多go数据库的3个条目进行3次富集分析，非常耗时。
57 run_go(gene_up,gene_down,pro='npc_VS_normal')
58}
59
60
61# 综合显示图
62if (F) {
63 go <- enrichGO(gene_up, OrgDb = "org.Hs.eg.db", ont="all")
64 library(ggplot2)
65 library(stringr)
66 barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")
67 barplot(go, split="ONTOLOGY",font.size =10)+
68 facet_grid(ONTOLOGY~., scale="free") +
69 scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
70 ggsave('pic/gene_up_GO_all_barplot.png')
71
72 go <- enrichGO(gene_down, OrgDb = "org.Hs.eg.db", ont="all")
73 barplot(go, split="ONTOLOGY",font.size =10)+
74 facet_grid(ONTOLOGY~., scale="free") +
75 scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
76 ggsave('pic/gene_down_GO_all_barplot.png')
77}
detailed plot
npc_VS_normal_kegg_up_down.png（kegg结果还可以接受）
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
step5-anno-GSEA
GSEA数据下载网页：https://www.gsea-msigdb.org/gsea/downloads.jsp
msigdb.v7.0.symbols.gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.symbols.gmt
msigdb.v7.0.entrez.gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.entrez.gmt
所有打包gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb_v7.0_files_to_download_locally.zip
1# 载入数据和R包
2# DEG为limma得到的差异分析结果
3if (T) {
4 rm(list = ls())
5 options(stringsAsFactors = F)
6 load(file = 'data/anno_DEG.Rdata')
7 library(ggplot2)
8 library(clusterProfiler)
9 library(org.Hs.eg.db)
10}
11
12
13### 对 MigDB中的全部基因集 做GSEA分析。
14### 按照FC的值对差异基因进行排序
15# http://www.bio-info-trainee.com/2105.html
16# http://www.bio-info-trainee.com/2102.html
17# 自行修改存放gmt文件路径d
18# GSEA每个gene set的具体结果保存在gsea_results这个list中
19# 而最终结果保存在gsea_results_df数据框中
20d='data/GSEA_gmt/msigdb_v7.0_files_to_download_locally/msigdb_v7.0_GMTs/symbols/'
21if(T){
22 geneList=DEG$logFC
23 names(geneList)=DEG$symbol
24 geneList=sort(geneList,decreasing = T)
25 #选择gmt文件（MigDB中的全部基因集）
26
27 gmts=list.files(d,pattern = 'all')
28 gmts
29
30 #GSEA分析
31 library(GSEABase) # BiocManager::install('GSEABase')
32 ## 下面使用lapply循环读取每个gmt文件，并且进行GSEA分析
33 ## 如果存在之前分析后保存的结果文件，就不需要重复进行GSEA分析。
34 f='data/gsea_results.Rdata'
35 if(!file.exists(f)){
36 gsea_results <- lapply(gmts, function(gmtfile){
37 # gmtfile=gmts[2]
38 filepath=paste0(d,gmtfile)
39 geneset <- read.gmt(filepath)
40 print(paste0('Now process the ',gmtfile))
41 egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
42 head(egmt)
43 # gseaplot(egmt, geneSetID = rownames(egmt[1,]))
44 return(egmt)
45 })
46 # 上面的代码耗时，所以保存结果到本地文件
47 save(gsea_results,file = f)
48 }
49 load(file = f)
50 #提取gsea结果，熟悉这个对象
51 gsea_results_list<- lapply(gsea_results, function(x){
52 cat(paste(dim(x@result)),'\n')
53 x@result
54 })
55}
56# 随便看几个结果图
57dev.off()
58gsea_results_df <- do.call(rbind, gsea_results_list)
59gseaplot(gsea_results[[2]],geneSetID = "NIKOLSKY_BREAST_CANCER_7P15_AMPLICON")
60gseaplot(gsea_results[[2]],'RICKMAN_HEAD_AND_NECK_CANCER_D')
step6-anno-GSVA
老实说，我对GSVA还不是很理解，但是知道在代码方面，做起来比较简单，有2个输入文件，一个是表达矩阵，一个是gmt文件即可。先把代码谢在这里，日后如果有需要，再仔细研究吧：
1### 对 MigDB中的全部基因集 做GSVA分析。
2## 还有ssGSEA, PGSEA
3# 载入数据
4if(T){
5 rm(list = ls())
6 options(stringsAsFactors = F)
7
8 library(ggplot2)
9 library(clusterProfiler)
10 library(org.Hs.eg.db)
11 load(file = 'data/step1-output.Rdata')
12 # 每次都要检测数据
13 dat[1:4,1:4]
14}
15
16# GSVA分析
17# 存放gene set的文件路径需要具体修改
18d='data/GSEA_gmt/msigdb_v7.0_files_to_download_locally/msigdb_v7.0_GMTs/symbols/'
19if (T) {
20 X=dat
21 table(group_list)
22 ## Molecular Signatures Database (MSigDb)
23 gmts=list.files(d,pattern = 'all')
24 gmts
25 library(GSVA) # BiocManager::install('GSVA')
26 if(!file.exists('data/gsva_msigdb.Rdata')){
27 es_max <- lapply(gmts, function(gmtfile){
28 # gmtfile=gmts[8];gmtfile
29 geneset <- read.gmt(file.path(d,gmtfile))
30 es.max <- gsva(X, geneset,
31 mx.diff=FALSE, verbose=FALSE,
32 parallel.sz=1)
33 return(es.max)
34 })
35 adjPvalueCutoff <- 0.001
36 logFCcutoff <- log2(2)
37 es_deg <- lapply(es_max, function(es.max){
38 # es.max=es_max[[1]]
39 table(group_list)
40 dim(es.max)
41 design <- model.matrix(~0+factor(group_list))
42 colnames(design)=levels(factor(group_list))
43 rownames(design)=colnames(es.max)
44 design
45 library(limma)
46 # contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
47 contrast.matrix<-makeContrasts("npc-normal",
48 levels = design)
49
50 contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
51
52 deg = function(es.max,design,contrast.matrix){
53 ##step1
54 fit <- lmFit(es.max,design)
55 ##step2
56 fit2 <- contrasts.fit(fit, contrast.matrix)
57 ##这一步很重要，大家可以自行看看效果
58
59 fit2 <- eBayes(fit2) ## default no trend !!!
60 ##eBayes() with trend=TRUE
61 ##step3
62 res <- decideTests(fit2, p.value=adjPvalueCutoff)
63 summary(res)
64 tempOutput = topTable(fit2, coef=1, n=Inf)
65 nrDEG = na.omit(tempOutput)
66 #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
67 head(nrDEG)
68 return(nrDEG)
69 }
70
71 re = deg(es.max,design,contrast.matrix)
72 nrDEG=re
73 head(nrDEG)
74 return(nrDEG)
75 })
76 gmts
77 save(es_max,es_deg,file='data/gsva_msigdb.Rdata')
78 }
79}
80
81# 画图展示，结果存放在pic/下
82if (T) {
83 load(file='data/gsva_msigdb.Rdata')
84 library(pheatmap)
85 lapply(1:length(es_deg), function(i){
86 # i=8
87 print(i)
88 dat=es_max[[i]]
89 df=es_deg[[i]]
90 df=df[df$P.Value<0.01 & abs(df$logFC) > 0.3,]
91 print(dim(df))
92 if(nrow(df)>5){
93 n=rownames(df)
94 dat=dat[match(n,rownames(dat)),]
95 ac=data.frame(g=group_list)
96 rownames(ac)=colnames(dat)
97 rownames(dat)=substring(rownames(dat),1,50)
98 pheatmap::pheatmap(dat,
99 fontsize_row = 8,height = 11,
100 annotation_col = ac,show_colnames = F,
101 filename = paste0('[pic/gsva_',strsplit(gmts[i],'[.]')[[1]][1],'.pdf'))
102
103 }
104 })
105
106 adjPvalueCutoff <- 0.001
107 logFCcutoff <- log2(2)
108 df=do.call(rbind ,es_deg)
109 es_matrix=do.call(rbind ,es_max)
110 df=df[df$P.Value<0.01 & abs(df$logFC) > 0.5,]
111 write.csv(df,file = 'data/GSVA_DEG.csv')
112}
因为用到的这个样本用GSVA没有得到显著性结果，所以没有图出来，具体也没有深究，有需要日后再仔细研究吧
step7-visualization
http://www.360doc.com/conteant/18/0309/18/33459258_735717104.shtml
1### 主要是关于KEGG方面的扩展图
2### 主要是关于KEGG方面的扩展图
3
4# 载入数据
5if (T) {
6 rm(list = ls())
7 options(stringsAsFactors = F)
8
9 library(ggplot2)
10 library(clusterProfiler)
11 library(org.Hs.eg.db)
12 load(file = 'data/anno_DEG.Rdata')
13
14 head(DEG)
15}
16
17# pre-process data
18if (T) {
19 gene_up= DEG[DEG$g == 'UP','ENTREZID']
20 gene_down=DEG[DEG$g == 'DOWN','ENTREZID']
21 gene_diff=c(gene_up,gene_down)
22 gene_all=as.character(DEG[ ,'ENTREZID'] )
23 # 制作差异基因list L
24 boxplot(DEG$logFC)
25
26 geneList=DEG$logFC
27 names(geneList)=DEG$ENTREZID
28 geneList=sort(geneList,decreasing = T)
29}
30
31
32# KEGG富集分析得到结果
33if (T) {
34 if (!file.exists("data/enrichkk.rdata")) {
35 gene_down
36 gene_up
37 enrichKK <- enrichKEGG(gene = gene_up,
38 organism = 'hsa',
39 #universe = gene_all,
40 pvalueCutoff = 0.1,
41 qvalueCutoff =0.1)
42 save(enrichKK,file = "data/enrichkk.rdata")
43 }
44 load(file = "data/enrichkk.rdata")
45 # 查看KEGG结果
46 head(enrichKK)[,1:6]
47 # 打开网页看相关KEGG通路图
48 browseKEGG(enrichKK, 'hsa05150')
49
50 # 将数据中的entrz-id变成symbol
51 # 更为易读
52 enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
53 enrichKK
54}
55
56
57## 可视化
58#条带图
59if (T) {
60 # par(mfrow=c(2,1))
61 barplot(enrichKK,showCategory=20)
62 ggsave("pic/barplot.png")
63}
64
65#气泡图
66if (T) {
67 dotplot(enrichKK)
68 ggsave("pic/dotplot.png")
69}
70
71#下面的图需要映射颜色，设置和示例数据一样的geneList
72
73# 展示top5通路的共同基因，要放大看。
74#Gene-Concept Network
75if (T) {
76 cnetplot(enrichKK, foldChange=geneList,colorEdge = TRUE, circular = F)
77 ggsave("pic/cnetplot.png")
78 cnetplot(enrichKK, foldChange=geneList, colorEdge = TRUE, circular = T)
79 ggsave("pic/cnetplot_circular.png")
80}
81
82
83#Enrichment Map
84if (T) {
85 emapplot(enrichKK)
86 ggsave("pic/Enrichment_Map.png")
87}
88
89#(4)展示通路关系,仅仅是针对于GO数据库结果。
90# goplot(enrichKK)
91#(5)Heatmap-like functional classification
92if (T) {
93 heatplot(enrichKK,foldChange = geneList)
94 ggsave("pic/Enrichment_Heatmap.png")
95}
条带图：
气泡图：
Gene-Concept Network图：每一个小蓝圈表示一个基因，其颜色表示FC值，每个KEGGterm圈的大小由里面包含基因的数目决定。
成环：更炫酷了，但是感觉图形展示不方便
不成环：信息展示更有力吧
Enrichment Map图：这里和上面的图类似，只不过不再显示具体的基因，而是直接画出每个term和term之间的关系，每个圆圈代表着一个term，圆圈大小代表着有多少个基因，颜色表示p值。
如果term和term之间有共同的基因，那么就会连接起来，聚在一起。
Heatmap-like functional classification：
和我们常规的热图不太像，这里纵轴是每个KEGG通路，横轴是涉及到的基因。颜色表示FC值。
上面所有的代码都来自生信技能树曾老板jimmy的帮助，同时我在测试运行的过程中又进行了部分改进和增补。
就用曾老板亲自编辑的感谢词来感谢吧：
Fat Yang thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes !