GEO再學(xué)習(xí)

GEO數(shù)據(jù)挖掘已經(jīng)成為生信學(xué)者必備技能,我以為自己會(huì)跑代碼了就是會(huì)了,其實(shí)呢,細(xì)細(xì)去領(lǐng)會(huì)每個(gè)細(xì)節(jié),還是會(huì)發(fā)現(xiàn)很多盲區(qū).聽了曉娟老師的課,發(fā)現(xiàn)很多很多細(xì)節(jié)我并不甚清,一定要明白,才能融會(huì)貫通.同時(shí),曉娟老師又把思路捋地非常清,練習(xí)實(shí)操幾次,一定可以像別人家的孩子那樣,做到數(shù)據(jù)挖掘,復(fù)現(xiàn)文章美圖.同時(shí),如果自己有了數(shù)據(jù),也可以這樣分析了!

GEO數(shù)據(jù)挖掘思路

• 1.下載數(shù)據(jù)
• 2.數(shù)據(jù)預(yù)處理
• 3.差異分析
• 4.差異分析后的可視化-熱圖紫岩、火山圖、富集分析

1.下載數(shù)據(jù)

(1)用GSE號(hào)下載

從文章中搜索GSE號(hào),然后去GEO網(wǎng)址:https://www.ncbi.nlm.nih.gov/geo/下載,以``GSE42872`為例

Platform (GPL) :芯片或測(cè)序平臺(tái)的相關(guān)信息
Sample (GSM) :單個(gè)樣本的數(shù)據(jù)
Series (GSE) :一系列GSM整合后的數(shù)據(jù)
DataSets (GDS):由GEO的工作人員整合的樣品數(shù)據(jù)

如果上面這些發(fā)懵,記不住,那看下圖

用GEOquery包的下載方法

平臺(tái)與對(duì)應(yīng)的基因注釋包

代碼如下

library(GEOquery)
f<-'GSE42872.Rdata'
####getGPL獲得平臺(tái)的注釋信息，但下載速度會(huì)慢很多
####而且注釋文件格式大多不如bioconductor包好用
if(!file.exists(f)){
  gset<-getGEO('GSE42872',destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=f)
}
#數(shù)據(jù)提取
load('GSE42872.Rdata')

(2)下載原始數(shù)據(jù)-txt.gz

某些情況下,下載的series Matrix File(s)文件并不是整理好的,這時(shí)可能需要下載原始文件,就是RAW.tar.參考學(xué)徒數(shù)據(jù)挖掘干擾MYC-WWP1通路重新激活PTEN的抑癌活性——3步搞定GSEA分析一樣可以整理成想要的表達(dá)矩陣信息.文章中的GSE號(hào)為gse126789

下載RAW data

下載后的文件是txt.gz結(jié)尾的

image-20190804040439655

下載數(shù)據(jù)后,處理的代碼如下

rm(list = ls())  
options(stringsAsFactors = F)
# 下載GSE126789_RAW.tar
# 批量讀取文件處理為表達(dá)矩陣
dir <- "/Users/mengmeng/Downloads/GSE126789_RAW/"#上一步下載的文件路徑

files <- list.files(path = dir, pattern = "*wwp1.*.txt.gz", recursive = T) #找到wwp1的WT和敲除樣本,也就是文件名,用的list.files函數(shù)僅僅是文件的名字,并不包括里面的數(shù)據(jù)
expr <- lapply(files,
               function(x){
                 # 只讀取基因名和count
                 expr <- read.table(file = file.path(dir,x), sep = "\t", header = T, stringsAsFactors = F)[,c(1,7)]#取出來(lái)第一列和第七列,即基因名和count計(jì)數(shù),矩陣中取不相鄰的兩列就可以[,c(1,7)]這么取,也可以cbind(b$Geneid,b$Count)取,但啰嗦
                 return(expr)#注意此時(shí)的files僅僅是文件名字,并不是解壓后文件含有的矩陣信息,而這時(shí)的files是character,是字符型,要清楚的是,lapply不僅可以作用于列表list,同樣可以應(yīng)用在向量上.向量包括數(shù)值、字符赞弥、邏輯型向量
               }) #這個(gè)循環(huán)非常棒,要理解function(x)這個(gè)里面的x是什么.files是什么,x就是什么,files是文件名,那么x也是文件名,所以才有file = file.path(dir,x),這個(gè)file.path(dir,x)就相當(dāng)于read.table后面括號(hào)中的文件名
