中文字幕理论片,69视频免费在线观看,亚洲成人app,国产1级毛片,刘涛最大尺度戏视频,欧美亚洲美女视频,2021韩国美女仙女屋vip视频

內(nèi)容目錄

step0-install_packagesstep1-download_data1.背景知識(shí)2.關(guān)于GEO中的幾個(gè)文件說明A. 2種family file(s)B. Series Matrix File(s)C. GSE64634_RAW.tar3.關(guān)于下載鏡像問題4.關(guān)于探針id轉(zhuǎn)換idmap1包——針對(duì)bioconductor包附加小功能——對(duì)基因名字進(jìn)行注釋（`annoGene`）idmap2包——萬(wàn)能芯片探針I(yè)D注釋平臺(tái)包（提取soft文件）idmap3包——idmap1和idmap2都無法注釋的平臺(tái)AnnoProbe包5. 整體代碼這里抄step2-checkstep3-DEGstep4-anno-go-kegg真正代碼detailed plot綜合顯示圖（更推薦）step5-anno-GSEAstep6-anno-GSVAstep7-visualization致謝

step0-install_packages

簡(jiǎn)單安裝R包和設(shè)置鏡像代碼：

step1-download_data

1.背景知識(shí)

這里通過使用GEOqueryR包來幫助我們下載數(shù)據(jù)，比手動(dòng)登錄GEO數(shù)據(jù)庫(kù)一個(gè)個(gè)點(diǎn)擊要方便很多！

而且我發(fā)現(xiàn)，最近NCBI更新后，GEO數(shù)據(jù)庫(kù)也更新了！

點(diǎn)擊后：

點(diǎn)擊后：

可以看到官網(wǎng)的代碼和我們目前用的代碼基本一致。

2.關(guān)于GEO中的幾個(gè)文件說明

下面一個(gè)個(gè)的來說吧！

A. 2種family file(s)

也就是SOFT formatted family file(s)和MINiML formatted family file(s)，通過字面，我們可以理解這是2種不同格式的探針說明文件。為什么我會(huì)這樣說呢？因?yàn)槲蚁螺d了soft文件來查看：

這時(shí)首行的內(nèi)容：

就是告訴你這個(gè)數(shù)據(jù)是GPL570這個(gè)平臺(tái)測(cè)序得到的：

中間有很多行關(guān)于GSM的，可能記錄的是用到這個(gè)平臺(tái)測(cè)序的sample：

接著就是一系列探針信息了：

我們可以通過在R種提取信息對(duì)我們得到的矩陣種探針做注釋。

B. Series Matrix File(s)

這里面包含了3部分資料：數(shù)據(jù)屬性+患者信息+表達(dá)矩陣

數(shù)據(jù)屬性在最前面的幾行，和患者信息之間有一個(gè)空行，但是他們都是以“!”開頭：

接下來是患者信息：

緊接著患者信息下一行會(huì)有個(gè)提示：!series_matrix_table_begin，然后下面就是表達(dá)矩陣信息了：

這個(gè)就是我們需要的表達(dá)矩陣信息了，矩陣中每一行是一個(gè)基因（也就是一個(gè)探針），每一列是一個(gè)樣本（GSM）

C. GSE64634_RAW.tar

這個(gè)我一開始不知道是什么東西，后來問了問jimmy，然后他給我發(fā)了個(gè)資料你要挖的公共數(shù)據(jù)集作者上傳了錯(cuò)誤的表達(dá)矩陣腫么辦。大致瞅了瞅，知道這個(gè)應(yīng)該是芯片的原始數(shù)據(jù)，然后也把這個(gè)文件給下載下來看了看：

根據(jù)常識(shí)猜了猜想：

X和Y應(yīng)該代表著芯片上的位置，這個(gè)可能和探針有對(duì)應(yīng)關(guān)系。

MEAN是信號(hào)的平均值，STDV是信號(hào)的標(biāo)準(zhǔn)差。

NPIXELS是像素點(diǎn)吧。

既然知道了這是原始數(shù)據(jù)，jimmy又給發(fā)了代碼，感覺不學(xué)習(xí)下有點(diǎn)對(duì)不起自己，那就跟著代碼過一遍吧：

 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文件數(shù)據(jù)
 7celFiles <- list.celfiles('./',listGzipped = T)
 8celFiles
 9affyRaw <- read.celfiles( celFiles )
10# 提取矩陣并做normalization
11eset <- rma(affyRaw)
12eset
13# 檢查數(shù)據(jù)
14dat=exprs(eset)
15dim(dat)
16boxplot(dat,las=2)

可以看到，整個(gè)過程非常的簡(jiǎn)單，就是利用了oligo這個(gè)R包而已。讀入所有的cel文件后，利用rma()函數(shù)將數(shù)據(jù)進(jìn)行了normalization，得到的是一個(gè)ExpressionSet對(duì)象??！

