结合TCGAbiolinks批量下载TCGA临床信息

TCGAbiolinks包常用来TCGA的下载和各种常用分析，本人在使用TCGAbiolinks下载和整理TCGA临床信息时碰到一些麻烦，具体如下：

TCGAbiolinks下载临床信息代码常规如下：

library(TCGAbiolinks)

projects <- TCGAbiolinks::getGDCprojects()$project_id
projects <- projects[grepl('^TCGA', projects, perl=TRUE)]

meta=list()
for(project in projects){  
  query <- GDCquery(project = project,
                    data.category = "Clinical", 
                    file.type = "xml")
  GDCdownload(query)
  clinical <- GDCprepare_clinic(query, clinical.info = "patient")
  meta[[project]]=clinical
}

meta=as.data.frame(rbind,meta)

通过以上代码就把TCGA的33种肿瘤的临床信息就下载完了。GDCprepare_clinic函数中关于临床信息的选项不多，就drug, admin, follow_up,radiation, patient, stage_event这么几个，而我们一般需要的是bcr_patient_barcode, vital_status, days_to_death, days_to_last_followup, days_to_birth, gender, stage_event这些信息。因此，还需要对这个代码进行一些修改，参考了生信技能树的代码。

另外一点就是TCGA-GBM, TCGA-LAML, TCGA-LGG, TCGA-PCPG, TCGA-SARC这几个 肿瘤是没有stage_event信息的，所以TCGA的33种肿瘤还不能一次性处理。

library(TCGAbiolinks)
library(XML)
library(dplyr)
library(stringr)
library(tibble)
projects <- getGDCprojects()$project_id
projects <- projects[grepl('^TCGA', projects, perl=TRUE)]
projects <- projects[order(projects)]

projects1=projects[-which(projects %in% c("TCGA-GBM", "TCGA-LAML", "TCGA-LGG", "TCGA-PCPG", "TCGA-SARC"))]

# No stage information was found in project2.
projects2=c("TCGA-GBM", "TCGA-LAML", "TCGA-LGG", "TCGA-PCPG", "TCGA-SARC")

clinical1=list()
for (project in projects1) {
  query <- GDCquery(project = project,
                    data.category = "Clinical",
                    file.type = "xml")
  GDCdownload(query)
  dir1="/harmonized/Clinical/Clinical_Supplement/"
  xmls=dir(paste0('GDCdata/',project,dir1),pattern = "*.xml$",recursive = T)
  
  cl=list()
  
  use_col = c('bcr_patient_barcode',
              'vital_status',
              'days_to_death',
              'days_to_last_followup',
              'days_to_birth',
              'gender' ,
              'stage_event')
  
  for(i in 1:length(xmls)){
    result = xmlParse(paste0('GDCdata/',project,dir1,xmls[[i]]))
    rootnode = xmlRoot(result)
    dat=xmlToDataFrame(rootnode[2])
    if(sum(use_col %in% colnames(dat))==7)
      cl[[i]]=dat[,use_col]
  }
  
  clinical1[[project]] = do.call(rbind,cl)
}

clinical2=list()
for (project in projects2) {
  query <- GDCquery(project = project,
                    data.category = "Clinical",
                    file.type = "xml")
  GDCdownload(query)
  dir1="/harmonized/Clinical/Clinical_Supplement/"
  xmls=dir(paste0('GDCdata/',project,dir1),pattern = "*.xml$",recursive = T)
  
  cl=list()
  
  use_col= c('bcr_patient_barcode',
              'vital_status',
              'days_to_death',
              'days_to_last_followup',
              'days_to_birth',
              'gender')
  
  for(i in 1:length(xmls)){
    result = xmlParse(paste0('GDCdata/',project,dir1,xmls[[i]]))
    rootnode = xmlRoot(result)
    dat=xmlToDataFrame(rootnode[2])
    if(sum(use_col %in% colnames(dat))==6)
      cl[[i]]=dat[,use_col]
  }
  
  clinical2[[project]] = do.call(rbind,cl)
}

