Merge ac87c9b into 379a442

.gitignore

@@ -69,4 +69,4 @@ vaccine-peptides-report.xlsx
 vaccine-peptides-report.json
 vaccine-peptides.csv
 neoepitope-report.xlsx
-
+all_variants.csv

run-vaxrank-b16-test-data.sh

@@ -15,6 +15,11 @@ vaxrank \
     --output-neoepitope-report neoepitope-report.xlsx \
     --output-json-file vaccine-peptides-report.json \
     --output-csv vaccine-peptides.csv \
+    --output-passing-variants-csv all_variants.csv \
+    --interferon-gamma-response-csv /Users/julia/code/gene-lists/interferon-gamma-response.csv \
+    --class1-mhc-presentation-pathway-csv /Users/julia/code/gene-lists/class1-mhc-presentation-pathway.csv \
+    --cancer-driver-genes-csv /Users/julia/code/gene-lists/cancer-driver-genes.csv \
+    --cancer-driver-variants-csv /Users/julia/code/gene-lists/cancer-driver-variants.csv \
     --output-reviewed-by "John Doe,Jane Doe" \
     --output-final-review "All the Does" \
     --output-patient-id "Test Patient"

test/data/b16.f10/b16.vcf

@@ -2,5 +2,6 @@
 #chr	pos	id	ref	alt	qual	filter	info
 chr4	45802539	.	G	C	.	.	.
 chr9	82927102	.	G	T	.	.	.
+chr11	101177505	.	T	C	.	.	.
 chr13	5864876	.	C	CG	.	.	.
-chrX	8125624	.	C	A	.	.	.
+chrX	8125624	.	C	A	.	.	.

test/test_mutant_protein_sequence.py

@@ -15,10 +15,11 @@
 from __future__ import absolute_import, print_function, division
 
 from nose.tools import eq_, assert_almost_equal
-from vaxrank.core_logic import ranked_vaccine_peptides
+from vaxrank.core_logic import VaxrankCoreLogic
 from mhctools import RandomBindingPredictor
 from isovar.cli.variant_sequences_args import make_variant_sequences_arg_parser
 from isovar.cli.rna_args import allele_reads_generator_from_args
+from varcode.cli import variant_collection_from_args
 
 from .testing_helpers import data_path
 
@@ -58,7 +59,9 @@ def test_mutant_amino_acids_in_mm10_chrX_8125624_refC_altA_pS460I():
         "--bam", data_path("b16.f10/b16.combined.sorted.bam"),
     ])
     reads_generator = allele_reads_generator_from_args(args)
-    ranked_list, _ = ranked_vaccine_peptides(
+    variants = variant_collection_from_args(args)
+    core_logic = VaxrankCoreLogic(
+        variants=variants,
         reads_generator=reads_generator,
         mhc_predictor=random_binding_predictor,
         vaccine_peptide_length=15,
@@ -67,6 +70,8 @@ def test_mutant_amino_acids_in_mm10_chrX_8125624_refC_altA_pS460I():
         min_alt_rna_reads=1,
         min_variant_sequence_coverage=1,
         variant_sequence_assembly=True)
+    core_logic.process_variants()
+    ranked_list = core_logic.ranked_vaccine_peptides()
 
     for variant, vaccine_peptides in ranked_list:
         eq_(
@@ -89,15 +94,19 @@ def test_mutant_amino_acids_in_mm10_chr9_82927102_refGT_altTG_pT441H():
         "--bam", data_path("b16.f10/b16.combined.sorted.bam"),
     ])
     reads_generator = allele_reads_generator_from_args(args)
-    ranked_list, _ = ranked_vaccine_peptides(
+    variants = variant_collection_from_args(args)
+    core_logic = VaxrankCoreLogic(
         reads_generator=reads_generator,
         mhc_predictor=random_binding_predictor,
+        variants=variants,
         vaccine_peptide_length=15,
         padding_around_mutation=5,
         min_alt_rna_reads=1,
         min_variant_sequence_coverage=1,
         variant_sequence_assembly=True,
         max_vaccine_peptides_per_variant=1)
+    core_logic.process_variants()
+    ranked_list = core_logic.ranked_vaccine_peptides()
 
     for variant, vaccine_peptides in ranked_list:
         vaccine_peptide = vaccine_peptides[0]
@@ -113,18 +122,21 @@ def test_keep_top_k_epitopes():
         "--bam", data_path("b16.f10/b16.combined.sorted.bam"),
     ])
     reads_generator = allele_reads_generator_from_args(args)
-
+    variants = variant_collection_from_args(args)
     keep_k_epitopes = 3
-    ranked_list, _ = ranked_vaccine_peptides(
+    core_logic = VaxrankCoreLogic(
         reads_generator=reads_generator,
         mhc_predictor=random_binding_predictor,
+        variants=variants,
         vaccine_peptide_length=15,
         padding_around_mutation=5,
         min_alt_rna_reads=1,
         min_variant_sequence_coverage=1,
         variant_sequence_assembly=True,
         max_vaccine_peptides_per_variant=1,
         num_mutant_epitopes_to_keep=keep_k_epitopes)
+    core_logic.process_variants()
+    ranked_list = core_logic.ranked_vaccine_peptides()
 
     for variant, vaccine_peptides in ranked_list:
         vaccine_peptide = vaccine_peptides[0]

test/test_shell_script.py

@@ -94,7 +94,7 @@ def test_pdf_report():
         assert getsize(f.name) > 0
 
 
-@patch('vaxrank.core_logic.vaccine_peptides_for_variant')
+@patch('vaxrank.core_logic.VaxrankCoreLogic.vaccine_peptides_for_variant')
 def test_report_no_peptides(mock_vaccine_peptides_for_variant):
     # simulate case where we have no epitopes for any variant
     mock_vaccine_peptides_for_variant.return_value = []

vaxrank/cli.py

@@ -28,10 +28,12 @@
     mhc_binding_predictor_from_args,
 )
 
