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包和設置鏡像代碼：

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包，或者說包裝了一個函數吧，如果不想了解具體原理，那么用法如下：

具體用法也非常簡單：

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”）

因為很多時候，用戶是找不到這些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./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}

具體用法也非常簡單：

關于我們得到的ExpressionSet對象，可以通過gset@annotation得到我們需要的注釋平臺信息。

前面說到了這個變量p2s_df，像我這么喜歡資源的人，當然也要保存一份在本地了。大家有興趣的可以自己write到本地保存。哈哈哈哈哈哈哈哈哈哈（jimmy應該不會打我吧）！

附加小功能——對基因名字進行注釋（`annoGene`）

同樣的，還是idmap1這個R包里的函數，如果安裝過這個R包就無須再安裝了，如果沒有安裝又想用這個功能，還真的沒有辦法，因為這些數據存在于這個包的自帶變量中（humanGTF、mouseGTF、ratGTF）:

 1# 安裝方法2(無效)：
 2source('https://raw./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./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./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步：

1. 用linux中的grep -v命令將表頭以“#”和“！”開頭的行去掉

1. 用R讀入數據

1tmp=read.table('data/GPL570.txt',header = T,sep = '\t',fill = T)

結果證明這樣有錯誤，但是，具體原因有空再去找吧。下面看正確的讀入方法——借助GEOqueryR包工具：

一般來說，大家關心的其實就是探針的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)

使用方法：

查看支持的平臺：

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(idmap3)
2ids=idmap3::get_pipe_IDs('GPL21827')

查看支持的平臺：

AnnoProbe包

芯片探針序列的基因注釋已經無需你自己親自做了——一個匯總

安裝方法：

https://codeload.github.com/jmzeng1314/AnnoProbe/legacy.tar.gz/master

1library(devtools)
2install_github('jmzeng1314/AnnoProbe')
3library(AnnoProbe)

使用方法：

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./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的數據，所以在這里也需要先轉置。

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)例如：

• 這里將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}

真正代碼

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/downloads.jsp

msigdb.v7.0.symbols.gmt：https://www./gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.symbols.gmt

msigdb.v7.0.entrez.gmt：https://www./gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.entrez.gmt

所有打包gmt：https://www./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./2105.html
16# http://www./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文件即可。先把代碼謝在這里，日后如果有需要，再仔細研究吧：

因為用到的這個樣本用GSVA沒有得到顯著性結果，所以沒有圖出來，具體也沒有深究，有需要日后再仔細研究吧