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

## # A tibble: 6 × 11
##   id           lost_to_followup disease_type days_to_lost_to_foll…¹ submitter_id
##   <chr>        <lgl>            <chr>        <lgl>                  <chr>       
## 1 58771370-50… NA               Myeloid Leu… NA                     TARGET-20-P…
## 2 28da5b52-29… NA               Myeloid Leu… NA                     TARGET-20-P…
## 3 28dae019-4d… NA               Myeloid Leu… NA                     TARGET-20-D…
## 4 28f7e68c-e0… NA               Myeloid Leu… NA                     TARGET-20-K…
## 5 28ffdfb9-f0… NA               Myeloid Leu… NA                     TARGET-20-P…
## 6 29312892-07… NA               Myeloid Leu… NA                     TARGET-20-P…
## # ℹ abbreviated name: ¹​days_to_lost_to_followup
## # ℹ 6 more variables: created_datetime <chr>, primary_site <chr>,
## #   updated_datetime <chr>, case_id <chr>, index_date <chr>, state <chr>

head(clindat[["diagnoses"]])

## # A tibble: 6 × 84
##   case_id       irs_stage iss_stage ajcc_pathologic_stage ann_arbor_clinical_s…¹
##   <chr>         <lgl>     <lgl>     <lgl>                 <lgl>                 
## 1 58771370-508… NA        NA        NA                    NA                    
## 2 28da5b52-29e… NA        NA        NA                    NA                    
## 3 28ffdfb9-f05… NA        NA        NA                    NA                    
## 4 29312892-078… NA        NA        NA                    NA                    
## 5 58aa0126-221… NA        NA        NA                    NA                    
## 6 29b8f9de-3ab… NA        NA        NA                    NA                    
## # ℹ abbreviated name: ¹​ann_arbor_clinical_stage
## # ℹ 79 more variables: created_datetime <dttm>, enneking_msts_stage <lgl>,
## #   inrg_stage <lgl>, enneking_msts_metastasis <lgl>,
## #   tissue_or_organ_of_origin <chr>, age_at_diagnosis <int>,
## #   esophageal_columnar_dysplasia_degree <lgl>, cog_liver_stage <lgl>,
## #   child_pugh_classification <lgl>, metastasis_at_diagnosis_site <lgl>,
## #   state <chr>, prior_treatment <lgl>, …

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)

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().
• 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"

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] 86

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 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] 10

presults = pquery |> results_all()
length(ids(presults))

## [1] 86

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

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] "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.

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

##    doc_count                           key
## 1     181726    annotated_somatic_mutation
## 2     171451                 aligned_reads
## 3     112687       simple_somatic_mutation
## 4     108323          structural_variation
## 5      73585           copy_number_segment
## 6      50608               gene_expression
## 7      39593   aggregated_somatic_mutation
## 8      39002          copy_number_estimate
## 9      35978              mirna_expression
## 10     33476                   slide_image
## 11     33146      masked_methylation_array
## 12     27344        biospecimen_supplement
## 13     24236    submitted_genotyping_array
## 14     23135     simple_germline_variation
## 15     21589       masked_somatic_mutation
## 16     16573        methylation_beta_value
## 17     14277           clinical_supplement
## 18     11324              pathology_report
## 19      7906            protein_expression
## 20      1426    submitted_expression_array
## 21       132 secondary_expression_analysis

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] 1027517

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     : '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-APOLLO

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     : '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] 858

Asking 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] TRUE

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 × 17
##   id         data_format access file_name submitter_id data_category acl_1 type 
##   <chr>      <chr>       <chr>  <chr>     <chr>        <chr>         <chr> <chr>
## 1 e10dc899-… TSV         contr… 984c9d22… 8a237d1d-8a… Transcriptom… phs0… gene…
## 2 13af947c-… TSV         open   984c9d22… 962e286d-54… Transcriptom… open  gene…
## 3 6ab604cb-… TSV         open   b4388ff2… 5d2257c4-0a… Transcriptom… open  gene…
## 4 0da0d708-… TSV         open   95493f7f… dd17894d-43… Transcriptom… open  gene…
## 5 400dc107-… TSV         open   cfdb5780… 0dad997d-7a… Transcriptom… open  gene…
## 6 3fadac67-… TSV         open   25ff83ac… 657b1f82-36… Transcriptom… open  gene…
## # ℹ 9 more variables: platform <chr>, file_size <int>, created_datetime <chr>,
## #   md5sum <chr>, updated_datetime <chr>, file_id <chr>, data_type <chr>,
## #   state <chr>, experimental_strategy <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' &
                            access == "open" &
                            analysis.workflow_type == 'STAR - Counts')
manifest_df = qfiles |> manifest()
nrow(manifest_df)

## [1] 429

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.

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.

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

BAM slicing

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-S2

library(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))

How many cases are included in all TARGET projects?

cases() |> filter(~ project.program.name=='TARGET') |> count()

## [1] 6543

How many cases are included in all TCGA projects?

cases() |> filter(~ project.program.name=='TCGA') |> count()

## [1] 11428

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_type

##   doc_count                  key
## 1      1098        primary tumor
## 2      1023 blood derived normal
## 3       162  solid tissue normal
## 4         7           metastatic

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] 162

res = 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"

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"

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"

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"

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"

How many of each type of file are available?

res = files() |> facet('type') |> aggregations()
res$type

##    doc_count                           key
## 1     181726    annotated_somatic_mutation
## 2     171451                 aligned_reads
## 3     112687       simple_somatic_mutation
## 4     108323          structural_variation
## 5      73585           copy_number_segment
## 6      50608               gene_expression
## 7      39593   aggregated_somatic_mutation
## 8      39002          copy_number_estimate
## 9      35978              mirna_expression
## 10     33476                   slide_image
## 11     33146      masked_methylation_array
## 12     27344        biospecimen_supplement
## 13     24236    submitted_genotyping_array
## 14     23135     simple_germline_variation
## 15     21589       masked_somatic_mutation
## 16     16573        methylation_beta_value
## 17     14277           clinical_supplement
## 18     11324              pathology_report
## 19      7906            protein_expression
## 20      1426    submitted_expression_array
## 21       132 secondary_expression_analysis

ggplot(res$type,aes(x = key,y = doc_count)) + geom_bar(stat='identity') +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

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

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)