TCGAbiolinks数据下载

微点

修改于 2019-10-30 17:24:44

1.2K0

修改于 2019-10-30 17:24:44

文章被收录于专栏：sci666sci666

TCGAbiolinks数据下载TCGA数据

下载TCGA数据的方法有很多，但比较好用的包我认为就是TCGAbiolinks,TCGAbiolinks是一个可用于检索，下载，并准备TCGA数据用于下游分析的R包，该包发表在著名的 核酸研究 杂志上。
TCGAbiolinks的优点在于具备一体化的下载整合，无需再使用复杂的方法对下载的单个数据重新进行整合，新手及临床医生尤其适合，我们的目的就是分析数据，没有必要去做些非必须的事。
UCSC xena 本来是一个很好的工具，但有一个问题是它存储的数据是log2(x+1) 转换的 RSEM nomalized data，关于这个数据有文章去使用，但我不确定能否直接用limma去做后续分析。笔者曾试图搞清楚这个数据究竟如何使用，发现众说纷纭，似乎UCSCxena本身也没有给出一个很清晰的答案，因此就表达谱数据而言，UCSCxena并不是最好的选择，当然也有说把log2转换的数据再转换回来的，总觉得并不合适。

load package

1if (!requireNamespace("BiocManager", quietly = TRUE))
2    install.packages("BiocManager")
3if (!requireNamespace("BiocManager", quietly = TRUE))
4    install.packages("TCGAbiolinks")
5library(TCGAbiolinks)
6library(DT)
7library(dplyr)
8library(SummarizedExperiment)

数据来源-根据TCGAbiolinks的官方说明

不同的数据来源 Legacy 与 harmonized
GDC Legacy Archive:以前在CGHUB和TCGA数据门户中存储的数据的原始数据，由TCGA数据协调中心（DCC）托管，在该门户中用GRCH37（HG19）和GRCH36（HG18）作为参考基因组
GDC harmonized database:可用数据与grch38（hg38）使用gdc生物信息学流程进行协调，该流程提供了生物标本和临床数据标准化的方法，简单讲就是对数据进行了一定标准化处理。harmonized数据库包括转录谱数据，甲基化数据，miRNA数据，但缺少芯片数据

数据检索

 1query.met <- GDCquery(project = "TCGA-COAD",
 2                      data.category = "DNA Methylation",
 3                      legacy = FALSE,
 4                      platform = c("Illumina Human Methylation 450"))
 5query.exp <- GDCquery(project = "TCGA-COAD",
 6                      data.category = "Transcriptome Profiling",
 7                      data.type = "Gene Expression Quantification", 
 8                      workflow.type = "HTSeq - FPKM-UQ")
 9query.exp[1:5,1:5]
10query.met[1:5,1:5]
11datatable(getResults(query.exp, cols = c("data_type","cases")),
12          filter = 'top',
13          options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
14          rownames = FALSE)

mRNA Expression

 1query <- GDCquery(project = "TCGA-BRCA",# Cancer type
 2                  data.category = "Transcriptome Profiling",
 3                  data.type = "Gene Expression Quantification", 
 4                  workflow.type = "HTSeq - Counts"# raw count 
 5                  )
 6GDCdownload(query,
 7            directory = "./project",#
 8            method = "api",
 9            files.per.chunk = 100)# 
10data <- GDCprepare(query)
11count_data=assay(data)
12count_data[1:5,1:5]
13dim(count_data)##56537 1222
14
15## clinical information
16colData(data)[1:5,1:5]
17
18## save data
19#save(count_data,file = "BRCA_count.Rdata")

clinical data

 1clinical <- GDCquery_clinic(project = "TCGA-BRCA", type = "clinical")
 2clinical[1:5,1:5]
 3dim(clinical)## 1097 68
 4
 5## save
 6#save(clinical,file="BRCA_clinical.Rdata")
 7#write.csv(clinical,file="TCGAbiolinks-BRCA-clinical.csv")
 8
 9## clinical-2
10clinical_2<-colData(data)
11#write.csv(clinical,file="TCGAbiolinks-BRCA-clinical.csv")
12
13## 获取所有TCGA的临床信息
14library(data.table)
15library(dplyr)
16library(regexPipes)
17clinical <- TCGAbiolinks:::getGDCprojects()$project_id %>% 
18    regexPipes::grep("TCGA",value=T) %>% ## TCGA
19    sort %>% 
20    plyr::alply(1,GDCquery_clinic, .progress = "text") %>% 
21    rbindlist
22dim(clinical)
23clinical[1:5,1:5]
24#readr::write_csv(clinical,path = paste0("BRCA_clin_indexed.csv"))

miRNA data

 1query <- GDCquery(project = "TCGA-BRCA", 
 2                    data.category = "Transcriptome Profiling", 
 3                    data.type = "miRNA Expression Quantification", 
 4                    workflow.type = "BCGSC miRNA Profiling")
 5## 检束结果
 6results<-getResults(query)
 7dim(results)
 8results[1:5,1:5]
 9colnames(results)
10
11## download data
12GDCdownload(query, 
13            method = "api", 
14            files.per.chunk = 20)# 减少下载失败风险
15mir_exp<- GDCprepare(query = query,
16                     summarizedExperiment=F)# set F
17
18## 数据结构
19dim(mir_exp)
20mir_exp[1:5,1:5]
21
22## save data
23#save(mir_exp,file="BRCA_miRNA_raw.Rdata")