使用R包RTCGA把TCGA数据本地化

生信技能树

发布于 2022-07-26 10:06:35

5880

发布于 2022-07-26 10:06:35

文章被收录于专栏：生信技能树

前面我们介绍了MSKCC和Broad研究所的网页工具可以帮助我们探索TCGA数据库的多个癌症的多组学数据，见：

虽然他们都有超级棒的网页工具，但是我们是生信工程师，还是习惯了自己写代码来批量操作，所以就有配套的R包：cgdsr 和 RTCGAToolbox，但是它们有一个弊端是需要实时联网在线根据自己的需求去下载数据，对网络环境不好的小伙伴来说是一个考验。

所以我们也推送了两个离线解决方案吧，首先是使用R包RTCGA把TCGA数据本地化。最新的数据版本是2016-01-28的，本来是可以选择以下的包（但是这些包都被移除了，并不在bioconductor官网啦）：

RTCGA.mutations.20160128
RTCGA.rnaseq.20160128
RTCGA.clinical.20160128
RTCGA.mRNA.20160128
RTCGA.miRNASeq.20160128
RTCGA.RPPA.20160128
RTCGA.CNV.20160128
RTCGA.methylation.20160128

本来是说旧版本已经可以考虑弃用了，但是因为上面的新版数据包都不存在了，我们反而是只能说退而求其次选择旧版数据。不过TCGA数据其实是2005就开始了，哪怕是截止到 2015-11-01 版本也不会过时很多。

下面是基于 2015-11-01 版本的 TCGA 数据，如果你使用默认安装方法，通常是安装的旧版数据：

RTCGA.mutations
RTCGA.rnaseq # 测序表达量矩阵
RTCGA.clinical
RTCGA.PANCAN12
RTCGA.mRNA # 芯片表达量矩阵
RTCGA.miRNASeq
RTCGA.RPPA
RTCGA.CNV
RTCGA.methylation

每一个包都是正规的bioconductor包，使用标准安装代码即可：

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("RTCGA")
BiocManager::install("RTCGA.clinical")
BiocManager::install("RTCGA.mRNA")

这个包不仅仅是提供数据整理和下载，也有几个好用的打包的函数，包括：

expressionsTCGA Gather Expressions for TCGA Datasets
mutationsTCGA Gather Mutations for TCGA Datasets

heatmapTCGA Create Heatmaps for TCGA Datasets 
kmTCGA Plot Kaplan-Meier Estimates of Survival Curves for Survival Data
pcaTCGA Plot Two Main Components of Principal Component Analysis

其中expressionsTCGA和mutationsTCGA是从本地R包里面获取表达量信息，另外的3个函数（heatmapTCGA，kmTCGA，pcaTCGA）出图都还凑合，后面我们会一一介绍。

支持的数据集

最主要的包RTCGA有一个函数infoTCGA可以查询他们团队整理的数据集，其实就是 https://gdac.broadinstitute.org/ 网页里面的数据集：

library(RTCGA) 
all_TCGA_cancers=infoTCGA()
DT::datatable(all_TCGA_cancers)

如下所示：

支持的数据集

既然也是 https://gdac.broadinstitute.org/，其实就跟通过R包RTCGAToolbox链接FireBrowse来探索TCGA等公共数据 类似了，因为都是broad研究所的整理。

首先需要使用expressionsTCGA和mutationsTCGA从本地R包里面获取表达量信息，然后使用前面提到了有3个函数（heatmapTCGA，kmTCGA，pcaTCGA）可以进行统计可视化，我们就来逐一介绍。

expressionsTCGA

首先我们需要对任意基因从任意癌症里面获取芯片表达数据，这里我们拿下面3种癌症做示范：

Breast invasive carcinoma (BRCA)
Ovarian serous cystadenocarcinoma (OV)
Lung squamous cell carcinoma (LUSC)

我们拿一些免疫基因作为例子，其中要注意的是mRNA并不是rnaseq两者不太一样，具体样本数量，可以看最前面的表格。

library(RTCGA)
library(RTCGA.mRNA)
cg=c('CD3D','CD3G','CD247','IFNG','IL2RG','IRF1','IRF4','LCK','OAS2','STAT1')
cg
expr <- expressionsTCGA(BRCA.mRNA, OV.mRNA, LUSC.mRNA,
                        extract.cols =  cg ) 
expr
nb_samples <- table(expr$dataset)
nb_samples
library(ggpubr)
ggboxplot(expr, x = "dataset", y = "CD3D",
          title = "CD3D", ylab = "Expression",
          color = "dataset", palette = "jco",
          rotate = TRUE)

