信息時(shí)代就是好，沒有自己的數(shù)據(jù)咱可以從免費(fèi)的數(shù)據(jù)庫找數(shù)據(jù)呀！目前 TCGA 和 GEO 數(shù)據(jù)庫簡直就是專為科研人員準(zhǔn)備的免費(fèi)的寶庫。前面一直講的空間轉(zhuǎn)錄組，今天來穿插著講一節(jié) TCGA & GEO 數(shù)據(jù)分析換換口味吧。
今天主要給大家介紹一下對(duì)轉(zhuǎn)錄組 TCGA & GEO 的數(shù)據(jù)進(jìn)行分析，包括差異分析、富集分析、單基因展示以及其它個(gè)性化做圖。除了想利用公共數(shù)據(jù)庫挖掘轉(zhuǎn)錄組數(shù)據(jù)的發(fā)文章的同學(xué)之外，這篇教程其實(shí)特別適合自己做了少量轉(zhuǎn)錄組樣本又不想做后期驗(yàn)證，想直接用公共數(shù)據(jù)庫的數(shù)據(jù)來驗(yàn)證自己數(shù)據(jù)結(jié)果的朋友。
目前各種 TCGA & GEO 數(shù)據(jù)下載的教程太多了，所以今天直接省略這部分部分，直接從后面的分析開始。廢話不多說，直接上干貨！
一、準(zhǔn)備文件
基因表達(dá)矩陣文件，可以是 TCGA 下載的轉(zhuǎn)錄組數(shù)據(jù)，也可以是 GEO 下載的轉(zhuǎn)錄組數(shù)據(jù)，格式如下（行是基因，列是樣本）：
樣本分組文件，格式如下（type 是分組）：
二、差異分析
一般 TCGA 或 GEO 上大量樣本數(shù)據(jù)的差異分析 limma 用的比較多，所以這里也使用 limma 包來分析差異基因，注意需要先安裝 limma R 包。
安裝 limma 包：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("limma")
分析差異：
## 讀取基因表達(dá)件，取 log2
geneexp = read.table("gene_exp_BRCA.txt",header=T,row.names=1,sep="\t")
geneexp = log2(geneexp+1)
## 讀取分組文件
group_file = read.table("sample_group_BRCA.txt",header=T,row.names=1,sep="\t",as.is =TRUE)
rownames(group_file) = gsub('-','.',rownames(group_file))
geneexp = geneexp[,rownames(group_file)]
### 差異分組設(shè)置
samps<-factor(group_file$type)
design <- model.matrix(~0+samps) ;
colnames(design) <- levels(samps)
### 模型擬合
fit <- lmFit(geneexp, design)
cont.matrix<-makeContrasts(Basal-Normal,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
final<-topTable(fit2, coef=1, number=dim(geneexp)[1], adjust.method="BH")
> head(final)
logFC AveExpr t P.Value adj.P.Val B
SDPR -4.358861 2.679133 -46.10828 8.888505e-122 4.738017e-117 267.4862
TPX2 4.030685 3.807811 45.39822 2.571983e-120 6.854979e-116 264.1584
UBE2C 4.089727 3.270516 43.72542 8.428935e-117 1.497681e-112 256.1483
CDC20 3.971916 3.450304 42.85468 6.270780e-115 8.356598e-111 251.8811
FIGF -3.220704 1.583802 -40.93757 1.053838e-110 1.123497e-106 242.2399
FAXDC2 -2.179867 1.637909 -39.05494 2.084335e-106 1.851758e-102 232.4283
三、繪制差異基因散點(diǎn)圖和火山圖
這一步需要用到 ggplot2 繪圖，沒有安裝的話需要先安裝 ggplot2 包。
# 繪制差異基因散點(diǎn)圖和火山圖
library(ggplot2)
g1 = "Normal"
g2 = "Basal"
g1_exp = geneexp[rownames(final),rownames(group_file)[which(group_file$type==g1)]]
g2_exp = geneexp[rownames(final),rownames(group_file)[which(group_file$type==g2)]]
g1_mean = apply(g1_exp,1,mean)
g2_mean = apply(g2_exp,1,mean)
type=rep('No',length(g1_mean))
type[which(final$logFC> 1 & final$adj.P.Val <0.05)] = "Up"
type[which(final$logFC < -1 & final$adj.P.Val < 0.05)] = "Down"
datam = data.frame(g1_mean,g2_mean,logFC=final$logFC,FDR=final$adj.P.Val,type,stringsAsFactors=FALSE)
## 散點(diǎn)圖
ggplot(datam,aes(g1_mean,g2_mean,colour=type))+
geom_point(stat="identity",size=1)+theme(legend.title=element_blank())+scale_color_manual(values =c("Down"='blue',"No"='grey',"Up"='orange'))+
labs(x=paste(g1,'Log2(FPKM+1)'),y=paste(g2,'Log2(FPKM+1)'),title=paste(g2,'VS',g1,sep=""))+
coord_cartesian(ylim=c(0,10),xlim=c(0,10))+geom_segment(aes(x = 0, y = 0, xend = 10, yend = 10),size=1,colour="#999999",linetype="dotted")+theme(plot.title = element_text(hjust = 0.5),title=element_text(face="bold",size=15,colour="black"),axis.title=element_text(face="bold",size=13,colour="black"),axis.text.x=element_text(face="bold",size=12,colour="black"),axis.text.y=element_text(face="bold",size=12,colour="black"),legend.text=element_text(face="bold",size=13,colour="black"))
散點(diǎn)圖結(jié)果如下：
## 繪制火山圖
ggplot(datam,aes(logFC,-log10(FDR),colour=type))+
geom_point(stat="identity",size=1.2)+theme(legend.title=element_blank())+scale_color_manual(values =c("Down"='blue',"No"='grey',"Up"='orange'))+
labs(x="Log2 (FC)",y="-Lg10 (FDR)",title=paste(g2,'VS',g1,sep=""))+coord_cartesian(xlim=c(-5,5))+
geom_hline(aes(yintercept=1.3),colour="white",size=1.1)+
geom_vline(aes(xintercept =-1),colour="white",size=1.1)+geom_vline(aes(xintercept =1),colour="white",size=1.1)+
theme(axis.title.y = element_text(vjust=-0.1),axis.title.x = element_text(vjust=-0.6),title = element_text(vjust=0.8))+theme(plot.title = element_text(hjust = 0.5),title=element_text(face="bold",size=15,colour="black"),axis.title=element_text(face="bold",size=13,colour="black"),axis.text.x=element_text(face="bold",size=10,colour="black"),axis.text.y=element_text(face="bold",size=10,colour="black"),legend.text=element_text(face="bold",size=12,colour="black"))
火山圖結(jié)果如下：
四、差異基因富集分析
這里我們介紹一下怎么用 clusterProfiler R 包來做差異基因富集分析。需要先安裝好 clusterProfiler，org.Hs.eg.db 兩個(gè) R 包。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
GO 富集分析
library(clusterProfiler)
library(org.Hs.eg.db)
### 提取差異基因 list
diffgenes <- final[(final[,"adj.P.Val"]<0.05 & abs(final[,"logFC"])>=1),]
genelist <- diffgenes[,"logFC"]
names(genelist) = rownames(diffgenes)
##id 轉(zhuǎn)換，將 SYMBOL 轉(zhuǎn)換成 ENTREZID
gene = bitr(rownames(diffgenes), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")
#####################GO 富集
ego <- enrichGO(
gene = gene$ENTREZID,
keyType = "ENTREZID",
OrgDb = "org.Hs.eg.db",
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE)
# 再用 setReadable 函數(shù)將基因 ID 映射到基因 Symbol
ego2 <- setReadable(ego, OrgDb = "org.Hs.eg.db", 'ENTREZID')
繪制柱狀圖：
barplot(ego2,showCategory = 20)
繪制氣泡圖：
dotplot(ego2, showCategory = 20)
繪制網(wǎng)絡(luò)圖 1：
cnetplot(ego2, showCategory = 3,foldChange=genelist, circular = TRUE, colorEdge = TRUE)
繪制網(wǎng)絡(luò)圖 2：
emapplot(ego2, foldChange=genelist, showCategory = 20)
KEGG 富集分析
注意這一步需要連接網(wǎng)絡(luò)，因?yàn)?clusterProfiler 是在線抓取新的 pathway 數(shù)據(jù)庫的。當(dāng)然也有用 kegg.db R 包直接內(nèi)置數(shù)據(jù)庫的情況，這里不做介紹。
######################KEGG 富集
kegg <- enrichKEGG(
gene = gene$ENTREZID,
keyType = "kegg",
organism = 'hsa',
pvalueCutoff = 1,
pAdjustMethod = "BH",
qvalueCutoff = 1,
use_internal_data = FALSE)
kegg2 <- setReadable(kegg, OrgDb = "org.Hs.eg.db", 'ENTREZID')
繪制柱狀圖：
barplot(kegg2,showCategory = 20)
繪制氣泡圖：
dotplot(kegg2, showCategory = 20)
繪制網(wǎng)絡(luò)圖 1：
cnetplot(kegg2, showCategory = 10,foldChange=genelist, circular = TRUE, colorEdge = TRUE)
繪制網(wǎng)絡(luò)圖 2：
emapplot(kegg2, foldChange=genelist, showCategory = 20)
五、差異基因熱圖標(biāo)記基因
當(dāng)我們已經(jīng)用自己的數(shù)據(jù)做過分析得到了差異基因，或者是已經(jīng)有自己關(guān)注的基因，想用 TCGA & GEO 數(shù)據(jù)來驗(yàn)證自己的數(shù)據(jù)或結(jié)論的時(shí)候，用 TCGA & GEO 數(shù)據(jù)差異基因熱圖同時(shí)標(biāo)記特定基因來展示結(jié)果就特別合適。
###### 差異基因熱圖標(biāo)記關(guān)注基因
library("ComplexHeatmap")
library("circlize")
diff_exp = geneexp[rownames(diffgenes),]
# 進(jìn)行 zscore 歸一化
diff_exp_scaled = t(apply(diff_exp, 1, scale))
colnames(diff_exp_scaled) = colnames(diff_exp)
# 需要標(biāo)記的基因
genes = c('CDC20','THRB','FAM13A','CCNA2','PRR15','CHI3L1','CLCA2','ABCA10','A2ML1','LCN2')
## 設(shè)置位置
gene_pos = which(rownames(diff_exp) %in% genes)
ha = rowAnnotation(foo = anno_mark(at = gene_pos, labels = genes))
m=round(max(abs(diff_exp_scaled)))
## 畫熱圖
Heatmap(diff_exp_scaled, colorRamp2(c(-m,0,m),c("blue", "#EEEEEE", "red")),
right_annotation = ha,name = "Z-score",show_row_names =FALSE,show_column_names =FALSE)
繪圖結(jié)果如下：
六、單個(gè)基因繪圖
要在 TCGA & GEO 數(shù)據(jù)中驗(yàn)證自己的關(guān)注的基因的差異情況，除了前面說的差異基因熱圖標(biāo)記特定基因之外，也可以對(duì)單個(gè)基因直接進(jìn)行繪圖。這里我們提供三種展示方式，選擇自己喜歡的一種展示形式就行。
用箱線圖展示：
## 單個(gè)基因箱線圖
exp2 = data.frame(t(diff_exp))
exp2$type = group_file$type
ggplot(exp2, aes(x=type, y=CDC20,fill=type)) +
geom_boxplot()+geom_jitter(alpha = .3, width =0.2,size=1)+labs(title = "") + ylab("log2 (FPKM)")+ xlab("")+theme_bw()+
theme(title=element_text(face="bold",size=16),axis.title=element_text(face="bold",size=15),axis.text.x=element_text(face="bold",angle=80,size=13,hjust=1),
axis.text.y=element_text(face="bold",size=12),legend.text=element_text(face="bold",size=13),legend.title=element_text(face="bold",size=13))
結(jié)果如下：
小提琴圖展示：
### 單個(gè)基因小提琴圖
ggplot(exp2, aes(x=type, y=CDC20,fill=type)) +
geom_violin(alpha=0.5) + geom_boxplot(alpha = .5,fill="white", width= .2)+labs(title = "") + ylab("log2 (FPKM)")+ xlab("")+theme_bw()+
theme(title=element_text(face="bold",size=16),axis.title=element_text(face="bold",size=15),axis.text.x=element_text(face="bold",angle=80,size=13,hjust=1),
axis.text.y=element_text(face="bold",size=12),legend.text=element_text(face="bold",size=13),legend.title=element_text(face="bold",size=13))
結(jié)果如下：
散點(diǎn)圖展示：
### 單個(gè)基因散點(diǎn)圖
ggplot(exp2, aes(x=type, y=CDC20)) +
geom_dotplot(binaxis='y', stackdir='center',binwidth = 0.1)+stat_summary(fun.y=median, geom="point", shape=18,size=3, color="red")+
theme(axis.title=element_text(face="bold",size=15),axis.text.x=element_text(face="bold",angle=80,size=13,hjust=1),axis.text.y=element_text(face="bold",size=12))
結(jié)果如下：
這里只展示一個(gè)基因的情況，如果是有多個(gè)基因需要繪圖的，把代碼里的基因名替換一下就好了。
好啦，今天的分享就到次為止啦，下次再繼續(xù)哦！
更多伯豪生物人工服務(wù)：