后面會(huì)比較我們利用rawdata得到的結(jié)果和我們直接下載的Series Matrix File(s)文件之間的區(qū)別！

3.關(guān)于下載鏡像問題

jimmy曾經(jīng)發(fā)表過GEO數(shù)據(jù)庫(kù)中國(guó)區(qū)鏡像橫空出世推文，因?yàn)?strong>我每次下載都是用vpn，但是怕以后萬(wàn)一沒有vpn了，所以還是學(xué)習(xí)下，以防萬(wàn)一嘛！再說了，可以學(xué)習(xí)下新知識(shí)，何樂而不為呢？

通過學(xué)習(xí)，發(fā)現(xiàn)實(shí)際上jimmy是開發(fā)了一個(gè)R包，或者說包裝了一個(gè)函數(shù)吧，如果不想了解具體原理，那么用法如下：

具體用法也非常簡(jiǎn)單：

1eSet=geoChina('GSE2345') 

通過安裝方法3（源碼）我們可以看出這個(gè)包的原理：

通過訪問下載 http://49.235.27.111/GEOmirror/GSE2nnn/GSE2546_eSet.Rdata

所以可以猜到，應(yīng)該是jimmy事先用循環(huán)的方式幫我們下載好了很多GEO數(shù)據(jù)，并做成了Rdata格式的文件！非常的良苦用心了??！

4.關(guān)于探針id轉(zhuǎn)換

因?yàn)镚EO數(shù)據(jù)矩陣中，橫坐標(biāo)都是探針的代號(hào)，如下圖：

我們只看這些代號(hào)并不能理解具體是什么基因，于是這就需要我們做id轉(zhuǎn)換了：將探針的代號(hào)轉(zhuǎn)為基因symbol。

idmap1包——針對(duì)bioconductor包

這里又要提到j(luò)immy“開發(fā)”的一個(gè)R包了：

第一個(gè)萬(wàn)能芯片探針I(yè)D注釋平臺(tái)R包——37個(gè)芯片平臺(tái)

老規(guī)矩，還是先學(xué)習(xí)下：

• 一般重要的芯片在R的bioconductor里面都是有包的。但是需要我們單獨(dú)下載對(duì)應(yīng)的包。具體關(guān)系如下：

（具體的R包名稱是bioc_package后加“.db”）

因?yàn)楹芏鄷r(shí)候，用戶是找不到這些GPL平臺(tái)對(duì)應(yīng)的R包，或者下載安裝困難，其實(shí)僅僅是需要探針與基因?qū)?yīng)關(guān)系，沒有必要去下載安裝這幾十個(gè)M的包。于是jimmy就開發(fā)了idmap1這個(gè)R包來幫我們：

安裝方法（這里好像只能用方法1，因?yàn)樵谳d入包的同時(shí)附帶有一個(gè)變量p2s_df加載，如果用方法2和3沒辦法得到這個(gè)變量）：

也可以直接下載：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}

具體用法也非常簡(jiǎn)單：

關(guān)于我們得到的ExpressionSet對(duì)象，可以通過gset@annotation得到我們需要的注釋平臺(tái)信息。

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

附加小功能——對(duì)基因名字進(jìn)行注釋（`annoGene`）

同樣的，還是idmap1這個(gè)R包里的函數(shù)，如果安裝過這個(gè)R包就無須再安裝了，如果沒有安裝又想用這個(gè)功能，還真的沒有辦法，因?yàn)?strong>這些數(shù)據(jù)存在于這個(gè)包的自帶變量中（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}

這里補(bǔ)充學(xué)習(xí)了2個(gè)函數(shù)：

R語(yǔ)言parse函數(shù)與eval函數(shù)的字符串轉(zhuǎn)命令行及執(zhí)行操作：

parse()函數(shù)能將字符串轉(zhuǎn)換為表達(dá)式expression；eval()函數(shù)能對(duì)表達(dá)式求解

idmap2包——萬(wàn)能芯片探針I(yè)D注釋平臺(tái)包（提取soft文件）

第二個(gè)萬(wàn)能芯片探針I(yè)D注釋平臺(tái)R包——有122個(gè)GPL

