In MSKCC-Epi-Bio/gnomeR: Wrangle and analyze IMPACT and TCGA mutation data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

In this vignette we will walk through a data example to show available {gnomeR} data processing, visualization and analysis functions. We will also outline some {gnomeR} helper functions to use when you encounter common pitfalls and format inconsistencies when working with mutation, CNA or structural variant data.

Setting up

Make sure {gnomeR} is installed & loaded. We will use dplyr for the purposes of this vignette as well.

library(gnomeR)
library(dplyr)

To demonstrate {gnomeR} functions, we will be using a random sample of 200 patients from a publicly available prostate cancer study retrieved from cBioPortal. Data on mutations, CNA, and structural variants for these patients are available in the gnomeR package (gnomeR::mutations, gnomeR::cna, gnomeR::sv).

Note: To access data from cBioPortal, you can use the {cbioportalR} package: https://github.com/karissawhiting/cbioportalR.

Data Formats

Mutation, CNA, or structural variant (also called fusion) data may be formatted differently depending on where you source it. Below we outline some differences you may encounter when downloading data from the cBioPortal website versus pulling it via the cBioPortal API, and review the important/required columns for each. Example data in this package was pulled using the API.

See cBioPortal documentation for more details on different file formats supported on cBioPortal, and their data schema and coding.

Mutation Data

The most common mutation data format is the Mutation Annotation Format (MAF) created as part of The Cancer Genome Atlas (TCGA) project. Each row in a MAF file represents a unique gene for a specified sample, therefore there are usually several rows per sample. To use MAF files, you need at minimum a sample ID column and a hugo symbol column, though often additional information like mutation type or location are also necessary.

MAF formats are fairly consistent across sources, however if you download the raw data from a study on cBioPortal using the interactive download button you might notice some differences in the variables names in comparison to data imported from the API. For instance, below web-downloaded MAF names are on the left and API-downloaded MAF names are on the right:

* `Tumor_Sample_Barcode` is called `sampleId` 
* `Hugo_symbol` is called `hugoGeneSymbol` 
* `Variant_Classification` is called `mutationType`
* `HGVSp_Short` is called `proteinChange`
* `Chromosome` is called 'chr'

Some of the other variables are named differently as well but those differences are more intuitive. You can refer to gnomeR::names_df for more information on possible MAF variable names.

Luckily, most {gnomeR} functions use the data dictionary in gnomeR::names_df to automatically recognize the most common MAF variable names and turn them into clean, snakecase names in resulting dataframes. For example:

gnomeR::mutations %>% names()

rename_columns(gnomeR::mutations) %>% names()

As you can see, some variables, such as linkMsa, were not transformed because they are not used in {gnomeR} functions.

CNA data

The discrete copy number data from cBioPortal contains values that would be derived from copy-number analysis algorithms like GISTIC 2.0 or RAE. CNA data is often presented in a long or wide format:

1. Long-format - Each row is a CNA event for a given gene and sample, therefore samples often have multiple rows. This is most common format you will receive when downloading data using the API.

gnomeR::cna[1:6, ]

1. Wide-format - Organized such that there is one column per sample and one row per gene. Thus, each sample's events are contained within one column. This is most common format you will receive when downloading data from the cBioPortal web browser.

gnomeR::cna_wide[1:6, 1:6]

{gnomeR} features two helper functions to easily pivot from wide- to long-format and vice-versa.

pivot_cna_wider(rename_columns(gnomeR::cna))
pivot_cna_longer(gnomeR::cna_wide)

These functions will also relabel CNA levels (numeric values) to characters as shown below:

allowed_cna_levels <- tibble::tribble(
               ~detailed_coding, ~numeric_coding,   ~simplified_coding,
                      "neutral",             "0",       "neutral",
          "homozygous deletion",            "-2",      "deletion",
                          "loh",          "-1.5",      "deletion",
          "hemizygous deletion",            "-1",      "deletion",
                         "gain",             "1", "amplification",
      "high level amplification",            "2", "amplification")

allowed_cna_levels %>% knitr::kable()

{gnomeR} automatically checks CNA data labels and recodes as needed within functions. You can also use the recode_cna() function to do it yourself if pivoting is unnecessary:

