前兩次教程“GEO數(shù)據(jù)挖掘（一）簡(jiǎn)單快速下載GEO數(shù)據(jù)”、“GEO數(shù)據(jù)挖掘（一）下載SRA庫(kù)原始測(cè)序數(shù)據(jù)”已經(jīng)分享過(guò)如何下載數(shù)據(jù)，本次將分享如何將geo數(shù)據(jù)內(nèi)的差異基因提取出來(lái)，進(jìn)行差異化分析，找到上調(diào)下調(diào)基因，再進(jìn)行可視化展示。首先，我們得知道為什么需要對(duì)基因進(jìn)行差異表達(dá)分析。因?yàn)樵谘芯磕[瘤等疾病的發(fā)病過(guò)程中，一些原本沉默的基因可能高表達(dá)，而一些原本正常表達(dá)的基因，其表達(dá)量可能就會(huì)下降，而恰恰就是這些相較于正常樣本基因表達(dá)量發(fā)生變化的基因，它們控制或者影響著腫瘤的發(fā)生。因此，如果要研究腫瘤發(fā)生的機(jī)制，需要設(shè)立case（實(shí)驗(yàn)組）和control（對(duì)照組）進(jìn)行差異差分析，探究表征實(shí)驗(yàn)組和對(duì)照組的差異表達(dá)基因，構(gòu)建基因富集通路研究腫瘤發(fā)病機(jī)制，然后重點(diǎn)解讀這些富集通路，并探究這些通路與樣本的表型之間有著怎樣的關(guān)系或者內(nèi)在機(jī)制。
那么接下來(lái)將分四步為大家?guī)?lái)基因差異分析及可視化的全套代碼分享：
第一步：提取表達(dá)矩陣、臨床信息、芯片編號(hào)
使用列表取子集的方法提取eSet列表（我們下載好的原始數(shù)據(jù)）里的第一個(gè)元素：eSet[[1]]；并使用exprs函數(shù)把它轉(zhuǎn)化成矩陣；通常使用：head（100），判斷數(shù)據(jù)是否需要轉(zhuǎn)換為log型。
#(1)提取表達(dá)矩陣exp----
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
#exp = log2(exp+1) #判斷是否需要log2
#(2)提取臨床信息----
使用pData()函數(shù)可以得到每個(gè)樣本的臨床信息，通常數(shù)據(jù)框的來(lái)源列（source）或者特征列（characterizes）會(huì)描述哪些樣本是control或者treatment。
pd <- pData(eSet[[1]])
#(3)提取芯片平臺(tái)編號(hào)----
數(shù)據(jù)通常使用的是不同的芯片探針，后續(xù)根據(jù)GPL平臺(tái)編號(hào)轉(zhuǎn)化成entrez ID或symbol ID；
gpl <- eSet[[1]]@annotation
p = identical(rownames(pd),colnames(exp));p
save(gse,exp,pd,gpl,file = "step1_output.Rdata")#儲(chǔ)存結(jié)果
第二步 構(gòu)建分組、芯片注釋
1.構(gòu)建分組信息
根據(jù)所GSE內(nèi)的樣本來(lái)源信息進(jìn)行分組。
rm(list = ls())
load("step1_output.Rdata")
library(stringr)
table(colnames(exp))
group_list = ifelse(str_detect(pd$source_name_ch1," tumor"),"test","normal")
group_list=factor(group_list,
? ? ? ? ? ? ? ? ? levels = c("test","normal"))
group_list
table(group_list)
2.芯片注釋
if(T){
? a = getGEO(gpl,destdir = ".")
? b = a@dataTable@table
? colnames(b)
? ids2 = b[,c("ID","GENE_SYMBOL")]
? colnames(ids2) = c("probe_id","symbol")
? ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"http:///"),]
}
save(group_list,ids2,exp,pd,file = "step2_output.Rdata")
第三步 數(shù)據(jù)檢查---PCA
然后利用FactoMineR，factoextra包進(jìn)行主成分分析，查看處理和對(duì)照是否有明顯分組。
繪制主成分分析圖
rm(list = ls())#清空環(huán)境
load("step1_output.Rdata")
load("step2_output.Rdata")
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra)
# pca的統(tǒng)一操作
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
? ? ? ? ? ? ? ? ? ? ? ? geom.ind = "point", # show points only (nbut not "text")
? ? ? ? ? ? ? ? ? ? ? ? col.ind = group_list, # color by groups
? ? ? ? ? ? ? ? ? ? ? ? palette = c("#00AFBB", "#E7B800"),
? ? ? ? ? ? ? ? ? ? ? ? addEllipses = TRUE, # Concentration ellipses
? ? ? ? ? ? ? ? ? ? ? ? legend.title = "Groups"
)
print(pca_plot)
ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png"))
save(pca_plot,file = "pca_plot.Rdata")
第四步 差異化分析及可視化展示—火山圖、熱圖
導(dǎo)入前兩步的Rdata，使用limma包進(jìn)行差異分析。簡(jiǎn)單來(lái)說(shuō)，對(duì)基因差異分析就是對(duì)每個(gè)基因都進(jìn)行差異檢驗(yàn)，檢驗(yàn)基因的logFC值、平均表達(dá)量、P.value是否顯著等。差異分析的輸入文件：1.表達(dá)矩陣；2.分組信息。Limma包做差異分析分為三個(gè)步驟：lmFit,、eBayes,、topTable。使用這個(gè)包需要三個(gè)數(shù)據(jù)：1.表達(dá)矩陣；2.差異比較矩陣；3.差異分析矩陣，分析結(jié)果需要重點(diǎn)看logFC值和P值。
rm(list = ls())
load("step1_output.Rdata")
load("step2_output.Rdata")
1.差異分析
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
head(deg)
2.為deg數(shù)據(jù)框添加幾列
#1.加probe_id列，把行名變成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
#tibble::rownames_to_column()
head(deg)
#merge將兩個(gè)表整合起來(lái)
table(deg$probe_id %in% ids$probe_id)
#deg <- inner_join(deg,ids,by="probe_id")
deg <- merge(x = deg,y = ids2, by="probe_id")
deg <- deg[!duplicated(deg$symbol),]
dim(deg)
#2.加change列：上調(diào)或下調(diào)，火山圖要用
#logFC_t=mean(deg$logFC)+2*sd(deg$logFC)
logFC_t=1.5
change=ifelse(deg$P.Value>0.05,'stable',
? ? ? ? ? ? ? ifelse( deg$logFC >logFC_t,'up',
? ? ? ? ? ? ? ? ? ? ? ifelse( deg$logFC < -logFC_t,'down','stable') )
)
deg <- mutate(deg,change)
head(deg)
table(deg$change)
3.加ENTREZID列，后面富集分析要用
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(unique(deg$symbol), fromType = "SYMBOL",
? ? ? ? ? ? toType = c( "ENTREZID"),
? ? ? ? ? ? OrgDb = org.Hs.eg.db)
head(s2e)
head(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
head(deg)
save(logFC_t,deg,file = "step4_output.Rdata")
數(shù)據(jù)可視化-----火山圖、熱圖
rm(list = ls())
load("step1_output.Rdata")
load("step2_output.Rdata")
load("step4_output.Rdata")
library(dplyr)
1.火山圖
dat <- mutate(deg,v=-log10(P.Value))
if(T){
? for_label <- dat%>%
? ? filter(symbol %in% c("RUNX2","FN1"))
}
if(F){
? for_label <- dat %>% head(10)
}
if(F) {
? x1 = dat %>%
? ? filter(change == "up") %>%
? ? head(3)
? x2 = dat %>%
? ? filter(change == "down") %>%
? ? head(3)
? for_label = rbind(x1,x2)
}
p <- ggplot(data = dat,
? ? ? ? ? ? aes(x = logFC,
? ? ? ? ? ? ? ? y = v)) +
? geom_point(alpha=0.4, size=3.5,
? ? ? ? ? ? aes(color=change)) +
? ylab("-log10(Pvalue)")+
? scale_color_manual(values=c("blue", "grey","red"))+
? geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
? geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) +
? theme_bw()
p
volcano_plot <- p +
? geom_point(size = 3, shape = 1, data = for_label) +
? ggrepel::geom_label_repel(
? ? aes(label = symbol),
? ? data = for_label,
? ? color="black"
? )
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png"))
2.繪制熱圖
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(pheatmap)
heatmap_plot <- pheatmap(n,
? ? ? ? show_colnames =F,
? ? ? ? show_rownames = F,
? ? ? ? annotation_col=annotation_col,
? ? ? ? scale = "row")
#保存結(jié)果
library(ggplot2)
png(file = paste0(gse,"heatmap.png"))
ggsave(plot = heatmap_plot,filename = paste0(gse,"heatmap.png"))
dev.off()
目前我們已經(jīng)完成了GEO數(shù)據(jù)挖掘的基因差異分析及可視化部分，已經(jīng)拿到了上調(diào)下調(diào)的顯著差異表達(dá)基因，但是疾病的發(fā)生往往不單單受某個(gè)基因的影響，如果想繼續(xù)了解疾病癌癥的發(fā)生機(jī)制，還需要大家繼續(xù)關(guān)注后續(xù)即將分享的“基因富集通路”方面教程。
請(qǐng)持續(xù)關(guān)注“GEO數(shù)據(jù)挖掘”系列文章，每周一個(gè)實(shí)用干貨帶您上手生信分析。
