SomaticCancerAlterations
Table of Contents
1 Motivation
Over the last years, large efforts have been taken to characterize the somatic landscape of cancers. Many of the conducted studies make their results publicly available, providing a valuable resource for investigating beyond the level of individual cohorts. The SomaticCancerAlterations package collects mutational data of several tumor types, currently focusing on the TCGA calls sets, and aims for a tight integration with and workflows. In the following, we will illustrate how to access this data and give examples for use cases.
2 Data Sets
The Cancer Genome Atlas (TCGA)1 is a consortium effort to analyze a variety of tumor types, including gene expression, methylation, copy number changes, and somatic mutations2. With the SomaticCancerAlterations package, we provide the callsets of somatic mutations for all publically available TCGA studies. Over time, more studies will be added, as they become available and unrestriced in their usage.
To get started, we get a list of all available data sets and access the metadata associated with each study.
all_datasets = scaListDatasets() print(all_datasets)
## [1] "gbm_tcga" "hnsc_tcga" "kirc_tcga" "luad_tcga" "lusc_tcga" "ov_tcga" ## [7] "skcm_tcga" "thca_tcga"
meta_data = scaMetadata() print(meta_data)
## Cancer_Type Center NCBI_Build Sequence_Source ## gbm_tcga GBM broad.mit.edu 37 WXS ## hnsc_tcga HNSC broad.mit.edu 37 Capture ## kirc_tcga KIRC broad.mit.edu 37 Capture ## luad_tcga LUAD broad.mit.edu 37 WXS ## lusc_tcga LUSC broad.mit.edu 37 WXS ## ov_tcga OV broad.mit.edu 37 WXS ## skcm_tcga SKCM broad.mit.edu 37 Capture ## thca_tcga THCA broad.mit.edu 37 WXS ## Sequencing_Phase Sequencer Number_Samples Number_Patients ## gbm_tcga Phase_I Illumina GAIIx 291 291 ## hnsc_tcga Phase_I Illumina GAIIx 319 319 ## kirc_tcga Phase_I Illumina GAIIx 297 293 ## luad_tcga Phase_I Illumina GAIIx 538 519 ## lusc_tcga Phase_I Illumina GAIIx 178 178 ## ov_tcga Phase_I Illumina GAIIx 142 142 ## skcm_tcga Phase_I Illumina GAIIx 266 264 ## thca_tcga Phase_I Illumina GAIIx 406 403 ## Cancer_Name ## gbm_tcga Glioblastoma multiforme ## hnsc_tcga Head and Neck squamous cell carcinoma ## kirc_tcga Kidney Chromophobe ## luad_tcga Lung adenocarcinoma ## lusc_tcga Lung squamous cell carcinoma ## ov_tcga Ovarian serous cystadenocarcinoma ## skcm_tcga Skin Cutaneous Melanoma ## thca_tcga Thyroid carcinoma
Next, we load a single dataset with the scaLoadDataset function.
ov = scaLoadDatasets("ov_tcga", merge = TRUE)
3 Exploring Mutational Data
The somatic variants of each study are represented as a object, ordered by genomic positions. Additional columns describe properties of the variant and relate it the the affected gene, sample, and patient.
head(ov, 3)
## GRanges object with 3 ranges and 14 metadata columns: ## seqnames ranges strand | Hugo_Symbol Entrez_Gene_Id ## <Rle> <IRanges> <Rle> | <factor> <integer> ## ov_tcga 1 [1334552, 1334552] * | CCNL2 81669 ## ov_tcga 1 [1961652, 1961652] * | GABRD 2563 ## ov_tcga 1 [2420688, 2420688] * | PLCH2 9651 ## Variant_Classification Variant_Type Reference_Allele ## <factor> <factor> <factor> ## ov_tcga Silent SNP C ## ov_tcga Silent SNP C ## ov_tcga Missense_Mutation SNP C ## Tumor_Seq_Allele1 Tumor_Seq_Allele2 Verification_Status ## <factor> <factor> <factor> ## ov_tcga C T Unknown ## ov_tcga C T Unknown ## ov_tcga C G Unknown ## Validation_Status Mutation_Status Patient_ID ## <factor> <factor> <factor> ## ov_tcga Valid Somatic TCGA-24-2262 ## ov_tcga Valid Somatic TCGA-24-1552 ## ov_tcga Valid Somatic TCGA-13-1484 ## Sample_ID index Dataset ## <factor> <integer> <factor> ## ov_tcga TCGA-24-2262-01A-01W-0799-08 3901 ov_tcga ## ov_tcga TCGA-24-1552-01A-01W-0551-08 3414 ov_tcga ## ov_tcga TCGA-13-1484-01A-01W-0545-08 1567 ov_tcga ## ------- ## seqinfo: 86 sequences from an unspecified genome
with(mcols(ov), table(Variant_Classification, Variant_Type))
## Variant_Type ## Variant_Classification DEL INS SNP ## 3'UTR 0 0 3 ## 5'Flank 0 0 1 ## 5'UTR 0 0 1 ## Frame_Shift_Del 79 0 0 ## Frame_Shift_Ins 0 16 0 ## IGR 0 0 5 ## In_Frame_Del 26 0 0 ## In_Frame_Ins 0 1 0 ## Intron 0 0 34 ## Missense_Mutation 0 0 4299 ## Nonsense_Mutation 0 0 285 ## Nonstop_Mutation 0 0 6 ## RNA 0 0 1 ## Silent 0 0 1417 ## Splice_Site 9 2 121 ## Translation_Start_Site 1 0 1
With such data at hand, we can identify the samples and genes haboring the most mutations.
head(sort(table(ov$Sample_ID), decreasing = TRUE))
## ## TCGA-09-2049-01D-01W-0799-08 TCGA-13-0923-01A-01W-0420-08 ## 119 118 ## TCGA-09-2050-01A-01W-0799-08 TCGA-25-1326-01A-01W-0492-08 ## 111 110 ## TCGA-25-1313-01A-01W-0492-08 TCGA-23-1110-01A-01D-0428-08 ## 104 102
head(sort(table(ov$Hugo_Symbol), decreasing = TRUE), 10)
## ## TP53 TTN PCDHAC2 MUC16 MUC17 PCDHGC5 USH2A CSMD3 CD163L1 ## 118 30 14 12 9 9 9 8 7 ## DYNC1H1 ## 7
4 Exploring Multiple Studies
Instead of focusing on an individual study, we can also import several at
once. The results are stored as a GRangesList in which each
element corresponds to a single study. This can be merged into a single GRanges
object with merge = TRUE.
three_studies = scaLoadDatasets(all_datasets[1:3]) print(elementLengths(three_studies))
## gbm_tcga hnsc_tcga kirc_tcga ## 22166 73766 26265
class(three_studies)
## [1] "SimpleGenomicRangesList" ## attr(,"package") ## [1] "GenomicRanges"
merged_studies = scaLoadDatasets(all_datasets[1:3], merge = TRUE) class(merged_studies)
## [1] "GRanges" ## attr(,"package") ## [1] "GenomicRanges"
We then compute the number of mutations per gene and study:
gene_study_count = with(mcols(merged_studies), table(Hugo_Symbol, Dataset)) gene_study_count = gene_study_count[order(apply(gene_study_count, 1, sum), decreasing = TRUE), ] gene_study_count = addmargins(gene_study_count) head(gene_study_count)
## Dataset ## Hugo_Symbol gbm_tcga hnsc_tcga kirc_tcga Sum ## Unknown 29 899 630 1558 ## TTN 121 401 125 647 ## TP53 101 323 8 432 ## MUC16 68 155 46 269 ## ADAM6 0 173 63 236 ## MUC4 17 32 130 179
Further, we can subset the data by regions of interests, and compute descriptive statistics only on the subset.
tp53_region = GRanges("17", IRanges(7571720, 7590863)) tp53_studies = subsetByOverlaps(merged_studies, tp53_region)
For example, we can investigate which type of somatic variants can be found in TP53 throughout the studies.
addmargins(table(tp53_studies$Variant_Classification, tp53_studies$Dataset))
## ## gbm_tcga hnsc_tcga kirc_tcga Sum ## Frame_Shift_Del 6 41 0 47 ## Frame_Shift_Ins 1 11 0 12 ## In_Frame_Del 2 7 0 9 ## In_Frame_Ins 0 2 0 2 ## Missense_Mutation 81 183 6 270 ## Nonsense_Mutation 4 54 0 58 ## Nonstop_Mutation 0 0 0 0 ## Silent 1 6 1 8 ## Splice_Site 6 19 1 26 ## Translation_Start_Site 0 0 0 0 ## RNA 0 0 0 0 ## Sum 101 323 8 432
To go further, how many patients have mutations in TP53 for each cancer type?
fraction_mutated_region = function(y, region) { s = subsetByOverlaps(y, region) m = length(unique(s$Patient_ID)) / metadata(s)$Number_Patients return(m) } mutated_fraction = sapply(three_studies, fraction_mutated_region, tp53_region) mutated_fraction = data.frame(name = names(three_studies), fraction = mutated_fraction)
library(ggplot2) p = ggplot(mutated_fraction) + ggplot2::geom_bar(aes(x = name, y = fraction, fill = name), stat = "identity") + ylim(0, 1) + xlab("Study") + ylab("Ratio") + theme_bw() print(p)
5 Data Provenance
5.1 TCGA Data
When importing the mutation data from the TCGA servers, we checked the data for consistency and fix common ambiguities in the annotation.
5.1.1 Processing
- Selection of the most recent somatic variant calls for each study. These were
stored as
*.maffiles in the TCGA data directory3. If both manually curated and automatically generated variant calls were available, the curated version was chosen. - Importing of the
*.maffiles into and checking for consistency with the TCGA MAF specifications4. Please note that these guidelines are currently only suggestions and most TCGA files violate some of these. - Transformation of the imported variants into a GRanges object, with one row for each reported variant. Only columns related to the genomic origin of the somatic variant were stored, additional columns describing higher-level effects, such as mutational consequences and alterations at the protein level, were dropped. The seqlevels information defining the chromosomal ranges were taken from the 1000genomes phase 2 reference assembly5.
- The patient barcode was extracted from the sample barcode.
- Metadata describing the design and analysis of the study was extracted.
- The processed variants were written to disk, with one file for each study. The metadata for all studies were stored as a single, separate object.
5.1.2 Selection Criteria of Data Sets
We included data sets in the package that were
- conducted by the Broad Institute.
- cleared for unrestricted access and usage6.
- sequenced with Illumina platforms.
5.1.3 Consistency Check
According to the TCGA specifications for the MAF files, we screened and
corrected for common artifacts in the data regarding annotation. This included:
- Transfering of all genomic coordinates to the NCBI 37 reference notation (with the chromosome always depicted as 'MT')
- Checking of the entries against all allowed values for this field (currently
for the columns
Hugo_Symbol,Chromosome,Strand,Variant_Classification,Variant_Type,Reference_Allele,Tumor_Seq_Allele1,Tumor_Seq_Allele2,Verification_Status,Validation_Status,Sequencer).
6 Alternatives
The TCGA data sets can be accessed in different ways. First, the TCGA itself offers access to certain types of its collected data7. Another approach has been taken by the cBioPortal for Cancer Genomics8 which has performed high-level analyses of several TCGA data sources, such as gene expression and copy number changes. This summarized data can be queried through an interface9.
7 Session Info
## R version 3.2.0 (2015-04-16) ## Platform: x86_64-unknown-linux-gnu (64-bit) ## Running under: Ubuntu 14.04.2 LTS ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 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 ## ## attached base packages: ## [1] stats4 parallel stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] ggbio_1.16.0 ggplot2_1.0.1 ## [3] GenomicRanges_1.20.1 GenomeInfoDb_1.4.0 ## [5] IRanges_2.2.0 S4Vectors_0.6.0 ## [7] BiocGenerics_0.14.0 SomaticCancerAlterations_1.4.0 ## ## loaded via a namespace (and not attached): ## [1] VariantAnnotation_1.14.0 reshape2_1.4.1 ## [3] splines_3.2.0 lattice_0.20-31 ## [5] colorspace_1.2-6 rtracklayer_1.28.0 ## [7] GenomicFeatures_1.20.0 RBGL_1.44.0 ## [9] survival_2.38-1 XML_3.98-1.1 ## [11] foreign_0.8-63 DBI_0.3.1 ## [13] BiocParallel_1.2.0 RColorBrewer_1.1-2 ## [15] lambda.r_1.1.7 plyr_1.8.1 ## [17] stringr_0.6.2 zlibbioc_1.14.0 ## [19] Biostrings_2.36.0 munsell_0.4.2 ## [21] gtable_0.1.2 futile.logger_1.4 ## [23] OrganismDbi_1.10.0 evaluate_0.6 ## [25] labeling_0.3 latticeExtra_0.6-26 ## [27] Biobase_2.28.0 knitr_1.9 ## [29] GGally_0.5.0 biomaRt_2.24.0 ## [31] AnnotationDbi_1.30.0 highr_0.4.1 ## [33] proto_0.3-10 Rcpp_0.11.5 ## [35] acepack_1.3-3.3 BSgenome_1.36.0 ## [37] scales_0.2.4 formatR_1.1 ## [39] graph_1.46.0 Hmisc_3.15-0 ## [41] XVector_0.8.0 Rsamtools_1.20.0 ## [43] gridExtra_0.9.1 digest_0.6.8 ## [45] biovizBase_1.16.0 grid_3.2.0 ## [47] tools_3.2.0 bitops_1.0-6 ## [49] exomeCopy_1.14.0 RCurl_1.95-4.5 ## [51] RSQLite_1.0.0 dichromat_2.0-0 ## [53] Formula_1.2-1 cluster_2.0.1 ## [55] futile.options_1.0.0 MASS_7.3-40 ## [57] reshape_0.8.5 rpart_4.1-9 ## [59] GenomicAlignments_1.4.0 nnet_7.3-9