TCGA Plot Caller Comparisons
-C. Savonen for ALSF CCDL
-2019
- -Purpose: After creating a TCGA consensus file with 02-merge_callers.R, this notebook plots these data in a series of plots to compare the consensus calls to the individual callers.
Summary of Findings:
-Individual Caller’s Reports
-Outline of analyses completed:
--
-
- Upset graph
-
- - VAF of combinations of callers -
- Base changes by caller -
- VAF distributions by caller -
- VAF correlations between callers -
Usage
-To run this from the command line, use:
-Rscript -e "rmarkdown::render('analyses/snv-callers/compare_snv_callers_plots-tcga.Rmd',
- clean = TRUE)"This assumes you are in the top directory of the repository.
-Setup
-Packages and functions
-Read in set up script.
- - - -rr upset_plot <- function(callers_per_mutation_df) { callers_per_mutation_df %>% dplyr::group_by(caller_combo) %>% dplyr::tally() %>% ggplot2::ggplot(ggplot2::aes(x = reorder(caller_combo, -n), y = n)) + ggplot2::geom_bar( position = , stat = ) + ggplot2::theme_classic() + ggupset::scale_x_mergelist(sep = -) + ggupset::axis_combmatrix(sep = -) + ggplot2::xlab(\) + ggplot2::ylab(\) }
Set up directories.
- - - -rr # Start up connection con <- DBI::dbConnect(RSQLite::SQLite(), file.path(scratch_dir, _snv_db.sqlite))
rr join_cols = c(, _Position, _Allele, , _Sample_Barcode, _Classification)
Connect to TCGA database
-Connect to SQLite database.
- - - -rr strelka <- dplyr::tbl(con, ) %>% dplyr::select(join_cols, ) %>% as.data.frame()
lancet <- dplyr::tbl(con, ) %>% dplyr::select(join_cols, ) %>% as.data.frame()
-mutect <- dplyr::tbl(con, ) %>% dplyr::select(join_cols, ) %>% as.data.frame()
-Note what columns we will join by.
- - - -rr # Bring out the data.frame all_caller_df <- strelka %>% dplyr::full_join(mutect, by = join_cols, suffix = c(_strelka, _mutect)) %>% dplyr::full_join(lancet, by = join_cols) %>% dplyr::rename(VAF_lancet = VAF)
Take a peek at what this looks like
-head(all_caller_df)
-Set up tables from the database but only call the columns we need.
- - - -strelka <- dplyr::tbl(con, "strelka") %>%
- dplyr::select(join_cols, "VAF") %>%
- as.data.frame()
-
-lancet <- dplyr::tbl(con, "lancet") %>%
- dplyr::select(join_cols, "VAF") %>%
- as.data.frame()
-
-mutect <- dplyr::tbl(con, "mutect") %>%
- dplyr::select(join_cols, "VAF") %>%
- as.data.frame()Because DBI does not support full_join, we had to use a series of left_join and union_all calls in order to get a full join of all the callers.
rr # Take out the consensus VAF for this combination long_caller_df <- long_vaf_df %>% dplyr::filter(caller != )
Determine the combination of callers for each mutation
-callers_per_mutation <- tapply( long_caller_df\(caller, long_caller_df\)index, paste0, collapse = -
-) %>% as.data.frame() %>% tibble::rownames_to_column()
Determine Median VAF for each mutation
-vaf_med <- tapply( long_caller_df\(vaf, long_caller_df\)index, median ) %>% as.data.frame() %>% tibble::rownames_to_column()
-Join the median VAF and the callers that call that mutation into one data.frame
-callers_per_mutation <- callers_per_mutation %>% dplyr::inner_join(vaf_med, by = ) %>% dplyr::left_join(long_caller_df %>% dplyr::select(Variant_Classification, index), by = ) %>% # Make column names more sensible dplyr::rename(caller_combo = ..x, median_vaf = ..y)
-Column `index` joining character vector and factor, coercing into character vectorrr upsettr_plot <- callers_per_mutation %>% upset_plot()
Print this out here
-upsettr_plot
-
Get median vaf for each mutation with each caller combination.
- - - -rr # We can save the plot like a normal ggplot png(file = file.path(plots_dir, -upset_plot.png)) upsettr_plot dev.off()
null device
- 1 Upset graph
-Make the upset plot.
- - - -rr callers_per_mutation %>% dplyr::filter(!(Variant_Classification %in% c(\5’Flank, \3’Flank, , ))) %>% upset_plot()
Save this plot to a png
- - - -rr callers_per_mutation %>% dplyr::filter(Variant_Classification %in% c(\5’Flank, \3’Flank, , )) %>% upset_plot()
Transcripts only upset plot.
- - - -rr # Make this plot callers_per_mutation %>% ggplot2::ggplot(ggplot2::aes(x = caller_combo, y = median_vaf)) + ggplot2::geom_violin() + ggplot2::theme_classic() + ggupset::scale_x_mergelist(sep = -) + ggupset::axis_combmatrix(sep = -) + ggplot2::xlab(\) + ggplot2::ylab(VAF Across Callers)
Non transcript upset plot.
- - - -rr # Save this plot ggplot2::ggsave(file.path(plots_dir, -upset_median_vaf_plot.png))
Saving 7 x 7 in imageVAF for combination sets of callers
-Graph mutations by combinations of callers and plot the VAF density. Get list of callers per mutation and the median vaf for each. Graph the median VAF for each combination of callers.
- - - -rr # Summarize by the number of times each base change shows up in each category. perc_change_df <- all_caller_df %>% # Make change variable dplyr::mutate(base_change = paste0(Allele, >, Reference_Allele)) %>% dplyr::mutate( # From the base_change variable, summarize insertions, deletions, and # changes that are more than one base into their own groups. change = dplyr::case_when( grepl(^-, base_change) ~ , grepl(-$, base_change) ~ , nchar(base_change) > 3 ~ change, TRUE ~ base_change ) ) %>% # Whittle down to necessary columns dplyr::select(change, dplyr::starts_with()) %>% # Make this long form tidyr::gather(key = , value = , -change) %>% # Drop the prefix dplyr::mutate(caller = gsub(_, \, caller)) %>% # Get rid of NA rows dplyr::filter(!is.na(vaf)) %>% # Summarize the number of mutations per caller dplyr::count(caller, change, name = ) %>% dplyr::add_count(caller, wt = count) %>% # Calculate the percent of each dplyr::mutate(percent = count / n) %>% # Drop nonsensical change dplyr::filter(!grepl(>|>N|C>C|T>T|A>A|G>G, change), !is.na(change)) %>% dplyr::mutate( change = as.factor(change), # Change factor level order so ins and del are at the end change = forcats::fct_relevel(change, , , _change, after = Inf) )
rr perc_change_df %>% ggplot2::ggplot(ggplot2::aes(x = change, y = percent)) + ggplot2::theme_classic() + ggplot2::geom_bar( position = , stat = , ggplot2::aes(fill = caller) ) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + ggplot2::xlab(\) + ggplot2::ylab(of callers’ mutations) + colorblindr::scale_fill_OkabeIto()
Base Changes
-Summarize base change information per caller.
- - - -rr # Save this plot ggplot2::ggsave(file.path(plots_dir, -perc_base_change_plot.png))
Saving 7 x 7 in imageMake a barplot illustrating the percent of the mutations for each caller that that are each type of change.
- - - -rr long_vaf_df %>% ggplot2::ggplot(ggplot2::aes(x = caller, y = vaf, color = caller)) + ggplot2::geom_violin() + ggplot2::theme_classic() + colorblindr::scale_color_OkabeIto()
rr # Save this plot ggplot2::ggsave(file.path(plots_dir, -vaf_violin_plot.png))
Saving 7 x 7 in imageVAF distributions
-Make a violin plot of each caller.
- - - -rr GGally::ggpairs(vaf_df %>% dplyr::select(-VAF_consensus, -Variant_Classification), mapping = ggplot2::aes(alpha = 0.05)) + ggplot2::theme_classic()
rr # Save this plot ggplot2::ggsave(file.path(plots_dir, -vaf_correlations_plot.png))
Saving 7 x 7 in imageVAF correlations
-Make the plot with ggpairs.
- - - -rr # Summarize by the number of times each base change shows up in each category. perc_var_df <- long_vaf_df %>% # Summarize the number of mutations per caller dplyr::count(caller, Variant_Classification, name = ) %>% dplyr::add_count(caller, wt = count) %>% # Calculate the percent of each dplyr::mutate(percent = count / n)
rr perc_var_df %>% ggplot2::ggplot(ggplot2::aes(x = caller, y = percent)) + ggplot2::theme_classic() + ggplot2::geom_bar( position = , stat = , ggplot2::aes(fill = Variant_Classification) ) + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + ggplot2::xlab(\) + ggplot2::ylab(of callers’ mutations)
Mutation Region barplot
-Summarize Variant_Classification information per caller.
rr # Save this plot ggplot2::ggsave(file.path(plots_dir, -variant_classification_plot.png))
Saving 7 x 7 in imageMake a barplot illustrating the percent of the mutations for each caller that that are each type of change.
- - - -rr sessionInfo()
R version 3.6.0 (2019-04-26)
-Platform: x86_64-pc-linux-gnu (64-bit)
-Running under: Debian GNU/Linux 9 (stretch)
-
-Matrix products: default
-BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
-
-locale:
- [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
- [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
- [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
- [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] stats graphics grDevices utils datasets methods base
-
-other attached packages:
-[1] optparse_1.6.2
-
-loaded via a namespace (and not attached):
- [1] tidyr_0.8.3 bit64_0.9-7 jsonlite_1.6
- [4] R.utils_2.9.0 assertthat_0.2.1 stats4_3.6.0
- [7] blob_1.1.1 GenomeInfoDbData_1.2.2 yaml_2.2.0
-[10] pillar_1.4.2 RSQLite_2.1.1 backports_1.1.4
-[13] glue_1.3.1 digest_0.6.20 polyclip_1.10-0
-[16] GenomicRanges_1.38.0 RColorBrewer_1.1-2 XVector_0.26.0
-[19] colorspace_1.4-1 plyr_1.8.4 cowplot_0.9.4
-[22] htmltools_0.3.6 R.oo_1.22.0 pkgconfig_2.0.2
-[25] GetoptLong_0.1.7 zlibbioc_1.32.0 purrr_0.3.2
-[28] scales_1.0.0 tweenr_1.0.1 getopt_1.20.3
-[31] ggforce_0.2.2 tibble_2.1.3 ggupset_0.1.0.9000
-[34] farver_1.1.0 IRanges_2.20.1 ggplot2_3.2.0
-[37] BiocGenerics_0.32.0 lazyeval_0.2.2 cli_1.1.0
-[40] magrittr_1.5 crayon_1.3.4 memoise_1.1.0
-[43] evaluate_0.14 GGally_1.4.0 R.methodsS3_1.7.1
-[46] fansi_0.4.0 MASS_7.3-51.4 forcats_0.4.0
-[49] rsconnect_0.8.13 tools_3.6.0 data.table_1.12.2
-[52] hms_0.4.2 GlobalOptions_0.1.1 ComplexHeatmap_2.2.0
-[55] stringr_1.4.0 S4Vectors_0.24.1 munsell_0.5.0
-[58] cluster_2.1.0 compiler_3.6.0 GenomeInfoDb_1.22.0
-[61] rlang_0.4.0 grid_3.6.0 RCurl_1.95-4.12
-[64] rstudioapi_0.10 rjson_0.2.20 circlize_0.4.8
-[67] colorblindr_0.1.0 bitops_1.0-6 base64enc_0.1-3
-[70] labeling_0.3 rmarkdown_1.13 gtable_0.3.0
-[73] reshape_0.8.8 DBI_1.0.0 reshape2_1.4.3
-[76] R6_2.4.0 knitr_1.23 dplyr_0.8.3
-[79] bit_1.1-14 utf8_1.1.4 zeallot_0.1.0
-[82] clue_0.3-57 rprojroot_1.3-2 shape_1.4.4
-[85] readr_1.3.1 stringi_1.4.3 parallel_3.6.0
-[88] Rcpp_1.0.1 vctrs_0.1.0 png_0.1-7
-[91] dbplyr_1.4.2 tidyselect_0.2.5 xfun_0.8 # Save this plot
-ggplot2::ggsave(file.path(plots_dir, "tcga-variant_classification_plot.png"))Where are the unique Lancet calls?
- - - -source(file.path("..", "chromosomal-instability", "util", "chr-break-plot.R"))
-source(file.path("..", "chromosomal-instability", "util", "chr-break-calculate.R"))# Set up Chr sizes
-chr_sizes <- readr::read_tsv(file.path(data_dir, "WGS.hg38.strelka2.unpadded.bed"),
- col_names = FALSE
-) %>%
- # Reformat the chromosome variable to drop the "chr"
- dplyr::mutate(X1 = factor(gsub("chr", "", X1),
- levels = c(1:22, "X", "Y", "M")
- )) %>%
- # Remove sex chromosomes
- dplyr::filter(!(X1 %in% c("X", "Y", "M")))Parsed with column specification:
-cols(
- X1 = [31mcol_character()[39m,
- X2 = [32mcol_double()[39m,
- X3 = [32mcol_double()[39m
-)# Make chromosome size named vector
-chr_sizes_vector <- chr_sizes$X3
-names(chr_sizes_vector) <- chr_sizes$X1only_lancet <- all_caller_df %>%
- dplyr::filter(!is.na(VAF_lancet) & is.na(VAF_strelka) & is.na(VAF_mutect)) %>%
- dplyr::mutate(Chromosome = gsub("chr", "", Chromosome))
-
-lancet_densities <- break_density(breaks_df = only_lancet,
- sample_id = "all",
- samples_col = "Tumor_Sample_Barcode",
- chrom_col = "Chromosome",
- start_col = "Start_Position",
- end_col = "Start_Position",
- window_size = 1e7,
- chr_sizes_vector = chr_sizes_vector)
-
-map_breaks_plot(lancet_densities,
- y_val = "total_counts",
- color = "blue",
- y_lab = "Total Num of Calls",
- main_title = "Lancet Only Calls"
- )
-
Session Info
- - - -sessionInfo()