clinical1=as.data.frame(do.call(rbind,clinical1))
clinical2=as.data.frame(do.call(rbind,clinical2))
clinical2$stage_event=NA
identical(colnames(clinical1),colnames(clinical2))
meta=rbind(clinical1,clinical2)%>%
  tibble::add_column(Cancer = stringr::str_split(rownames(.),'-',simplify = T)[,2], .before="bcr_patient_barcode")
meta$Cancer=stringr::str_split(meta$Cancer,'[.]',simplify = T)[,1]

meta=meta[!duplicated(meta$bcr_patient_barcode),]
meta=meta[order(meta$Cancer,meta$bcr_patient_barcode),]
rownames(meta) <- meta$bcr_patient_barcode

# rename colname
colnames(meta)=c('Cancer','ID','event','death','last_followup','age','gender','stage')
# define missing value as NA
meta[meta==""]=NA

查看一下现在的临床信息：

图片

接下来就是计算生存期，整理分期信息了。

# calculate survival time
table(meta$event)
# Alive  Dead
# 8405  2748
meta$time = ifelse(meta$event=="Alive",
                   meta$last_followup,
                   meta$death)
meta$time = round(as.numeric(meta$time)/30,2)
meta$age = round(-(as.numeric(meta$age)/365),0)

meta$gender=ifelse(meta$gender=="MALE","M","F")

# define event，Alive=0，Dead=1
meta$event=ifelse(meta$event=='Alive',
                  0,
                  1)
table(meta$event)
# 0     1
# 8405  2748

k1 = meta$time>=0.1;table(k1)
# FALSE  TRUE
# 598 10477
k2 = !(is.na(meta$time)|is.na(meta$event));table(k2)
# FALSE  TRUE
# 92 11075
meta = meta[k1&k2,]
nrow(meta)
# [1] 10477

# stage
meta$stage=str_split(meta$stage,"Stage ",simplify = T)[,2]
a = str_extract_all(str_sub(meta$stage,1,4),"I|V");head(a)
b = sapply(a,paste,collapse = "")
table(b)
#      I   II  III   IV
# 2562 2487 2206 2259  963

meta$stage = b
# set missing value as NA
table(meta$stage)
#      I   II  III   IV
# 2562 2487 2206 2259  963

meta[meta=="" | meta=="not reported"]=as.character(NA)
table(meta$stage,useNA = "always")
# I    II   III   IV  <NA>
# 2487 2206 2259  963 2562

meta=meta[,c('Cancer','event','time','age','gender','stage')]

整理好的临床信息如下：

图片

看看每个肿瘤病例数：

table(meta$Cancer)

 # ACC BLCA BRCA CESC CHOL COAD DLBC ESCA  GBM HNSC KICH KIRC KIRP LAML  LGG LIHC LUAD LUSC MESO   OV PAAD PCPG PRAD READ 
 # 91  393 1025  283   47  340   44  150  592  518  109  521  260  173  504  341  478  470   85  581  179  178  498  127 
# SARC SKCM STAD TGCT THCA THYM UCEC  UCS  UVM 
# 251  435  384  132  498  122  536   54   78

最后就是和表达矩阵进行匹配了，log2(tpm+1)形式的表达矩阵长这样：

图片

匹配临床信息和表达矩阵：

# match or merge meta and expression
exprSet=tpm%>%
  filter(Group=="Tumor")%>%
  tibble::add_column(ID = stringr::str_sub(rownames(.),1,12),
                     .before="Cancer") %>%
  filter(!duplicated(ID)) %>%
  remove_rownames %>%
  tibble::column_to_rownames("ID")
exprSet=exprSet[,-(1:2)]

s = intersect(rownames(meta),rownames(exprSet));length(s)
# [1] 9678
exprSet = exprSet[s,]
meta = meta[s,]
dim(exprSet)
# [1]  9678 19571
dim(meta)
#[1] 9678    6
identical(rownames(meta),rownames(exprSet))
#[1] TRUE

如果要合并临床信息和表达矩阵的话，可以执行 以下代码：

meta$ID=rownames(meta)
exprSet$ID=rownames(exprSet)
meta=inner_join(meta,exprSet)%>%
tibble::column_to_rownames("ID")

image.png

致谢：
TCGAbiolinks
生信技能树