Update MAF fields for v19 SNV consensus

.circleci/config.yml

@@ -225,21 +225,21 @@ jobs:
 #          name: RNA-Seq composition
 #          command: ./scripts/run_in_ci.sh Rscript -e "rmarkdown::render('analyses/rna-seq-composition/rna-seq-composition.Rmd', clean = TRUE)"
 
-    #   - run:
-    #       name: TCGA SNV Caller Analysis
-    #       command: ./scripts/run_in_ci.sh bash analyses/snv-callers/run_caller_consensus_analysis-tcga.sh
+      - 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-pbta.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-pbta.sh
 
-    #  - run:
-    #      name: Tumor mutation burden with TCGA
-    #      command: ./scripts/run_in_ci.sh Rscript -e "rmarkdown::render('analyses/tmb-compare/compare-tcga-pbta.Rmd', clean = TRUE)"
+      - run:
+          name: Tumor mutation burden with TCGA
+          command: ./scripts/run_in_ci.sh Rscript -e "rmarkdown::render('analyses/tmb-compare/compare-tcga-pbta.Rmd', clean = TRUE)"
 
-    #  - run:
-    #      name: Exploration of nonsynonymous filter
-    #      command: ./scripts/run_in_ci.sh bash analyses/snv-callers/explore_variant_classifications/run_explorations.sh
+      - run:
+          name: Exploration of nonsynonymous filter
+          command: ./scripts/run_in_ci.sh bash analyses/snv-callers/explore_variant_classifications/run_explorations.sh
 
       # This analysis was used to explore the TCGA PBTA data when the BED files used to calculate TCGA
       # were incorrect https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/568
@@ -275,9 +275,9 @@ jobs:
     #     name: d3b TMB code
     #     command: ./scripts/run_in_ci.sh bash analyses/tmb-compare/TMB_d3b_code/run_tmb_d3b.sh
 
-    #  - run:
-    #      name: Compare TMB calculations
-    #      command: ./scripts/run_in_ci.sh Rscript -e "rmarkdown::render('analyses/tmb-compare/compare-tmb-calculations.Rmd', clean = TRUE)"
+      - run:
+          name: Compare TMB calculations
+          command: ./scripts/run_in_ci.sh Rscript -e "rmarkdown::render('analyses/tmb-compare/compare-tmb-calculations.Rmd', clean = TRUE)"
 
       - run:
           name: Run survival plots

analyses/snv-callers/run_caller_consensus_analysis-pbta.sh

@@ -26,6 +26,7 @@ vaf_cutoff=${OPENPBTA_VAF_CUTOFF:-0}
 run_plots_nb=${OPENPBTA_PLOTS:-0}
 
 ################################ Set Up Database ################################
+echo "Setting up Database"
 python3 analyses/snv-callers/scripts/01-setup_db.py \
   --db-file $dbfile \
   --strelka-file data/pbta-snv-strelka2.vep.maf.gz \
@@ -35,19 +36,22 @@ python3 analyses/snv-callers/scripts/01-setup_db.py \
   --meta-file data/pbta-histologies.tsv
 
 ##################### Merge callers' files into total files ####################
+echo "Merging callers"
 Rscript analyses/snv-callers/scripts/02-merge_callers.R \
   --db_file $dbfile \
   --output_file $consensus_file \
   --vaf_filter $vaf_cutoff \
   --overwrite
 
 ########################## Add consensus to db ################################
+echo "Adding consensus to database"
 python3 analyses/snv-callers/scripts/01-setup_db.py \
   --db-file $dbfile \
   --consensus-file $consensus_file
 
 ############# Create intersection BED files for TMB calculations ###############
 # Make All mutations BED files
+echo "Making intersection bed files"
 bedtools intersect \
   -a data/WGS.hg38.strelka2.unpadded.bed \
   -b data/WGS.hg38.mutect2.vardict.unpadded.bed \
@@ -56,6 +60,7 @@ bedtools intersect \
 #################### Make coding regions file
 # Convert GTF to BED file for use in bedtools
 # Here we are only extracting lines with as a CDS i.e. are coded in protein
+echo "Making CDS bed file"
 gunzip -c data/gencode.v27.primary_assembly.annotation.gtf.gz \
   | awk '$3 ~ /CDS/' \
   | convert2bed --do-not-sort --input=gtf - \
@@ -64,6 +69,7 @@ gunzip -c data/gencode.v27.primary_assembly.annotation.gtf.gz \
   > $cds_file
 
 ######################### Calculate consensus TMB ##############################
+echo "Calculating TMB"
 Rscript analyses/snv-callers/scripts/03-calculate_tmb.R \
   --db_file $dbfile \
   --output analyses/snv-callers/results/consensus \
@@ -79,5 +85,6 @@ gzip $consensus_file
 ############################# Comparison Plots #################################
 if [ "$run_plots_nb" -gt "0" ]
 then
+ echo "Making comparison plots"
  Rscript -e "rmarkdown::render('analyses/snv-callers/compare_snv_callers_plots.Rmd', clean = TRUE)"
 fi

