fix: create tnscope mnvs

BALSAMIC/constants/cluster_analysis.json

@@ -388,7 +388,7 @@
 		"time": "01:00:00",
 		"n": 4
 	},
-	"modify_tnscope_infofield": {
+	"post_process_tnscope": {
 		"time": "01:00:00",
 		"n": 4
 	},

BALSAMIC/constants/rules.py

@@ -200,8 +200,9 @@
     "vep_annotate_germlineVAR_tumor",
     "vep_annotate_germlineVAR_normal",
     # SNVs
-    "modify_tnscope_infofield",
-    "modify_tnscope_infofield_umi",
+    "bcftools_split_tnscope_variants",
+    "sentieon_tnscope_umi",
+    "sentieon_tnscope_umi_tn",
     "gatk_update_vcf_sequence_dictionary",
     "bcftools_filter_tnscope_clinical_tumor_only",
     "bcftools_filter_tnscope_clinical_tumor_normal",

BALSAMIC/constants/variant_filters.py

@@ -329,6 +329,9 @@ class WgsSNVFilters(BaseSNVFilters):
 
 class TgaSNVFilters(BaseSNVFilters):
     research = [
+        VCFFilter(
+            filter_name="MERGED", Description="SNV Merged with neighboring variants"
+        ),
         VCFFilter(tag_value=0.01, filter_name="SWEGENAF", field="INFO"),
         VCFFilter(tag_value=0.005, filter_name="balsamic_high_pop_freq", field="INFO"),
     ]

BALSAMIC/constants/workflow_params.py

@@ -115,7 +115,7 @@
     },
 }
 
