TCGA临床数据（化疗数据、用药反应等）和生存信息（4种临床结局）整理

GDC官网中有3种来源的TCGA临床信息：

• XML files: 原始信息，包含所有临床数据，包括化疗信息、使用的药物名字等，而且会包含每次更新的信息（比如随访时间、生存信息等）；
• BCR Biotab: 从XML文件中解析得到的TSV文件（tab分割文件）
• indexed clinical: 也是从XML文件中提取的，但是是经过整理的（比如这里面包含的生存信息就是只保留最新的）

indexed-clinical和XML中的临床信息主要有两个区别:

• XML有包含更多信息，比如化疗信息、药物信息、随访数据（会更新）、生物样本等，indexed-clinical可以看做是XML的子集
• indexed-clinical中的信息是经过整理的，里面的数据都是最新的（比如随访数据）。如果某个患者在第一次随访时是活着的，但是第二次随访是死亡了，那么在indexed-clinical中这个人的生存状态就是死亡（最新的），但是在XML文件中会有两个记录，一个是第一次的活着，还有一个是第二次的死亡。

TCGAbiolinks下载的SummarisedExperiment对象中的临床数据是indexed clinical data，但是和直接使用GDCquery_clinic得到的临床信息(这个也是indexed clinical data)不太一样，并且还会添加一些从相关文献中得到的信息（这部分信息会添加paper_前缀）。

从SE对象中提取的临床信息

我们以TCGA-COAD为例进行演示。

下面是从TCGA-COAD原始的SummarisedExperiment对象（通过easyTCGA::getmrnaexpr("TCGA-COAD")下载的）提取的临床信息：

load("G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-COAD_clinical.rdata")

dim(clin_info)

## [1] 524 109

#colnames(clin_info) 109列

一共109列临床信息，还是挺多的，但是很多信息平常都用不到。而且并不是所有的癌症都是109列，因为每个癌症的信息是不一样的，比如乳腺癌就会有分型信息（PAM50分型）。

这里面并没有用药信息、化疗信息、用药反应等。

下面我们通过TCGAbiolinks下载临床数据。

library(TCGAbiolinks)

indexed clinical data

通过以下方式可以下载indexed clinical data：

clinical.indexed <- GDCquery_clinic(project = "TCGA-COAD", type = "clinical")

dim(clinical.indexed)

## [1] 461  70

#colnames(clinical.indexed)

一共70列，虽然列名中有一些和化疗药物、放疗有关，但是并没有用药信息（可以自己点开看看）。

这个indexed clinical data是从XML整理出来的，虽然没有XML全，但是好在信息都是最新的，而且里面的病理分期是分开的（TNM）是3列数据，XML中的病理分期是在一起的。所以通常来说这个临床信息是最好用的(除了没有药物等信息基本没啥缺点).

所以这个数据其实可以和XML中的数据结合使用，提取里面好用的即可。

这种方法得到的临床信息和从SummarisedExperiment对象中提取的临床信息也不是完全一样（如下所示），但是比较关键的信息都是有的，比如病理分期、生存状态、生存时间等：

# 取个交集看看
intersect(colnames(clin_info), colnames(clinical.indexed))

 [1] "submitter_id"                "state"                      
 [3] "synchronous_malignancy"      "ajcc_pathologic_stage"      
 [5] "days_to_diagnosis"           "last_known_disease_status"  
 [7] "tissue_or_organ_of_origin"   "days_to_last_follow_up"     
 [9] "age_at_diagnosis"            "primary_diagnosis"          