df <- do.call(cbind, expr)#得到的每一個(gè)矩陣橫向疊加(cbind就是列相加)
rownames(df) <- df[,1]
df <- df[,seq(2,ncol(df),by=2)] #去掉重復(fù)讀取的基因名#seq函數(shù):從第二列開始,直到最后一列,by間隔為2來(lái)取
colnames(df) <- substr(files,1,10) #列名為樣本名,用substr提取字符串,僅僅保留樣本名如GSM3613325
df <- df[apply(df, 1, sum)!=0,] #去掉在所有樣本中count都為0的基因,取出來(lái)不等于0的矩陣
grouplist <- c(rep("KO",3),rep("WT",3)) #記錄下分組信息,新建一個(gè)grouplist向量,即'KO'重復(fù)3次,'WT'重復(fù)3次
expr <- df
save(expr,grouplist,file = "./expr_group.Rdata")

(3)下載原始數(shù)據(jù)-CEL.gz

同上面一樣,如果下載的RAW.tar數(shù)據(jù)是CEL.gz,那就去看jimmy老師的文章吧,用affy包讀取affymetix的基因表達(dá)芯片數(shù)據(jù)-CEL格式數(shù)據(jù)

2.數(shù)據(jù)預(yù)處理

本次對(duì)表達(dá)矩陣進(jìn)行差異分析時(shí),用的是limma包.通常使用limma進(jìn)行差異分析處理時(shí)莹痢，需要經(jīng)過(guò)log2后的矩陣作為表達(dá)矩陣輸入堕战。根據(jù)log2FC的定義默刚，這個(gè)數(shù)字表示變化倍數(shù)經(jīng)過(guò)log2后的一個(gè)值，比如log2FC=1,則變化為2倍；log2FC=2,則變化為4倍筋蓖。這是常用的一種表述方法卸耘。在使用limma函數(shù)計(jì)算時(shí)，如果輸入的矩陣沒(méi)有經(jīng)過(guò)log2處理粘咖，則會(huì)把FC當(dāng)成log2FC輸入蚣抗，這或許是因?yàn)閘imma默認(rèn)輸入的是log2后的表達(dá)式。這里有必要提到log的一個(gè)運(yùn)算瓮下，即翰铡，可見(jiàn)對(duì)于已經(jīng)log2后的數(shù)據(jù)，計(jì)算log2FC = log2(A/B)只需要直接使用log2A-log2B讽坏。所以如果給出的是一個(gè)未經(jīng)log2的數(shù)值锭魔，函數(shù)也會(huì)直接相減以得到log2FC，這就導(dǎo)致計(jì)算出來(lái)的差異表達(dá)高達(dá)幾百甚至上千路呜。在判斷counts是否需要重新標(biāo)準(zhǔn)化以及是否需要log2時(shí)赂毯，可以根據(jù)數(shù)值大小粗略估計(jì)。如果表達(dá)豐度的數(shù)值在50以內(nèi)拣宰，通常是經(jīng)過(guò)log2轉(zhuǎn)化的。如果數(shù)字在幾百幾千烦感，則是未經(jīng)轉(zhuǎn)化的巡社。因?yàn)?的幾十次方已經(jīng)非常巨大，如果2的幾百次方手趣，則不符合實(shí)際情況晌该。參考

1.https://blog.csdn.net/tuanzide5233/article/details/88542805,

2.https://blog.csdn.net/tuanzide5233/article/details/88542558

3.http://www.bio-info-trainee.com/1209.html

用代碼查看矩陣是否需要log

qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
ex <- log2(ex+1)

• 用箱線圖看數(shù)據(jù)的分布

boxplot(ex,las=2,cex.axis=0.6,main='data check')

image-20190804071930092

• 用PCA看分組,即檢驗(yàn)現(xiàn)在的表達(dá)矩陣中的樣本信息所對(duì)應(yīng)的分組信息是否有以下的情況:
  • 是否有離群樣本
  • 實(shí)驗(yàn)組和對(duì)照組是否正確(有沒(méi)有標(biāo)反)
  • 有沒(méi)有批次效應(yīng)