可以看到在这3个癌症里面CD3D基因的表达量范围是不一样的：

3个癌症里面CD3D基因的表达量范围是不一样

不过，这个数据是mRNA并不是rnaseq，而且看起来是被zscore了的，这样的值在不同数据集里面的对比起来是有问题的。

大家可以安装另外的包，去看看转录组测序数据里面的这些基因的表达量。而且测序数据里面的基因还比较麻烦，不仅仅是要symbol还要entrez的ID，具体需要看 https://wiki.nci.nih.gov/display/TCGA/RNASeq+Version+2 的解释。

mutationsTCGA

也是同样的安装：

> BiocManager::install('RTCGA.mutations')
# RTCGA.mutations_20151101.22.0.tar.gz'
# Content type 'application/x-gzip' length 108858817 bytes (103.8 MB)

可以看到仍然是 2015-11-01的数据，文件比较大，所以仍然是考验网络速度。

library(RTCGA.mutations)  
library(dplyr)
mutationsTCGA(BRCA.mutations, OV.mutations) %>%
  filter(Hugo_Symbol == 'TP53') %>%
  filter(substr(bcr_patient_barcode, 14, 15) ==
           "01") %>% # cancer tissue
  mutate(bcr_patient_barcode =
           substr(bcr_patient_barcode, 1, 12)) ->
  BRCA_OV.mutations
table(BRCA_OV.mutations$dataset)

可以看到，我们很轻松的就从OV和BRCA两个癌症里面，拿到了有TP53突变的病人的ID ，可以看到：

> table(BRCA_OV.mutations$dataset)

BRCA.mutations   OV.mutations 
           306            279 
> as.data.frame(tail(sort(table(BRCA_OV.mutations$Variant_Classification))))
               Var1 Freq
1       Splice_Site   20
2   Frame_Shift_Ins   23
3   Splice_Site_SNP   32
4   Frame_Shift_Del   64
5 Nonsense_Mutation   70
6 Missense_Mutation  352

两个癌症数据集里面的都有几百个病人有TP53基因的突变，而且绝大部分突变类型都是Nonsense_Mutation和Missense_Mutation。

这样的突变信息，就可以对病人进行分组后，看生存差异，我在生信技能树多次分享过生存分析的细节；

heatmapTCGA

我感觉完全没有必要使用它的函数，因为热图基本上大家都是可以自定义，再不济至少可以pheatmap一下。比如前面的3个癌症数据集的表达量矩阵。简单的代码如下所示：

ar=data.frame(group=expr$dataset)
rownames(ar)=expr$bcr_patient_barcode
head(ar)
ct=expr[,3:ncol(expr)]
rownames(ct)=expr$bcr_patient_barcode
pheatmap::pheatmap(ct,
                   show_rownames = F,
                   annotation_row = ar)

可以看到，虽然前面3个癌症里面CD3D基因的表达量范围有差异，但是这些免疫集联合起来就没有那么强烈的癌症特异性了，说明免疫这个变量在每个癌症内部都是很具有异质性的，所以不同癌症很难根据免疫进行区分。

免疫这个变量在每个癌症内部都是很具有异质性

pcaTCGA

前面我们介绍了仅仅是根据指定的基因列表，就可以筛选表达量矩阵，并且进行合理的分组，见：

好像也没有必要使用作者的pcaTCGA函数。

kmTCGA

我们前面筛选到了OV和BRCA两个癌症的TP53突变信息，现在就可以结合起来临床信息，使用它的kmTCGA函数，进行生存分析：

library(RTCGA.clinical)
survivalTCGA(
  BRCA.clinical,
  OV.clinical,
  extract.cols = "admin.disease_code"
) %>%
  dplyr::rename(disease = admin.disease_code) ->
  BRCA_OV.clinical

BRCA_OV.clinical %>%
  left_join(
    BRCA_OV.mutations,
    by = "bcr_patient_barcode"
  ) %>%
  mutate(TP53 =
           ifelse(!is.na(Variant_Classification), "Mut","WILDorNOINFO")) ->
  BRCA_OV.clinical_mutations

BRCA_OV.clinical_mutations %>%
  select(times, patient.vital_status, disease, TP53) -> BRCA_OV.2plot

kmTCGA(
  BRCA_OV.2plot,
  explanatory.names = c("TP53", "disease"),
  break.time.by = 400,
  xlim = c(0,2000),
  pval = TRUE) -> km_plot

km_plot

如下所示：