[11] "prior_malignancy"            "year_of_diagnosis"          
[13] "prior_treatment"             "ajcc_staging_system_edition"
[15] "ajcc_pathologic_t"           "morphology"                 
[17] "ajcc_pathologic_n"           "ajcc_pathologic_m"          
[19] "classification_of_tumor"     "diagnosis_id"               
[21] "icd_10_code"                 "site_of_resection_or_biopsy"
[23] "tumor_grade"                 "progression_or_recurrence"  
[25] "alcohol_history"             "exposure_id"                
[27] "race"                        "gender"                     
[29] "ethnicity"                   "vital_status"               
[31] "age_at_index"                "days_to_birth"              
[33] "year_of_birth"               "demographic_id"             
[35] "days_to_death"               "bcr_patient_barcode"        
[37] "year_of_death"   

XML临床信息

可以通过官网下载XML格式的临床信息,点点点即可,如下图所示,官网的选择结果：

当然也可以通过TCGAbiolinks下载XML格式的临床数据,但是要注意，由于一个病人可能有多个临床信息（比如有多次化疗信息），所以一次只能解析一个表格，需要通过clinical.info参数指定要解析的数据：

使用TCGAbiolinks下载官网的XML临床信息，注意这种方式其中的几个参数和官网都是一一对应的，所以数量啥的都是和官网完全一致的：

# 下载XML临床数据
query1 <- GDCquery(
      project = "TCGA-COAD",
      data.category = "Clinical",
      data.type = "Clinical Supplement",
      data.format = "bcr xml"
    )

## --------------------------------------
## o GDCquery: Searching in GDC database
## --------------------------------------
## Genome of reference: hg38
## --------------------------------------------
## oo Accessing GDC. This might take a while...
## --------------------------------------------
## ooo Project: TCGA-COAD
## --------------------
## oo Filtering results
## --------------------
## ooo By data.format
## ooo By data.type
## ----------------
## oo Checking data
## ----------------
## ooo Checking if there are duplicated cases
## ooo Checking if there are results for the query
## -------------------
## o Preparing output
## -------------------

GDCdownload(query1)

## Downloading data for project TCGA-COAD
## GDCdownload will download 459 files. A total of 17.981211 MB
## Downloading as: Wed_Jan_10_20_02_01_2024.tar.gz

# 解析patient信息,指定clinical.info
clinical.patient.xml <- GDCprepare_clinic(query1, clinical.info = "patient")

## To get the following information please change the clinical.info argument
## => new_tumor_events: new_tumor_event 
## => drugs: drug 
## => follow_ups: follow_up 
## => radiations: radiation
## Parsing follow up version: follow_up_v1.0
## Adding stage event information
## Updating days_to_last_followup and vital_status from follow_up information using last entry
## Parsing follow up version: follow_up_v1.0

dim(clinical.patient.xml) # 也是459，和上面的图片一样
## [1] 459  76

#colnames(clinical.patient.xml)

这里面的很多信息和上面两种是重复的，毕竟上面两种都是从XML中提取的。而且这个信息看似很全,但是很多信息都是NA,比如其中的生存信息就是不全的,大部分都是NA.

解析药物信息，这里面就有各种化疗药物及用药反应信息：

clinical.drug.xml <- GDCprepare_clinic(query1, clinical.info = "drug")

colnames(clinical.drug.xml)

##  [1] "bcr_patient_barcode"                "tx_on_clinical_trial"              
##  [3] "regimen_number"                     "bcr_drug_barcode"                  
##  [5] "bcr_drug_uuid"                      "total_dose"                        
##  [7] "total_dose_units"                   "prescribed_dose"                   
##  [9] "prescribed_dose_units"              "number_cycles"                     
## [11] "days_to_drug_therapy_start"         "days_to_drug_therapy_end"          
## [13] "therapy_types"                      "drug_name"                         
## [15] "clinical_trail_drug_classification" "regimen_indication"                
## [17] "regimen_indication_notes"           "route_of_administrations"          
## [19] "therapy_ongoing"                    "measure_of_response"               
## [21] "day_of_form_completion"             "month_of_form_completion"          
## [23] "year_of_form_completion"            "project"

dim(clinical.drug.xml)

## [1] 593  24

# 用药信息和药物反应,CP,PR等
clinical.drug.xml[156:160,c("drug_name","measure_of_response")]

