There are two main classes of operations when working with the NCI GDC.
Both classes of operation are reviewed in detail in the following sections.
Querying metadata
Vast amounts of metadata about cases (patients, basically), files, projects, and so-called annotations are available via the NCI GDC API. Typically, one will want to query metadata to either focus in on a set of files for download or transfer or to perform so-called aggregations (pivot-tables, facets, similar to the R table() functionality).
Querying metadata starts with creating a “blank” query. One will often then want to filter the query to limit results prior to retrieving results. The GenomicDataCommons package has helper functions for listing fields that are available for filtering.
In addition to fetching results, the GDC API allows faceting, or aggregating,, useful for compiling reports, generating dashboards, or building user interfaces to GDC data (see GDC web query interface for a non-R-based example).
Creating a query
A query of the GDC starts its life in R. Queries follow the four metadata endpoints available at the GDC. In particular, there are four convenience functions that each create GDCQuery objects (actually, specific subclasses of GDCQuery):
projects()cases()files()annotations()
pquery = projects()
The pquery object is now an object of (S3) class, GDCQuery (and gdc_projects and list). The object contains the following elements:
- fields: This is a character vector of the fields that will be returned when we retrieve data. If no fields are specified to, for example, the
projects() function, the default fields from the GDC are used (see default_fields())
- filters: This will contain results after calling the
filter() method and will be used to filter results on retrieval.
- facets: A character vector of field names that will be used for aggregating data in a call to
aggregations().
- archive: One of either “default” or “legacy”.
- token: A character(1) token from the GDC. See the authentication section for details, but note that, in general, the token is not necessary for metadata query and retrieval, only for actual data download.
Looking at the actual object (get used to using str()!), note that the query contains no results.
str(pquery)
## List of 5
## $ fields : chr [1:10] "dbgap_accession_number" "disease_type" "intended_release_date" "name" ...
## $ filters: NULL
## $ facets : NULL
## $ legacy : logi FALSE
## $ expand : NULL
## - attr(*, "class")= chr [1:3] "gdc_projects" "GDCQuery" "list"
Retrieving results
[ GDC pagination documentation ]
[ GDC sorting documentation ]
With a query object available, the next step is to retrieve results from the GDC. The GenomicDataCommons package. The most basic type of results we can get is a simple count() of records available that satisfy the filter criteria. Note that we have not set any filters, so a count() here will represent all the project records publicly available at the GDC in the “default” archive"
pcount = count(pquery)
# or
pcount = pquery %>% count()
pcount
## [1] 39
The results() method will fetch actual results.
presults = pquery %>% results()
These results are returned from the GDC in JSON format and converted into a (potentially nested) list in R. The str() method is useful for taking a quick glimpse of the data.
str(presults)
## List of 8
## $ dbgap_accession_number: chr [1:10] NA NA NA NA ...
## $ disease_type :List of 10
## ..$ TCGA-CHOL : chr "Cholangiocarcinoma"
## ..$ TCGA-OV : chr "Ovarian Serous Cystadenocarcinoma"
## ..$ TCGA-LIHC : chr "Liver Hepatocellular Carcinoma"
## ..$ TCGA-HNSC : chr "Head and Neck Squamous Cell Carcinoma"
## ..$ TCGA-COAD : chr "Colon Adenocarcinoma"
## ..$ TCGA-ACC : chr "Adrenocortical Carcinoma"
## ..$ TCGA-MESO : chr "Mesothelioma"
## ..$ TARGET-AML: chr "Acute Myeloid Leukemia"
## ..$ TARGET-NBL: chr "Neuroblastoma"
## ..$ TCGA-LGG : chr "Brain Lower Grade Glioma"
## $ released : logi [1:10] TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ state : chr [1:10] "legacy" "legacy" "legacy" "legacy" ...
## $ primary_site :List of 10
## ..$ TCGA-CHOL : chr "Bile Duct"
## ..$ TCGA-OV : chr "Ovary"
## ..$ TCGA-LIHC : chr "Liver"
## ..$ TCGA-HNSC : chr "Head and Neck"
## ..$ TCGA-COAD : chr "Colorectal"
## ..$ TCGA-ACC : chr "Adrenal Gland"
## ..$ TCGA-MESO : chr "Pleura"
## ..$ TARGET-AML: chr "Blood"
## ..$ TARGET-NBL: chr "Nervous System"
## ..$ TCGA-LGG : chr "Brain"
## $ project_id : chr [1:10] "TCGA-CHOL" "TCGA-OV" "TCGA-LIHC" "TCGA-HNSC" ...
## $ id : chr [1:10] "TCGA-CHOL" "TCGA-OV" "TCGA-LIHC" "TCGA-HNSC" ...
## $ name : chr [1:10] "Cholangiocarcinoma" "Ovarian Serous Cystadenocarcinoma" "Liver Hepatocellular Carcinoma" "Head and Neck Squamous Cell Carcinoma" ...
## - attr(*, "row.names")= int [1:10] 1 2 3 4 5 6 7 8 9 10
## - attr(*, "class")= chr [1:3] "GDCprojectsResults" "GDCResults" "list"
A default of only 10 records are returned. We can use the size and from arguments to results() to either page through results or to change the number of results. Finally, there is a convenience method, results_all() that will simply fetch all the available results given a query. Note that results_all() may take a long time and return HUGE result sets if not used carefully. Use of a combination of count() and results() to get a sense of the expected data size is probably warranted before calling results_all()
length(ids(presults))
## [1] 10
presults = pquery %>% results_all()
length(ids(presults))
## [1] 39
# includes all records
length(ids(presults)) == count(pquery)
## [1] TRUE
Extracting subsets of results or manipulating the results into a more conventional R data structure is not easily generalizable. However, the purrr, rlist, and data.tree packages are all potentially of interest for manipulating complex, nested list structures. For viewing the results in an interactive viewer, consider the listviewer package.
In the case of the projects entity, the default results (using default fields, that is) can be simplified easily with as.data.frame.
head(as.data.frame(presults))
## project_id dbgap_accession_number released state id
## 1 TARGET-AML phs000465 TRUE legacy TARGET-AML
## 2 TARGET-CCSK phs000466 TRUE legacy TARGET-CCSK
## 3 TARGET-NBL phs000467 TRUE legacy TARGET-NBL
## 4 TARGET-OS phs000468 TRUE legacy TARGET-OS
## 5 TARGET-RT phs000470 TRUE legacy TARGET-RT
## 6 TARGET-WT phs000471 TRUE legacy TARGET-WT
## name disease_type primary_site
## 1 Acute Myeloid Leukemia Acute My.... Blood
## 2 Clear Cell Sarcoma of the Kidney Clear Ce.... Kidney
## 3 Neuroblastoma Neurobla.... Nervous ....
## 4 Osteosarcoma Osteosarcoma Bone
## 5 Rhabdoid Tumor Rhabdoid.... Kidney
## 6 High-Risk Wilms Tumor High-Ris.... Kidney
Fields and Values
[ GDC fields documentation ]
Central to querying and retrieving data from the GDC is the ability to specify which fields to return, filtering by fields and values, and faceting or aggregating. The GenomicDataCommons package includes two simple functions, available_fields() and default_fields(). Each can operate on a character(1) endpoint name (“cases”, “files”, “annotations”, or “projects”) or a GDCQuery object.
default_fields('files')
## [1] "access" "acl"
## [3] "created_datetime" "data_category"
## [5] "data_format" "data_type"
## [7] "error_type" "experimental_strategy"
## [9] "file_autocomplete" "file_id"
## [11] "file_name" "file_size"
## [13] "file_state" "md5sum"
## [15] "origin" "platform"
## [17] "revision" "state"
## [19] "state_comment" "submitter_id"
## [21] "tags" "type"
## [23] "updated_datetime"
# The number of fields available for files endpoint
length(available_fields('files'))
## [1] 592
# The first few fields available for files endpoint
head(available_fields('files'))
## [1] "access" "acl"
## [3] "analysis.analysis_id" "analysis.analysis_type"
## [5] "analysis.created_datetime" "analysis.input_files.access"
The fields to be returned by a query can be specified following a similar paradigm to that of the dplyr package. The select() function is a verb that resets the fields slot of a GDCQuery; note that this is not quite analogous to the dplyr select() verb that limits from already-present fields. We completely replace the fields when using select() on a GDCQuery.
# Default fields here
qcases = cases()
qcases$fields
## [1] "aliquot_ids" "analyte_ids"
## [3] "case_autocomplete" "case_id"
## [5] "created_datetime" "days_to_index"
## [7] "days_to_lost_to_followup" "disease_type"
## [9] "index_date" "lost_to_followup"
## [11] "portion_ids" "primary_site"
## [13] "sample_ids" "slide_ids"
## [15] "state" "submitter_aliquot_ids"
## [17] "submitter_analyte_ids" "submitter_id"
## [19] "submitter_portion_ids" "submitter_sample_ids"
## [21] "submitter_slide_ids" "updated_datetime"
# set up query to use ALL available fields
# Note that checking of fields is done by select()
qcases = cases() %>% GenomicDataCommons::select(available_fields('cases'))
head(qcases$fields)
## [1] "case_id" "aliquot_ids"
## [3] "analyte_ids" "annotations.annotation_id"
## [5] "annotations.case_id" "annotations.case_submitter_id"
Finding fields of interest is such a common operation that the GenomicDataCommons includes the grep_fields() function and the field_picker() widget. See the appropriate help pages for details.
Facets and aggregation
[ GDC facet documentation ]
The GDC API offers a feature known as aggregation or faceting. By specifying one or more fields (of appropriate type), the GDC can return to us a count of the number of records matching each potential value. This is similar to the R table method. Multiple fields can be returned at once, but the GDC API does not have a cross-tabulation feature; all aggregations are only on one field at a time. Results of aggregation() calls come back as a list of data.frames (actually, tibbles).
# total number of files of a specific type
res = files() %>% facet(c('type','data_type')) %>% aggregations()
res$type
## key doc_count
## 1 aligned_reads 45988
## 2 annotated_somatic_mutation 45577
## 3 simple_somatic_mutation 45577
## 4 copy_number_segment 44752
## 5 gene_expression 34722
## 6 mirna_expression 22976
## 7 methylation_beta_value 12359
## 8 biospecimen_supplement 11328
## 9 clinical_supplement 11169
## 10 aggregated_somatic_mutation 144
## 11 masked_somatic_mutation 132
Using aggregations() is an also easy way to learn the contents of individual fields and forms the basis for faceted search pages.
Filtering
[ GDC filtering documentation ]
The GenomicDataCommons package uses a form of non-standard evaluation to specify R-like queries that are then translated into an R list. That R list is, upon calling a method that fetches results from the GDC API, translated into the appropriate JSON string. The R expression uses the formula interface as suggested by Hadley Wickham in his vignette on non-standard evaluation
It’s best to use a formula because a formula captures both the expression to evaluate and the environment where the evaluation occurs. This is important if the expression is a mixture of variables in a data frame and objects in the local environment [for example].
For the user, these details will not be too important except to note that a filter expression must begin with a “~”.
qfiles = files()
qfiles %>% count() # all files
## [1] 274724
To limit the file type, we can refer back to the section on faceting to see the possible values for the file field “type”. For example, to filter file results to only “gene_expression” files, we simply specify a filter.
qfiles = files() %>% filter(~ type == 'gene_expression')
# here is what the filter looks like after translation
str(get_filter(qfiles))
## List of 2
## $ op :Classes 'scalar', 'character' chr "="
## $ content:List of 2
## ..$ field: chr "type"
## ..$ value: chr "gene_expression"
What if we want to create a filter based on the project (‘TCGA-OVCA’, for example)? Well, we have a couple of possible ways to discover available fields. The first is based on base R functionality and some intuition.
grep('pro',available_fields('files'),value=TRUE)
## [1] "cases.diagnoses.progression_free_survival"
## [2] "cases.diagnoses.progression_free_survival_event"
## [3] "cases.diagnoses.progression_or_recurrence"
## [4] "cases.project.dbgap_accession_number"
## [5] "cases.project.disease_type"
## [6] "cases.project.intended_release_date"
## [7] "cases.project.name"
## [8] "cases.project.primary_site"
## [9] "cases.project.program.dbgap_accession_number"
## [10] "cases.project.program.name"
## [11] "cases.project.program.program_id"
## [12] "cases.project.project_id"
## [13] "cases.project.releasable"
## [14] "cases.project.released"
## [15] "cases.project.state"
## [16] "cases.samples.days_to_sample_procurement"
## [17] "cases.samples.method_of_sample_procurement"
## [18] "cases.samples.portions.slides.number_proliferating_cells"
## [19] "cases.tissue_source_site.project"
Interestingly, the project information is “nested” inside the case. We don’t need to know that detail other than to know that we now have a few potential guesses for where our information might be in the files records. We need to know where because we need to construct the appropriate filter.
files() %>% facet('cases.project.project_id') %>% aggregations()
## $cases.project.project_id
## key doc_count
## 1 TCGA-BRCA 27207
## 2 TCGA-LUAD 14804
## 3 TCGA-UCEC 13604
## 4 TCGA-LUSC 13124
## 5 TCGA-HNSC 12895
## 6 TCGA-LGG 12603
## 7 TCGA-THCA 12703
## 8 TCGA-OV 13054
## 9 TCGA-PRAD 12568
## 10 TCGA-COAD 11824
## 11 TCGA-SKCM 11265
## 12 TCGA-KIRC 12272
## 13 TCGA-STAD 10731
## 14 TCGA-BLCA 10193
## 15 TCGA-GBM 9657
## 16 TCGA-LIHC 9511
## 17 TCGA-CESC 7349
## 18 TCGA-KIRP 7368
## 19 TCGA-SARC 6282
## 20 TCGA-ESCA 4473
## 21 TCGA-PAAD 4433
## 22 TCGA-PCPG 4422
## 23 TCGA-READ 4012
## 24 TCGA-LAML 3954
## 25 TCGA-TGCT 3636
## 26 TARGET-NBL 2806
## 27 TCGA-THYM 2974
## 28 TARGET-AML 1873
## 29 TCGA-ACC 2108
## 30 TARGET-WT 1324
## 31 TCGA-MESO 2050
## 32 TCGA-UVM 1928
## 33 TCGA-KICH 1853
## 34 TCGA-UCS 1364
## 35 TCGA-CHOL 1157
## 36 TCGA-DLBC 1163
## 37 TARGET-OS 4
## 38 TARGET-RT 174
## 39 TARGET-CCSK 2
We note that cases.project.project_id looks like it is a good fit. We also note that TCGA-OV is the correct project_id, not TCGA-OVCA. Note that unlike with dplyr and friends, the filter() method here replaces the filter and does not build on any previous filters.
qfiles = files() %>% filter( ~ cases.project.project_id == 'TCGA-OV' & type == 'gene_expression')
str(get_filter(qfiles))
## List of 2
## $ op :Classes 'scalar', 'character' chr "and"
## $ content:List of 2
## ..$ :List of 2
## .. ..$ op :Classes 'scalar', 'character' chr "="
## .. ..$ content:List of 2
## .. .. ..$ field: chr "cases.project.project_id"
## .. .. ..$ value: chr "TCGA-OV"
## ..$ :List of 2
## .. ..$ op :Classes 'scalar', 'character' chr "="
## .. ..$ content:List of 2
## .. .. ..$ field: chr "type"
## .. .. ..$ value: chr "gene_expression"
qfiles %>% count()
## [1] 1137
Asking for a count() of results given these new filter criteria gives 1137 results. Generating a manifest for bulk downloads is as simple as asking for the manifest from the current query.
manifest_df = qfiles %>% manifest()
head(manifest_df)
## # A tibble: 6 x 5
## id
## <chr>
## 1 33976290-14b3-4c8a-9c91-8049b440137c
## 2 c663d732-facc-4eb7-9cb7-6ac5bd59da28
## 3 a39645aa-6277-4035-9a19-923dc294cc73
## 4 179f0ecf-d1b1-4968-8c78-15938b62a70e
## 5 f6dc8e84-b797-478b-a573-353dde598f94
## 6 b6e15e57-e69f-4a79-b98e-e786a9a11cad
## # ... with 4 more variables: filename <chr>, md5 <chr>, size <int>,
## # state <chr>
Note that we might still not be quite there. Looking at filenames, there are suspiciously named files that might include “FPKM”, “FPKM-UQ”, or “counts”. Another round of grep and available_fields, looking for “type” turned up that the field “analysis.workflow_type” has the appropriate filter criteria.
qfiles = files() %>% filter( ~ cases.project.project_id == 'TCGA-OV' &
type == 'gene_expression' &
analysis.workflow_type == 'HTSeq - Counts')
manifest_df = qfiles %>% manifest()
nrow(manifest_df)
## [1] 379
The GDC Data Transfer Tool can be used (from R, transfer() or from the command-line) to orchestrate high-performance, restartable transfers of all the files in the manifest. See the bulk downloads section for details.
Authentication
[ GDC authentication documentation ]
The GDC offers both “controlled-access” and “open” data. As of this writing, only data stored as files is “controlled-access”; that is, metadata accessible via the GDC is all “open” data and some files are “open” and some are “controlled-access”. Controlled-access data are only available after going through the process of obtaining access.
After controlled-access to one or more datasets has been granted, logging into the GDC web portal will allow you to access a GDC authentication token, which can be downloaded and then used to access available controlled-access data via the GenomicDataCommons package.
The GenomicDataCommons uses authentication tokens only for downloading data (see transfer and gdcdata documentation). The package includes a helper function, gdc_token, that looks for the token to be stored in one of three ways (resolved in this order):
- As a string stored in the environment variable,
GDC_TOKEN
- As a file, stored in the file named by the environment variable,
GDC_TOKEN_FILE
- In a file in the user home directory, called
.gdc_token
As a concrete example:
token = gdc_token()
transfer(...,token=token)
# or
transfer(...,token=get_token())