gnomeR::cna %>%
  mutate(alteration_recoded = recode_cna(alteration))%>%
  select(hugoGeneSymbol, sampleId, alteration, alteration_recoded)%>%
  head()

Preparing Data For Analysis

Process Data with create_gene_binary()

Often the first step to analyzing genomic data is organizing it in an event matrix. This matrix will have one row for each sample in your cohort and one column for each type of genomic event. Each cell will take a value of 0 (no event on that gene/sample), 1 (event on that gene/sample) or NA (missing data or gene not tested on panel). The create_gene_binary() function helps you process your data into this format for use in downstream analysis.

You can create_gene_binary() from any single type of data (mutation, CNA or fusion):

create_gene_binary(mutation = gnomeR::mutations)[1:6, 1:6]

or you can process several types of alterations into a single matrix. Supported data types are:

• mutations
• copy number amplifications
• copy number deletions
• gene fusions

All datasets should be in long-format.

When processing multiple types of alteration data, by default there will be a separate column for each type of alteration on that gene. For example, if the TP53 gene could have up to 4 columns: TP53 (mutation), TP53.Amp (amplification), TP53.Del (deletion), and TP53.fus (structural variant). Further, if no events are observed within the data set for a type of alteration (let's say no TP53 mutations but some TP53 amplifications), columns with all zeros will be excluded. For example, in colnames(all_bin) below, there is at least one ERG.fus event, but no ERG mutation events.

Note the use of the samples argument. This allows you to specify exactly which samples are in your resulting data frame. This argument allows you to retain samples that have no genetic events as a row of 0 and NA. If you do not specify such samples within your sample cohort, rows with no alterations will be excluded from the final matrix.

samples <- unique(gnomeR::mutations$sampleId)[1:10]

all_bin <- create_gene_binary(
    samples = samples,
    mutation = gnomeR::mutations,
    cna = gnomeR::cna,
    fusion = gnomeR::sv
)

all_bin[1:6, 1:6]

colnames(all_bin)

Notes on some helpful create_gene_binary() arguments:

• mut_type- by default, any germline mutations will be omitted because data is often incomplete, but you can choose to leave them in if needed.
• specify_panel- If you are working across a set of samples that was sequenced on several different gene panels, this argument will insert NAs for the genes that weren't tested for any given sample. You can pass a string "impact" indicating automatically guessing panels and processsing IMPACT samples based on ID, or you can pass a data frame with columns sample_id and gene_panel for more fine grained control of NA annotation.
• recode_aliases - Sometimes genes have several accepted names or change names over time. This can be an issue if genes are coded under multiple names in studies, or if you are working across studies. By default, this function will search for aliases for genes in your data set and resolved them to their current most common name.

Collapse Data with summarize_by_gene()

If the type of alteration event (mutation, amplification, deletion, structural variant) does not matter for your analysis, and you want to see if any event occurred for a gene, pipe your create_gene_binary() object through the summarize_by_gene() function. As you can see, this compresses all alteration types of the same gene into one column. So, where in all_bin there was an ERG.fus column but no ERG column, now summarize_by_gene() only has an ERG column with a 1 for any type of event.

dim(all_bin)

all_bin_summary <- all_bin %>% 
  summarize_by_gene()

all_bin_summary[1:6, 1:6]

colnames(all_bin_summary)

Analyzing Data

Once you have processed the data into a binary format, you may want to visualize and summarize it with the following helper functions:

Summarize Alterations with subset_by_frequency() and tbl_genomic()

You can use the subset_by_frequency() function along with tbl_genomic() function to easily display highest prevalence alterations in your data set.

First, subset your data set by top genes at a give threshold (t). The threshold is a value between 0 and 1 to indicate % prevalence cutoff to keep. You can subset at this threshold on the alteration level:

top_prev_alts <- all_bin %>%
  subset_by_frequency(t = .1)

colnames(top_prev_alts)

Or gene level:

top_prev_genes <- all_bin_summary %>%
  subset_by_frequency(t = .1)

colnames(top_prev_genes)

Next, we can pass our subset of genes or alterations to the tbl_genomic() function, which can be used to display summary tables of alterations. It is built off the {gtsummary} package and therefore you can use most customizations available in that package to alter the look of your tables.

The gene_binary argument expects binary data as generated by the create_gene_binary(), summarize_by_gene() or subset_by_frequency() functions. Example below shows the summary using gene data for ten samples.

In the example below for tb1, if the gene is altered in at least 15% of samples, the gene is included in the summary table.

samples <- unique(mutations$sampleId)[1:10]

gene_binary <- create_gene_binary(
  samples = samples,
  mutation = mutations,
  cna = cna,
  mut_type = "somatic_only", snp_only = FALSE,
  specify_panel = "no"
)

tbl1 <- gene_binary %>% 
  subset_by_frequency(t = .15) %>%
  tbl_genomic() 

We can add additional customizations to the table with the following gtsummary functions:

tbl1 %>%
  gtsummary::bold_labels() 

Annotate Gene Pathways with add_pathways()

The add_pathways() function allows you add columns to your gene binary matrix that annotate custom gene pathways, or oncogenic signaling pathways (add citation).

The function expects a binary matrix as obtained from the gene_binary() function and will return a gene binary with additional columns added for specified pathways.

There are a set of default pathways available in the package that can be viewed using gnomeR::pathways (Sanchez-Vega, F et al., 2018). This new data frame will include columns for mutations, CNAs, structural variants, and pathways. You can subset to only the pathways if you choose.

# available pathways
names(gnomeR::pathways)

pathways <- add_pathways(gene_binary, pathways = c("Notch", "p53")) %>%
  select(sample_id, pathway_Notch, pathway_p53)

head(pathways)

Data Visualizations

The mutation_viz functions allows you to visualize data for the variables related to variant classification, variant type, SNV class as well as top variant genes.

mutation_viz(mutations)

Customizing your colors

{gnomeR} comes with 3 distinct palettes that can be useful for plotting high dimensional genomic data: pancan, main, and sunset. pancan and main offer a wide range of colors that can be useful to map to discrete scales. sunset has fewer colors but offers a color spectrum useful for interpolation for continuous variables. You can view hex codes of all colors with gnomer_colors, and the 3 distinct palettes with gnomer_palettes$pancan, gnomer_palettes$main, gnomer_palettes$sunset. The gnomer_palette() is used to set/subset specific palettes, create a continuous palette, and/or plot palettes to show the user the specific colors they chose.

The code below will show the first 4 colors from the pancan palette for a discrete palette.

gnomer_palette(
  name = "pancan",
  n = 4,
  type = "discrete",
  plot_col = TRUE,
  reverse = FALSE
)

If you wanted to make a continuous palette you can change the type= option and specify how many colors your want to use to create a gradient palette. type = continuous uses grDevices::colorRampPalette() as the engine to create the gradient palette. Additionally, gnomer_palette() accepts additional arguments to be passed to colorRampPalette() with the parameter .... The example below uses 20 colors.

gnomer_palette(
  name = "sunset",
  n = 20,
  type = "continuous",
  plot_col = TRUE,
  reverse = FALSE
)

Examples how to use gnomer_palette():

library(ggplot2)

ggplot(iris, aes(Sepal.Width, Sepal.Length, color = Species)) +
  geom_point(size = 4) +
  scale_color_manual(values = gnomer_palette("pancan"))

# use a continuous color scale - interpolates between colors
ggplot(iris, aes(Sepal.Width, Sepal.Length, color = Sepal.Length)) +
  geom_point(size = 4, alpha = .6) +
  scale_color_gradientn(colors = 
                          gnomer_palette("sunset", type = "continuous"))

Set Color Palettes Globally

If you do not want to constantly have to set the palette for each plot you can set a palette theme globally for your entire session using set_gnomer_palette(). Additionally, you can set the discrete and gradient (appropriate for continuous variables) palette independently. With the examples below you can see the default colors change for each plot without having to call a scale_* function.

set_gnomer_palette(palette = "main", gradient = "sunset")

ggplot(mtcars, aes(wt, mpg, color = factor(cyl))) +
  geom_point()

ggplot(mtcars, aes(wt, mpg, color = cyl)) +
  geom_point()

You can reset the palettes back to default ggplot color palettes with reset_gnomer_palette().

reset_gnomer_palette()

ggplot(mtcars, aes(wt, mpg, color = cyl)) +
  geom_point()

MSKCC-Epi-Bio/gnomeR documentation built on Oct. 17, 2024, 3:39 p.m.