這次是根據(jù)[soft文件](#####A. 2種family file(s))進(jìn)行提取信息得到的！

如果是我自己來處理這樣的文件，我應(yīng)該會(huì)分2步：

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

1. 用R讀入數(shù)據(jù)

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

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

一般來說，大家關(guān)心的其實(shí)就是探針的ID，以及基因的symbol列。有了這個(gè)變量后，就可以按照R語(yǔ)言基本操作來提取我們需要的信息了。

注意：我檢查了得到的結(jié)果，里面存在有的探針I(yè)D對(duì)應(yīng)2個(gè)基因，如下圖：

雖然不知道這些代表著什么意思，但是，我將這個(gè)數(shù)據(jù)和bioconductor包里的hgu133plus2.db數(shù)據(jù)做了比較，結(jié)果是這樣的：自己提取的結(jié)果中如果是一個(gè)ID對(duì)應(yīng)2個(gè)基因，那么這個(gè)探針在bioconductor包里基本上找不到數(shù)據(jù)。而其他一個(gè)ID對(duì)應(yīng)2個(gè)基因的結(jié)果均和bioconductor包一致。

當(dāng)然了，我這只是單獨(dú)來看一個(gè)平臺(tái)的探針，而在idmap2jimmy已經(jīng)幫我們整理好了，直接用就行了??！

安裝方法：

（比較慢）

也可以直接下載：https://codeload.github.com/jmzeng1314/idmap2/legacy.tar.gz/master

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

使用方法：

查看支持的平臺(tái)：

1gpl_list[,1:4] # gpl_list是R包自帶的變量

idmap3包——idmap1和idmap2都無法注釋的平臺(tái)

第三個(gè)萬(wàn)能芯片探針I(yè)D注釋平臺(tái)R包——idmap1和idmap2都無法注釋的平臺(tái)

如果我們拿到的soft注釋文件中是序列信息，那么我們?cè)撛趺醋瞿兀?/p>

應(yīng)該是先將序列比對(duì)到參考基因組上，然后通過提取基因注釋文件中的數(shù)據(jù)得到基因symbol！

而在idmap3包中，jimmy已經(jīng)幫我們做好了??！說他寵粉也是真的了??！我都懶得做，他居然還寫了個(gè)循環(huán)來完成了這個(gè)事。

安裝方法：

（比較慢）

也可以直接下載：https://codeload.github.com/jmzeng1314/idmap3/legacy.tar.gz/master

使用方法：

1library(idmap3)
2ids=idmap3::get_pipe_IDs('GPL21827')

查看支持的平臺(tái)：

AnnoProbe包

芯片探針序列的基因注釋已經(jīng)無需你自己親自做了——一個(gè)匯總

安裝方法：

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# 下載的文件會(huì)自動(dòng)保存到./data/下
 5# 得到的gset是一個(gè)ExpressionSet對(duì)象
 6
 7
 8# 只需要修改file的名字，即可下載得到相應(yīng)的geo數(shù)據(jù)
 9# 獲取患者信息，需要自行修改
10file='GSE64634'
11# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64634
12
13if (T) {
14  # 數(shù)據(jù)下載
15  if (T) {
16    library(GEOquery)
17    # 這個(gè)包需要注意兩個(gè)配置，一般來說自動(dòng)化的配置是足夠的。
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)       ## 平臺(tái)文件
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)# 根據(jù)disease state:ch1一列得知分組信息
36    group_list=c(rep('normal',4),rep('npc',12))# nasopharyngeal carcinoma NPC 鼻咽癌
37    table(group_list)
38  }
39
40  # 對(duì)數(shù)據(jù)進(jìn)行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這個(gè)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,]# 調(diào)整順序，讓dat的順序和ids中的一致
83
84    ids$median=apply(dat,1,median)
85    ids=ids[order(ids$symbol,ids$median,decreasing = T),]# 按照基因名、中位數(shù)大小排序
86    ids=ids[!duplicated(ids$symbol),]# 只保留相同symbol中中位數(shù)最大的探針
87    dat=dat[ids$probe_id,]# 調(diào)整順序，讓dat的順序和ids中的一致
88    rownames(dat)=ids$symbol# id轉(zhuǎn)換
89    dat[1:4,1:4]
90  }
91}
92
93save(dat,group_list,file = 'data/step1-output.Rdata')

step2-check

• 畫PCA圖時(shí)要求是行名時(shí)樣本名，列名時(shí)探針名，因此此時(shí)需要轉(zhuǎn)置。格式要求data.frame。

• 關(guān)于scale：在基因表達(dá)量的數(shù)據(jù)中，有些基因表達(dá)量極低，有些基因表達(dá)量極高，因此把每個(gè)基因在不同處理和重復(fù)中的數(shù)據(jù)轉(zhuǎn)換為平均值為0，方差為1的數(shù)據(jù)，所以在這里也需要先轉(zhuǎn)置。

all_samples_PCA.png

heatmap_top1000_sd.png

cor_top500_mad.png

step3-DEG

關(guān)于差異分析是否需要比較矩陣

差異分析的統(tǒng)計(jì)學(xué)本質(zhì)探究

• 原始版火山圖：plot(logFC,-log10(P.Value))