主成分分析（PCA)是一種降維分析，將多個(gè)維度的數(shù)據(jù)簡(jiǎn)化降維绿渣，突出主要成分朝群，所以叫做主成分分析。PCA是一種數(shù)據(jù)降維技巧中符，它能將大量相關(guān)變量轉(zhuǎn)化為一組很少的不相關(guān)變量姜胖，這些無(wú)關(guān)變量稱為主成分。例如淀散，使用PCA可將30個(gè)相關(guān)（很可能冗余）的環(huán)境變量轉(zhuǎn)化為5個(gè)無(wú)關(guān)的成分變量右莱，并且盡可能地保留原始數(shù)據(jù)集的信息。

用一個(gè)通俗的例子來(lái)舉例:

image-20190804071242809

這個(gè)是大英帝國(guó)四個(gè)地區(qū)的17種食物的消費(fèi)數(shù)據(jù)（單位克/人/周）档插，你能看出來(lái)四個(gè)地方的消費(fèi)有何不同嗎慢蜓？或者說(shuō)，四個(gè)地方哪幾個(gè)更為接近郭膛？誰(shuí)更為特殊晨抡？

如果光肉眼看，肯定是很費(fèi)勁的，而且也不一定準(zhǔn)耘柱？

從數(shù)據(jù)的結(jié)構(gòu)上來(lái)看如捅，以英國(guó)4個(gè)地區(qū)為自變量的話，這么這個(gè)數(shù)據(jù)實(shí)際上含有17個(gè)維度帆谍，每一個(gè)維度都不是完全一樣伪朽，但是每一個(gè)維度的近似程度或者疏遠(yuǎn)程度是有差別的，也就是說(shuō)汛蝙，每一個(gè)維度對(duì)整體變異度的貢獻(xiàn)不一樣烈涮。

通俗點(diǎn)說(shuō)，每種食物的消耗量的差別對(duì)“四個(gè)地區(qū)的飲食區(qū)別”這一結(jié)果的影響程度是不一樣的窖剑。

那么坚洽，我們就需要把哪些對(duì)“四個(gè)地區(qū)的飲食區(qū)別”影響最大的食物給篩選出來(lái)，然后排序西土，并依據(jù)相關(guān)公式進(jìn)行整合計(jì)算讶舰，計(jì)算出一個(gè)新的參數(shù)——也就是主成分。排名第一的主成分對(duì)整體變異度的貢獻(xiàn)最大需了。

對(duì)上述數(shù)據(jù)進(jìn)行PCA分析跳昼，結(jié)果如下：

image-20190804071350641

這是第一主成分的區(qū)分情況，我們?cè)偌由系诙€(gè)維度：PC2

image-20190804071408923

PC1對(duì)四個(gè)地區(qū)的區(qū)分結(jié)果最好肋乍，PC2就看不出有什么區(qū)分結(jié)果了鹅颊。當(dāng)然，數(shù)據(jù)不同，PC1和PC2的貢獻(xiàn)度也不同，有時(shí)候PC2對(duì)數(shù)據(jù)的區(qū)分力度也很大驼仪，要看具體的數(shù)據(jù).

比如有8種樣本，6個(gè)是近似病理類型的帝雇，另外兩個(gè)是另一大類的疾病類型的樣本。這時(shí)候?qū)?個(gè)樣本的蛋白質(zhì)組學(xué)數(shù)據(jù)進(jìn)行PCA分析蛉拙，那么尸闸，鑒定到的數(shù)千個(gè)蛋白的豐度信息經(jīng)過(guò)PCA的計(jì)算，其結(jié)果一定是某6個(gè)聚集在一起刘离，另外2個(gè)聚集在一起室叉。

假設(shè)樣本中平均鑒定到3000個(gè)蛋白，那么也就意味著本數(shù)據(jù)有3000個(gè)維度硫惕，每個(gè)蛋白（維度）的相對(duì)定量信息的變化在這8個(gè)樣本中的分布都是不一樣的茧痕，有些蛋白在8個(gè)樣本中變異很小，有些很大恼除，有些居中踪旷。那么曼氛，應(yīng)用PCA進(jìn)行降維，將3000個(gè)維度降至2-3個(gè)維度令野，從而簡(jiǎn)化了數(shù)據(jù)舀患，突出了主要矛盾。參考:

http://bbs.foodmate.net/thread-1071830-1-1.html

