TCGA臨床數據（化療數據、用藥反應等）和生存信息（4種臨床結局）

專注R語言在生物醫學中的使用

設為“星标”，精彩不錯過

GDC官網中有3種來源的TCGA臨床信息：

-和XML中的臨床信息主要有兩個區别:

下載的對象中的臨床數據是 data，但是和直接使用得到的臨床信息(這個也是 data)不太一樣，并且還會添加一些從相關文獻中得到的信息（這部分信息會添加前綴）。

從SE對象中提取的臨床信息

我們以TCGA-COAD為例進行演示。

下面是從TCGA-COAD原始的對象（通過::("TCGA-COAD")下載的）提取的臨床信息：

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

dim(clin_info)

##[1]524109

#colnames(clin_info)109列

一共109列臨床信息，還是挺多的，但是很多信息平常都用不到。而且并不是所有的癌症都是109列，因為每個癌症的信息是不一樣的，比如乳腺癌就會有分型信息（PAM50分型）。

這裡面并沒有用藥信息、化療信息、用藥反應等。

下面我們通過下載臨床數據。

library(TCGAbiolinks)

data

通過以下方式可以下載 data：

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

dim(clinical.indexed)

##[1]46170

#colnames(clinical.indexed)

一共70列，雖然列名中有一些和化療藥物、放療有關，但是并沒有用藥信息（可以自己點開看看）。

這個 data是從XML整理出來的，雖然沒有XML全，但是好在信息都是最新的，而且裡面的病理分期是分開的（TNM）是3列數據，XML中的病理分期是在一起的。所以通常來說這個臨床信息是最好用的(除了沒有藥物等信息基本沒啥缺點).

所以這個數據其實可以和XML中的數據結合使用，提取裡面好用的即可。

這種方法得到的臨床信息和從對象中提取的臨床信息也不是完全一樣（如下所示），但是比較關鍵的信息都是有的，比如病理分期、生存狀态、生存時間等：

#取個交集看看
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格式的臨床信息,點點點即可,如下圖所示,官網的選擇結果：

當然也可以通過下載XML格式的臨床數據,但是要注意，由于一個病人可能有多個臨床信息（比如有多次化療信息），所以一次隻能解析一個表格，需要通過參數指定要解析的數據：

使用下載官網的XML臨床信息，注意這種方式其中的幾個參數和官網都是一一對應的，所以數量啥的都是和官網完全一緻的：

#下載XML臨床數據
query1<-GDCquery(
project="TCGA-COAD",
data.category="Clinical",
data.type="ClinicalSupplement",
data.format="bcrxml"
)

##--------------------------------------
##oGDCquery:SearchinginGDCdatabase
##--------------------------------------
##Genomeofreference:hg38
##--------------------------------------------
##ooAccessingGDC.Thismighttakeawhile...
##--------------------------------------------
##oooProject:TCGA-COAD
##--------------------
##ooFilteringresults
##--------------------
##oooBydata.format
##oooBydata.type
##----------------
##ooCheckingdata
##----------------
##oooCheckingifthereareduplicatedcases
##oooCheckingifthereareresultsforthequery
##-------------------
##oPreparingoutput
##-------------------

GDCdownload(query1)

##DownloadingdataforprojectTCGA-COAD
##GDCdownloadwilldownload459files.Atotalof17.981211MB
##Downloadingas:Wed_Jan_10_20_02_01_2024.tar.gz

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

##Togetthefollowinginformationpleasechangetheclinical.infoargument
##=>new_tumor_events:new_tumor_event
##=>drugs:drug
##=>follow_ups:follow_up
##=>radiations:radiation
##Parsingfollowupversion:follow_up_v1.0
##Addingstageeventinformation
##Updatingdays_to_last_followupandvital_statusfromfollow_upinformationusinglastentry
##Parsingfollowupversion:follow_up_v1.0

dim(clinical.patient.xml)#也是459，和上面的圖片一樣
##[1]45976

#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]59324

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

##drug_namemeasure_of_response
##156LeucovorinCompleteResponse
##157FluorouracilCompleteResponse
##158FluorouracilCompleteResponse
##159FolinicacidPartialResponse
##1605-FluorouracilPartialResponse

解析信息,這個信息和無病生存期有關(但是這個我們有更好的選擇,用TCGA-CDR的數據):

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

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

##[1]11810

解析随訪信息，這裡面的随訪數據是最新的數據，但其實也是很久沒更新過了，-中的随訪信息就是從這裡更新.但是這裡面有很多重複的信息,比如某個人有多次随訪,都會記錄在這裡,如果你自己整理的話就很浪費時間.

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

##Parsingfollowupversion:follow_up_v1.0

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

##[1]57718

這裡面也有、藥物反應（CR/PR等）等信息。