• 原始版MA圖：plot(AveExpr,logFC)

  1# 載入數(shù)據(jù),檢查數(shù)據(jù)
  2if (T) {
  3  rm(list = ls()) 
  4  options(stringsAsFactors = F)
  5  library(ggpubr)
  6  load(file = 'data/step1-output.Rdata')
  7  # 每次都要檢測(cè)數(shù)據(jù)
  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# 利用定義好的函數(shù)檢查數(shù)據(jù)
 23bplot(dat[1,])
 24bplot(dat[2,])
 25bplot(dat[3,])
 26bplot(dat[4,])
 27
 28
 29# limma
 30library(limma)
 31# 方法1：不制作比較矩陣，簡(jiǎn)單
 32# 但是做不到隨心所欲的指定任意兩組進(jìn)行比較
 33if (T) {
 34  design=model.matrix(~factor( group_list ))
 35  fit=lmFit(dat,design)
 36  fit2=eBayes(fit)
 37  ## 上面是limma包用法的一種方式 
 38  options(digits = 4) #設(shè)置全局的數(shù)字有效位數(shù)為4
 39  topTable(fit2,coef=2,adjust='BH') 
 40}
 41
 42# 方法2：制作比較矩陣
 43# 可以隨心所欲的指定任意兩組進(jìn)行比較
 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  # 這個(gè)矩陣聲明，我們要把 npc 組跟 Normal 進(jìn)行差異分析比較
 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是把所有結(jié)果都打印出來
 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# 結(jié)果存放在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  # 每次都要檢測(cè)數(shù)據(jù)
129  dat[1:4,1:4]
130  table(group_list)
131  x=deg$logFC
132  names(x)=rownames(deg)
133
134  # cg中存放著變化上升和下降的前100個(gè)基因名
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

• 關(guān)于基因名和ID之間的各種轉(zhuǎn)換：bitr(geneID, fromType, toType, OrgDb, drop = TRUE)例如：

• 這里將KEGG和GO富集分析函數(shù)自定義為3個(gè)函數(shù)——run_kegg、run_go、kegg_plot

  每次運(yùn)行前先運(yùn)行下面代碼：

  1# KEGG pathway analysis
  2# 具體結(jié)果數(shù)據(jù)在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個(gè)基因集分開做超幾何分布檢驗(yàn)
  9  # 首先是上調(diào)基因集。
 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  # 首先是下調(diào)基因集。
 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  # 最后是上下調(diào)合并后的基因集。
 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  #畫圖設(shè)置, 這個(gè)圖很丑，大家可以自行修改。
 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算法跟上面的使用差異基因集做超幾何分布檢驗(yàn)不一樣。
 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    # 這里找不到顯著下調(diào)的通路，可以選擇調(diào)整閾值，或者其它。
 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# 具體結(jié)果數(shù)據(jù)在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  # 因?yàn)間o數(shù)據(jù)庫(kù)非常多基因集，所以運(yùn)行速度會(huì)很慢。
 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  # 只有第一次需要運(yùn)行，就保存結(jié)果哈，下次需要探索結(jié)果，就載入即可。
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結(jié)果，放在一幅圖里展示
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結(jié)果還可以接受）

GO系列結(jié)果過于冗余：

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數(shù)據(jù)下載網(wǎng)頁(yè)：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# 載入數(shù)據(jù)和R包
 2# DEG為limma得到的差異分析結(jié)果
 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### 對(duì) MigDB中的全部基因集 做GSEA分析。
14### 按照FC的值對(duì)差異基因進(jìn)行排序
15# http://www.bio-info-trainee.com/2105.html
16# http://www.bio-info-trainee.com/2102.html 
17# 自行修改存放gmt文件路徑d
18# GSEA每個(gè)gene set的具體結(jié)果保存在gsea_results這個(gè)list中
19# 而最終結(jié)果保存在gsea_results_df數(shù)據(jù)框中
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循環(huán)讀取每個(gè)gmt文件，并且進(jìn)行GSEA分析
33  ## 如果存在之前分析后保存的結(jié)果文件，就不需要重復(fù)進(jìn)行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    # 上面的代碼耗時(shí)，所以保存結(jié)果到本地文件
47    save(gsea_results,file = f)
48  }
49  load(file = f)
50  #提取gsea結(jié)果，熟悉這個(gè)對(duì)象
51  gsea_results_list<- lapply(gsea_results, function(x){
52    cat(paste(dim(x@result)),'\n')
53    x@result
54  })
55}
56# 隨便看幾個(gè)結(jié)果圖
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

老實(shí)說，我對(duì)GSVA還不是很理解，但是知道在代碼方面，做起來比較簡(jiǎn)單，有2個(gè)輸入文件，一個(gè)是表達(dá)矩陣，一個(gè)是gmt文件即可。先把代碼謝在這里，日后如果有需要，再仔細(xì)研究吧：

因?yàn)橛玫降倪@個(gè)樣本用GSVA沒有得到顯著性結(jié)果，所以沒有圖出來，具體也沒有深究，有需要日后再仔細(xì)研究吧