那么PC1和PC2是什么呢?上面舉例中,比如3000個(gè)蛋白,3000個(gè)維度,PC1是多少個(gè)維度的組合呢?是哪些個(gè)維度的組合呢?下面這張R語(yǔ)言實(shí)戰(zhàn)中的圖片可以回答,并不是某一個(gè)維度,而是

image-20190804072619643

PCA的目標(biāo)是用一組較少的不相關(guān)變量代替大量相關(guān)變量气破，同時(shí)盡可能保留初始變量的信息聊浅，這些推導(dǎo)所得的變量稱為主成分，它們是觀測(cè)變量的線性組合现使。如第一主成分為

image-20190804073137978

它是k個(gè)觀測(cè)變量的加權(quán)組合低匙，對(duì)初始變量集的方差解釋性最大。第二主成分也是初始變量的線性組合碳锈，對(duì)方差的解釋性排第二顽冶，同時(shí)與第一主成分正交（不相關(guān)）。后面每一個(gè)主成分都最大化它對(duì)方差的解釋程度售碳，同時(shí)與之前所有的主成分都正交.目標(biāo)是能用較少的主成分來(lái)解釋全部變量强重。

####PCA看分組情況
library("FactoMineR")
library("factoextra") 
df.pca <- PCA(t(ex), graph = FALSE)
fviz_pca_ind(df.pca,
             geom.ind = "point",
             col.ind = group_list, 
             addEllipses = TRUE, 
             legend.title = "Groups"
)
save(ex,pd,group_list,file='ex_pd.Rdata')

當(dāng)處理組和未處理組分得很開時(shí),說(shuō)明實(shí)驗(yàn)分組

處理組和未處理組分得很開

3.差異分析

下面是GSE5282的數(shù)據(jù)集,整理好所得到的表達(dá)矩陣和分組矩陣.

image-20190804073913072

image-20190804074017364

表達(dá)矩陣中的樣本名的順序要和分組信息中的順序一一對(duì)應(yīng)上,非常重要

