TCGAbiolinks: Downloading and preparing files for analysis
28 October 2020
TCGAbiolinks has provided a few functions to download and prepare data from GDC for analysis. This section starts by explaining the different downloads methods and the SummarizedExperiment object, which is the default data structure used in TCGAbiolinks, followed by some examples.
Downloading and preparing data for analysis
Data download: Methods differences
There are two methods to download GDC data using TCGAbiolinks:
- client: this method creates a MANIFEST file and download the data using GDC Data Transfer Tool this method is more reliable but it might be slower compared to the api method.
- api: this methods used the GDC Application Programming Interface (API) to download the data. This will create a MANIFEST file and the data downloaded will be compressed into a tar.gz file. If the size and the number of the files are too big this tar.gz will be too big which might have a high probability of download failure. To solve that we created the
files.per.chunkargument which will split the files into small chunks, for example, if chunks.per.download is equal to 10 we will download only 10 files inside each tar.gz.
Data prepared: SummarizedExperiment object
A SummarizedExperiment object has three main matrices that can be accessed using the SummarizedExperiment package):
- Sample matrix information is accessed via
colData(data): stores sample information. TCGAbiolinks will add indexed clinical data and subtype information from marker TCGA papers. - Assay matrix information is accessed via
assay(data): stores molecular data - Feature matrix information (gene information) is accessed via
rowRanges(data): stores metadata about the features, including their genomic ranges
Summarized Experiment: annotation information
When using the function GDCprepare there is an argument called SummarizedExperiment which defines the output type a Summarized Experiment (default option) or a data frame. To create a summarized Experiment object we annotate the data with genomic positions with last patch release version of the genome available. For legacy data (data aligned to hg19) TCGAbiolinks is using GRCh37.p13 and for harmonized data (data aligned to hg38) now it is using GRCh38.p7 (May 2017).
Unfortunately, some of the updates changes/remove gene symbols, change coordinates, etc. Which might introduce some loss of data. For example, if the gene was removed we cannot map it anymore and that information will be lost in the SummarizedExperiment.
If you set SummarizedExperiment to FALSE, you will get the data unmodified just as they are in the files and ad your own annotation.
Also, there are no updated for DNA methylation data. But the last metadata available can be found here: http://zwdzwd.github.io/InfiniumAnnotation
Arguments
GDCdownload
| Argument | Description |
|---|---|
| query | A query for GDCquery function |
| token.file | Token file to download controlled data (only for method = “client”) |
| method | Uses the API (POST method) or gdc client tool. Options “api”, “client”. API is faster, but the data might get corrupted in the download, and it might need to be executed again |
| directory | Directory/Folder where the data was downloaded. Default: GDCdata |
| files.per.chunk | This will make the API method only download n (files.per.chunk) files at a time. This may reduce the download problems when the data size is too large. Expected a integer number (example files.per.chunk = 6) |
GDCprepare
| Argument | Description |
|---|---|
| query | A query for GDCquery function |
| save | Save result as RData object? |
| save.filename | Name of the file to be save if empty an automatic will be created |
| directory | Directory/Folder where the data was downloaded. Default: GDCdata |
| summarizedExperiment | Create a summarizedExperiment? Default TRUE (if possible) |
| remove.files.prepared | Remove the files read? Default: FALSE This argument will be considered only if save argument is set to true |
| add.gistic2.mut | If a list of genes (gene symbol) is given, columns with gistic2 results from GDAC firehose (hg19) and a column indicating if there is or not mutation in that gene (hg38) (TRUE or FALSE - use the MAF file for more information) will be added to the sample matrix in the summarized Experiment object. |
| mut.pipeline | If add.gistic2.mut is not NULL this field will be taken in consideration. Four separate variant calling pipelines are implemented for GDC data harmonization. Options: muse, varscan2, somaticsniper, MuTect2. For more information: https://gdc-docs.nci.nih.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline/ |
| mutant_variant_classification | List of mutant_variant_classification that will be consider a sample mutant or not. Default: “Frame_Shift_Del”, “Frame_Shift_Ins”, “Missense_Mutation”, “Nonsense_Mutation”, “Splice_Site”, “In_Frame_Del”, “In_Frame_Ins”, “Translation_Start_Site”, “Nonstop_Mutation” |
Search and download data from legacy database using GDC api method
In this example we will download gene expression data from legacy database (data aligned against genome of reference hg19) using GDC api method and we will show object data and metadata.
query <- GDCquery(project = "TCGA-GBM",
data.category = "Gene expression",
data.type = "Gene expression quantification",
platform = "Illumina HiSeq",
file.type = "normalized_results",
experimental.strategy = "RNA-Seq",
barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"),
legacy = TRUE)
GDCdownload(query, method = "api", files.per.chunk = 10)
data <- GDCprepare(query)# Gene expression aligned against hg19.
datatable(as.data.frame(colData(data)),
options = list(scrollX = TRUE, keys = TRUE, pageLength = 5),
rownames = FALSE)# Only first 100 to make render faster
datatable(assay(data)[1:100,],
options = list(scrollX = TRUE, keys = TRUE, pageLength = 5),
rownames = TRUE)## GRanges object with 100 ranges and 3 metadata columns:
## seqnames ranges strand | gene_id entrezgene
## <Rle> <IRanges> <Rle> | <character> <numeric>
## A1BG chr19 58856544-58864865 - | A1BG 1
## A2M chr12 9220260-9268825 - | A2M 2
## NAT1 chr8 18027986-18081198 + | NAT1 9
## NAT2 chr8 18248755-18258728 + | NAT2 10
## RP11-986E7.7 chr14 95058395-95090983 + | RP11-986E7.7 12
## ... ... ... ... . ... ...
## ADORA1 chr1 203059782-203136533 + | ADORA1 134
## ADORA2A chr22 24813847-24838328 + | ADORA2A 135
## ADORA2B chr17 15848231-15879060 + | ADORA2B 136
## ADORA3 chr1 112025970-112106584 - | ADORA3 140
## ADPRH chr3 119298115-119308792 + | ADPRH 141
## ensembl_gene_id
## <character>
## A1BG ENSG00000121410
## A2M ENSG00000175899
## NAT1 ENSG00000171428
## NAT2 ENSG00000156006
## RP11-986E7.7 ENSG00000273259
## ... ...
## ADORA1 ENSG00000163485
## ADORA2A ENSG00000128271
## ADORA2B ENSG00000170425
## ADORA3 ENSG00000121933
## ADPRH ENSG00000144843
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengthsSearch and download data for two samples from database
In this example we will download gene expression quantification from harmonized database (data aligned against genome of reference hg38). Also, it shows the object data and metadata.
# Gene expression aligned against hg38
query <- GDCquery(project = "TCGA-GBM",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - FPKM-UQ",
barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"))
GDCdownload(query)
data <- GDCprepare(query)datatable(as.data.frame(colData(data)),
options = list(scrollX = TRUE, keys = TRUE, pageLength = 5),
rownames = FALSE)datatable(assay(data)[1:100,],
options = list(scrollX = TRUE, keys = TRUE, pageLength = 5),
rownames = TRUE)## GRanges object with 100 ranges and 3 metadata columns:
## seqnames ranges strand | ensembl_gene_id
## <Rle> <IRanges> <Rle> | <character>
## ENSG00000000003 chrX 100627109-100639991 - | ENSG00000000003
## ENSG00000000005 chrX 100584802-100599885 + | ENSG00000000005
## ENSG00000000419 chr20 50934867-50958555 - | ENSG00000000419
## ENSG00000000457 chr1 169849631-169894267 - | ENSG00000000457
## ENSG00000000460 chr1 169662007-169854080 + | ENSG00000000460
## ... ... ... ... . ...
## ENSG00000005421 chr7 95297676-95324707 - | ENSG00000005421
## ENSG00000005436 chr2 75652000-75710989 - | ENSG00000005436
## ENSG00000005448 chr2 74421678-74425755 + | ENSG00000005448
## ENSG00000005469 chr7 87345681-87399795 + | ENSG00000005469
## ENSG00000005471 chr7 87401697-87480435 - | ENSG00000005471
## external_gene_name original_ensembl_gene_id
## <character> <character>
## ENSG00000000003 TSPAN6 ENSG00000000003.13
## ENSG00000000005 TNMD ENSG00000000005.5
## ENSG00000000419 DPM1 ENSG00000000419.11
## ENSG00000000457 SCYL3 ENSG00000000457.12
## ENSG00000000460 C1orf112 ENSG00000000460.15
## ... ... ...
## ENSG00000005421 PON1 ENSG00000005421.7
## ENSG00000005436 GCFC2 ENSG00000005436.12
## ENSG00000005448 WDR54 ENSG00000005448.15
## ENSG00000005469 CROT ENSG00000005469.10
## ENSG00000005471 ABCB4 ENSG00000005471.14
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengthsGDCprepare: Outputs
This function is still under development, it is not working for all cases. See the tables below with the status. Examples of query, download, prepare can be found in this gist.
Harmonized data
| Data.category | Data.type | Workflow Type | Status |
|---|---|---|---|
| Transcriptome Profiling | Gene Expression Quantification | HTSeq - Counts | Data frame or SE (losing 5% of information when mapping to genomic regions) |
| HTSeq - FPKM-UQ | Returning only a (losing 5% of information when mapping to genomic regions) | ||
| HTSeq - FPKM | Returning only a (losing 5% of information when mapping to genomic regions) | ||
| Isoform Expression Quantification | Not needed | ||
| miRNA Expression Quantification | Not needed | Returning only a dataframe for the moment | |
| Copy number variation | Copy Number Segment | Returning only a dataframe for the moment | |
| Masked Copy Number Segment | Returning only a dataframe for the moment | ||
| Gene Level Copy Number Scores | Returning only a dataframe for the moment | ||
| Simple Nucleotide Variation | |||
| Raw Sequencing Data | |||
| Biospecimen | Slide Image | ||
| Biospecimen | Biospecimen Supplement | ||
| Clinical |
Legacy data
| Data.category | Data.type | Platform | file.type | Status |
|---|---|---|---|---|
| Transcriptome Profiling | ||||
| Copy number variation | - | Affymetrix SNP Array 6.0 | nocnv_hg18.seg | Working |
| - | Affymetrix SNP Array 6.0 | hg18.seg | Working | |
| - | Affymetrix SNP Array 6.0 | nocnv_hg19.seg | Working | |
| - | Affymetrix SNP Array 6.0 | hg19.seg | Working | |
| - | Illumina HiSeq | Several | Working | |
| Simple Nucleotide Variation | Simple somatic mutation | |||
| Raw Sequencing Data | ||||
| Biospecimen | ||||
| Clinical | ||||
| Protein expression | MDA RPPA Core | - | Working | |
| Gene expression | Gene expression quantification | Illumina HiSeq | normalized_results | Working |
| Illumina HiSeq | results | Working | ||
| HT_HG-U133A | - | Working | ||
| AgilentG4502A_07_2 | - | Data frame only | ||
| AgilentG4502A_07_1 | - | Data frame only | ||
| HuEx-1_0-st-v2 | FIRMA.txt | Not Preparing | ||
| gene.txt | Not Preparing | |||
| Isoform expression quantification | ||||
| miRNA gene quantification | ||||
| Exon junction quantification | ||||
| Exon quantification | ||||
| miRNA isoform quantification | ||||
| DNA methylation | Illumina Human Methylation 450 | Not used | Working | |
| Illumina Human Methylation 27 | Not used | Working | ||
| Illumina DNA Methylation OMA003 CPI | Not used | Working | ||
| Illumina DNA Methylation OMA002 CPI | Not used | Working | ||
| Illumina Hi Seq | Not working | |||
| Raw Microarray Data | ||||
| Structural Rearrangement | ||||
| Other |
Examples
Legacy archive
DNA methylation: Get all TCGA IDAT files
#-------------------------------------------------------
# Example to idat files from TCGA projects
#-------------------------------------------------------
projects <- TCGAbiolinks:::getGDCprojects()$project_id
projects <- projects[grepl('^TCGA',projects,perl=T)]
match.file.cases.all <- NULL
for(proj in projects){
print(proj)
query <- GDCquery(project = proj,
data.category = "Raw microarray data",
data.type = "Raw intensities",
experimental.strategy = "Methylation array",
legacy = TRUE,
file.type = ".idat",
platform = "Illumina Human Methylation 450")
match.file.cases <- getResults(query,cols=c("cases","file_name"))
match.file.cases$project <- proj
match.file.cases.all <- rbind(match.file.cases.all,match.file.cases)
tryCatch(GDCdownload(query, method = "api", files.per.chunk = 20),
error = function(e) GDCdownload(query, method = "client"))
}
# This will create a map between idat file name, cases (barcode) and project
readr::write_tsv(match.file.cases.all, path = "idat_filename_case.txt")
# code to move all files to local folder
for(file in dir(".",pattern = ".idat", recursive = T)){
TCGAbiolinks::move(file,basename(file))
}DNA methylation: aligned against hg19
Protein expression
Gene expression: aligned against hg19
# Aligned against Hg19
query.exp.hg19 <- GDCquery(project = "TCGA-GBM",
data.category = "Gene expression",
data.type = "Gene expression quantification",
platform = "Illumina HiSeq",
file.type = "normalized_results",
experimental.strategy = "RNA-Seq",
barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"),
legacy = TRUE)
GDCdownload(query.exp.hg19)
data <- GDCprepare(query.exp.hg19)Harmonized database
Copy Number
GISTIC2
Gene expression: aligned against hg38
For more examples, please check: http://rpubs.com/tiagochst/TCGAbiolinks_RNA-seq_new_projects
# mRNA pipeline: https://gdc-docs.nci.nih.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/
query.exp.hg38 <- GDCquery(project = "TCGA-GBM",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - FPKM-UQ",
barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"))
GDCdownload(query.exp.hg38)
expdat <- GDCprepare(query = query.exp.hg38,
save = TRUE,
save.filename = "exp.rda")miRNA
library(TCGAbiolinks)
query.mirna <- GDCquery(project = "TARGET-AML",
experimental.strategy = "miRNA-Seq",
data.category = "Transcriptome Profiling",
barcode = c("TARGET-20-PATDNN","TARGET-20-PAPUNR"),
data.type = "miRNA Expression Quantification")
GDCdownload(query.mirna)
mirna <- GDCprepare(query = query.mirna,
save = TRUE,
save.filename = "mirna.rda")
query.isoform <- GDCquery(project = "TARGET-AML",
experimental.strategy = "miRNA-Seq",
data.category = "Transcriptome Profiling",
barcode = c("TARGET-20-PATDNN","TARGET-20-PAPUNR"),
data.type = "Isoform Expression Quantification")
GDCdownload(query.isoform)
isoform <- GDCprepare(query = query.isoform,
save = TRUE,
save.filename = "mirna-isoform.rda")DNA methylation: aligned against hg38
#--------------------------------------
# DNA methylation data
#--------------------------------------
# DNA methylation aligned to hg38
query_met.hg38 <- GDCquery(project= "TCGA-LGG",
data.category = "DNA Methylation",
platform = "Illumina Human Methylation 450",
barcode = c("TCGA-HT-8111-01A-11D-2399-05","TCGA-HT-A5R5-01A-11D-A28N-05"))
GDCdownload(query_met.hg38)
data.hg38 <- GDCprepare(query_met.hg38)DNA methylation IDAT
# Using sesame http://bioconductor.org/packages/sesame/
# Please cite 10.1093/nar/gky691 and doi: 10.1093/nar/gkt090.
library(TCGAbiolinks)
proj <- "TCGA-ACC"
query <- GDCquery(project = proj,
data.category = "Raw microarray data",
data.type = "Raw intensities",
experimental.strategy = "Methylation array",
legacy = TRUE,
barcode = c("TCGA-OR-A5JT","CGA-OR-A5LG","TCGA-OR-A5JX"),
file.type = ".idat",
platform = "Illumina Human Methylation 450")
tryCatch(GDCdownload(query, method = "api", files.per.chunk = 20),
error = function(e) GDCdownload(query, method = "client"))
betas <- GDCprepare(query)