2020-01-07-TCGA之RTCGA

2020-01-07 本文已影响0人程凉皮儿

TCGA数据的另一个重要的包：RTCGA
参考学习资料：
TCGA之miRNA预测临床表现
 RTCGA（1）数据概况与数据下载
 RTCGA（2）数据分析与可视化
 https://www.bioconductor.org/packages/release/bioc/vignettes/RTCGA/inst/doc/RTCGA_Workflow.html
https://www.bioconductor.org/packages/release/bioc/manuals/RTCGA/man/RTCGA.pdf
TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析

RTCGA包安装

if(!require(BiocManager)){
  install.packages("BiocManager")
}
BiocManager::install("RTCGA")
library('RTCGA')
library('magrittr')

infoTCGA函数用于查看所有癌症类型：

> infoTCGA()%>%
+   rownames()%>%
+   sub('-counts','',.) -> cohort
> cohort
 [1] "ACC"      "BLCA"     "BRCA"     "CESC"     "CHOL"     "COAD"     "COADREAD"
 [8] "DLBC"     "ESCA"     "FPPP"     "GBM"      "GBMLGG"   "HNSC"     "KICH"    
[15] "KIPAN"    "KIRC"     "KIRP"     "LAML"     "LGG"      "LIHC"     "LUAD"    
[22] "LUSC"     "MESO"     "OV"       "PAAD"     "PCPG"     "PRAD"     "READ"    
[29] "SARC"     "SKCM"     "STAD"     "STES"     "TGCT"     "THCA"     "THYM"    
[36] "UCEC"     "UCS"      "UVM"

共收录了38种癌症
checkTCGA用于查询数据集的名字及数据集的公开日期：
checkTCGA查询的信息可以用于downTCGA以下载特定的TCGA数据。

checkTCGA('DataSets', 'BRCA' ) #BRCA的TCGA数据集
checkTCGA('Dates') #所有的TCGA数据的公布日期，最新的为2016-01-28

数据下载与导入

有两种方式下载：downTCGA下载数据、R包的形式下载数据，推荐以R包的形式下载数据。
第一种
下载ACC的miRNA基因表达数据

dir.create( 'hre')
downloadTCGA( cancerTypes = 'ACC', dataSet = 'miR_gene_expression',
              destDir = 'hre', date = tail( checkTCGA('Dates'), 2 )[1] )

下载的数据需要使用readTCGA的方式去读取：

dir.create('data')
# 下载临床数据
# dataset = "clinical" 是默认的参数，可以忽略
downloadTCGA( cancerTypes = c('BRCA', 'OV'),
              destDir = 'data' )
# 读取数据集
sapply( c('BRCA', 'OV'), function( element ){
  folder <- grep( paste0( '(_',element,'\\.', '|','_',element,'-FFPE)', '.*Clinical'),
                  list.files('data/'),value = TRUE)
  path <- paste0( 'data/', folder, '/', element, '.clin.merged.txt')
  assign( value = readTCGA( path, 'clinical' ),
          x = paste0(element, '.clin.data'), envir = .GlobalEnv)
})

第一种方法比较繁琐，不是大众首推的方式。
第二种方式：通过R包的形式下载数据

RTCGA将TCGA分类汇总后制成了bioconductor包来管理，现在在用的有两个版本2015-11-01与2016-01-28，建议使用最新版，数据是最新的。

RTCGA.mutations
RTCGA.rnaseq
RTCGA.clinical
RTCGA.PANCAN12
RTCGA.mRNA
RTCGA.miRNASeq
RTCGA.RPPA
RTCGA.CNV
RTCGA.methylation
在这里面，数据包RTCGA.rnaseq, RTCGA.mRNA, RTCGA.RPPA, RTCGA.miRNASeq, RTCGA.methylation的内容通过expressionsTCGA函数获取，数据包RTCGA.clinical的内容通过survivalTCGA函数获取，数据包RTCGA.mutations的内容通过mutationsTCGA函数获取。

以RTCGA.rnaseq包为例，安装之后可以通过expressionsTCGA函数获取相应的TCGA数据，其包含的数据是以“癌症缩写+rnaseq”格式的数据框，如下图：

library("BiocManager")
BiocManager::install("RTCGA.rnaseq")
library("RTCGA.rnaseq")
#获取BRCA、OV及HNSC中的VENTX基因的表达量
expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq,
                extract.cols = "VENTX|27287")

如下图，获取的结果中，第一列为病人barcode，不同的病人的不同样本有不同的barcode用于区分，第二列为dataset，代表来源于哪个数据框，本例中为BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq，第三列及以后列（如有）为自己选择的感兴趣的基因表达量。

> expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq,
+                 extract.cols = "VENTX|27287")
# A tibble: 2,085 x 3
   bcr_patient_barcode          dataset     `VENTX|27287`
   <chr>                        <chr>               <dbl>
 1 TCGA-3C-AAAU-01A-11R-A41B-07 BRCA.rnaseq          6.55
 2 TCGA-3C-AALI-01A-11R-A41B-07 BRCA.rnaseq          8.70
 3 TCGA-3C-AALJ-01A-31R-A41B-07 BRCA.rnaseq         13.6 
 4 TCGA-3C-AALK-01A-11R-A41B-07 BRCA.rnaseq          6.21
 5 TCGA-4H-AAAK-01A-12R-A41B-07 BRCA.rnaseq         34.5 
 6 TCGA-5L-AAT0-01A-12R-A41B-07 BRCA.rnaseq         26.4 
 7 TCGA-5L-AAT1-01A-12R-A41B-07 BRCA.rnaseq         39.9 
 8 TCGA-5T-A9QA-01A-11R-A41B-07 BRCA.rnaseq          6.28
 9 TCGA-A1-A0SB-01A-11R-A144-07 BRCA.rnaseq         14.4 
10 TCGA-A1-A0SD-01A-11R-A115-07 BRCA.rnaseq         12.9 
# ... with 2,075 more rows

另外，这些bioconductor数据包还可以通过convertTCGA转换成相应的bioconductor对象，如RTCGA.rnaseq,

RTCGA.RPPA, RTCGA.PANCAN12, mRNA, RTCGA.methylation可以转换成ExpressionSet对象，RTCGA.CNV可以转换成GRanges对象。

生存分析

RTCGA.clinical包中分癌症类型，有38个数据表，如BRCA.clinical, OV.clinical等等，癌症的编号详见http://gdac.broadinstitute.org/。每个数据表包括"admin.bcr"、"admin.day_of_dcc_upload"、"admin.disease_code" 等共计3703列。如下所示：

#BiocManager::install("RTCGA.clinical")
library("RTCGA.clinical")
library("tidyverse")
...
> BRCA.clinical%>%colnames()%>%head
[1] "admin.bcr"                          "admin.day_of_dcc_upload"           
[3] "admin.disease_code"                 "admin.file_uuid"                   
[5] "admin.month_of_dcc_upload"          "admin.patient_withdrawal.withdrawn"
> BRCA.clinical%>%ncol()
[1] 3703

survivalTCGA()可以从RTCGA.clinical中获取临床数据，如果不注明extract.cols参数，那么结果只有三列："times 、bcr_patient_barcode、patient.vital_status"，分别是生存时间、barcode及病人生存状态。

kmTCGA()用于绘制生存曲线，参数explanatory.names代表分组：

# 获取生存数据，不同的疾病
survivalTCGA(BRCA.clinical, OV.clinical, extract.cols = "admin.disease_code") -> BRCAOV.survInfo
# 绘图
kmTCGA(BRCAOV.survInfo, explanatory.names = "admin.disease_code", pval = TRUE)
kmTCGA(BRCAOV.survInfo, explanatory.names = "admin.disease_code", main = "",xlim = c(0,4000))

生存分析

截断x轴

# 获取生存数据，不同的治疗方式
survivalTCGA(BRCA.clinical,
             extract.cols = c("patient.drugs.drug.therapy_types.therapy_type")) %>%
  filter(patient.drugs.drug.therapy_types.therapy_type %in%
           c("chemotherapy", "hormone therapy")) %>%
  rename(therapy = patient.drugs.drug.therapy_types.therapy_type) -> BRCA.survInfo.chemo
# 绘图
kmTCGA(BRCA.survInfo.chemo, explanatory.names = "therapy", xlim = c(0, 3000), conf.int = FALSE)

image.png

表达分析

以RTCGA.rnaseq包为例，expressionsTCGA ()用于获取表达数据。

expressionsTCGA()如指定参数extract.cols，则返回特定基因在各个样本的表达量，如不指定则返回全部基因，其值的形式为“Gene symbol|Gene ID”，如 "VENTX|27287"。
PCA

library(RTCGA.rnaseq)
# 获取全部基因的表达情况
expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq) %>%
  rename(cohort = dataset) %>%
  filter(substr(bcr_patient_barcode, 14, 15) == "01") -> BRCA.OV.HNSC.rnaseq.cancer

