GEPIA2详解（中国智造-肿瘤数据库）
记得某一次讲座上，听到北京大学张泽民教授演讲时提到了他们实验室开发的GEPIA(Gene Expression Profiling Interactive Analysis)，基因表达谱数据动态分析网页工具，不知道为什么获得了几百的引用，毕竟他老人家引以为傲的工作应该是各大癌症的单细胞CNS文章啊！ 我这里斗胆猜测一下，我们生信技能树果子3年前的推文：中国制造:碉堡的TCGA可视化网站GEPIA 当居首功！有趣的是，开发这个工具的第一作者唐泽方已经入职IBM，与生物信息学渐行渐远，果子呢，也差不多退出生信技能树，自立门户了，他的果子学生信在一小撮人那边的牌子也很响亮。 三年过去了，GEPIA早就更新到了第二版，很多粉丝后台留言希望我们再介绍一波，但是网页工具的介绍，我的确写不出来，恰好有学徒的学徒投稿，所以加急分享出来！ 下面是正文：
图片来源：GEPIA2
作者：李瑞萌
审校：Jimmy
GEPIA2 是北京大学张泽民老师实验室开发的一个网站，能够对TCGA和GTEx项目共9736个肿瘤样本、8587个正常样本的RNA－seq表达数据进行分析。目前该网站已经有两篇文章发表。
参考文献：
Tang, Z. et al. (2017) GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res, 10.1093/nar/gkx247.
Tang, Z. et al. (2019) GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res, 10.1093/nar/gkz430.
GEPIA2的两大功能：
Expression Analysis
Custom Data Analysis
1
Expression Analysis
细分为8个功能
General
Differential Genes
Expression DIY
Survival Analysis
Isoform Details
Correlation Analysis
Similar Genes Detection
Dimensionality Reduction
1
General
在搜索框内输入感兴趣的gene symbol或者Ensembl ID，可以得到该gene及其isoform的详细信息，并且以body map、bar plot、dot plot形式表示其在肿瘤样本和正常样本中的表达情况。当然也可以输入Isoform symbol或者Isoform ID。
General
其他7个功能的实现也很简单：输入基因名称（或者isoform、gene signature)，选择癌症数据集，设置一些参数，即可得到基因列表或者可视化结果。以‘Expression DIY’为例，如下图：
这个网站用起来非常友好，哪里不会点哪里。不理解参数，点击’help‘；忘记肿瘤名称了，点击Cancer name。而点击‘example’，会弹出一个新网页 “Examples for GEPIA2 Usage“；这个网页提供了一些用于可视化的Rscript代码。
2
Differential Genes
在某一肿瘤/正常组织中差异表达的基因或者isoform，并且显示这些基因在染色体上的位置分布。
差异基因在染色体上的位置分布
3
Expression DIY
可以画四种图
(1) profile：用dot plot分析一个基因或者isoform在不同肿瘤样本（和正常样本）中的表达情况。
FOXD1基因的表达情况
(2) Box Plot：分析一个基因或者isoform或者a multi-gene signature在不同肿瘤样本和正常样本中的表达情况。也可以对其在某一肿瘤不同亚型中的表达情况进行研究，如下图。
CD163基因在BRCA三种亚型的表达情况
(3) Stage Plot：用violin plot分析一个基因或者isoform在肿瘤不同阶段的表达情况。
(4) Multiple Genes Comparison：用heatmap分析多个基因在不同肿瘤样本（和正常样本）中的表达情况。
4
Survival Analysis
(1) Survival Analysis: 一个基因、isoform或者a multi-gene signature在任意癌症中的OS或者DFS。也可以对其在任意肿瘤不同亚型中的OS或者DFS进行研究。
FOXD1基因在BRCA的OS
(2) Most Differential Survival Genes: 获得在某种癌症中，与生存相关的基因/isoform列表。
(3) Survival Map: 用heatmap表示多个基因或者isoform在多种癌症中的生存分析结果。
The heat map shows the hazard ratios in logarithmic scale (log10) for different genes. The red and blue blocks denote higher and lower risks, respectively. The rectangles with frames mean the significant unfavorable and favorable results in prognostic analyses.
5
Isoform Details
(1) Isoform Usage: 结合violin plot和bar plot，可以找到肿瘤特异性的isoform以及在某一肿瘤类型中发生的isoform 'switch' 事件。
violin plot表示 the expression level (log2(TPM + 1)) of each isoform in a certain gene。
bar plot表示 the isoform usage (from 0% to 100%) distribution。
(2) isoform protein domain structure plots based on the prediction of Pfam.
6
Correlation Analysis
对两个基因、两个isoform或者两个signature在任意肿瘤中的相关性进行分析。
7
Similar Genes Detection
在 "TCGA Tumor", "TCGA Normal" 或者 "GTEx"样本中，搜索具有相似表达特征的基因、isoform或者signature。
8
Dimensionality Reduction
输入要研究的基因列表，选择感兴趣的"TCGA Tumor", "TCGA Normal" 或者 "GTEx"样本集以及其他参数，就会得到2D plot、3D plot以及每个主成分解释方差的比例(bar plot)。
2
Custom Data Analysis
该网站支持用户上传自己的肿瘤RNA-seq数据，以识别molecular subtype, TCGA immune subtype以及pan-cancer subtype。并且，基因和isoform的表达谱数据也可以与TCGA、GTEx数据进行比较。
接下来通过代码重复GEPIA2网站的可视化结果，以Stage plot和Correlation Analysis为例。
01
下载数据并读入R中
(1) 从UCSC xena下载表达量数据和临床信息
首先，下载UCSC Toil RNA-seq Recompute数据。
这个数据包含10535个样本，数据量比较大、耗内存，电脑配置不高的话可能容易死机；文件中的数据是log2(tmp+0.001)，使用的注释文件是gencode v23版本。
接着，下载临床信息，后面的分析需要用到“ajcc_pathologic_tumor_stage”等信息。
临床信息
(2) 下载注释文件
注释文件
# 读入表达量数据library(data.table)
exprSet=fread('./tcga_RSEM_gene_tpm.gz')
# 读入注释文件
library("rtracklayer")
gtf1=rtracklayer::import('/home/guest/annotation/gencode.v23.annotation.gtf.gz')
gtf_df=as.data.frame(gtf1)
02
整理表达矩阵
# 根据注释文件筛选基因：过滤得到protein_coding基因，并找到其ensembl id对应的gene symbol
mRNA_exprSet[1:5,1:2]
# gene_id TCGA-19-1787-01
# OR4F5 | ENSG00000186092.4 | protein_coding -9.9658
#FO538757.3 | ENSG00000279928.1 | protein_coding -4.6082
#FO538757.2 | ENSG00000279457.3 | protein_coding 4.3793
# OR4F29 | ENSG00000278566.1 | protein_coding -9.9658
# OR4F16 | ENSG00000273547.1 | protein_coding -9.9658
# 读入临床信息文件
clinical=fread('./Survival_SupplementalTable_S1_20171025_xena_sp.gz')
clinical_part=as.data.frame(clinical[,c(1,3,7)])#提取sample、cancer type abbreviation、ajcc_pathologic_tumor_stage信息
# 筛选样本：筛选出BRCA样本，并将癌旁组织剔除掉
clinical_BRCA #无癌旁组织的BRCA临床信息
# sample cancer type abbreviation ajcc_pathologic_tumor_stage group
#529 TCGA-3C-AAAU-01 BRCA Stage X tumor
#530 TCGA-3C-AALI-01 BRCA Stage IIB tumor
#531 TCGA-3C-AALJ-01 BRCA Stage IIB tumor
#532 TCGA-3C-AALK-01 BRCA Stage IA tumor
aaa=match(clinical_BRCA$sample,colnames(mRNA_exprSet))
#这里的mRNA_exprSet是所有肿瘤类型的表达矩阵，match()后得到mRNA_exprSet中BRCA样本的index
aaa=na.omit(aaa)
mRNA_exprSet=mRNA_exprSet[,aaa]#删除对应不上的样本
03
可视化
1.Stage plot（小提琴图）
# log2(tpm+0.001)逆向得到tpm
mRNA_exprSet=2^mRNA_exprSet-0.001 #tpm
mRNA_exprSet_log=log2(mRNA_exprSet+1) #log2(tpm+1)
# 将4个基因的表达量数据与临床信息整合 (使用mRNA_exprSet_log数据)
frame[1:3,1:8]
# sample PARP1 MMP9 PDCD1 CD274
#1 TCGA-3C-AAAU-01 6.815453 4.140807 0.2868778 1.3334116
#2 TCGA-3C-AALI-01 6.568070 5.979131 1.9107394 1.8678698
#3 TCGA-3C-AALJ-01 5.740120 8.044370 1.6368965 0.8560014
# cancer type abbreviation ajcc_pathologic_tumor_stage group
#1 BRCA Stage X tumor
#2 BRCA Stage IIB tumor
#3 BRCA Stage IIB tumor
# 去除tumor stage信息不明确的样本，便于后面stage plot
frame_new=frame[grep('Stage',frame$ajcc_pathologic_tumor_stage),]
# violin plot
library(vioplot)
gname=levels(factor(frame_new$ajcc_pathologic_tumor_stage))
gname
# [1] "Stage I" "Stage IA" "Stage IB" "Stage II" "Stage IIA"
# "Stage IIB" "Stage III" "Stage IIIA"
# [9] "Stage IIIB" "Stage IIIC" "Stage IV" "Stage X"
sample_id=frame_new[1]
frame_new=frame_new[,-1]
frame_new[8]=sample_id
gene=paste(c('PARP1','MMP9','PD1','PDL1'),'_gr',sep='')
for (x in 1:4){
assign(gene[x],lapply(1:12,function(i){
bb=gname[i]
y=frame_new[frame_new$ajcc_pathologic_tumor_stage== bb,]
y1=y[x]
}))
}#根据ajcc_pathologic_tumor_stage，对所有样本分组
pdf('vioplot_Result1.pdf',width = 2 + length(gname),height = 5)
for (i in c('PARP1','MMP9','PD1','PDL1')){
aov_model = aov(get(i) ~ ajcc_pathologic_tumor_stage,data = frame_new)
stat_result = summary(aov_model)[[1]][1,][4:5]#方差分析
dist = max(range(frame_new[i]))/10
ylim = range(frame_new[i])
ylim[2] = ylim[2]+dist
boxplot(0,0 , xlim=c(0.5,length(gname) + 0.5), ylim=ylim, xaxt="n", col="yellow")
vioplot(get(paste(i,'_gr',sep='')),col='blue',cex.axis = 1,font.axis =2,tcl=-0.3,add=T)
title(ylab = paste('log2(TPM+1) ',i,sep=''), cex.lab=1.5)
text(x = length(gname) + 0.5 ,y=max(frame_new[i])*1.05,
labels = paste("F value = ",signif(stat_result[1,1],3),"\nPr(>F) = ",
signif(stat_result[1,2],3),sep=""),cex = 1,col = "black",pos = 2)
axis(1, at = 1:length(gname), labels = gname, tick = TRUE)
}
dev.off()
代码运行结果：
PARP1
MMP9
PD1
PDL1
GEPIA2网站结果(与上面的顺序一致）：
PARP1
MMP9
PD1
PDL1
2.Correlation Analysis
df[1:3,]
# df记录四个基因的tpm信息，并未log化
# PARP1 MMP9 PDCD1 CD274
#TCGA-3C-AAAU-01 111.63047 16.64034 0.2199972 1.5199788
#TCGA-3C-AALI-01 93.88252 62.08087 2.7600174 2.6499326
#TCGA-3C-AALJ-01 52.45006 262.99554 2.1099611 0.8100147
library(pspearman)
method = "spearman"
# 计算使用未log化的数据，可视化使用log化的数据
pdf('scatterplot_Result.pdf',width = 6,height = 5.5)
for (j in 1:3)
{
for (i in (j+1):4) {
rvalue = signif(cor(x = df[,i], y = df[,j],method = method),2)
cpvalue = signif(as.numeric(spearman.test(x = df[,i], y = df[,j],approximation="t-distribution")[3]),2)
plot(x = log2(df[,i] + 1),y = log2(df[,j] + 1),main = NULL,cex.lab = 1.5,xlab = paste("log2(",colnames(df)[i]," TPM)",sep = ""),
ylab = paste("log2(",colnames(df)[j]," TPM)",sep = ""), pch = 19,cex=0.5)
range_y = range(log2(df[,j] + 1))
text(x = min(log2(df[,i] + 1)) ,y = max(log2(df[,j] + 1)) - (range_y[2] - range_y[1])/15,labels = paste("p-value = ",as.character(cpvalue),"\nR = ",as.character(rvalue),sep=""),cex = 1.3,col = "black",pos = 4)
}}
dev.off()
代码运行结果：
GEPIA2网站结果(与上面结果顺序一致）：
比较与总结
1.数据源：TCGA数据有多种下载方式，最开始我重复这些图的时候，使用的数据是从GDC下载的，代码运行的结果与原图有些差异；后来搜到了这个网站，发现它使用的是UCSC xena项目的数据，并且从TCGA文献中收集不同肿瘤亚型的信息。
2.比较GEPIA2网站与'自己写代码' 的可视化结果：它们的p-value、F value等还是有些差异，可能是因为我们的数据不太一样，对基因或者样本的过滤标准也不同。比较奇怪的是，PARP1 vs PD1以及PARP1 vs mmp9的散点图与网站的可视化图形比较相似，但是p-value和R值很不一样；而PARP1 vs PDL1的结果与网站结果比较一致。这个问题先留下，以后再研究研究。
3.我觉得GEPIA2网站使用起来太友好了，哪里不会点哪里，上手特别快。如果你对R语言不太熟悉，推荐使用GEPIA2网站；如果你是生信新手，想写代码得到GEPIA2的可视化结果，推荐看看生信技能树视频，以及参考GEPIA2提供的代码。
4.“使用R语言写代码”可以通过设置一些参数，向图片上添加拟合线或者其他内容；也可以实现count、FPKM与TPM之间的转换并得到相应的可视化结果。但GEPIA2提供的plot参数比较少；并且使用的是TPM值，设置参数时可以选择是否log2(TPM + 1) ，并不提供count、FPKM值的可视化结果。
参考：
哔哩哔哩【生信技能树】-- TCGA肿瘤数据库知识图谱