analyses/snv-callers/scripts/01-setup_db.py

@@ -101,19 +101,6 @@
     ('Matched_Norm_Sample_Barcode', 'TEXT'),
     ('Match_Norm_Seq_Allele1', 'TEXT'),
     ('Match_Norm_Seq_Allele2', 'TEXT'),
-    ('Tumor_Validation_Allele1', 'TEXT'),
-    ('Tumor_Validation_Allele2', 'TEXT'),
-    ('Match_Norm_Validation_Allele1', 'TEXT'),
-    ('Match_Norm_Validation_Allele2', 'TEXT'),
-    ('Verification_Status', 'TEXT'),
-    ('Validation_Status', 'TEXT'),
-    ('Mutation_Status', 'TEXT'),
-    ('Sequencing_Phase', 'TEXT'),
-    ('Sequence_Source', 'TEXT'),
-    ('Validation_Method', 'TEXT'),
-    ('Score', 'TEXT'),
-    ('BAM_File', 'TEXT'),
-    ('Sequencer', 'TEXT'),
     ('Tumor_Sample_UUID', 'TEXT'),
     ('Matched_Norm_Sample_UUID', 'TEXT'),
     ('HGVSc', 'TEXT'),
@@ -161,7 +148,6 @@
     ('AF', 'TEXT'),
     ('AFR_AF', 'TEXT'),
     ('AMR_AF', 'TEXT'),
-    ('ASN_AF', 'TEXT'),
     ('EAS_AF', 'TEXT'),
     ('EUR_AF', 'TEXT'),
     ('SAS_AF', 'TEXT'),
@@ -180,31 +166,11 @@
     ('TSL', 'TEXT'),
     ('HGVS_OFFSET', 'TEXT'),
     ('PHENO', 'TEXT'),
-    ('MINIMISED', 'TEXT'),
-    ('ExAC_AF', 'TEXT'),
-    ('ExAC_AF_AFR', 'TEXT'),
-    ('ExAC_AF_AMR', 'TEXT'),
-    ('ExAC_AF_EAS', 'TEXT'),
-    ('ExAC_AF_FIN', 'TEXT'),
-    ('ExAC_AF_NFE', 'TEXT'),
-    ('ExAC_AF_OTH', 'TEXT'),
-    ('ExAC_AF_SAS', 'TEXT'),
     ('GENE_PHENO', 'TEXT'),
     ('FILTER', 'TEXT'),
     ('flanking_bps', 'TEXT'),
     ('vcf_id', 'TEXT'),
     ('vcf_qual', 'REAL'),
-    ('ExAC_AF_Adj', 'TEXT'),
-    ('ExAC_AC_AN_Adj', 'TEXT'),
-    ('ExAC_AC_AN', 'TEXT'),
-    ('ExAC_AC_AN_AFR', 'TEXT'),
-    ('ExAC_AC_AN_AMR', 'TEXT'),
-    ('ExAC_AC_AN_EAS', 'TEXT'),
-    ('ExAC_AC_AN_FIN', 'TEXT'),
-    ('ExAC_AC_AN_NFE', 'TEXT'),
-    ('ExAC_AC_AN_OTH', 'TEXT'),
-    ('ExAC_AC_AN_SAS', 'TEXT'),
-    ('ExAC_FILTER', 'TEXT'),
     ('gnomAD_AF', 'TEXT'),
     ('gnomAD_AFR_AF', 'TEXT'),
     ('gnomAD_AMR_AF', 'TEXT'),
@@ -215,9 +181,12 @@
     ('gnomAD_OTH_AF', 'TEXT'),
     ('gnomAD_SAS_AF', 'TEXT'),
     ('vcf_pos', 'INTEGER'),
+    ('HotSpotAllele', 'INTEGER'),
     ('VAF', 'REAL')
 ]
 
+common_cols = [col for col, type in maf_types]
+
 needed_cols = [
     'Hugo_Symbol',
     'Entrez_Gene_Id',
@@ -304,7 +273,9 @@
         # process the chunk
         chunk['VAF'] = (chunk['t_alt_count'] /
                         (chunk['t_ref_count'] + chunk['t_alt_count']))
-        if table_name not in ('strelka', 'lancet', 'consensus'):
+        if table_name in ('strelka', 'lancet', 'consensus'):
+            chunk = chunk[common_cols]
+        else:
             chunk = chunk[needed_cols]
         chunk.to_sql(table_name, con, if_exists='append')
     # create indexes

analyses/snv-callers/scripts/03-calculate_tmb.R

@@ -284,8 +284,8 @@ strelka_mutect_maf_df <- strelka_mutect_maf_df %>%
     ),
   by = "Tumor_Sample_Barcode"
   ) %>%
-  # Remove samples if they are labeled "Panel"
-  dplyr::filter(experimental_strategy != "Panel")
+  # Remove samples if they are not WGS or WXS
+  dplyr::filter(experimental_strategy %in% c("WGS", "WXS"))
 
 ############################# Set Up BED Files #################################
 # Make a data.frame of the unique BED file paths and their names