-SLEEP_BEFORE_START = 600
+SLEEP_BEFORE_START = 800
 
 WORKFLOW_PARAMS = {
     "bam_post_processing": {

BALSAMIC/containers/varcall_py3/varcall_py3.yaml

@@ -104,6 +104,7 @@ dependencies:
   - pthread-stubs=0.4
   - pycosat=0.6.4
   - pycparser=2.20
+  - pyfaidx=0.8.1.3
   - pyopenssl=20.0.1
   - pysam=0.19.1
   - pysocks=1.7.1

BALSAMIC/snakemake_rules/umi/modify_tnscope_infofield_umi.rule

@@ -1,15 +1,14 @@
 
 rule modify_tnscope_infofield_umi:
     input:
-        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.preprocess.vcf.gz",
+        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
     output:
-        vcf_tnscope_umi = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
+        vcf_tnscope_umi = vcf_dir + "sentieon_tnscope/SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.post_process.vcf.gz",
     benchmark:
         Path(benchmark_dir,'modify_tnscope_infofield_umi_' + config[ "analysis" ][ "case_id" ] + ".tsv").as_posix()
     singularity:
         Path(singularity_image, config["bioinfo_tools"].get("bcftools") + ".sif").as_posix()
     params:
-        housekeeper_id = {"id": config["analysis"]["case_id"], "tags": "research"},
         modify_tnscope_infofield = get_script_path("modify_tnscope_infofield.py"),
         tmpdir = tempfile.mkdtemp(prefix=tmp_dir),
         case_name = config["analysis"]["case_id"]

BALSAMIC/snakemake_rules/umi/sentieon_varcall_tnscope.rule

@@ -12,12 +12,13 @@ rule sentieon_tnscope_umi:
         bed = config["panel"]["capture_kit"],
         dbsnp = config["reference"]["dbsnp"]
     output:
-        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.preprocess.vcf.gz",
+        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
         namemap = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.sample_name_map"
     benchmark:
         Path(benchmark_dir, "sentieon_tnscope_umi_" + config["analysis"]["case_id"] + ".tsv").as_posix()
     params:
         tmpdir = tempfile.mkdtemp(prefix=tmp_dir),
+        housekeeper_id = {"id": config["analysis"]["case_id"],"tags": "research"},
         sentieon_exec = config_model.sentieon.sentieon_exec,
         sentieon_lic = config_model.sentieon.sentieon_license,
         tumor_af = params.tnscope_umi.filter_tumor_af,

BALSAMIC/snakemake_rules/umi/sentieon_varcall_tnscope_tn.rule

@@ -12,12 +12,13 @@ rule sentieon_tnscope_umi_tn:
         bed = config["panel"]["capture_kit"],
         dbsnp = config["reference"]["dbsnp"]
     output:
-        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.preprocess.vcf.gz",
+        vcf_tnscope_umi = vcf_dir +  "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
         namemap = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.sample_name_map"
     benchmark:
         Path(benchmark_dir, "sentieon_tnscope_umi_" + config["analysis"]["case_id"] + ".tsv").as_posix()
     params:
         tmpdir = tempfile.mkdtemp(prefix=tmp_dir),
+        housekeeper_id = {"id": config["analysis"]["case_id"],"tags": "research"},
         sentieon_exec = config_model.sentieon.sentieon_exec,
         sentieon_lic = config_model.sentieon.sentieon_license,
         tumor_af = params.tnscope_umi.filter_tumor_af,

BALSAMIC/snakemake_rules/variant_calling/snv_quality_filter.rule

@@ -86,7 +86,7 @@ if config["analysis"]["sequencing_type"] == 'targeted' and config["analysis"]["a
     input:
       vcf = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.vcf.gz",
     output:
-      vcf_filtered = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.vcf.gz",
+      vcf_filtered = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.pre_process.vcf.gz",
     benchmark:
       Path(benchmark_dir,'bcftools_quality_filter_vardict_tumor_only_' + config["analysis"]["case_id"] + ".tsv").as_posix()
     singularity:
@@ -155,7 +155,7 @@ elif config["analysis"]["sequencing_type"] == 'targeted' and config["analysis"][
     input:
       vcf = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.vcf.gz",
     output:
-      vcf_filtered = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.vcf.gz",
+      vcf_filtered = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.pre_process.vcf.gz",
     benchmark:
       Path(benchmark_dir,'bcftools_quality_filter_vardict_tumor_normal_' + config["analysis"]["case_id"] + ".tsv").as_posix()
     singularity:
@@ -223,7 +223,7 @@ elif config["analysis"]["sequencing_type"] == 'targeted' and config["analysis"][
 if config_model.analysis.analysis_workflow == AnalysisWorkflow.BALSAMIC_UMI and config["analysis"]["analysis_type"] == 'paired':
   rule bcftools_quality_filter_TNscope_umi_tumor_normal:
     input:
-      vcf = vcf_dir +  "SNV.somatic."+ config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
+      vcf = vcf_dir + "sentieon_tnscope/SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.post_process.vcf.gz",
     output:
       vcf_filtered = vcf_dir +  "SNV.somatic."+ config["analysis"]["case_id"] + ".tnscope_umi.research.vcf.gz",
     benchmark:
@@ -251,9 +251,9 @@ if config_model.analysis.analysis_workflow == AnalysisWorkflow.BALSAMIC_UMI and
 elif config_model.analysis.analysis_workflow == AnalysisWorkflow.BALSAMIC_UMI and config["analysis"]["analysis_type"] == 'single':
   rule bcftools_quality_filter_TNscope_umi_tumor_only:
     input:
-      vcf=vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.vcf.gz",
+      vcf = vcf_dir + "sentieon_tnscope/SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.post_process.vcf.gz",
     output:
-      vcf_filtered=vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope_umi.research.vcf.gz"
+      vcf_filtered = vcf_dir +  "SNV.somatic."+ config["analysis"]["case_id"] + ".tnscope_umi.research.vcf.gz"
     benchmark:
       Path(benchmark_dir,'bcftools_quality_filter_TNscope_umi_tumor_only' + config["analysis"][
         "case_id"] + ".tsv").as_posix()

BALSAMIC/snakemake_rules/variant_calling/tnscope_post_process.rule

@@ -5,16 +5,16 @@
 
 rule bcftools_split_tnscope_variants:
     input:
-        ref = config["reference"]["reference_genome"],
         vcf = vcf_dir + "sentieon_tnscope/ALL.somatic." + config["analysis"]["case_id"] + ".tnscope.vcf.gz",
     output:
-        vcf_tnscope = vcf_dir + "sentieon_tnscope/SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.preprocess.vcf",
+        vcf_tnscope = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.vcf.gz",
         vcf_tnscope_sv = vcf_dir + "SV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.vcf.gz",
     benchmark:
         Path(benchmark_dir,'bcftools_split_tnscope_variants_' + config[ "analysis" ][ "case_id" ] + ".tsv").as_posix()
     singularity:
         Path(singularity_image, config["bioinfo_tools"].get("bcftools") + ".sif").as_posix()
     params:
+        housekeeper_id = {"id": config["analysis"]["case_id"], "tags": "research"},
         tmpdir = tempfile.mkdtemp(prefix=tmp_dir),
         case_name = config["analysis"]["case_id"]
     threads:
@@ -26,40 +26,52 @@ rule bcftools_split_tnscope_variants:
 export TMPDIR={params.tmpdir};
 mkdir -p {params.tmpdir};
 
-bcftools view --include 'INFO/SVTYPE=="."' -o {output.vcf_tnscope} {input.vcf} ; 
+bcftools view --include 'INFO/SVTYPE=="."' -O z -o {output.vcf_tnscope} {input.vcf} ; 
 bcftools view --include 'INFO/SVTYPE!="."' -O z -o {output.vcf_tnscope_sv} {input.vcf}; 
 tabix -p vcf -f {output.vcf_tnscope_sv};
+tabix -p vcf -f {output.vcf_tnscope};
         """
 
-rule modify_tnscope_infofield:
+rule post_process_tnscope:
     input:
-        vcf_tnscope = vcf_dir + "sentieon_tnscope/SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.preprocess.vcf",
+        vcf_tnscope = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.pre_process.vcf.gz",
+        ref = config["reference"]["reference_genome"],
     output:
-        vcf_tnscope = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.vcf.gz",
+        vcf_tnscope = vcf_dir + "SNV.somatic." + config["analysis"]["case_id"] + ".tnscope.research.vcf.gz",
     benchmark:
-        Path(benchmark_dir,'modify_tnscope_infofield_' + config[ "analysis" ][ "case_id" ] + ".tsv").as_posix()
+        Path(benchmark_dir,'post_process_tnscope_' + config[ "analysis" ][ "case_id" ] + ".tsv").as_posix()
     singularity:
         Path(singularity_image, config["bioinfo_tools"].get("bcftools") + ".sif").as_posix()
     params:
-        housekeeper_id = {"id": config["analysis"]["case_id"], "tags": "research"},
+        merge_mnvs = get_script_path("merge_mnp.py"),
         modify_tnscope_infofield = get_script_path("modify_tnscope_infofield.py"),
+        edit_vcf_script= get_script_path("edit_vcf_info.py"),
         tmpdir = tempfile.mkdtemp(prefix=tmp_dir),
         case_name = config["analysis"]["case_id"],
-        edit_vcf_script = get_script_path("edit_vcf_info.py"),
+        sentieon_exec = config_model.sentieon.sentieon_exec,
+        sentieon_lic = config_model.sentieon.sentieon_license,
+        matched_normal_filternames = ",".join(BaseSNVFilters.MATCHED_NORMAL_FILTER_NAMES),
         variant_caller= "tnscope"
     threads:
-        get_threads(cluster_config, 'modify_tnscope_infofield')
+        get_threads(cluster_config, 'post_process_tnscope')
     message:
-        "Add DP and AF tumor sample info and FOUND_IN to INFO field for case: {params.case_name}"
+        "Merge TNscope SNVs with same phaseID to MNVs."
+        "Add DP and AF tumor sample info and FOUND_IN to INFO field: {params.case_name}"
     shell:
         """
-export TMPDIR={params.tmpdir};
-mkdir -p {params.tmpdir};
+export SENTIEON_TMPDIR={params.tmpdir};
+export SENTIEON_LICENSE={params.sentieon_lic};
+
+{params.sentieon_exec} pyexec {params.merge_mnvs} --preserve_filters {params.matched_normal_filternames} --max_distance 5 {input.vcf_tnscope} {input.ref} > {params.tmpdir}/tnscope.research.mnv.vcf ;
+
+python {params.modify_tnscope_infofield} {params.tmpdir}/tnscope.research.mnv.vcf {params.tmpdir}/tnscope.research.mnv.add_info_fields.vcf ;
 
-python {params.modify_tnscope_infofield} {input.vcf_tnscope} {params.tmpdir}/vcf_tnscope_snvs_modified.vcf ;
-python {params.edit_vcf_script} -i {params.tmpdir}/vcf_tnscope_snvs_modified.vcf -o {params.tmpdir}/vcf_tnscope_snvs_modified_found_in_added.vcf -c {params.variant_caller};
-bgzip {params.tmpdir}/vcf_tnscope_snvs_modified_found_in_added.vcf ;
-mv {params.tmpdir}/vcf_tnscope_snvs_modified_found_in_added.vcf.gz {output.vcf_tnscope} ;
+python {params.edit_vcf_script} -i {params.tmpdir}/tnscope.research.mnv.add_info_fields.vcf -o {params.tmpdir}/tnscope.research.mnv.add_info_fields.added_found_in.vcf -c {params.variant_caller};
+
+bgzip {params.tmpdir}/tnscope.research.mnv.add_info_fields.added_found_in.vcf ;
+
+mv {params.tmpdir}/tnscope.research.mnv.add_info_fields.added_found_in.vcf.gz {output.vcf_tnscope} ;
 tabix -p vcf -f {output.vcf_tnscope} ;
         """
 
+

BALSAMIC/workflows/balsamic.smk

@@ -20,6 +20,7 @@ from BALSAMIC.constants.analysis import (
 from BALSAMIC.constants.paths import BALSAMIC_DIR
 from BALSAMIC.constants.rules import SNAKEMAKE_RULES
 from BALSAMIC.constants.variant_filters import (
+    BaseSNVFilters,
     SVDB_FILTER_SETTINGS,
     MANTA_FILTER_SETTINGS,
     WgsSNVFilters,

CHANGELOG.rst

@@ -1,28 +1,32 @@
 [X.X.X]
--------
+--------
 
 Added:
 ^^^^^^
 * Added option to disable hard filter of variants in matched normal https://github.com/Clinical-Genomics/BALSAMIC/pull/1509
+* Added check to verify sample sex for all workflows https://github.com/Clinical-Genomics/BALSAMIC/pull/1516
+
 
 Changed:
 ^^^^^^^^
+* Reworked bcftools filters https://github.com/Clinical-Genomics/BALSAMIC/pull/1509
+* Renamed high_normal_tumor_af_frac to in_normal https://github.com/Clinical-Genomics/BALSAMIC/pull/1509
+* check to verify sample sex for all workflows https://github.com/Clinical-Genomics/BALSAMIC/pull/1516
+* Merging SNVs into MNVs in TNscope TGA https://github.com/Clinical-Genomics/BALSAMIC/pull/1524
+* Change raw delivery SNV file for TGA to before any post-processing https://github.com/Clinical-Genomics/BALSAMIC/pull/1524
 
 
 Removed:
 ^^^^^^^^
 * Remove WGS-level GC-bias metric from TGA workflow https://github.com/Clinical-Genomics/BALSAMIC/pull/1521
 
 
-Changed:
-^^^^^^^^
-* Reworked bcftools filters https://github.com/Clinical-Genomics/BALSAMIC/pull/1509
-* Renamed high_normal_tumor_af_frac to in_normal https://github.com/Clinical-Genomics/BALSAMIC/pull/1509
-* check to verify sample sex for all workflows https://github.com/Clinical-Genomics/BALSAMIC/pull/1516
+Fixed:
+^^^^^^
 
 
 [16.0.0]
--------
+--------
 
 Added:
 ^^^^^^

docs/balsamic_filters.rst

@@ -245,6 +245,7 @@ The `TNscope <https://www.biorxiv.org/content/10.1101/250647v1.abstract>`_ algor
 *min_base_qual*: Minimal base quality to consider in calling
 
 ::
+
     min_base_qual = 15
 
 *min_tumor_allele_frac*: Set the minimum tumor AF to be considered as potential variant site.
@@ -313,6 +314,33 @@ The `TNscope <https://www.biorxiv.org/content/10.1101/250647v1.abstract>`_ algor
     marks variant with soft-filter `in_normal` variant if: AF(normal) / AF(tumor) > 0.3
 
 
+**Post-processing of TNscope variants**
+
+After quality-filtering TNscope variants and before merging with VarDict variants the phased SNVs and InDels from TNscope are merged together to MNVs using a slightly modified script from `Sentieon-scripts <https://github.com/Sentieon/sentieon-scripts/blob/master/merge_mnp/merge_mnp.py>`_ which can be found in ``BALSAMIC/assets/scripts/merge_mnp.py``
+
+This was done to avoid multiple representations of the same variant as VarDict already outputs these types of variants as MNVs, and because VEP isn't coded to handle phased SNVs in the interpretation of protein effect.
+
+In the merging of phased SNVs to MNV we need to handle how to consolidate information from multiple variants into a single metric, and importantly also for the FILTER column.
+
+An example is a MNV created by merging a phased germline SNV with a somatic SNV. This has been solved as follows:
+
+- `MNV_CONFLICTING_FILTERS`: Is a filter given to MNVs with constituent variants with different filters (such as `in_normal` and `PASS`)
+
+.. note::
+
+    However, as we may have multiple filters which means similar things, such as germline_risk and in_normal, MNVs constituted by variants with only these filters set aren't exactly "conflicting".
+
+Therefore the logic for setting `MNV_CONFLICTING_FILTERS` has been made a bit more complex, and in summary there are 3 possible outcomes for filters when merging SNVs/InDels into MNVs:
+
+1. Single filter such as PASS, when all constituting variants all have the same filter and no other.
+2. Multiple filters, such as in_normal,germline_risk, when all constituting variants have at least 1 of the matched normal filters.
+3. `MNV_CONFLICTING_FILTERS` when the merged variants have conflicting filters, and they don't all contain matched normal filters.
+
+.. note::
+
+    In addition to this a few more fields are added to the INFO field of the created MNVs containing comma-separated lists of AD, AF, and FILTER from its constituting variants.
+
+
 **Post-call Observation database Filters**
 
 

docs/balsamic_pon.rst

@@ -8,7 +8,7 @@ Currently two PON-methods are implemented in BALSAMIC to correct for biases and
 - To produce normalised CN-profiles for WGS cases visualised in ``GENS``.
 
 Sharing PON for publications
-======================
+============================
 
 If a PON has been used for the analysis of samples in a research project and a publication requires that the PON is uploaded to some database, a request can be made to Clinical Genomics, and depending on the status of the consent of the individuals from which the samples used in the construction of the PON has been derived it may or may not be possible to share the PON.
 

docs/balsamic_sv_cnv.rst

@@ -58,7 +58,7 @@ It is mandatory to provide the gender of the sample from BALSAMIC version >= 10.
 Further details about a specific caller can be found in the links for the repositories containing the documentation for SV and CNV callers along with the links for the articles are listed in `bioinfo softwares <https://balsamic.readthedocs.io/en/latest/bioinfo_softwares.html>`_.
 
 **Difficult to detect clinically relevant SVs**
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 **IGH::DUX4 rearrangements**
 

docs/user_guide.rst

@@ -5,7 +5,7 @@ Short tutorial
 Here a short tutorial is provided for BALSAMIC (**version** = 16.0.0).
 
 Regarding fastq-inputs
----------------------
+---------------------------
 
 Previous versions of BALSAMIC only accepted one fastq-pair per sample, which required concatenation of fastq-pairs if multiple existed.