The GenomicDataCommons Package
Sunday, October 13, 2024
Abstract
The National Cancer Institute (NCI) has established the Genomic Data Commons (GDC). The GDC provides the cancer research community with an open and unified repository for sharing and accessing data across numerous cancer studies and projects via a high-performance data transfer and query infrastructure. The GenomicDataCommons Bioconductor package provides basic infrastructure for querying, accessing, and mining genomic datasets available from the GDC. We expect that the Bioconductor developer and the larger bioinformatics communities will build on the GenomicDataCommons package to add higher-level functionality and expose cancer genomics data to the plethora of state-of-the-art bioinformatics methods available in Bioconductor.
1 What is the GDC?
From the Genomic Data Commons (GDC) website:
The National Cancer Institute’s (NCI’s) Genomic Data Commons (GDC) is a data sharing platform that promotes precision medicine in oncology. It is not just a database or a tool; it is an expandable knowledge network supporting the import and standardization of genomic and clinical data from cancer research programs. The GDC contains NCI-generated data from some of the largest and most comprehensive cancer genomic datasets, including The Cancer Genome Atlas (TCGA) and Therapeutically Applicable Research to Generate Effective Therapies (TARGET). For the first time, these datasets have been harmonized using a common set of bioinformatics pipelines, so that the data can be directly compared. As a growing knowledge system for cancer, the GDC also enables researchers to submit data, and harmonizes these data for import into the GDC. As more researchers add clinical and genomic data to the GDC, it will become an even more powerful tool for making discoveries about the molecular basis of cancer that may lead to better care for patients.
The data model for the GDC is complex, but it worth a quick overview and a graphical representation is included here.
The data model is encoded as a so-called property graph. Nodes represent entities such as Projects, Cases, Diagnoses, Files (various kinds), and Annotations. The relationships between these entities are maintained as edges. Both nodes and edges may have Properties that supply instance details.
The GDC API exposes these nodes and edges in a somewhat simplified set of RESTful endpoints.
2 Quickstart
This quickstart section is just meant to show basic functionality. More details of functionality are included further on in this vignette and in function-specific help.
This software is available at Bioconductor.org and can be downloaded via
BiocManager::install.
To report bugs or problems, either
submit a new issue
or submit a bug.report(package='GenomicDataCommons') from within R (which
will redirect you to the new issue on GitHub).
2.1 Installation
Installation can be achieved via Bioconductor’s BiocManager package.
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install('GenomicDataCommons')library(GenomicDataCommons)2.2 Check connectivity and status
The GenomicDataCommons package relies on having network
connectivity. In addition, the NCI GDC API must also be operational
and not under maintenance. Checking status can be used to check this
connectivity and functionality.
GenomicDataCommons::status()## $commit
## [1] "7c8ffd0436bb0bb4dafed2d191586309ba6618bf"
##
## $data_release
## [1] "Data Release 41.0 - August 28, 2024"
##
## $data_release_version
## $data_release_version$major
## [1] 41
##
## $data_release_version$minor
## [1] 0
##
## $data_release_version$release_date
## [1] "2024-08-28"
##
##
## $status
## [1] "OK"
##
## $tag
## [1] "7.5.1"
##
## $version
## [1] 1And to check the status in code:
stopifnot(GenomicDataCommons::status()$status=="OK")2.3 Find data
The following code builds a manifest that can be used to guide the
download of raw data. Here, filtering finds gene expression files
quantified as raw counts using STAR from ovarian cancer patients.
ge_manifest <- files() %>%
filter( cases.project.project_id == 'TCGA-OV') %>%
filter( type == 'gene_expression' ) %>%
filter( analysis.workflow_type == 'STAR - Counts') %>%
manifest()
head(ge_manifest)2.4 Download data
After the 858 gene expression files
specified in the query above. Using multiple processes to do the download very
significantly speeds up the transfer in many cases. On a standard 1Gb
connection, the following completes in about 30 seconds. The first time the
data are downloaded, R will ask to create a cache directory (see ?gdc_cache
for details of setting and interacting with the cache). Resulting
downloaded files will be stored in the cache directory. Future access to
the same files will be directly from the cache, alleviating multiple downloads.
fnames <- lapply(ge_manifest$id[1:20], gdcdata)If the download had included controlled-access data, the download above would
have needed to include a token. Details are available in
the authentication section below.
2.5 Metadata queries
2.5.1 Clinical data
Accessing clinical data is a very common task. Given a set of case_ids,
the gdc_clinical() function will return a list of four tibbles.
- demographic
- diagnoses
- exposures
- main
case_ids = cases() %>% results(size=10) %>% ids()
clindat = gdc_clinical(case_ids)
names(clindat)## [1] "demographic" "diagnoses" "exposures" "follow_ups" "main"head(clindat[["main"]])head(clindat[["diagnoses"]])2.5.2 General metadata queries
The GenomicDataCommons package can access the significant
clinical, demographic, biospecimen, and annotation information
contained in the NCI GDC. The gdc_clinical() function will often
be all that is needed, but the API and GenomicDataCommons package
make much flexibility if fine-tuning is required.
expands = c("diagnoses","annotations",
"demographic","exposures")
clinResults = cases() %>%
GenomicDataCommons::select(NULL) %>%
GenomicDataCommons::expand(expands) %>%
results(size=50)
str(clinResults[[1]],list.len=6)## chr [1:50] "58771370-5082-485e-ac68-13edfbd9ef0c" ...# or listviewer::jsonedit(clinResults)3 Basic design
This package design is meant to have some similarities to the “hadleyverse” approach of dplyr. Roughly, the functionality for finding and accessing files and metadata can be divided into:
- Simple query constructors based on GDC API endpoints.
- A set of verbs that when applied, adjust filtering, field selection, and faceting (fields for aggregation) and result in a new query object (an endomorphism)
- A set of verbs that take a query and return results from the GDC
In addition, there are exhiliary functions for asking the GDC API for information about available and default fields, slicing BAM files, and downloading actual data files. Here is an overview of functionality1 See individual function and methods documentation for specific details..
- Creating a query
projects()cases()files()annotations()
- Manipulating a query
filter()facet()select()
- Introspection on the GDC API fields
mapping()available_fields()default_fields()grep_fields()available_values()available_expand()
- Executing an API call to retrieve query results
results()count()response()
- Raw data file downloads
gdcdata()transfer()gdc_client()
- Summarizing and aggregating field values (faceting)
aggregations()
- Authentication
gdc_token()
- BAM file slicing
slicing()
4 Usage
There are two main classes of operations when working with the NCI GDC.
- Querying metadata and finding data files (e.g., finding all gene expression quantifications data files for all colon cancer patients).
- Transferring raw or processed data from the GDC to another computer (e.g., downloading raw or processed data)
Both classes of operation are reviewed in detail in the following sections.
4.1 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).
4.1.1 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 (seedefault_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 4
## $ fields : chr [1:10] "dbgap_accession_number" "disease_type" "intended_release_date" "name" ...
## $ filters: NULL
## $ facets : NULL
## $ expand : NULL
## - attr(*, "class")= chr [1:3] "gdc_projects" "GDCQuery" "list"4.1.2 Retrieving results
[ GDC pagination 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] 86The 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 9
## $ id : chr [1:10] "TARGET-AML" "MATCH-Z1I" "HCMI-CMDC" "MATCH-W" ...
## $ primary_site :List of 10
## ..$ TARGET-AML: chr [1:2] "Unknown" "Hematopoietic and reticuloendothelial systems"
## ..$ MATCH-Z1I : chr [1:12] "Bronchus and lung" "Gallbladder" "Pancreas" "Unknown" ...
## ..$ HCMI-CMDC : chr [1:24] "Breast" "Rectum" "Nasal cavity and middle ear" "Bronchus and lung" ...
## ..$ MATCH-W : chr [1:18] "Breast" "Renal pelvis" "Corpus uteri" "Bladder" ...
## ..$ MATCH-Z1D : chr [1:15] "Breast" "Corpus uteri" "Bones, joints and articular cartilage of other and unspecified sites" "Prostate gland" ...
## ..$ MATCH-Z1A : chr [1:14] "Corpus uteri" "Bladder" "Rectum" "Ovary" ...
## ..$ MATCH-U : chr [1:11] "Meninges" "Ovary" "Liver and intrahepatic bile ducts" "Kidney" ...
## ..$ MATCH-Q : chr [1:13] "Corpus uteri" "Rectum" "Ovary" "Parotid gland" ...
## ..$ TCGA-PCPG : chr [1:7] "Other endocrine glands and related structures" "Heart, mediastinum, and pleura" "Connective, subcutaneous and other soft tissues" "Spinal cord, cranial nerves, and other parts of central nervous system" ...
## ..$ MATCH-H : chr [1:11] "Ovary" "Liver and intrahepatic bile ducts" "Unknown" "Anus and anal canal" ...
## $ dbgap_accession_number: chr [1:10] "phs000465" "phs002058" NA "phs001948" ...
## $ project_id : chr [1:10] "TARGET-AML" "MATCH-Z1I" "HCMI-CMDC" "MATCH-W" ...
## $ disease_type :List of 10
## ..$ TARGET-AML: chr [1:2] "Not Applicable" "Myeloid Leukemias"
## ..$ MATCH-Z1I : chr [1:6] "Squamous Cell Neoplasms" "Epithelial Neoplasms, NOS" "Nevi and Melanomas" "Ductal and Lobular Neoplasms" ...
## ..$ HCMI-CMDC : chr [1:16] "Cystic, Mucinous and Serous Neoplasms" "Ductal and Lobular Neoplasms" "Adenomas and Adenocarcinomas" "Complex Mixed and Stromal Neoplasms" ...
## ..$ MATCH-W : chr [1:8] "Cystic, Mucinous and Serous Neoplasms" "Ductal and Lobular Neoplasms" "Adenomas and Adenocarcinomas" "Neoplasms, NOS" ...
## ..$ MATCH-Z1D : chr [1:8] "Cystic, Mucinous and Serous Neoplasms" "Ductal and Lobular Neoplasms" "Adenomas and Adenocarcinomas" "Complex Mixed and Stromal Neoplasms" ...
## ..$ MATCH-Z1A : chr [1:6] "Adenomas and Adenocarcinomas" "Neoplasms, NOS" "Mesothelial Neoplasms" "Squamous Cell Neoplasms" ...
## ..$ MATCH-U : chr [1:7] "Adenomas and Adenocarcinomas" "Neoplasms, NOS" "Nerve Sheath Tumors" "Mesothelial Neoplasms" ...
## ..$ MATCH-Q : chr [1:5] "Cystic, Mucinous and Serous Neoplasms" "Adenomas and Adenocarcinomas" "Neoplasms, NOS" "Squamous Cell Neoplasms" ...
## ..$ TCGA-PCPG : chr "Paragangliomas and Glomus Tumors"
## ..$ MATCH-H : chr [1:4] "Epithelial Neoplasms, NOS" "Adenomas and Adenocarcinomas" "Gliomas" "Neoplasms, NOS"
## $ name : chr [1:10] "Acute Myeloid Leukemia" "Genomic Characterization CS-MATCH-0007 Arm Z1I" "NCI Cancer Model Development for the Human Cancer Model Initiative" "Genomic Characterization CS-MATCH-0007 Arm W" ...
## $ releasable : logi [1:10] TRUE FALSE TRUE FALSE FALSE FALSE ...
## $ state : chr [1:10] "open" "open" "open" "open" ...
## $ released : logi [1:10] TRUE TRUE TRUE TRUE TRUE TRUE ...
## - 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] 10presults = pquery %>% results_all()
length(ids(presults))## [1] 86# includes all records
length(ids(presults)) == count(pquery)## [1] TRUEExtracting 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.
4.1.3 Fields and Values
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] "average_base_quality" "average_insert_size"
## [5] "average_read_length" "cancer_dna_fraction"
## [7] "channel" "chip_id"
## [9] "chip_position" "contamination"
## [11] "contamination_error" "created_datetime"
## [13] "data_category" "data_format"
## [15] "data_type" "error_type"
## [17] "experimental_strategy" "file_autocomplete"
## [19] "file_id" "file_name"
## [21] "file_size" "genome_doubling"
## [23] "imaging_date" "magnification"
## [25] "md5sum" "mean_coverage"
## [27] "msi_score" "msi_status"
## [29] "pairs_on_diff_chr" "plate_name"
## [31] "plate_well" "platform"
## [33] "proc_internal" "proportion_base_mismatch"
## [35] "proportion_coverage_10X" "proportion_coverage_10x"
## [37] "proportion_coverage_30X" "proportion_coverage_30x"
## [39] "proportion_reads_duplicated" "proportion_reads_mapped"
## [41] "proportion_targets_no_coverage" "read_pair_number"
## [43] "revision" "stain_type"
## [45] "state" "state_comment"
## [47] "subclonal_genome_fraction" "submitter_id"
## [49] "tags" "total_reads"
## [51] "tumor_ploidy" "tumor_purity"
## [53] "type" "updated_datetime"
## [55] "wgs_coverage"# The number of fields available for files endpoint
length(available_fields('files'))## [1] 1230# 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] "consent_type" "created_datetime"
## [7] "days_to_consent" "days_to_lost_to_followup"
## [9] "diagnosis_ids" "disease_type"
## [11] "index_date" "lost_to_followup"
## [13] "portion_ids" "primary_site"
## [15] "sample_ids" "slide_ids"
## [17] "state" "submitter_aliquot_ids"
## [19] "submitter_analyte_ids" "submitter_diagnosis_ids"
## [21] "submitter_id" "submitter_portion_ids"
## [23] "submitter_sample_ids" "submitter_slide_ids"
## [25] "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.
See the appropriate help pages for details.
4.1.4 Facets and aggregation
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$typeUsing aggregations() is an also easy way to learn the contents of individual
fields and forms the basis for faceted search pages.
4.1.5 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] 1027517To 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 : 'scalar' 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) %>%
head()## [1] "analysis.input_files.proc_internal"
## [2] "analysis.input_files.proportion_base_mismatch"
## [3] "analysis.input_files.proportion_coverage_10X"
## [4] "analysis.input_files.proportion_coverage_10x"
## [5] "analysis.input_files.proportion_coverage_30X"
## [6] "analysis.input_files.proportion_coverage_30x"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() %>%
head()## $cases.project.project_id
## doc_count key
## 1 54096 FM-AD
## 2 61173 TCGA-BRCA
## 3 72791 CPTAC-3
## 4 51339 TARGET-AML
## 5 48755 MP2PRT-ALL
## 6 32026 TCGA-LUAD
## 7 28417 TCGA-UCEC
## 8 29489 TCGA-HNSC
## 9 29021 TCGA-OV
## 10 29352 TCGA-KIRC
## 11 28245 TCGA-THCA
## 12 29334 TCGA-LUSC
## 13 28573 TCGA-LGG
## 14 27793 TCGA-PRAD
## 15 24039 TCGA-GBM
## 16 25726 TCGA-COAD
## 17 25210 TCGA-STAD
## 18 24540 TCGA-SKCM
## 19 23394 TCGA-BLCA
## 20 20820 TCGA-LIHC
## 21 27014 MMRF-COMMPASS
## 22 17584 TARGET-ALL-P2
## 23 16418 TCGA-CESC
## 24 16431 TCGA-KIRP
## 25 20761 HCMI-CMDC
## 26 14184 TCGA-SARC
## 27 16794 BEATAML1.0-COHORT
## 28 13367 TARGET-NBL
## 29 10654 CGCI-BLGSP
## 30 14954 REBC-THYR
## 31 10894 TCGA-PAAD
## 32 10234 TCGA-ESCA
## 33 10457 TCGA-PCPG
## 34 10290 TCGA-TGCT
## 35 8887 TCGA-READ
## 36 8839 TCGA-LAML
## 37 9244 CPTAC-2
## 38 7098 TCGA-THYM
## 39 6415 TARGET-WT
## 40 5856 CGCI-HTMCP-CC
## 41 5134 TCGA-ACC
## 42 4759 TCGA-KICH
## 43 4874 TCGA-MESO
## 44 4549 TCGA-UVM
## 45 5454 CMI-MBC
## 46 4114 TARGET-OS
## 47 5286 NCICCR-DLBCL
## 48 3323 TCGA-UCS
## 49 3736 TARGET-ALL-P3
## 50 2726 TCGA-DLBC
## 51 2738 TCGA-CHOL
## 52 1712 CGCI-HTMCP-DLBCL
## 53 1826 EXCEPTIONAL_RESPONDERS-ER
## 54 1796 CDDP_EAGLE-1
## 55 1628 OHSU-CNL
## 56 1571 MP2PRT-WT
## 57 1455 APOLLO-LUAD
## 58 1419 MATCH-I
## 59 1036 TARGET-RT
## 60 1305 CMI-MPC
## 61 1093 WCDT-MCRPC
## 62 1091 MATCH-W
## 63 1090 MATCH-Z1A
## 64 896 CGCI-HTMCP-LC
## 65 980 MATCH-S1
## 66 896 ORGANOID-PANCREATIC
## 67 891 MATCH-Z1D
## 68 852 MATCH-Q
## 69 806 CMI-ASC
## 70 810 MATCH-B
## 71 783 MATCH-Y
## 72 700 MATCH-R
## 73 694 MATCH-Z1B
## 74 671 MATCH-P
## 75 660 MATCH-Z1I
## 76 553 CTSP-DLBCL1
## 77 547 BEATAML1.0-CRENOLANIB
## 78 545 MATCH-U
## 79 510 MATCH-N
## 80 509 MATCH-H
## 81 339 TRIO-CRU
## 82 263 MATCH-C1
## 83 185 TARGET-CCSK
## 84 101 TARGET-ALL-P1
## 85 61 MATCH-S2
## 86 42 VAREPOP-APOLLOWe 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 : 'scalar' chr "and"
## $ content:List of 2
## ..$ :List of 2
## .. ..$ op : 'scalar' chr "="
## .. ..$ content:List of 2
## .. .. ..$ field: chr "cases.project.project_id"
## .. .. ..$ value: chr "TCGA-OV"
## ..$ :List of 2
## .. ..$ op : 'scalar' chr "="
## .. ..$ content:List of 2
## .. .. ..$ field: chr "type"
## .. .. ..$ value: chr "gene_expression"qfiles %>% count()## [1] 858Asking for a count() of results given these new filter criteria gives r qfiles %>% count() results. Filters can be chained (or nested) to
accomplish the same effect as multiple & conditionals. The count()
below is equivalent to the & filtering done above.
qfiles2 = files() %>%
filter( cases.project.project_id == 'TCGA-OV') %>%
filter( type == 'gene_expression')
qfiles2 %>% count()## [1] 858(qfiles %>% count()) == (qfiles2 %>% count()) #TRUE## [1] TRUEGenerating a manifest for bulk downloads is as simple as asking for the manifest from the current query.
manifest_df = qfiles %>% manifest()
head(manifest_df)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' &
access == "open" &
analysis.workflow_type == 'STAR - Counts')
manifest_df = qfiles %>% manifest()
nrow(manifest_df)## [1] 429The 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.
4.2 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())4.3 Datafile access and download
4.3.1 Data downloads via the GDC API
The gdcdata function takes a character vector of one or more file
ids. A simple way of producing such a vector is to produce a
manifest data frame and then pass in the first column, which will
contain file ids.
fnames = gdcdata(manifest_df$id[1:2],progress=FALSE)Note that for controlled-access data, a
GDC authentication token is required. Using the
BiocParallel package may be useful for downloading in parallel,
particularly for large numbers of smallish files.
4.3.2 Bulk downloads
The bulk download functionality is only efficient (as of v1.2.0 of the GDC Data Transfer Tool) for relatively large files, so use this approach only when transferring BAM files or larger VCF files, for example. Otherwise, consider using the approach shown above, perhaps in parallel.
# Requires gcd_client command-line utility to be isntalled
# separately.
fnames = gdcdata(manifest_df$id[3:10], access_method = 'client')4.3.3 BAM slicing
5 Use Cases
5.1 Cases
5.1.1 How many cases are there per project_id?
res = cases() %>% facet("project.project_id") %>% aggregations()
head(res)## $project.project_id
## doc_count key
## 1 18004 FM-AD
## 2 2492 TARGET-AML
## 3 1587 TARGET-ALL-P2
## 4 1510 MP2PRT-ALL
## 5 1345 CPTAC-3
## 6 1132 TARGET-NBL
## 7 1098 TCGA-BRCA
## 8 995 MMRF-COMMPASS
## 9 826 BEATAML1.0-COHORT
## 10 652 TARGET-WT
## 11 617 TCGA-GBM
## 12 608 TCGA-OV
## 13 585 TCGA-LUAD
## 14 560 TCGA-UCEC
## 15 537 TCGA-KIRC
## 16 528 TCGA-HNSC
## 17 516 TCGA-LGG
## 18 507 TCGA-THCA
## 19 504 TCGA-LUSC
## 20 500 TCGA-PRAD
## 21 489 NCICCR-DLBCL
## 22 470 TCGA-SKCM
## 23 461 TCGA-COAD
## 24 449 REBC-THYR
## 25 443 TCGA-STAD
## 26 412 TCGA-BLCA
## 27 383 TARGET-OS
## 28 377 TCGA-LIHC
## 29 342 CPTAC-2
## 30 339 TRIO-CRU
## 31 324 CGCI-BLGSP
## 32 307 TCGA-CESC
## 33 291 TCGA-KIRP
## 34 278 HCMI-CMDC
## 35 263 TCGA-TGCT
## 36 261 TCGA-SARC
## 37 212 CGCI-HTMCP-CC
## 38 200 CMI-MBC
## 39 200 TCGA-LAML
## 40 191 TARGET-ALL-P3
## 41 185 TCGA-ESCA
## 42 185 TCGA-PAAD
## 43 179 TCGA-PCPG
## 44 176 OHSU-CNL
## 45 172 TCGA-READ
## 46 124 TCGA-THYM
## 47 113 TCGA-KICH
## 48 101 WCDT-MCRPC
## 49 92 TCGA-ACC
## 50 87 APOLLO-LUAD
## 51 87 TCGA-MESO
## 52 84 EXCEPTIONAL_RESPONDERS-ER
## 53 80 TCGA-UVM
## 54 70 CGCI-HTMCP-DLBCL
## 55 70 ORGANOID-PANCREATIC
## 56 69 TARGET-RT
## 57 63 CMI-MPC
## 58 60 MATCH-I
## 59 58 TCGA-DLBC
## 60 57 TCGA-UCS
## 61 56 BEATAML1.0-CRENOLANIB
## 62 52 MP2PRT-WT
## 63 51 TCGA-CHOL
## 64 50 CDDP_EAGLE-1
## 65 45 CTSP-DLBCL1
## 66 45 MATCH-W
## 67 45 MATCH-Z1A
## 68 41 MATCH-S1
## 69 39 CGCI-HTMCP-LC
## 70 36 CMI-ASC
## 71 36 MATCH-Z1D
## 72 35 MATCH-Q
## 73 33 MATCH-B
## 74 31 MATCH-Y
## 75 29 MATCH-Z1B
## 76 28 MATCH-P
## 77 28 MATCH-R
## 78 26 MATCH-Z1I
## 79 24 TARGET-ALL-P1
## 80 23 MATCH-U
## 81 21 MATCH-H
## 82 21 MATCH-N
## 83 13 TARGET-CCSK
## 84 11 MATCH-C1
## 85 7 VAREPOP-APOLLO
## 86 3 MATCH-S2library(ggplot2)
ggplot(res$project.project_id,aes(x = key, y = doc_count)) +
geom_bar(stat='identity') +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
5.1.2 How many cases are included in all TARGET projects?
cases() %>% filter(~ project.program.name=='TARGET') %>% count()## [1] 65435.1.3 How many cases are included in all TCGA projects?
cases() %>% filter(~ project.program.name=='TCGA') %>% count()## [1] 114285.1.4 What is the breakdown of sample types in TCGA-BRCA?
# The need to do the "&" here is a requirement of the
# current version of the GDC API. I have filed a feature
# request to remove this requirement.
resp = cases() %>% filter(~ project.project_id=='TCGA-BRCA' &
project.project_id=='TCGA-BRCA' ) %>%
facet('samples.sample_type') %>% aggregations()
resp$samples.sample_type5.1.5 Fetch all samples in TCGA-BRCA that use “Solid Tissue” as a normal.
# The need to do the "&" here is a requirement of the
# current version of the GDC API. I have filed a feature
# request to remove this requirement.
resp = cases() %>% filter(~ project.project_id=='TCGA-BRCA' &
samples.sample_type=='Solid Tissue Normal') %>%
GenomicDataCommons::select(c(default_fields(cases()),'samples.sample_type')) %>%
response_all()
count(resp)## [1] 162res = resp %>% results()
str(res[1],list.len=6)## List of 1
## $ id: chr [1:162] "3d676bba-154b-4d22-ab59-d4d4da051b94" "ac075bc0-1b59-4557-beea-541694faee03" "b2aac45b-2073-4c7a-adb9-769a4fdcc111" "b3264748-947a-43aa-b227-b294fbcc8447" ...head(ids(resp))## [1] "3d676bba-154b-4d22-ab59-d4d4da051b94"
## [2] "ac075bc0-1b59-4557-beea-541694faee03"
## [3] "b2aac45b-2073-4c7a-adb9-769a4fdcc111"
## [4] "b3264748-947a-43aa-b227-b294fbcc8447"
## [5] "b6b1dc9a-91f4-4b0a-afd5-62c9a90c0d5e"
## [6] "17c1d42c-cb84-4655-a4cd-b54bae17ecaf"5.1.6 Get all TCGA case ids that are female
cases() %>%
GenomicDataCommons::filter(~ project.program.name == 'TCGA' &
"cases.demographic.gender" %in% "female") %>%
GenomicDataCommons::results(size = 4) %>%
ids()## [1] "85a85a11-7200-4e96-97af-6ba26d680d59"
## [2] "7922df77-f09a-488c-a1be-58646ceb9b3e"
## [3] "8727855e-120a-4216-a803-8cc6cd1159be"
## [4] "82e96c6c-a88c-4e52-be56-7f24f6c7b835"5.1.7 Get all TCGA-COAD case ids that are NOT female
cases() %>%
GenomicDataCommons::filter(~ project.project_id == 'TCGA-COAD' &
"cases.demographic.gender" %exclude% "female") %>%
GenomicDataCommons::results(size = 4) %>%
ids()## [1] "733d8b6a-ca9d-4a69-8c9c-1f88733e8b68"
## [2] "ad440651-a2de-4bb1-90da-1e5e8975ab59"
## [3] "65bb7520-f055-43a8-b735-1152fa2c9e04"
## [4] "13eff2e5-e33a-485f-9ba4-8a7ccb3c7528"5.1.8 Get all TCGA cases that are missing gender
cases() %>%
GenomicDataCommons::filter(~ project.program.name == 'TCGA' &
missing("cases.demographic.gender")) %>%
GenomicDataCommons::results(size = 4) %>%
ids()## [1] "81d11171-5d9e-4950-98f8-b0ab0f2d7908"
## [2] "7f36eff2-aa69-4c4d-8101-8801c812a36b"
## [3] "dde3c31f-51cb-4236-9aa9-0c9eda6105eb"
## [4] "f88560c8-3475-47f2-9d48-4b8311bc1132"5.1.9 Get all TCGA cases that are NOT missing gender
cases() %>%
GenomicDataCommons::filter(~ project.program.name == 'TCGA' &
!missing("cases.demographic.gender")) %>%
GenomicDataCommons::results(size = 4) %>%
ids()## [1] "85a85a11-7200-4e96-97af-6ba26d680d59"
## [2] "7922df77-f09a-488c-a1be-58646ceb9b3e"
## [3] "8727855e-120a-4216-a803-8cc6cd1159be"
## [4] "82e96c6c-a88c-4e52-be56-7f24f6c7b835"5.2 Files
5.2.1 How many of each type of file are available?
res = files() %>% facet('type') %>% aggregations()
res$typeggplot(res$type,aes(x = key,y = doc_count)) + geom_bar(stat='identity') +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
5.2.2 Find gene-level RNA-seq quantification files for GBM
q = files() %>%
GenomicDataCommons::select(available_fields('files')) %>%
filter(~ cases.project.project_id=='TCGA-GBM' &
data_type=='Gene Expression Quantification')
q %>% facet('analysis.workflow_type') %>% aggregations()## list()# so need to add another filter
file_ids = q %>% filter(~ cases.project.project_id=='TCGA-GBM' &
data_type=='Gene Expression Quantification' &
analysis.workflow_type == 'STAR - Counts') %>%
GenomicDataCommons::select('file_id') %>%
response_all() %>%
ids()5.3 Slicing
5.3.1 Get all BAM file ids from TCGA-GBM
I need to figure out how to do slicing reproducibly in a testing environment and for vignette building.
q = files() %>%
GenomicDataCommons::select(available_fields('files')) %>%
filter(~ cases.project.project_id == 'TCGA-GBM' &
data_type == 'Aligned Reads' &
experimental_strategy == 'RNA-Seq' &
data_format == 'BAM')
file_ids = q %>% response_all() %>% ids()bamfile = slicing(file_ids[1],regions="chr12:6534405-6538375",token=gdc_token())
library(GenomicAlignments)
aligns = readGAlignments(bamfile)6 Troubleshooting
6.1 SSL connection errors
- Symptom: Trying to connect to the API results in:
Error in curl::curl_fetch_memory(url, handle = handle) :
SSL connect error- Possible solutions: The issue
is that the GDC supports only recent security Transport Layer Security (TLS),
so the only known fix is to upgrade the system
opensslto version 1.0.1 or later.- [Mac OS],
- [Ubuntu]
- [Centos/RHEL].
After upgrading
openssl, reinstall the Rcurlandhttrpackages.
7 sessionInfo()
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.5.1 GenomicDataCommons_1.28.2
## [3] magrittr_2.0.3 knitr_1.48
## [5] BiocStyle_2.32.1
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 sass_0.4.9 utf8_1.2.4
## [4] generics_0.1.3 tidyr_1.3.1 xml2_1.3.6
## [7] hms_1.1.3 digest_0.6.37 evaluate_1.0.1
## [10] grid_4.4.1 bookdown_0.40 fastmap_1.2.0
## [13] jsonlite_1.8.9 GenomeInfoDb_1.40.1 tinytex_0.53
## [16] BiocManager_1.30.25 httr_1.4.7 purrr_1.0.2
## [19] fansi_1.0.6 scales_1.3.0 UCSC.utils_1.0.0
## [22] jquerylib_0.1.4 cli_3.6.3 rlang_1.1.4
## [25] crayon_1.5.3 XVector_0.44.0 munsell_0.5.1
## [28] withr_3.0.1 cachem_1.1.0 yaml_2.3.10
## [31] tools_4.4.1 tzdb_0.4.0 dplyr_1.1.4
## [34] colorspace_2.1-1 GenomeInfoDbData_1.2.12 BiocGenerics_0.50.0
## [37] curl_5.2.3 vctrs_0.6.5 R6_2.5.1
## [40] magick_2.8.5 stats4_4.4.1 lifecycle_1.0.4
## [43] zlibbioc_1.50.0 S4Vectors_0.42.1 IRanges_2.38.1
## [46] pkgconfig_2.0.3 gtable_0.3.5 pillar_1.9.0
## [49] bslib_0.8.0 Rcpp_1.0.13 glue_1.8.0
## [52] highr_0.11 xfun_0.48 tibble_3.2.1
## [55] GenomicRanges_1.56.2 tidyselect_1.2.1 farver_2.1.2
## [58] htmltools_0.5.8.1 labeling_0.4.3 rmarkdown_2.28
## [61] readr_2.1.5 compiler_4.4.18 Developer notes
- The
S3object-oriented programming paradigm is used. - We have adopted a functional programming style with functions and methods that often take an “object” as the first argument. This style lends itself to pipeline-style programming.
- The GenomicDataCommons package uses the alternative request format (POST) to allow very large request bodies.