+import pandas as pd
 import serializable
 from varcode.cli import variant_collection_from_args
 
-from .core_logic import ranked_vaccine_peptides
+from .core_logic import VaxrankCoreLogic
+from .gene_pathway_check import GenePathwayCheck
 from .report import (
     make_ascii_report,
     make_html_report,
@@ -59,6 +61,7 @@ def new_run_arg_parser():
     add_output_args(arg_parser)
     add_optional_output_args(arg_parser)
     add_supplemental_report_args(arg_parser)
+    add_gene_list_args(arg_parser)
     return arg_parser
 
 
@@ -76,6 +79,7 @@ def cached_run_arg_parser():
     add_output_args(arg_parser)
     add_optional_output_args(arg_parser)
     add_supplemental_report_args(arg_parser)
+    add_gene_list_args(arg_parser)
     return arg_parser
 
 
@@ -180,6 +184,11 @@ def add_output_args(arg_parser):
         type=int,
         help="Number of mutations to report")
 
+    output_args_group.add_argument(
+        "--output-passing-variants-csv",
+        default="",
+        help="Path to CSV file containing some metadata about every passing variant")
+
 
 def add_vaccine_peptide_args(arg_parser):
     vaccine_peptide_group = arg_parser.add_argument_group("Vaccine peptide options")
@@ -218,6 +227,24 @@ def add_supplemental_report_args(arg_parser):
         default="",
         help="Local path to COSMIC vcf")
 
+def add_gene_list_args(arg_parser):
+    gene_list_args_group = arg_parser.add_argument_group("Interesting gene list file paths")
+    gene_list_args_group.add_argument(
+        "--interferon-gamma-response-csv",
+        default="",
+        help="Local path to interferon-gamma response CSV file")
+    gene_list_args_group.add_argument(
+        "--class1-mhc-presentation-pathway-csv",
+        default="",
+        help="Local path to MHC class I presentation pathway CSV file")
+    gene_list_args_group.add_argument(
+        "--cancer-driver-genes-csv",
+        default="",
+        help="Local path to cancer driver genes CSV file")
+    gene_list_args_group.add_argument(
+        "--cancer-driver-variants-csv",
+        default="",
+        help="Local path to cancer driver variants CSV file")
 
 def check_args(args):
     if not (args.output_csv or
@@ -272,7 +299,15 @@ def ranked_variant_list_with_metadata(args):
     reads_generator = allele_reads_generator_from_args(args)
     mhc_predictor = mhc_binding_predictor_from_args(args)
 
-    ranked_list, variants_count_dict = ranked_vaccine_peptides(
+    gene_pathway_check = GenePathwayCheck(
+        interferon_gamma_response_csv=args.interferon_gamma_response_csv,
+        class1_mhc_presentation_pathway_csv=args.class1_mhc_presentation_pathway_csv,
+        cancer_driver_genes_csv=args.cancer_driver_genes_csv,
+        cancer_driver_variants_csv=args.cancer_driver_variants_csv
+    )
+
+    core_logic = VaxrankCoreLogic(
+        variants=variants,
         reads_generator=reads_generator,
         mhc_predictor=mhc_predictor,
         vaccine_peptide_length=args.vaccine_peptide_length,
@@ -282,21 +317,32 @@ def ranked_variant_list_with_metadata(args):
         min_variant_sequence_coverage=args.min_variant_sequence_coverage,
         min_epitope_score=args.min_epitope_score,
         num_mutant_epitopes_to_keep=args.num_epitopes_per_peptide,
-        variant_sequence_assembly=args.variant_sequence_assembly)
+        variant_sequence_assembly=args.variant_sequence_assembly,
+        gene_pathway_check=gene_pathway_check
+    )
+    core_logic.process_variants()
 
-    ranked_list_for_report = ranked_list[:args.max_mutations_in_report]
+    variants_count_dict = core_logic.variant_counts()
+    assert len(variants) == variants_count_dict['num_total_variants'], \
+        "Len(variants) is %d but variants_count_dict came back with %d" % (len(variants), variants_count_dict['num_total_variants'])
 
+    if args.output_passing_variants_csv:
+        variant_metadata_dicts = core_logic.variant_properties()
+        df = pd.DataFrame(variant_metadata_dicts)
+        df.to_csv(args.output_passing_variants_csv, index=False)
+
+    ranked_list = core_logic.ranked_vaccine_peptides()
+    ranked_list_for_report = ranked_list[:args.max_mutations_in_report]
     patient_info = PatientInfo(
         patient_id=args.output_patient_id,
         vcf_paths=variants.sources,
         bam_path=args.bam,
         mhc_alleles=mhc_alleles,
-        num_somatic_variants=len(variants),
+        num_somatic_variants=variants_count_dict['num_total_variants'],
         num_coding_effect_variants=variants_count_dict['num_coding_effect_variants'],
         num_variants_with_rna_support=variants_count_dict['num_variants_with_rna_support'],
         num_variants_with_vaccine_peptides=variants_count_dict['num_variants_with_vaccine_peptides']
     )
-
     # return variants, patient info, and command-line args
     data = {
         'variants': ranked_list_for_report,