TCGA Consensus and Comparison Revised (1 of 2)

.circleci/config.yml

@@ -172,9 +172,13 @@ jobs:
          ################################
 
 
+      - run:
+          name: TCGA SNV Caller Analysis 
+          command: ./scripts/run_in_ci.sh bash analyses/snv-callers/run_caller_consensus_analysis-tcga.sh     
+
       - run:
           name: SNV Caller Analysis 
-          command: OPENPBTA_VAF_CUTOFF=0.5 ./scripts/run_in_ci.sh bash analyses/snv-callers/run_caller_consensus_analysis.sh
+          command: OPENPBTA_VAF_CUTOFF=0.5 ./scripts/run_in_ci.sh bash analyses/snv-callers/run_caller_consensus_analysis-pbta.sh
 
       - run:
           name: Lancet WXS vs WGS test

analyses/README.md

@@ -39,7 +39,7 @@ Note that _nearly all_ modules use the harmonized clinical data file (`pbta-hist
 | [`sample-distribution-analysis`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/sample-distribution-analysis) | `pbta-histologies.tsv` | Produces plots and tables that illustrate the distribution of different histologies in the PBTA data | N/A
 | [`selection-strategy-comparison`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/selection-strategy-comparison) | `pbta-gene-expression-rsem-fpkm.polya.rds` <br> `pbta-gene-expression-rsem-fpkm.stranded.rds` | Comparison of RNA-seq data from different selection strategies | N/A
 | [`sex-prediction-from-RNASeq`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/sex-prediction-from-RNASeq) | `pbta-gene-expression-kallisto.stranded.rds` <br> `pbta-histologies.tsv` | *In progress*; predicts genetic sex using RNA-seq data ([#84](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/6)) | N/A
-| [`snv-callers`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/snv-callers) | `pbta-snv-lancet.vep.maf.gz` <br> `pbta-snv-mutect2.vep.maf.gz` <br> `pbta-snv-strelka2.vep.maf.gz` <br> `pbta-snv-vardict.vep.maf.gz` | Generates consensus SNV and indel calls; calculates tumor mutation burden using the consensus calls | `results/consensus/pbta-snv-consensus-mutation.maf.tsv.gz` <br> `results/consensus/pbta-snv-consensus-mutation-tmb.tsv` (included in data download; too large for tracking via GitHub)
+| [`snv-callers`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/snv-callers) | `pbta-snv-lancet.vep.maf.gz` <br> `pbta-snv-mutect2.vep.maf.gz` <br> `pbta-snv-strelka2.vep.maf.gz` <br> `pbta-snv-vardict.vep.maf.gz` <br> `tcga-snv-lancet.vep.maf.gz` <br> `tcga-snv-mutect2.vep.maf.gz` <br> `tcga-snv-strelka2.vep.maf.gz` | Generates consensus SNV and indel calls for PBTA and TCGA data; calculates tumor mutation burden using the consensus calls | `results/consensus/pbta-snv-consensus-mutation.maf.tsv.gz` <br> `results/consensus/pbta-snv-consensus-mutation-tmb.tsv` <br> `results/consensus/pbta-snv-consensus-mutation-tmb-coding.tsv`(included in data download; too large for tracking via GitHub) <br> `results/consensus/tcga-snv-consensus-mutation.maf.tsv.gz` <br> `results/consensus/tcga-snv-consensus-mutation-tmb.tsv` <br> `results/consensus/tcga-snv-consensus-mutation-tmb-coding.tsv`
 | [`ssgsea-hallmark`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/ssgsea-hallmark) | `pbta-gene-counts-rsem-expected_count.stranded.rds` | *Deprecated*; performs GSVA using Hallmark gene sets | N/A
 | [`survival-analysis`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/survival-analysis) | TBD | *In progress*; will eventually contain functions for various types of survival analysis ([#18](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/18)) | N/A
 | [`sv-analysis`](https://github.com/AlexsLemonade/OpenPBTA-analysis/tree/master/analyses/sv-analysis) | `pbta-sv-manta.tsv.gz` <br> `independent-specimens.wgs.primary-plus.tsv` | *In progress*; chromothripsis analysis per [#27](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/27)| N/A

analyses/snv-callers/.gitignore

@@ -1,3 +1,4 @@
 # ignore folders with big results files
 results
-ref_files
+ref_files/*
+!ref_files/gencode.v19.basic.exome.hg38liftover.bed

analyses/snv-callers/README.md

@@ -30,15 +30,22 @@ To run the evaluations and comparisons of all the SNV callers, call the bash scr
 ```
 bash run_caller_analysis.sh
 ```
+
+For the TCGA data, it has its own script to run the same methods:
+
+```
+bash run_caller_analysis-tcga.sh
+```
+
 This bash script will return:
 
 - Comparison plots in a notebook: [`compare_snv_callers_plots.nb.html`](https://cansavvy.github.io/openpbta-notebook-concept/snv-callers/compare_snv_callers_plots.nb.html).
 - A zip file containing:
-  - `pbta-snv-consensus-mutation.maf.tsv` - is  [MAF-like file](#consensus-mutation-call) that contains the snvs that were called by all three of these callers for a given sample are saved to this file.
+  - `pbta/tcga-snv-consensus-mutation.maf.tsv` - is  [MAF-like file](#consensus-mutation-call) that contains the snvs that were called by all three of these callers for a given sample are saved to this file.
   These files combine the [MAF file data](https://docs.gdc.cancer.gov/Data/File_Formats/MAF_Format/) from 3 different SNV callers: [Mutect2](https://software.broadinstitute.org/cancer/cga/mutect), [Strelka2](https://github.com/Illumina/strelka), and [Lancet](https://github.com/nygenome/lancet).
   See the methods on the callers' settings [here](https://github.com/AlexsLemonade/OpenPBTA-manuscript/blob/master/content/03.methods.md#somatic-single-nucleotide-variant-calling) and see [the methods of this caller analysis and comparison below](#summary-of-methods).  
-  - `pbta-snv-consensus-mutation-tmb-coding.tsv` - Tumor Mutation burden calculations using *coding only* mutations use the consensus of Lancet, Mutect2, and Strelka2.
-  - `pbta-snv-consensus-mutation-tmb-all.tsv` - Tumor Mutation burden calculations using *all* mutations use the consensus of Mutect2, and Strelka2. (Lancet was excluded because it has a [coding region bias in the way it was run](https://github.com/AlexsLemonade/OpenPBTA-manuscript/blob/master/content/03.methods.md#snv-and-indel-calling)).
+  - `pbta/tcga-snv-consensus-mutation-tmb-coding.tsv` - Tumor Mutation burden calculations using *coding only* mutations use the consensus of Lancet, Mutect2, and Strelka2.
+  - `pbta/tcga-snv-consensus-mutation-tmb-all.tsv` - Tumor Mutation burden calculations using *all* mutations use the consensus of Mutect2, and Strelka2. (Lancet was excluded because it has a [coding region bias in the way it was run](https://github.com/AlexsLemonade/OpenPBTA-manuscript/blob/master/content/03.methods.md#snv-and-indel-calling)).
 
 ## Summary of Methods
 
@@ -62,7 +69,7 @@ As Strelka2 does not call multinucleotide variants (MNV), but instead calls each
 ### Tumor Mutation Burden Calculation
 
 For each experimental strategy and TMB calculation, the intersection of the genomic regions effectively being surveyed are used.
-These genomic regions are used for first filtering mutations to these regions and then for using the size in bp of the genomic regions surveyed as the TMB denominator. 
+These genomic regions are used for first filtering mutations to these regions and then for using the size in bp of the genomic regions surveyed as the TMB denominator.
 
 #### All mutations TMB
 
@@ -84,7 +91,7 @@ SNVs outside of these coding sequences are filtered out before being summed and
 ```
 WGS_coding_only_TMB = (total # coding sequence snvs called by all three of Strelka, Lancet, and Mutect2 ) / intersection_strelka_lancet_mutect_CDS_genome_size
 ```
-Because the same WXS BED file applies to all callers, that file is intersected with the coding sequences for filtering and for determining the denominator. 
+Because the same WXS BED file applies to all callers, that file is intersected with the coding sequences for filtering and for determining the denominator.
 ```
 WXS_coding_only_TMB = (total # coding sequence snvs called by all three of Strelka, Lancet, and Mutect2 ) /
 intersection_wxs_CDS_genome_size

analyses/snv-callers/compare_snv_callers_plots-tcga.Rmd

@@ -93,7 +93,7 @@ Connect to SQLite database.
 ```{r}
 # Start up connection
 con <- DBI::dbConnect(RSQLite::SQLite(), 
-                      file.path(scratch_dir, "tcga_snv_db.sqlite"))
+                      file.path(scratch_dir, "tcga_v2_snv_db.sqlite"))
 ```
 
 Note what columns we will join by.
@@ -389,52 +389,6 @@ perc_var_df %>%
 ggplot2::ggsave(file.path(plots_dir, "tcga-variant_classification_plot.png"))
 ```
 
-## Where are the unique Lancet calls? 
-
-```{r}
-source(file.path("..", "chromosomal-instability", "util", "chr-break-plot.R"))
-source(file.path("..", "chromosomal-instability", "util", "chr-break-calculate.R"))
-```
-
-```{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")))
-# Make chromosome size named vector
-chr_sizes_vector <- chr_sizes$X3
-names(chr_sizes_vector) <- chr_sizes$X1
-```
-
-```{r}
-only_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
 
 ```{r}