通过R包cgdsr链接cbioportal来探索TCGA等公共数据

通过R包cgdsr链接cbioportal来探索TCGA等公共数据 by 生信技能树

众所周知，TCGA数据库是目前最综合全面的癌症病人相关组学数据库，包括的多组学数据有：

DNA Sequencing (WGS/WES)
mRNA/miRNA Sequencing
Protein Expression Array
Array-based Expression
DNA Methylation Array
Copy Number Array

知名的肿瘤研究机构都有着自己的TCGA数据库探索工具，比如：

Broad Institute FireBrowse portal, The Broad Institute
cBioPortal for Cancer Genomics, Memorial Sloan-Kettering Cancer Center
TCGA Batch Effects, MD Anderson Cancer Center
Regulome Explorer, Institute for Systems Biology
Next-Generation Clustered Heat Maps, MD Anderson Cancer Center

其中MSKCC的就是cBioPortal，主页在：https://www.cbioportal.org/

这个网页工具可以很方便的查询感兴趣的目标基因在指定的癌症数据集里面，是否有突变，有表达量差异或者生存意义。包含的数据集列表在：https://www.cbioportal.org/datasets，值得注意的是这里面可不仅仅是TCGA数据集哦，也并不是说每个数据集都有多组学哦，如果我们按照样品数量排序，很容易看到：

cBioPortal网站文献列表

只要是ngs数据文献愿意遵守这个cBioPortal规则，就可以把其突变信息和表达量信息以及病人临床信息，按照固定的格式整理好，并且上传到cBioPortal网站，就可以使用cBioPortal的网页工具去探索它。

我们这里不纠结这个网页工具的用法了，那个很容易摸索，我们介绍生信工程师喜欢的编程语言来操作这个网页工具，就是 cgdsr: R-Based API for Accessing the MSKCC Cancer Genomics Data Server (CGDS)

认识cgdsr

既然是一个R包，首先就需要安装它，超级简单的代码：

install.packages('cgdsr')

接下来查看有多少个文献，多少个样品列表信息，每个文献都是有多少个数据集，分别对应的函数是：

getCancerStudies ，获取有哪些文献
getGeneticProfiles ，查询指定文献里面的数据集有哪些
getCaseLists，查询指定文献里面的样品列表信息有哪些

首先是文献信息，需要使用getCancerStudies 来获取有哪些文献：

library(cgdsr)
library(DT) 
# Get list of cancer studies at server ## 获取有哪些文献
mycgds <- CGDS("https://www.cbioportal.org/")
all_TCGA_studies <- getCancerStudies(mycgds)
all_TCGA_studies[1:3, 1:2]
#write.csv(all_TCGA_studies,paste0(Sys.time(),"all_TCGA_studies.csv"),row.names = F) 
DT::datatable(all_TCGA_studies)

这个代码其实就是实现了：https://www.cbioportal.org/datasets，网页里面的文献列表获取，可以看到也是340个文献。

因为每个文献都是突变数据和表达量数据，所以可以使用getGeneticProfiles ，查询指定文献里面的数据集有哪些。

cg_study <-  all_TCGA_studies[6,1]
cg_study # [1] "laml_tcga"
all_dataset <- getGeneticProfiles(mycgds, cg_study)
dim(all_dataset) # 11个数据集
cg_dataset <- all_dataset[9,1]
cg_dataset #[1] "laml_tcga_rna_seq_v2_mrna"

这个laml_tcga对应的文献提供了11种数据集，我们这里以 laml_tcga_rna_seq_v2_mrna 举例，它是 mRNA gene expression (RNA Seq V2 RSEM) 的简称。

使用 getCaseLists，查询指定文献里面的样品列表信息有哪些，代码如下所示：

all_tables <- getCaseLists(mycgds, cg_study)
all_tables[1:3, 1:2]
dim(all_tables) ## 共10种样本列表方式
# 这样的样品列表可以去获取其它属性信息
cg_table <-  all_tables[1,1]
cg_table # [1] "laml_tcga_all"

可以看到，虽然是同一个文献，前面的查询发现有11个数据集，但是只有10个样品列表，如下所示：

10个样品列表

这些样品列表可以用来获取临床信息，一般来说是使用全部的样品列表即可。

clinicaldata <- getClinicalData(mycgds, cg_table)

这些样品列表也可以用来获取组学信息，包括突变数据和表达量矩阵。

