全部流程來(lái)自：GEO數(shù)據(jù)庫(kù)挖掘—生信技能樹B站視頻，建議去看原文！

第一步：找到相關(guān)的GEO數(shù)據(jù)集（文獻(xiàn)/搜索），以胃癌gastric cancer為例

可去文獻(xiàn)中查找，用于練習(xí)

第二步：運(yùn)行R包GEOquery獲取數(shù)據(jù)（非常看網(wǎng)速，盡量下載下一點(diǎn)的包）

library(GEOquery)

eSet <- getGEO("GSE118916",

? ? ? ? ? ? ? destdir = '.',? ?#下載在當(dāng)前目錄

? ? ? ? ? ? ? getGPL = F)??#平臺(tái)信息不要

#查看下載得到的文件大小，與GEO網(wǎng)頁(yè)中標(biāo)識(shí)的是否一致

第三步：得到表達(dá)矩陣

#使用列表取子集的方法提取eSet列表里的第一個(gè)元素：eSet[[1]]；并使用exprs函數(shù)把它轉(zhuǎn)化成矩陣：exp <- exprs(eSet[[1]])

eSet[[1]]

exp <- exprs(eSet[[1]])

exp[1:4,1:4]#查看前四行

#代碼邏輯，把probe_id轉(zhuǎn)換為--->symbol ID/entrez ID

#分兩步走：過濾probe_id，得到每個(gè)基因所對(duì)應(yīng)的唯一的probe_id

得到probe_id與symbol ID這件的轉(zhuǎn)換關(guān)系

ID轉(zhuǎn)換的第一步必須要加載特定的R包，下載哪個(gè)包，需要根據(jù)GPL來(lái)定

eSet[[1]]@annotation #查看GPL平臺(tái)信息，即芯片信息

#Google直接搜索相關(guān)R包，安裝并調(diào)用

#如果有對(duì)應(yīng)的R包

library(R包)

#如果找不到R包

library(Biobase)

library(GEOquery)

gpl<-getGEO('GPL15297',destdir=".")

colnames(Table(gpl))

head(Table(gpl)[,c(1,15,19)])## you need to check this , which column do you need

probe2symbol=Table(gpl)[,c(1,15)]

write.csv(Table(gpl)[,c(1,15)],"GPL15207.csv")

tail(sort(table(ids$symbol)))#看高頻數(shù)的樣本

table(rownames(exp)%in%ids$probe_id)#查看exp與id是否一致

exp=exp[rownames(exp)%in%ids$probe_id,]#去除不一致的項(xiàng)

ids=ids[match(rownames(exp),ids$probe_id),]#使exp與id完全對(duì)應(yīng)

tmp=by(exp,ids$symbol,function(x)rownames(x)[which.max(rowMeans(x))])

probes=as.character(tmp)

dim(exp)

exp=exp[rownames(exp)%in%probes,]# 過濾有多個(gè)探針的基因

dim(exp)

rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]

head(exp)

第四步：獲取分組信息

pd <-pData(eSet[[1]])

library(stringr)

group_list=ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")group_list

str_detect(pd$title,"Control")

方法二：自己根據(jù)CSV里文件類型自己寫代碼

group_list=c(rep("control",times=3),rep("treat",times=3))

第五步：檢查表達(dá)矩陣

方法一：管家基因

exp['GAPDH',]#檢驗(yàn)用的管家基因

exp['ACTB',]#同上

boxplot(exp)

#使用ggplot2畫各個(gè)樣本表達(dá)量的boxplot圖

# 準(zhǔn)備畫圖所需數(shù)據(jù)exp_L

library(reshape2)

head(exp)

exp_L=melt(exp)

head(exp_L)

colnames(exp_L)=c('symbol','sample','value')

head(exp_L)

# 獲得分組信息

library(stringr)

group_list=ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")

group_list

exp_L$group=rep(group_list,each=nrow(exp))

head(exp_L)

# ggplot2畫圖 library(ggplot2)

p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()

print(p)

##boxplot圖精修版p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()

p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")

p=p+theme_set(theme_set(theme_bw(base_size=20)))

p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())

print(p)

#查案畫出的箱型圖，如果均值不在一條線上，說明樣本有批次效應(yīng)，需要人工校正

library(limma)

exp=normalizeBetweenArrays(exp)

#除了箱型圖還可以畫其他的圖

p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_violin()

print(p)

p=ggplot(exp_L,aes(value,fill=group))+geom_histogram(bins=200)+facet_wrap(~sample,nrow=4)

print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density()+facet_wrap(~sample,nrow=4)

print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density()

print(p)

第四步：檢查樣本信息，檢查樣本分組信息，一般看PCA圖，hclust圖

