R語言使用cgdsr包獲取TCGA資料範例詳解

2022-06-25 14:01:36

TCGA資料來源
TCGA資料庫探索工具
檢視任意資料集的樣本列表方式
選定資料形式及樣本列表後獲取感興趣基因的資訊,下載mRNA資料
選定樣本列表獲取臨床資訊
綜合性獲取
下載mRNA資料
獲取病例列表的臨床資料
從cBioPortal下載點突變資訊
從cBioPortal下載拷貝數變異資料
把拷貝數及點突變資訊結合畫熱圖

TCGA資料來源

眾所周知,TCGA資料庫是目前最綜合全面的癌症病人相關組學資料庫，包括的測序資料有：

DNA Sequencing

miRNA Sequencing

Protein Expression

mRNA Sequencing

Total RNA Sequencing

Array-based Expression

DNA Methylation

Copy Number

TCGA資料庫探索工具

知名的腫瘤研究機構都有著自己的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

其中cBioPortal更是被包裝到R包裡面

這裡介紹如何使用R語言的cgdsr包來獲取任意TCGA資料。

cgdsr包:R語言工具包，可以下載TCGA資料。
DT包:data.table包，簡稱DT包，是R語言中的資料視覺化工具包。DT包可以將Javascript中的方法運用到R中，也能將矩陣或者資料表在網頁中視覺化為表格，以及其它的一些功能。

> setwd("C:/Users/YLAB/Documents/R/win-library/4.1/")
> install.packages("R.methodsS3_1.8.1.zip",repos=NULL)#安裝
> install.packages("R.oo_1.24.0.zip",repos=NULL)#安裝
> install.packages("data.table")
> BiocManager::install("cgdsr", force = TRUE)#安裝
> library(cgdsr)
> library(DT)
#建立一個cgdsr物件
> mycgds <- CGDS("http://www.cbioportal.org/")
#檢查下載是否成功，如果是FAILED就是沒成功。
> test(mycgds)
getCancerStudies...  OK
getCaseLists (1/2) ...  OK
getCaseLists (2/2) ...  OK
getGeneticProfiles (1/2) ...  OK
getGeneticProfiles (2/2) ...  OK
getClinicalData (1/1) ...  OK
getProfileData (1/6) ...  OK
getProfileData (2/6) ...  OK
getProfileData (3/6) ...  OK
getProfileData (4/6) ...  OK
getProfileData (5/6) ...  OK
getProfileData (6/6) ...  OK
all_TCGA_studies <- getCancerStudies(mycgds)
> DT::datatable(all_TCGA_studies)

檢視任意資料集的樣本列表方式

上表的cancer_study_id其實就是資料集的名字，我們任意選擇一個資料集，比如stad_tcga_pub ，可以檢視它裡面有多少種樣本列表方式。

stad2014 <- "stad_tcga_pub"
## 獲取在stad2014資料集中有哪些表格（每個表格都是一個樣本列表）
all_tables <- getCaseLists(mycgds, stad2014)
dim(all_tables) ## 共6種樣本列表方式
[1] 6 5
DT::datatable(all_tables[,1:3])

檢視任意資料集的資料形式

## 而後獲取可以下載哪幾種資料，一般是mutation，CNV和表達量資料
all_dataset <- getGeneticProfiles(mycgds, stad2014)
DT::datatable(all_dataset,
                  extensions = 'FixedColumns',
                  options = list(                    #dom = 't',
                    scrollX = TRUE,
                    fixedColumns = TRUE
                  ))

一般來說，TCGA的一個專案資料就幾種，如下：

選定資料形式及樣本列表後獲取感興趣基因的資訊,下載mRNA資料

my_dataset <- 'stad_tcga_pub_rna_seq_v2_mrna'
my_table <- "stad_tcga_pub_rna_seq_v2_mrna" 
BRCA1 <- getProfileData(mycgds, "BRCA1", my_dataset, my_table)
dim(BRCA1)
[1] 265   1

樣本個數差異很大，不同癌症熱度不一樣。

選定樣本列表獲取臨床資訊

## 如果我們需要繪製survival curve，那麼需要獲取clinical資料
clinicaldata <- getClinicalData(mycgds, my_table)
DT::datatable(clinicaldata,
                  extensions = 'FixedColumns',
                  options = list(                    #dom = 't',
                    scrollX = TRUE,
                    fixedColumns = TRUE
                  ))

綜合性獲取

只需要根據癌症列表選擇自己感興趣的研究資料集即可，然後選擇好感興趣的資料形式及對應的樣本量。就可以獲取對應的資訊：

library(cgdsr)
library(DT)
 mycgds <- CGDS("http://www.cbioportal.org")   