BRCA.OV.HNSC.rnaseq.cancer %>%ncol()
    # 20533
# 由于基因数太多了，随机选1000个基因绘图
BRCA.OV.HNSC.rnaseq.cancer%>%select(sample(colnames(.),1000),-bcr_patient_barcode,cohort) %>%
  pcaTCGA("cohort")->pca.rnaseq
plot(pca.rnaseq)

分组
热图
hearmap用于绘制热图，主要有四个参数，第一个为包含所需变量的数据库，第二个和第三个分别是热图的X、Y轴分组，第四个变量为展示的变量。

# 获取MET、ZNF500、ZNF501三个基因在ACC、BLCA、BRCA、OV癌症中的表达
expressionsTCGA(ACC.rnaseq, BLCA.rnaseq, BRCA.rnaseq, OV.rnaseq,
                extract.cols = c("MET|4233", "ZNF500|26048", "ZNF501|115560")) %>%
  rename(cohort = dataset, MET = `MET|4233`) %>%  #cancer samples
  filter(substr(bcr_patient_barcode, 14, 15) == "01") %>%
  mutate(MET = cut(MET,
                   round(quantile(MET, probs = seq(0,1,0.25)), -2),
                   include.lowest = TRUE,
                   dig.lab = 5)) -> ACC_BLCA_BRCA_OV.rnaseq
# 以癌症为第一分组、MET基因表达量为第二分组，计算各分组ZNF500、ZNF501的基因表达中位数
ACC_BLCA_BRCA_OV.rnaseq %>%
  select(-bcr_patient_barcode) %>%
  group_by(cohort, MET) %>%
  summarise_all(median) %>%
  mutate(ZNF500 = round(`ZNF500|26048`),
         ZNF501 = round(`ZNF501|115560`)) -> ACC_BLCA_BRCA_OV.rnaseq.medians
# 绘制ZNF501热图
heatmapTCGA(ACC_BLCA_BRCA_OV.rnaseq.medians,
            "cohort", "MET", "ZNF501", title = "Heatmap of ZNF501 expression")

热图
箱线图
boxplotTCGA用于绘制箱线图，主要的参数有data、x、y、fill，分别是数据库、x轴、y轴，染色填充，两个有意思的参数coord.flip坐标轴翻转（默认翻转），facet.names分组（facet）变量，xlab和ylab用于定义坐标轴标题，legend.title定义图例，legend定义图例的位置。

# 获取RNAseq表达量数据
expressionsTCGA(ACC.rnaseq, BLCA.rnaseq, BRCA.rnaseq, OV.rnaseq,
                extract.cols = "MET|4233") %>%
  rename(cohort = dataset,
         MET = `MET|4233`) %>%
  #cancer samples
  filter(substr(bcr_patient_barcode, 14, 15) == "01") -> ACC_BLCA_BRCA_OV.rnaseq
h1 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "cohort", "MET") 
h2 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "cohort", "log1p(MET)") # 数据变换，可以压缩异常值的离群趋势
h3 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "reorder(cohort,log1p(MET), median)", "log1p(MET)") # 调整x的顺序
h4 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq,"cohort", "log1p(MET)",facet.names = "cohort")
library(patchwork)
h1 + h2 + h3 + h4 +plot_layout()

箱线图

RNAseq的数据演示结束

获取miRNA数据

有些miRNA并无表达，因此选择miRNA满足以下条件：其在至少十个病人中表达量均大于1。

rm(list = ls())
options(stringsAsFactors = F)
#BiocManager::install("RTCGA.miRNASeq")
library(RTCGA.miRNASeq)
rowName <- rownames(LUAD.miRNASeq)[seq(1,nrow(LUAD.miRNASeq),by=3)]
expr<-LUAD.miRNASeq[seq(1,nrow(LUAD.miRNASeq),by=3),3:ncol(LUAD.miRNASeq)]
expr<- apply(expr,2,as.numeric) 
expr<- na.omit(expr)
dim(expr)
expr <- expr[,apply(expr, 2,function(x){sum(x>1)>10})]
rownames(expr) <- rowName
dim(expr)

获取临床信息

library(RTCGA.clinical) 
meta <- LUAD.clinical
meta <- meta[,c('patient.bcr_patient_barcode','patient.vital_status',
                'patient.days_to_death','patient.days_to_last_followup',
                'patient.race',
                'patient.age_at_initial_pathologic_diagnosis',
                'patient.gender' ,
                'patient.stage_event.pathologic_stage')]

数据清洗