# 更改表達(dá)矩陣列名

head(exp)

colnames(exp)=paste(group_list,1:6,sep='')

head(exp)

# 定義nodePar

nodePar<-list(lab.cex=0.6,pch=c(NA,19),cex=0.7,col="blue")

# 聚類

hc=hclust(dist(t(exp)))

par(mar=c(5,5,5,10))

# 繪圖

plot(as.dendrogram(hc),nodePar=nodePar,horiz=TRUE)

#繪PCA圖

library(ggfortify)

# 互換行和列，再dim一下df=as.data.frame(t(exp))

# 不要view df，列太多，軟件會(huì)卡??；

dim(df)

dim(exp)

exp[1:6,1:6]

df[1:6,1:6]

df$group=group_list

autoplot(prcomp(df[,1:ncol(df)-1)]),data=df,colour='group')

#保存數(shù)據(jù)

save(exp,group_list,file="step2output.Rdata")

第五步：差異分析

#limma包需要以下三個(gè)數(shù)據(jù)，表達(dá)矩陣(exp)；分組矩陣(design)；差異比較矩陣(contrast.matrix)

#表達(dá)矩陣(exp)我們?cè)缇偷玫搅?，不用再制作了；我們也得到了存放分組信息的向量group_list，用它來(lái)制作我們的分組矩陣

rm(list=ls())## 魔幻操作，一鍵清空

options(stringsAsFactors=F)

load(file="step2output.Rdata")#上步得到的EXP分組信息

dim(exp)

library(limma)# 做分組矩陣?

design<-model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exp)

design#查看得到的分組矩陣

#輸入數(shù)據(jù)—差異比較矩陣

contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse="-"),levels=design)

contrast.matrix

#limma包做差異分析，只有三個(gè)步驟：lmFit，eBayes，topTable

##step1

fit<-lmFit(exp,design)

##step2

fit2<-contrasts.fit(fit,contrast.matrix)##這一步很重要，大家可以自行看看效果fit2<-eBayes(fit2) ##default no trend!!!

##eBayes()withtrend=TRUE

##step3

tempOutput=topTable(fit2,coef=1,n=Inf)nrDEG=na.omit(tempOutput)

#可以保存數(shù)據(jù) write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)

head(nrDEG)

save(nrDEG,file="DEGoutput.Rdata") #保存差異數(shù)據(jù)nrDEG

#檢查差異數(shù)據(jù)的準(zhǔn)確性，進(jìn)行可視化

#熱圖

rm(list=ls())## 魔幻操作，一鍵清空

options(stringsAsFactors=F)

load(file="DEGoutput.Rdata")

load(file="DEGinput.Rdata")

library(pheatmap)

choose_gene=head(rownames(nrDEG),25)

choose_matrix=exp[choose_gene,]

choose_matrix=t(scale(t(choose_matrix)))

pheatmap(choose_matrix)

#火山圖

rm(list=ls()) ## 魔幻操作，一鍵清空

options(stringsAsFactors=F)

load(file="DEGoutput.Rdata")

colnames(nrDEG)

plot(nrDEG$logFC,-log10(nrDEG$P.Value))

DEG=nrDEG

logFC_cutoff<-with(DEG,mean(abs(logFC))+2*sd(abs(logFC)))DEG$change=as.factor(ifelse(DEG$P.Value<0.05&abs(DEG$logFC)>logFC_cutoff,ifelse(DEG$logFC>logFC_cutoff,'UP','DOWN'),'NOT'))

head(DEG)

g=ggplot(data=DEG,aes(x=logFC,y=-log10(P.Value),color=change))+geom_point(alpha=0.4,size=1.75)+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'))

print(g)

第六步：富集分析

#KEGG富集分析

library(clusterProfiler)

gene_up=deg[deg$change=='up','ENTREZID']

gene_down=deg[deg$change=='down','ENTREZID']

gene_diff=c(gene_up,gene_down)

gene_all=deg[,'ENTREZID']

kk.up<-enrichKEGG(gene=gene_up,organism='hsa',universe=gene_all,pvalueCutoff=0.9,qvalueCutoff=0.9)

head(kk.up)[,1:6]

dim(kk.up)

kk.down<-enrichKEGG(gene=gene_down,organism='hsa',universe=gene_all,pvalueCutoff=0.9,qvalueCutoff=0.9)

head(kk.down)[,1:6]

dim(kk.down)

kk.diff<-enrichKEGG(gene=gene_diff,organism='hsa',pvalueCutoff=0.05)

head(kk.diff)[,1:6]

kegg_diff_dt<-kk.diff@result

down_kegg<-kk.down@result%>%