rm(list = ls()) 
options(stringsAsFactors = F)
load('ex_pd.Rdata')
ex <- ex[,1:6]
#####可以自己生成，group_list不可以有空格贸人；
group_list <- group_list[1:6]
group_list<- gsub(' ','_',trimws(group_list))#
###構(gòu)建實(shí)驗(yàn)設(shè)計(jì)矩陣
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(ex)
library(limma)#limma
####構(gòu)建比較矩陣,設(shè)置用來(lái)對(duì)比的組別
contrast.matrix<-makeContrasts(contrasts=c('EGF_4h-Control_4h'),levels = design)
######limma三部曲,只需要?dú)w一化后的數(shù)據(jù)间景、實(shí)驗(yàn)矩陣、比較矩陣的輸入
fit <- lmFit(ex,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等統(tǒng)計(jì)的結(jié)果艺智，topTable對(duì)p值校驗(yàn)拱燃，對(duì)基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
dim(tT)
####提取出我們需要的指標(biāo)
tT <- subset(tT, select=c('adj.P.Val',"P.Value","logFC"))
save(tT,file='DEG.Rdata')

>**多重檢驗(yàn)時(shí)，產(chǎn)生adjust.pvalue的需求**

接下來(lái),進(jìn)行id轉(zhuǎn)換

image-20190804214503093

同上面的一樣,在平臺(tái)-注釋包對(duì)應(yīng)內(nèi)找到GPL1355對(duì)應(yīng)的R包注釋集,是rat2302,下載的時(shí)候加上'db'.

BiocManager::install("rat2302.db",ask = F,update = F)
library('rat2302.db')
load('DEG.Rdata')
#####
ls("package:rat2302.db")
p2se<- select(rat2302.db,keys = keys(rat2302.db),columns = c('SYMBOL','ENTREZID'))#select這個(gè)函數(shù)可以直接將包中的'SYMBOL','ENTREZID'提取出來(lái),相比用to table 方便很多
tT$probe_id <- rownames(tT)
tT <- merge(tT,p2se,by.x='probe_id',by.y='PROBEID')#大小寫不同,可以用by.x,by.y,對(duì)應(yīng)上
save(tT,file='DEG_symbol.Rdata')

4.差異分析后的可視化-熱圖力惯、火山圖、富集分析

熱圖

load('ex_pd.Rdata')
load('DEG.Rdata')
ex <- ex[,1:6]
save(ex,pd,group_list,file='ex_pd.Rdata')
load('ex_pd.Rdata')
######pheatmap 這是取前100個(gè)gene進(jìn)行展示召嘶，可根據(jù)自己的需求進(jìn)行調(diào)整
choose_gene<-rownames(tT)[1:100]
#####用normalize后的數(shù)據(jù)進(jìn)行展示
data_matrix <-ex[choose_gene,]
#####熱圖更詳細(xì)的了解http://www.reibang.com/p/0e5ce59fa79e
#####scale對(duì)數(shù)據(jù)處理父晶，求得的結(jié)果為z-score，即數(shù)據(jù)與均值之間差幾個(gè)標(biāo)準(zhǔn)差 
####t是對(duì)數(shù)據(jù)做轉(zhuǎn)置處理弄跌，為了適應(yīng)scale的要求
data_matrix=t(scale(t(data_matrix)))
#####查看scale處理后數(shù)據(jù)的范圍
fivenum(data_matrix)
####目的是避免出現(xiàn)極大極小值影響可視化的效果
data_matrix[data_matrix>1.2]=1.2
data_matrix[data_matrix< -1.2] = -1.2
library(pheatmap)
pheatmap(data_matrix)
####調(diào)整下顏色甲喝，使正負(fù)值顏色的對(duì)比會(huì)更加鮮明
require(RColorBrewer)
bk = c(seq(-1.2,1.2, length=100))
annotation_col = data.frame(Group = rep(c("EGF_4h", "Control_4h"), c(3,3)))#annotation_col
rownames(annotation_col)<-colnames(ex)
####annotation_col和annotation_row的格式需為數(shù)據(jù)框
####breaks參數(shù)用于值匹配顏色
####看下，color和breaks的對(duì)應(yīng)铛只，進(jìn)行理解
pheatmap(data_matrix,
         breaks=bk,#從bk開始
         show_rownames = F,#默認(rèn)是展示
         annotation_col = annotation_col,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         filename = 'test.png')
####可以調(diào)整的內(nèi)容有是否聚類埠胖、是否分組、是否顯示行和列的內(nèi)容等

image-20190804220955224

test.png

火山圖

load(file = "DEG_symbol.Rdata")
####主要用兩個(gè)值淳玩，logFC和p.value
####可參考進(jìn)行l(wèi)ogFC_cutoff的計(jì)算
logFC_cutoff <- with(tT,mean(abs( logFC)) + 2*sd(abs( logFC)) )
###logFC_cutoff也可自行根據(jù)需求設(shè)置
####依據(jù)logFC和adj.P.val為火山圖的顏色分組做準(zhǔn)備
tT$change = ifelse(tT$adj.P.Val < 0.05 & abs(tT$logFC) > logFC_cutoff,
                         ifelse(tT$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
####記錄上調(diào)直撤、下調(diào)基因數(shù)目，作為title的內(nèi)容
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(tT[tT$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(tT[tT$change =='DOWN',]))

load(file = "DEG_symbol.Rdata")
####主要用兩個(gè)值蜕着，logFC和p.value
####可參考進(jìn)行l(wèi)ogFC_cutoff的計(jì)算
logFC_cutoff <- with(tT,mean(abs( logFC)) + 2*sd(abs( logFC)) )
###logFC_cutoff也可自行根據(jù)需求設(shè)置
####依據(jù)logFC和adj.P.val為火山圖的顏色分組做準(zhǔn)備
tT$change = ifelse(tT$adj.P.Val < 0.05 & abs(tT$logFC) > logFC_cutoff,
                         ifelse(tT$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
####記錄上調(diào)谋竖、下調(diào)基因數(shù)目红柱，作為title的內(nèi)容
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(tT[tT$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(tT[tT$change =='DOWN',]))

library(ggplot2)
g = ggplot(data=tT, 
           aes(x=logFC, y=-log10(P.Value), color=change)) +
    geom_point(alpha=0.4, size=1.75) + #出點(diǎn)狀圖
    theme_set(theme_set(theme_bw(base_size=20)))+ 
    xlab("log2 fold change") + ylab("-log10 p-value") + 
    ggtitle( this_tile ) + 
    theme(plot.title = element_text(size=15,hjust = 0.5))+
    scale_colour_manual(values = c('blue','black','red'))+
   annotate('text',x=tT$logFC[tT$logFC>5],y=-log10(tT$P.Value[tT$logFC>5]),label=tT$SYMBOL[tT$logFC>5])
print(g)
ggsave(g,filename = 'volcano.png')
save(tT,file='DEG_change.Rdata')

火山圖

load('DEG_change.Rdata')
tT <- na.omit(tT)
####富集前的基因處理
library(clusterProfiler)
gene_up= tT[tT$change == 'UP','ENTREZID'] 
gene_down=tT[tT$change == 'DOWN','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
gene_all = tT$ENTREZID
###############################KEGG富集分析#####################################################
## 1.KEGG pathway analysis
#上調(diào)、下調(diào)蓖乘、差異锤悄、所有基因
#clusterProfiler作kegg富集分析：
save(kk.down,file='down.Rdata')
kk.up<- enrichKEGG(gene         = gene_up,
                   organism     = 'rno',
                   universe     = gene_all,
                   pvalueCutoff = 1)
kk.down <- enrichKEGG(gene         =  gene_down,
                      organism     = 'rno',
                      universe     = gene_all,
                      pvalueCutoff = 1)
#######################下述ggplot代碼用于畫圖########################################
down_kegg <- kk.down@result
down_kegg <- down_kegg[down_kegg$pvalue<0.05,]
library(ggplot2)
ggplot(down_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")

image-20190804235105060

image-20190804234953521

KEGG注釋

#首先是KEGG的注釋
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load('DEG_change.Rdata')
tT <- na.omit(tT)
####富集前的基因處理
library(clusterProfiler)
gene_up= tT[tT$change == 'UP','ENTREZID'] 
gene_down=tT[tT$change == 'DOWN','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
gene_all = tT$ENTREZID
#######KEGG富集分析
## 1.KEGG pathway analysis
#上調(diào)嘉抒、下調(diào)零聚、差異、所有基因
#clusterProfiler作kegg富集分析：
save(kk.down,file='down.Rdata')
kk.up<- enrichKEGG(gene         = gene_up,
                   organism     = 'rno',
                   universe     = gene_all,
                   pvalueCutoff = 1)
kk.down <- enrichKEGG(gene         =  gene_down,
                      organism     = 'rno',
                      universe     = gene_all,
                      pvalueCutoff = 1)
#######################下述ggplot代碼用于畫圖########################################
down_kegg <- kk.down@result
down_kegg <- down_kegg[down_kegg$pvalue<0.05,]
#down_kegg
up_kegg<- kk.up@result
up_kegg <- up_kegg[up_kegg$pvalue<0.05,]
library(ggplot2)
ggplot(down_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")
#up_kegg
library(ggplot2)
ggplot(up_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")

image-20190805001223395

down_kegg

image-20190805001139289

up_kegg

富集程度看顏色深淺,而柱狀圖長(zhǎng)短代表Count數(shù),最后出圖就是上面富集分析(enrichKEGG函數(shù))可視化后的展示)

down_kegg-dotplot

GO注釋

主要包括這三個(gè)方面的注釋

image-20190805002532500

#代碼如下,可畫圖展示
library(org.Rn.eg.db)
ego_CC <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "CC",
                   pAdjustMethod = "BH")
#生物過(guò)程
ego_BP <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "BP",
                   pAdjustMethod = "BH")
#分子功能：
ego_MF <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "MF",
                   pAdjustMethod = "BH")
#可視化展示-Barplot
barplot(ego_CC, showCategory=12,title="bar_CC")
barplot(ego_BP, showCategory=12,title="bar_BP")
barplot(ego_MF, showCategory=12,title="bar_MF")
#可視化展示-Dotplot
dotplot(ego_CC,showCategory = 12,title="dot_CC")
dotplot(ego_BP,showCategory = 12,title="dot_BP")
dotplot(ego_MF,showCategory = 12,title="dot_MF")
save(kk.up,kk.down,ego_BP,ego_CC,ego_MF,file='annotation.Rdata')

最后友情宣傳生信技能樹

• 生物信息學(xué)“義診”

• 生物信息學(xué)"拍賣會(huì)"

• 全國(guó)巡講:R基礎(chǔ)些侍，Linux基礎(chǔ)和RNA-seq實(shí)戰(zhàn)演練 : 預(yù)告：12月28-30長(zhǎng)沙站

• 廣州珠江新城GEO數(shù)據(jù)挖掘滾動(dòng)開班

• DNA及RNA甲基化數(shù)據(jù)分析與課題設(shè)計(jì)