调整NA值、调整数据格式（数值变量、因子变量）。
生存时间拆分在两列中：死亡时间和随访时间，因此需要合并到一起。
上面获得的miRNA数据expr共有561行，而clinical数据meta共有522行，因此meta需要剔除部分不匹配的数据（下面的match函数部分）。

group_list <- ifelse(substr(rownames(expr),14,15)=='01','tumor','normal')
#table(group_list)
expr <- expr[group_list=='tumor',]

meta[,3][is.na(meta[,3])]=0
meta[,4][is.na(meta[,4])]=0
meta$days <- as.numeric(meta[,3])+as.numeric(meta[,4])
meta <- meta[,c(1:2,5:9)]
#colnames(meta)
colnames(meta)=c('ID','event','race','age','gender','stage',"days") 
meta$event=ifelse(meta$event=='alive',0,1)
meta$stage <- str_split(meta$stage,' ',simplify = T)[,2]
table(meta$stage)
meta$age <- as.numeric(meta$age)
meta$age_group <- ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')
meta$time <- meta$days/30
meta$race[is.na(meta$race)] <- "unknown"
meta$ID <- toupper(meta$ID) 
meta <- meta[match(substr(rownames(expr),1,12), meta$ID),]

meta[1:4,1:4]

挑选miRNA，构建cox回归模型

由于生存分析绘图常用survminer，因此使用survminer绘制cox图时有一个坑：绘图前需要手动指定好cox回归的模型。
cox回归还可以展示森林图。具体如下，虽然总p值为0.018，但是8个miRNA均没有达到显著水平（置信区间横跨无效线）。

library(survival)
library(survminer)

e <- expr[,c('hsa-mir-31','hsa-mir-196b','hsa-mir-766','hsa-mir-519a-1',
                 'hsa-mir-375','hsa-mir-187','hsa-mir-331','hsa-mir-101-1')]
e <- log2(e+1)
colnames(e) <- c('miR31','miR196b','miR766','miR519a1','miR375','miR187','miR331','miR101')
dat <- cbind(meta,e)

dat$gender <- factor(dat$gender)
dat$stage <- factor(dat$stage)

colnames(dat) 
s <- Surv(time, event) ~ miR31+miR196b+miR766+miR519a1+miR375+miR187+miR331+miR101
model <- coxph(s, data = dat )
summary(model,data=dat)
options(scipen=1)

fit <- survfit(model, data=dat)
fit$call$formula <- s #手动指定模型
ggsurvplot(fit, palette = "#2E9FDF", ggtheme = theme_minimal())

ggforest(model, data =dat, 
         main = "Hazard ratio", #标题
         cpositions = c(0.10, 0.22, 0.4), #前三列距离
         fontsize = 1.0, #字体大小
         refLabel = "1",#相对变量的数值标签
         noDigits = 4 #HR值及置信区间的有效数字位数
         )

生存曲线

森林图，目前还不是太明白，需要再学习

绘制风险因子关联图
跟森林图一样都是不熟悉的东西，需要进一步学习

new_dat=dat

fp <- predict(model,new_dat,type="risk")
fp <- predict(model,new_dat,type="expected") 
plot(fp,meta$days)
fp <- predict(model,new_dat) 
basehaz(model) 
library(Hmisc)
options(scipen=200)
with(new_dat,rcorr.cens(fp,Surv(time, event)))

library(cowplot)
library(pheatmap)
# https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
fp_dat=data.frame(s=1:length(fp),v=as.numeric(sort(fp )))
sur_dat=data.frame(s=1:length(fp),
                   t=meta[names(sort(fp )),'time'] ,
                   e=meta[names(sort(fp )),'event']  ) 
sur_dat$e=ifelse(sur_dat$e==0,'alive','death')
exp_dat=new_dat[names(sort(fp )),10:17]
plot.point=ggplot(fp_dat,aes(x=s,y=v))+geom_point()
plot.sur=ggplot(sur_dat,aes(x=s,y=t))+geom_point(aes(col=e))
mycolors <- colorRampPalette(c("blue", "white", "red"), bias = 1.2)(100)
tmp=t(scale(exp_dat))
tmp[tmp > 1] = 1
tmp[tmp < -1] = -1
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = F)
plot_grid(plot.point, plot.sur, plot.h$gtable,
          labels = c("A", "B","C"),
          align = 'v',ncol = 1)

虽然不知道怎么解释，但是觉得应该能看懂

2020-01-07-TCGA之RTCGA

RTCGA包安装

数据下载与导入

生存分析

表达分析

获取miRNA数据

猜你喜欢

热点阅读