? ? filter(pvalue<0.05) %>%mutate(group=-1)

up_kegg<-kk.up@result%>%

? ? filter(pvalue<0.05) %>%mutate(group=1)

kegg_plot<-function(up_kegg,down_kegg)

{

dat=rbind(up_kegg,down_kegg)colnames(dat)

dat$pvalue=-log10(dat$pvalue)

dat$pvalue=dat$pvalue*dat$group

dat=dat[order(dat$pvalue,decreasing=F),]

g_kegg<-ggplot(dat,aes(x=reorder(Description,order(pvalue,decreasing=F)),y=pvalue,fill=group))+geom_bar(stat="identity")+scale_fill_gradient(low="blue",high="red",guide=FALSE)+scale_x_discrete(name="Pathway names")+scale_y_continuous(name="log10P-value")+coord_flip()+theme_bw()+theme(plot.title=element_text(hjust=0.5))+ggtitle("Pathway Enrichment")

}

g_kegg<-kegg_plot(up_kegg,down_kegg)

g_kegg

ggsave(g_kegg,filename='kegg_up_down.png')#繪圖并保存

#GSEA是另一個(gè)常用的富集分析方法，目的是看看基因全局表達(dá)量的變化是否有某些特定的基因集合的傾向性。

data(geneList,package="DOSE")

head(geneList)

length(geneList)

names(geneList)

boxplot(geneList)

boxplot(deg$logFC)

geneList=deg$logFC

names(geneList)=deg$ENTREZID

geneList=sort(geneList,decreasing=T)

kk_gse<-gseKEGG(geneList=geneList,organism='hsa',nPerm=1000,minGSSize=120,pvalueCutoff=0.9,verbose=FALSE)

head(kk_gse)[,1:6]

gseaplot(kk_gse,geneSetID=rownames(kk_gse[1,]))

down_kegg<-kk_gse[kk_gse$pvalue<0.05&kk_gse$enrichmentScore<0,];down_kegg$group=-1

up_kegg<-kk_gse[kk_gse$pvalue<0.05&kk_gse$enrichmentScore>0,];up_kegg$group=1

gse_kegg=kegg_plot(up_kegg,down_kegg)

print(gse_kegg)

ggsave(gse_kegg,filename='kegg_up_down_gsea.png')

#GO富集分析

GO富集分析原理：有一個(gè)term注釋了100個(gè)差異表達(dá)基因參與了哪個(gè)過程，注釋完之后（模式生物都有現(xiàn)成的注釋包，不用我們自己注釋），計(jì)算相對(duì)于背景它是否顯著集中在某條通路、某一個(gè)細(xì)胞學(xué)定位、某一種生物學(xué)功能。

對(duì)GO database analysis一般從三個(gè)層面進(jìn)行：

Cellular component，CC? 細(xì)胞成分

Biological process， BP? 生物學(xué)過程

Molecular function，MF? 分子功能

這三個(gè)層面具體是指：

Cellular component解釋的是基因存在在哪里，在細(xì)胞質(zhì)還是在細(xì)胞核？如果存在細(xì)胞質(zhì)那在哪個(gè)細(xì)胞器上？如果是在線粒體中那是存在線粒體膜上還是在線粒體的基質(zhì)當(dāng)中？這些信息都叫Cellular component。

Biological process是在說明該基因參與了哪些生物學(xué)過程，比如，它參與了rRNA的加工或參與了DNA的復(fù)制，這些信息都叫Biological process

Molecular function在講該基因在分子層面的功能是什么？它是催化什么反應(yīng)的？

library(clusterProfiler)

gene_up=deg[deg$change=='up','ENTREZID']

gene_down=deg[deg$change=='down','ENTREZID']

gene_diff=c(gene_up,gene_down)

head(deg)

#細(xì)胞組分

ego_CC<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="CC",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

#生物過程，重要

ego_BP<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="BP",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

#分子功能，很重要

ego_MF<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="MF",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

save(ego_CC,ego_BP,ego_MF,file="ego_GPL6244.Rdata")

rm(list=ls())

load(file="ego_GPL6244.Rdata")

#第一種，條帶圖，按p從小到大排的#如果運(yùn)行了沒出圖，就dev.new()

barplot(ego_CC,showCategory=20,title="EnrichmentGO_CC")

barplot(ego_BP,showCategory=20,title="EnrichmentGO_CC")

#第二種，點(diǎn)圖，按富集數(shù)從大到小的dotplot(ego_CC,title="EnrichmentGO_BP_dot")

#保存

pdf(file="dotplot_GPL6244.pdf")

dotplot(ego_CC,title="EnrichmentGO_BP_dot")

dev.off()