如果需要最原始的XML信息，可以從XML文件中提取，如果要整理好的數據，可以使用-，這樣看這個就顯得不是那麼重要了。

下載BCR 臨床信息可以使用以下代碼：

query<-GDCquery(
project="TCGA-COAD",
data.category="Clinical",
data.type="ClinicalSupplement",
data.format="bcrbiotab"
)

GDCdownload(query)

clinical.BCRtab.all<-GDCprepare(query)#結果是7，也是和上面圖片一樣的

##Newnames:
##·`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中的也不是完全一樣哈。

簡單看下信息：

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中的信息無疑是最全的，但是也是最原始的，原始意味着對初學不友好，需要很多自己整理的過程，此時 data是更好的選擇，因為是經過整理過的。

我在包中添加了一個函數，可以1行代碼下載XML以及 data，并且會自動幫你提取其中的數據，并保存為rdata格式。

比如使用如下1行代碼：

getclinical("TCGA-COAD")

即可得到以下數據：

生存時間用哪個

生存時間用哪個？up還是？

可以看到當是Alive時，up是有時間的，此時是NA，當是Dead時，up是NA，此時是有時間的。

目前很多人都在用的做法是二者相加。但是也有人發現有些會同時有兩列時間，比如：

這個結果是從SE對象中提取出來的，我又看了下XML（ = ""）中的信息，發現此時隻有是有時間的：

所以我做了如下處理，如果生存狀态是Dead，就使用的時間，如果是Alive，就使用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]4584

head(coad_surv_gdc)

##sampleOS_PATIENTOS.time
##1TCGA-D5-6540-01A0TCGA-D5-6540491
##2TCGA-AA-3525-11A0TCGA-AA-35251
##3TCGA-AA-3525-01A0TCGA-AA-35251
##4TCGA-AA-3815-01A0TCGA-AA-38151005
##5TCGA-D5-6923-01A0TCGA-D5-6923378
##6TCGA-G4-6322-01A0TCGA-G4-6322792

這個可能也不是完全正确，但是我也找不到官方的說法到底用哪個，畢竟用啥的都有，我之前一直都是用的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]5394

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

##[1]4374

head(coad_surv_xena)

##sampleOS_PATIENTOS.time
##1TCGA-AA-3492-01A1TCGA-AA-34921
##2TCGA-AA-3492-11A1TCGA-AA-34921
##3TCGA-G4-6626-01A1TCGA-G4-66261
##4TCGA-AA-3525-01A0TCGA-AA-35251
##5TCGA-AA-3525-11A0TCGA-AA-35251
##6TCGA-AY-6196-01A0TCGA-AY-61966

取個交集看看：

#新發現：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]3815

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 TCGA Pan- Data to Drive High-

這篇文章提取出了4種TCGA的随訪結局（詳情自己讀文章）：

這個結果可以直接下載：

也可以通過XENA（）下載，但是兩者略有區别（不影響使用）。

我們下載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(`cancertypeabbreviation`=="COAD")%>%
select(1,2,26:33)

##sample_PATIENTOSOS.timeDSSDSS.timeDFIDFI.timePFI
##1TCGA-3L-AA1B-01TCGA-3L-AA1B0475047504750
##2TCGA-4N-A93T-01TCGA-4N-A93T01460146NANA0
##3TCGA-4T-AA8H-01TCGA-4T-AA8H0385038503850
##4TCGA-5M-AAT4-01TCGA-5M-AAT4149149NANA1
##5TCGA-5M-AAT6-01TCGA-5M-AAT612901290NANA1
##6TCGA-5M-AATE-01TCGA-5M-AATE0120001200NANA1
##7TCGA-A6-2670-01TCGA-A6-267007750775NANA0
##8TCGA-A6-2671-01TCGA-A6-26711133111331NANA1
##9TCGA-A6-2671-11TCGA-A6-26711133111331NANA1
##PFI.time
##1475
##2146
##3385
##449
##5219
##6810
##7775
##8535
##9535

結果直接給出了各種随訪結局和時間，其中OS和OS.time這兩個信息和GDC的差别不大（充分說明GDC官方每次更新的都不是這些内容）。

所以我把這個結果直接整合進中了，方便大家使用，這樣如果是OS，那麼可以用GDC整理的，也可以用這個，如果是DSS/DFI/PFI，就直接用這個就好了。

可通過以下方式加載：

library(easyTCGA)#安裝最新版本

data("tcga_cdr_surv")
dim(tcga_cdr_surv)

##[1]1259134

tcga_cdr_surv[1:6,26:33]

##OSOS.timeDSSDSS.timeDFIDFI.timePFIPFI.time
##1113551135517541754
##21167711677NANA1289
##30209102091153153
##414231423NANA1126
##513651365NANA150
##602703027030270302703

參考資料官方文檔：GDC官方文檔：#-data