TCGAbiolinks has provided a few functions to download and prepare data from GDC for analysis. This section starts by explaning 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 downlaod 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 whicih might have a high probability of download failure. To solve that we created the files.per.chunk argument 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 funciton 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 unmofied 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

Related discussions: issue 91, issue 50

Arguments

`GDCdownload`

Argument	Description
query	A query for GDCquery function
token.file	Token file to download controled 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)

rowRanges(data)

## GRanges object with 100 ranges and 3 metadata columns:
##                seqnames                 ranges strand |      gene_id
##                   <Rle>              <IRanges>  <Rle> |  <character>
##           A1BG    chr19   [58856544, 58864865]      - |         A1BG
##            A2M    chr12   [ 9220260,  9268825]      - |          A2M
##           NAT1     chr8   [18027986, 18081198]      + |         NAT1
##           NAT2     chr8   [18248755, 18258728]      + |         NAT2
##   RP11-986E7.7    chr14   [95058395, 95090983]      + | RP11-986E7.7
##            ...      ...                    ...    ... .          ...
##         ADORA1     chr1 [203059782, 203136533]      + |       ADORA1
##        ADORA2A    chr22 [ 24813847,  24838328]      + |      ADORA2A
##        ADORA2B    chr17 [ 15848231,  15879060]      + |      ADORA2B
##         ADORA3     chr1 [112025970, 112106584]      - |       ADORA3
##          ADPRH     chr3 [119298115, 119308792]      + |        ADPRH
##                entrezgene ensembl_gene_id
##                 <numeric>     <character>
##           A1BG          1 ENSG00000121410
##            A2M          2 ENSG00000175899
##           NAT1          9 ENSG00000171428
##           NAT2         10 ENSG00000156006
##   RP11-986E7.7         12 ENSG00000273259
##            ...        ...             ...
##         ADORA1        134 ENSG00000163485
##        ADORA2A        135 ENSG00000128271
##        ADORA2B        136 ENSG00000170425
##         ADORA3        140 ENSG00000121933
##          ADPRH        141 ENSG00000144843
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

Search 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)

rowRanges(data)

## 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 seqlengths

`GDCprepare`: 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
Simple Nucleotide Variation
Raw Sequencing Data
Biospecimen
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",chunks.per.download = 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

query_meth.hg19 <- 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"), 
                            legacy = TRUE)
GDCdownload(query_meth.hg19)
data.hg19 <- GDCprepare(query_meth.hg19)

Protein expression

query <- GDCquery(project = "TCGA-GBM",
                  data.category = "Protein expression",
                  legacy = TRUE, 
                  barcode = c("TCGA-OX-A56R-01A-21-A44T-20","TCGA-08-0357-01A-21-1898-20"))
GDCdownload(query)
data <- GDCprepare(query, save = TRUE, 
                   save.filename = "gbmProteinExpression.rda",
                   remove.files.prepared = TRUE)

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

query <- GDCquery(project = "TCGA-ACC", 
                  data.category = "Copy Number Variation",
                  data.type = "Copy Number Segment",
                  barcode = c( "TCGA-OR-A5KU-01A-11D-A29H-01", "TCGA-OR-A5JK-01A-11D-A29H-01"))
GDCdownload(query)
data <- GDCprepare(query)

Gene expression: aligned against hg38

# 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")

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)

TCGAbiolinks: Downloading and preparing files for analysis

2 February 2018

Downloading and preparing data for analysis

Arguments

`GDCdownload`

`GDCprepare`

Search and download data from legacy database using GDC api method

Search and download data for two samples from database

`GDCprepare`: Outputs

Harmonized data

Legacy data

Examples

Legacy archive

DNA methylation: Get all TCGA IDAT files

DNA methylation: aligned against hg19

Protein expression

Gene expression: aligned against hg19

Harmonized database

Copy Number

Gene expression: aligned against hg38

DNA methylation: aligned against hg38

TCGAbiolinks: Downloading and preparing files for analysis

2 February 2018

Downloading and preparing data for analysis

Arguments

GDCdownload

GDCprepare

Search and download data from legacy database using GDC api method

Search and download data for two samples from database

GDCprepare: Outputs

Harmonized data

Legacy data

Examples

Legacy archive

DNA methylation: Get all TCGA IDAT files

DNA methylation: aligned against hg19

Protein expression

Gene expression: aligned against hg19

Harmonized database

Copy Number

Gene expression: aligned against hg38

DNA methylation: aligned against hg38

`GDCdownload`

`GDCprepare`

`GDCprepare`: Outputs