TCGA数据库不仅仅存储了测序数据,更为珍贵的是TCGA数据库也存储了相关的临床信息,有了这些临床信息,后续可以进行多种分析,可以进行某个基因和临床预后的相关性分析,根据临床数据进行建模预测。总之,这一部分数据的深入挖掘一定是重头戏,通过这一部分数据和临床结合,能更好的贴近临床,贴近实际。
好了,依然是实践促进认知,直接以实例进行解读。
1. 使用TCGA官网数据进行整理
之前的数据下载中TCGA数据官网网页下载及gdc-client下载,从TCGA官网下载了这些文件
其中clinical文件下有有这两种文件
一种是TSV文件格式,另一种是json文件格式,所以后续可以有两种方法进行文件的整理。
tsv文件即Tab-separated values,指的是制表符分隔文件,与之类似的csv文件Comma-separated values,指的是逗号分隔符文件,后者更常见一些,所以使用R读取时,csv文件有自带的读取函数,tsv自行加载package再读取就会更方便一些。
1.1 读取json文件获得临床数据
1.1.1 从TCGA官网下载的json文件可以直接使用R读取,利用jsonlite 这个package中的fromJSON函数还是很便捷。
cli_kirc <- jsonlite::fromJSON("clinical.cart.2020-05-22.json")
1.1.2 读入之后,了解一下文件的特征
- 了解整个数据有什么
> class(cli_kirc)
[1] "data.frame"
> dim(cli_kirc)
[1] 530 4
> colnames(cli_kirc)
[1] "exposures" "demographic" "diagnoses" "case_id"
看似只有4列,其实这其中包含了很多的数据,从这当中也能深刻理解R数据存放的特点,数据框中可以放list,list里面可以接着存放data.frame。确实是很有趣的方面。
继续看看这四列还有什么信息。
- exposures中存放了什么
> class(cli_kirc$exposures)
[1] "list"
> length(cli_kirc$exposures)
[1] 530
> class(cli_kirc$exposures[[1]])
[1] "data.frame"
> dim(cli_kirc$exposures[[1]])
[1] 1 12
> colnames(cli_kirc$exposures[[1]])
[1] "created_datetime" "cigarettes_per_day" "exposure_id" "state" "updated_datetime"
[6] "weight" "alcohol_history" "years_smoked" "alcohol_intensity" "height"
[11] "bmi" "submitter_id"
在暴露因素中可以看到存储了12中信息,其中 "submitter_id" 是连结其他数据的节点。
- demographic 中存放的数据
> class(cli_kirc$demographic)
[1] "data.frame"
> dim(cli_kirc$demographic)
[1] 530 14
> colnames(cli_kirc$demographic)
[1] "ethnicity" "year_of_death" "year_of_birth" "age_at_index" "gender"
[6] "state" "updated_datetime" "demographic_id" "days_to_birth" "race"
[11] "vital_status" "submitter_id" "created_datetime" "days_to_death"
人口统计数据中,存放了14种信息,同样,"submitter_id" 是连结其他数据的节点。
- diagnoses 中存放的数据
> class(cli_kirc$diagnoses)
[1] "list"
> length(cli_kirc$diagnoses)
[1] 530
> class(cli_kirc$diagnoses[[1]])
[1] "data.frame"
> dim(cli_kirc$diagnoses[[1]])
[1] 1 29
> colnames(cli_kirc$diagnoses[[1]])
[1] "primary_diagnosis" "last_known_disease_status" "ajcc_pathologic_stage"
[4] "tissue_or_organ_of_origin" "state" "classification_of_tumor"
[7] "updated_datetime" "prior_malignancy" "days_to_last_known_disease_status"
[10] "progression_or_recurrence" "year_of_diagnosis" "prior_treatment"
[13] "morphology" "created_datetime" "ajcc_pathologic_m"
[16] "tumor_stage" "synchronous_malignancy" "ajcc_pathologic_n"
[19] "tumor_grade" "diagnosis_id" "age_at_diagnosis"
[22] "treatments" "days_to_diagnosis" "days_to_last_follow_up"
[25] "icd_10_code" "days_to_recurrence" "site_of_resection_or_biopsy"
[28] "submitter_id" "ajcc_pathologic_t"
诊断相关信息里面包含了29项信息,包含了我们非常关心的生存相关数据的收集,以及患者病理诊断的分级,TNM分期,病理诊断,治疗这些临床数据。那么有这些数据,下一步要怎么研究,那就要看自己的方向和兴趣了,这就体现了数据挖掘的魅力了,矿藏就在那里,就看用的人怎么用了,有些粗犷开发,有些细致提炼,这就是体现差距的地方了。
- case id主要还是所下载文件夹的名称
1.1.3 以上数据的整理
上面数据的处理可以分别将三个数据集建立单独的数据,然后根据需要进行合并,提取需要的数据。
demographic <- cli_kirc$demographic
exposure <- cli_kirc$exposures
diagnoses <- cli_kirc$diagnoses
> class(demographic);class(exposure);class(diagnoses)
[1] "data.frame"
[1] "list"
[1] "list"
1.1 利用tsv文件进行数据的整理
利用从TCGA官网下载的tsv文件,就显得很简单了,限速步骤只是对tsv文件的读入R即可。经过对比,readr 这个package比较好用。
library(readr)
cli_k <- read_tsv(file = 'clinical.tsv')
看看这个数据的大小和特征
> class(cli_k)
[1] "spec_tbl_df" "tbl_df" "tbl" "data.frame"
> dim(cli_k)
[1] 1060 149
要比json中读取的多很多,从列名看,因为当中有很多list中的内容没有展开,从观察项目,也就是行名来看,因为有些患者的治疗中有两种或几种治疗,在tsv文件中把这些都放在同等位置,所以相比于json文件中会多出这些。
实际上直接读入的tsv数据依然是比较复杂的类型,因为从数据类型中可以看出,本身就有四种类型,虽然这4种类型只是基于不同的数据读入模式,本质是有所类同的,可以说在具体分析中,具体对待。
接下来就是把相应的数据进行筛选,和测序数据进行合并,得到更多的挖掘信息。
2. 利用R package进行下载
针对于TCGA这样的宝库,肯定少不了有人做成R package,来解决TCGAS数据下载这些问题。虽然对于TCGA数据下载有好几种,但是仍然需要注意的是,这些R 包下载的数据的时效性是不能保证的,下载后最好还是能够看一看是否是最新的数据。TCGA官网对于临床数据的更新要比
2.1 利用RTCGA.clinical package进行下载
例子采用RTCGA.clinical package下载TCGA-BLCA数据的代码
library(RTCGA.clinical)
??RTCGA.clinical ##多多查看说明文档
BLCA_clinical <- BLCA.clinical
write.csv(BLCA_clinical,quote = F,"BLCA_clinical_Data.csv")
其实按照提取临床信息,上面的已经将临床数据下载到了,探索一下里面的数据有什么。
> class(BLCA_clinical)
[1] "data.frame"
> dim(BLCA_clinical)
[1] 412 2129
这个数据集列名确实很多,只需要将我们需要的列名挑选出来做先骨干的分析就可以,下面的是BLCA数据中,要做生存分析的数据。
BLCA_clinicaldata<-BLCA_clinical[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'
)]
str(BLCA_clinicaldata)
BLCA_clinicaldata[1:4,1:4]
rownames(BLCA_clinicaldata) <- NULL
BLCA_clinicaldata <- tibble::column_to_rownames(BLCA_clinicaldata, var = 'patient.bcr_patient_barcode')
##使用barcode,只是方便后续与测序数据衔接
看一看筛选后的数据
> str(BLCA_clinicaldata)
'data.frame': 412 obs. of 8 variables:
$ patient.bcr_patient_barcode : chr "tcga-2f-a9ko" "tcga-2f-a9kp" "tcga-2f-a9kq" "tcga-2f-a9kr" ...
$ patient.vital_status : chr "alive" "dead" "alive" "alive" ...
$ patient.days_to_death : chr NA "364" NA NA ...
$ patient.days_to_last_followup : chr "678" NA "2555" "3148" ...
$ patient.race : chr "white" "white" "white" NA ...
$ patient.age_at_initial_pathologic_diagnosis: chr "63" "66" "69" "59" ...
$ patient.gender : chr "male" "male" "male" "female" ...
$ patient.stage_event.pathologic_stage : chr "stage iv" "stage iv" "stage iii" "stage iii" ...
其实以上的信息包括了barcode,status:即就是生存状态,days_to_death,days_to_last_followup ,age_at_initial_pathologic_diagnosisrace:诊断时间,我觉得也可以理解成发病时间,可以用这个时间做一些分析,gender 性别
2.2 利用TCGAbiolinks package进行下载
TCGAbiolinks是一个分析处理TCGA数据的R包,通过GDC API来查询和下载TCGA的数据,同时提供了差异分析,生存分析,富集分析等常见的分析功能
在bioconductor有这个包的整个说明和代码,可以对照进行。还有关于这个包的使用的讲解文章生信修炼手册的使用TCGAbiolinks下载TCGA的数据,生信技能树的TCGA数据下载—TCGAbiolinks包参数详解。
有文章对TCGA-biolinks的数据,Xene上面下载的数据和TCGA官网数据进行了比较,同样,结论依然是很多二次数据库更新不一定能跟得上TCGA的官网数据。文章见UCSC Xena和TCGAbiolinks的HNSC临床数据不一致。
所以总体而言,还是直接从TCGA官网进行下载数据会更新,当然更权威。
网友评论