##mycancerstudy = getCancerStudies(mycgds)[25,1]
mycancerstudy = 'brca_tcga'   getCaseLists(mycgds,mycancerstudy)[,1]
##  [1] "brca_tcga_3way_complete"          "brca_tcga_all"                   
##  [3] "brca_tcga_protein_quantification" "brca_tcga_sequenced"             
##  [5] "brca_tcga_cna"                    "brca_tcga_methylation_hm27"      
##  [7] "brca_tcga_methylation_hm450"      "brca_tcga_mrna"                  
##  [9] "brca_tcga_rna_seq_v2_mrna"        "brca_tcga_rppa"                  
## [11] "brca_tcga_cnaseq"

getGeneticProfiles(mycgds,mycancerstudy)[,1]
##  [1] "brca_tcga_rppa"                          
##  [2] "brca_tcga_rppa_Zscores"                  
##  [3] "brca_tcga_protein_quantification"        
##  [4] "brca_tcga_protein_quantification_zscores"
##  [5] "brca_tcga_gistic"                        
##  [6] "brca_tcga_mrna"                          
##  [7] "brca_tcga_mrna_median_Zscores"           
##  [8] "brca_tcga_rna_seq_v2_mrna"               
##  [9] "brca_tcga_rna_seq_v2_mrna_median_Zscores"
## [10] "brca_tcga_linear_CNA"                    
## [11] "brca_tcga_methylation_hm450"             
## [12] "brca_tcga_mutations"

下載mRNA資料

mycaselist ='brca_tcga_rna_seq_v2_mrna' 
 mygeneticprofile = 'brca_tcga_rna_seq_v2_mrna'  
# Get data slices for a specified list of genes, genetic profile and case list
expr=getProfileData(mycgds,c('BRCA1','BRCA2'),mygeneticprofile,mycaselist)
DT::datatable(expr)

很簡單就得到了指定基因在指定癌症的表達量

獲取病例列表的臨床資料

myclinicaldata = getClinicalData(mycgds,mycaselist)
DT::datatable(myclinicaldata,
                  extensions = 'FixedColumns',
                  options = list(                    #dom = 't',
                    scrollX = TRUE,
                    fixedColumns = TRUE
                  ))
## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## http://rstudio.github.io/DT/server.html

從cBioPortal下載點突變資訊

#突變基因名稱集合
mutGene=c("EGFR", "PTEN", "TP53", "ATRX")
#檢索基因和遺傳圖譜的基因組圖譜資料
mut_df <- getProfileData(mycgds, 
                         caseList ="gbm_tcga_sequenced", 
                         geneticProfile = "gbm_tcga_mutations",
                         genes = mutGene
)
mut_df <- apply(mut_df,2,as.factor)
mut_df[mut_df == "NaN"] = ""
mut_df[is.na(mut_df)] = ""
mut_df[mut_df != ''] = "MUT" 
DT::datatable(mut_df)

從cBioPortal下載拷貝數變異資料

mutGene=c("TP53","UGT2B7","CYP3A4")
cna<-getProfileData(mycgds,mutGene,"gbm_tcga_gistic","gbm_tcga_sequenced")
cna<-apply(cna,2,function(x) as.character(factor(x,levels = c(-2:2),labels = c("HOMDEL","HETLOSS","DIPLOID","GAIN","AMP"))))
cna[is.na(cna)]=""
cna[cna=="DIPLOID"]=""
DT::datatable(cna)

把拷貝數及點突變資訊結合畫熱圖

下面的函數，主要是配色比較複雜，其實原理很簡單，就是一個熱圖。

library(ComplexHeatmap)
library(grid)
conb <- data.frame(matrix(paste(as.matrix(cna),as.matrix(mut_df),sep = ";"),                          nrow=nrow(cna),ncol=ncol(cna),                          dimnames=list(row.names(mut_df),colnames(cna))))
mat <- as.matrix(t(conb))
DT::datatable((mat)) 
alt <- apply(mat,1,function(x)strsplit(x,";"))
alt <- unique(unlist(alt))
alt <- alt[which(alt !="")]
alt <-c("background",alt)
alter_fun = list(  background = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="#CCCCCC",col=NA))  },  HOMDEL = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="blue3",col=NA))  },  HETLOSS = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="cadetblue1",col=NA))  },  GAIN = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="pink",col=NA))  },  AMP = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="red",col=NA))  },  MUT = function(x,y,w,h){    grid.rect(x,y,w-unit(0.5,"mm"),h-unit(0.5,"mm"),              gp=gpar(fill="#008000",col=NA))  })
col <- c("MUT"="#008000","AMP"="red","HOMDEL"="blue3",         "HETLOSS"="cadetblue1","GAIN"="pink")
alt = intersect(names(alter_fun),alt)
alt_fun_list <- alter_fun[alt]
col <- col[alt]
oncoPrint(mat=mat,alter_fun = alt_fun_list,          get_type = function(x) strsplit(x,";")[[1]],          col = col)