2020-01-07-TCGA之RTCGA

作者: 程凉皮儿 | 来源:发表于2020-01-07 22:45 被阅读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.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

        本文链接:https://www.haomeiwen.com/subject/yjqdactx.html