##          drug_name measure_of_response
## 156     Leucovorin   Complete Response
## 157   Fluorouracil   Complete Response
## 158   Fluorouracil   Complete Response
## 159   Folinic acid    Partial Response
## 160 5-Fluorouracil    Partial Response

解析new_tumor_event信息,这个信息和无病生存期有关(但是这个我们有更好的选择,用TCGA-CDR的数据):

clinical.nte.xml <- GDCprepare_clinic(query1, clinical.info = "new_tumor_event")

#colnames(clinical.nte.xml)
dim(clinical.nte.xml)

## [1] 118  10

解析随访信息，这里面的随访数据是最新的数据，但其实也是很久没更新过了，indexed-clinicla中的随访信息就是从这里更新.但是这里面有很多重复的信息,比如某个人有多次随访,都会记录在这里,如果你自己整理的话就很浪费时间.

clinical.followup.xml <- GDCprepare_clinic(query1, clinical.info = "follow_up")

## Parsing follow up version: follow_up_v1.0

#colnames(clinical.followup.xml)
dim(clinical.followup.xml)

## [1] 577  18

这里面也有new_tumor_event、药物反应（CR/PR等）等信息。

Biotab

如果需要最原始的XML信息，可以从XML文件中提取，如果要整理好的数据，可以使用indexed-clinical，这样看这个Biotab就显得不是那么重要了。

下载BCR Biotab临床信息可以使用以下代码：

query <- GDCquery(
    project = "TCGA-COAD", 
    data.category = "Clinical",
    data.type = "Clinical Supplement", 
    data.format = "bcr biotab"
)

GDCdownload(query)

clinical.BCRtab.all <- GDCprepare(query) # 结果是7，也是和上面图片一样的

## New names:
## • `history_other_malignancy` -> `history_other_malignancy...11`
## • `history_other_malignancy` -> `history_other_malignancy...44`

# 看看包含哪些信息
names(clinical.BCRtab.all)

## [1] "clinical_nte_coad"                "clinical_patient_coad"           
## [3] "clinical_omf_v4.0_coad"           "clinical_follow_up_v1.0_nte_coad"
## [5] "clinical_drug_coad"               "clinical_radiation_coad"         
## [7] "clinical_follow_up_v1.0_coad"

有随访信息，也有放化疗信息等，但是和XML中的也不是完全一样哈。

简单看下new_tumor_event信息：

colnames(clinical.BCRtab.all[["clinical_nte_coad"]])

##  [1] "bcr_patient_uuid"                                    
##  [2] "bcr_patient_barcode"                                 
##  [3] "new_tumor_event_dx_days_to"                          
##  [4] "new_tumor_event_site_surgery"                        
##  [5] "new_tumor_event_radiation_tx"                        
##  [6] "new_tumor_event_pharmaceutical_tx"                   
##  [7] "days_to_new_tumor_event_additional_surgery_procedure"
##  [8] "new_neoplasm_event_occurrence_anatomic_site"         
##  [9] "new_neoplasm_event_type"                             
## [10] "new_neoplasm_occurrence_anatomic_site_text"          
## [11] "new_tumor_event_additional_surgery_procedure"        
## [12] "progression_determined_by"                           
## [13] "residual_disease_post_new_tumor_event_margin_status"

通过以上简单的探索发现，XML中的信息无疑是最全的，但是也是最原始的，原始意味着对初学不友好，需要很多自己整理的过程，此时indexed clinical data是更好的选择，因为是经过整理过的。

我在easyTCGA包中添加了一个getclinical函数，可以1行代码下载XML以及indexed clinical data，并且会自动帮你提取其中的数据，并保存为rdata格式。

比如使用如下1行代码：

getclinical("TCGA-COAD")

即可得到以下数据：

生存时间用哪个

生存时间用哪个？days_to_last_follow_up还是days_to_death？