df <- getProfileData(mycgds, c('CD4','CD8A'), cg_dataset, cg_table)
dim(df) # 200个样品
df = na.omit(df)
dim(df) # 173个样品
cor(df[,1],df[,2])

因为我们使用全部的样品列表去获取表达量矩阵，而这个laml_tcga文献里面是200个病人，但并不是所有的病人都有转录组测序，其实就173个有，所以需要去除NA值。当然了，最好的方法是，获取表达量矩阵的时候，就使用对应的样品列表，而不是选择全部的病人样品列表。

这个 getProfileData 函数是万能的，后续如果我们要数量使用它来举例，基本上就是靠这一个函数即可。

玩转cgdsr

前面我们提到了通过R包cgdsr链接cbioportal来探索TCGA等公共数据，而我喜欢把TCGA数据库的应用划分为8个领域：

1、探索各类肿瘤不同临床特征（性别、年龄、种族、临床分期）的预后（生存曲线）
2、探索各类肿瘤与对照的单个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）水平的差异情况（箱线图）
3、探索各类肿瘤与对照的全局（mRNA，lncRNA，miRNA，甲基化，蛋白）水平的差异情况（差异分析流程）
4、探索各类肿瘤中两个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）水平相关性（散点图）
5、探索各类肿瘤中多个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）水平总结（热图）
6、探索各类肿瘤中单个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）与所有其它分子相关性并且排序
7、探索各类肿瘤中单个基因突变或者单个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）水平的预后（生存曲线）
8、探索各类肿瘤不同临床特征（性别、年龄、种族、临床分期）分组后的单个分子（mRNA，lncRNA，miRNA，甲基化，蛋白）特性的分布

下面我们就演示一下如何实现上面的需求：

性别这个临床特征是否有临床预后意义

> phe=clinicaldata[,c('OS_MONTHS','OS_STATUS','SEX')]
> head(phe)
                OS_MONTHS  OS_STATUS    SEX
TCGA.AB.2882.03     11.99 1:DECEASED Female
TCGA.AB.2991.03     59.99   0:LIVING Female
TCGA.AB.2847.03     19.97 1:DECEASED   Male
TCGA.AB.2979.03     22.04   0:LIVING Female
TCGA.AB.2818.03      9.95 1:DECEASED Female
TCGA.AB.2927.03      2.96 1:DECEASED Female

有了 OS_MONTHS OS_STATUS SEX，做一个简单的生存分析就太简单了。我在生信技能树多次分享过生存分析的细节；

单个基因的表达量和拷贝数变异的关系探索

同样的使用getProfileData函数，去指定的文献里面，根据指定的样品列表，去获取指定的数据集信息

df = getProfileData(mycgds, "NF1", c("gbm_tcga_gistic","gbm_tcga_mrna"), "gbm_tcga_all")
head(df)
df=na.omit(df)
table(df$gbm_tcga_gistic)
boxplot(df[,2] ~ df[,1], main="NF1 : CNA status vs mRNA expression", 
        xlab="CNA status", ylab="mRNA expression", outpch = NA)
stripchart(df[,2] ~ df[,1], vertical=T, add=T, method="jitter",pch=1,col='red')

如下所示：

单个基因的表达量和拷贝数变异的关系探索

两个基因的表达量相关性散点图

这个其实前面演示了：

df = getProfileData(mycgds, c("MDM2","MDM4"), "gbm_tcga_mrna", "gbm_tcga_all")
head(df)
df=na.omit(df)
plot(df, main="MDM2 and MDM4 mRNA expression", 
     xlab="MDM2 mRNA expression", ylab="MDM4 mRNA expression")

两个不同数据集的同一个基因表达量差异箱线图

df.pri = getProfileData(mycgds, "PTEN", "prad_mskcc_mrna_median_Zscores", "prad_mskcc_primary")
head(df.pri)
df.met = getProfileData(mycgds, "PTEN", "prad_mskcc_mrna_median_Zscores", "prad_mskcc_mets")
head(df.met)
boxplot(list(t(df.pri),t(df.met)), main="PTEN expression in primary and metastatic tumors", xlab="Tumor type", ylab="PTEN mRNA expression",names=c('primary','metastatic'), outpch = NA)
stripchart(list(t(df.pri),t(df.met)), vertical=T, add=T, method="jitter",pch=1,col='red')

如下所示：

两个不同数据集的同一个基因表达量差异箱线图

ixxmu / mp_duty