可以看到当vital_status是Alive时，days_to_last_follow_up是有时间的，此时days_to_death是NA，当vital_status是Dead时，days_to_last_follow_up是NA，此时days_to_death是有时间的。

目前很多人都在用的做法是二者相加。但是也有人发现有些sample会同时有两列时间，比如：

这个结果是从SE对象中提取出来的，我又看了下XML（clinical.info = "follow_up"）中的信息，发现此时只有days_to_death是有时间的：

所以我做了如下处理，如果生存状态是Dead，就使用days_to_death的时间，如果是Alive，就使用days_to_last_follow_up的时间。

我把这个结果和从XENA下载的结果比较了一下，发现是正确的：

整理一下SE对象中的生存信息：

library(tidyverse)

coad_surv_gdc <- clin_info %>% 
  select(sample_submitter_id, patient, vital_status,days_to_last_follow_up,
         days_to_death) %>% 
  mutate(OS_time = if_else(vital_status == "Dead", days_to_death,
                            days_to_last_follow_up),
         vital_status = if_else(vital_status == "Dead",1,0)) %>% 
  select(sample=sample_submitter_id, OS=vital_status,
         `_PATIENT`=patient,OS.time=OS_time
         )

rownames(coad_surv_gdc) <- NULL

coad_surv_gdc %>% distinct(`_PATIENT`,.keep_all = T) %>% dim()

## [1] 458   4

head(coad_surv_gdc)

##             sample OS     _PATIENT OS.time
## 1 TCGA-D5-6540-01A  0 TCGA-D5-6540     491
## 2 TCGA-AA-3525-11A  0 TCGA-AA-3525       1
## 3 TCGA-AA-3525-01A  0 TCGA-AA-3525       1
## 4 TCGA-AA-3815-01A  0 TCGA-AA-3815    1005
## 5 TCGA-D5-6923-01A  0 TCGA-D5-6923     378
## 6 TCGA-G4-6322-01A  0 TCGA-G4-6322     792

这个可能也不是完全正确，但是我也找不到官方的说法到底用哪个，毕竟用啥的都有，我之前一直都是用的days_to_last_follow_up......如果大家知道更好的方式，欢迎告诉我。

把从XENA下载的生存数据读取进来：

coad_surv_xena <- data.table::fread("../000files/TCGA-COAD.survival.tsv",
                                    data.table = F)
names(coad_surv_xena)

## [1] "sample"   "OS"       "_PATIENT" "OS.time"

dim(coad_surv_xena)

## [1] 539   4

coad_surv_xena %>% distinct(`_PATIENT`,.keep_all = T) %>% dim()

## [1] 437   4

head(coad_surv_xena)

##             sample OS     _PATIENT OS.time
## 1 TCGA-AA-3492-01A  1 TCGA-AA-3492       1
## 2 TCGA-AA-3492-11A  1 TCGA-AA-3492       1
## 3 TCGA-G4-6626-01A  1 TCGA-G4-6626       1
## 4 TCGA-AA-3525-01A  0 TCGA-AA-3525       1
## 5 TCGA-AA-3525-11A  0 TCGA-AA-3525       1
## 6 TCGA-AY-6196-01A  0 TCGA-AY-6196       6

取个交集看看：

# 新发现：inner_join不会自动去重，intersect会自动去重！
common_surv <- inner_join(coad_surv_gdc %>% 
                            distinct(`_PATIENT`,.keep_all = T),
                          coad_surv_xena %>% 
                            distinct(`_PATIENT`,.keep_all = T),
                          by = c("sample", "OS", "OS.time"))
dim(common_surv)

## [1] 381   5

381个相同的，效果还是不错的。

探索下不同的那些。

diff_surv1 <- setdiff(coad_surv_gdc %>% 
                        distinct(`_PATIENT`,.keep_all = T),
                      coad_surv_xena %>% 
                        distinct(`_PATIENT`,.keep_all = T)) %>% 
  drop_na() %>% 
  filter(OS.time > 0) #在gdc，不在xena

diff_surv2 <- setdiff(coad_surv_xena %>% distinct(`_PATIENT`,.keep_all = T),
                      coad_surv_gdc %>% distinct(`_PATIENT`,.keep_all = T)) #在xena，不在gdc

TCGA-PAN CDR

GDC官方发表在cell上的文章：An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics

这篇文章提取出了4种TCGA的随访结局（详情自己读文章）：

• OS：overall survival，总生存
• DSS：disease-specific survival，疾病特异生存
• DFI：disease-free interval，无病生存
• PFI：progression-free interval，无进展生存

这个结果可以直接下载：https://api.gdc.cancer.gov/data/1b5f413e-a8d1-4d10-92eb-7c4ae739ed81

也可以通过XENA（https://tcga-pancan-atlas-hub.s3.us-east-1.amazonaws.com/download/Survival_SupplementalTable_S1_20171025_xena_sp）下载，但是两者略有区别（不影响使用）。

我们下载XENA的，读取进来看看：

tcga_pan_surv <- data.table::fread("G:/bioinfo/000files/Survival_SupplementalTable_S1_20171025_xena_sp", data.table = F)

# 筛选TCGA-COAD的
tcga_pan_surv %>% 
  filter(`cancer type abbreviation` == "COAD") %>% 
  select(1,2,26:33)

##              sample     _PATIENT OS OS.time DSS DSS.time DFI DFI.time PFI
## 1   TCGA-3L-AA1B-01 TCGA-3L-AA1B  0     475   0      475   0      475   0
## 2   TCGA-4N-A93T-01 TCGA-4N-A93T  0     146   0      146  NA       NA   0
## 3   TCGA-4T-AA8H-01 TCGA-4T-AA8H  0     385   0      385   0      385   0
## 4   TCGA-5M-AAT4-01 TCGA-5M-AAT4  1      49   1       49  NA       NA   1
## 5   TCGA-5M-AAT6-01 TCGA-5M-AAT6  1     290   1      290  NA       NA   1
## 6   TCGA-5M-AATE-01 TCGA-5M-AATE  0    1200   0     1200  NA       NA   1
## 7   TCGA-A6-2670-01 TCGA-A6-2670  0     775   0      775  NA       NA   0
## 8   TCGA-A6-2671-01 TCGA-A6-2671  1    1331   1     1331  NA       NA   1
## 9   TCGA-A6-2671-11 TCGA-A6-2671  1    1331   1     1331  NA       NA   1
##     PFI.time
## 1        475
## 2        146
## 3        385
## 4         49
## 5        219
## 6        810
## 7        775
## 8        535
## 9        535

结果直接给出了各种随访结局和时间，其中OS和OS.time这两个信息和GDC的差别不大（充分说明GDC官方每次更新的都不是这些内容）。

所以我把这个结果直接整合进easyTCGA中了，方便大家使用，这样如果是OS，那么可以用GDC整理的，也可以用这个，如果是DSS/DFI/PFI，就直接用这个就好了。

可通过以下方式加载：

library(easyTCGA) # 安装最新版本

data("tcga_cdr_surv")
dim(tcga_cdr_surv)

## [1] 12591    34

tcga_cdr_surv[1:6,26:33]

##   OS OS.time DSS DSS.time DFI DFI.time PFI PFI.time
## 1  1    1355   1     1355   1      754   1      754
## 2  1    1677   1     1677  NA       NA   1      289
## 3  0    2091   0     2091   1       53   1       53
## 4  1     423   1      423  NA       NA   1      126
## 5  1     365   1      365  NA       NA   1       50
## 6  0    2703   0     2703   0     2703   0     2703

参考资料

1. TCGAbiolinks官方文档：https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/clinical.html
2. GDC官方文档：https://docs.gdc.cancer.gov/Encyclopedia/pages/Clinical_Data/#clinical-data