Merge pull request #1156 from maxplanck-ie/allelic_whatshap_dna

Allelic whatshap dna

Author: katsikora

conda-recipe/meta.yaml

@@ -1,4 +1,4 @@
-{% set version = "3.4.0" %}
+{% set version = "3.5.0" %}
 
 package:
   name: snakepipes

docs/content/News.rst

@@ -1,6 +1,13 @@
 snakePipes News
 ===============
 
+snakePipes 3.5.0
+----------------
+
+* Allelic-whatshap mode was added to DNAmapping. Requires --phasedVcf.
+* Allelic-whatshap mode now works with alignment-free mode and fromBAM in the mRNAseq workflow.
+
+
 snakePipes 3.4.0
 ----------------
 

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 name = "snakePipes"
 description = 'Snakemake workflows and wrappers for NGS data processing from the MPI-IE'
 readme = "README.md"
-version = "3.4.0"
+version = "3.5.0"
 keywords = [
     "DNAmapping",
     "ChIPSeq",

snakePipes/shared/rules/multiQC.snakefile

@@ -58,6 +58,9 @@ def multiqc_input_check(return_value):
             infiles.append( expand("allelic_bams/{sample}.filtered.SNPsplit_report.yaml", sample = samples) )
             infiles.append( expand("allelic_bams/{sample}.filtered.SNPsplit_sort.yaml", sample = samples) )
             indir += "allelic_bams"
+        if "allelic-whatshap" in mode:
+            infiles.append( expand("allelic_bams/{sample}.{suffix}.sorted.bam",sample=samples,suffix=['allele_flagged', 'genome1', 'genome2', 'unassigned']))
+            indir += "allelic_bams"
     elif pipeline=="rnaseq":
         # must be RNA-mapping, add files as per the mode
         if ( "alignment" in mode or "deepTools_qc" in mode or "three-prime-seq" in mode ) and not "allelic-mapping" in mode and not "allelic-counting" in mode and not "allelic-whatshap" in mode:

snakePipes/workflows/DNAmapping/DNAmapping.py

@@ -26,7 +26,8 @@ def parse_args(defaults={"verbose": False, "configFile": None,
                          "alignerOpts": "", "mateOrientation": "--fr",
                          "UMIDedup": False, "UMIDedupOpts": "",
                          "UMIDedupSep": "_", "UMIBarcode": False, "cutntag": False,
-                         "bcPattern": "NNNNCCCCCCCC", "aligner":"Bowtie2"}):
+                         "bcPattern": "NNNNCCCCCCCC", "aligner":"Bowtie2",
+                         "pvcf": None}):
     """
     Parse arguments from the command line.
     """
@@ -47,7 +48,7 @@ def parse_args(defaults={"verbose": False, "configFile": None,
     optional.add_argument("-m", "--mode",
                           dest="mode",
                           help="workflow running modes (available: 'mapping,"
-                          "allelic-mapping')(default: '%(default)s')",
+                          "allelic-mapping, allelic-whatshap')(default: '%(default)s')",
                           default=defaults["mode"])
 
     parserCommon.commonOptions(optional, defaults)
@@ -113,6 +114,11 @@ def parse_args(defaults={"verbose": False, "configFile": None,
                           choices=["Bowtie2","bwa","bwa-mem2"],
                           default=defaults["aligner"])
 
+    optional.add_argument("--phasedVcf",
+                          dest="pvcf",
+                          help="Phased vcf required for whatshap haplotagging. (default: '%(default)s')",
+                          default=defaults["pvcf"])
+
     return parser
 
 

snakePipes/workflows/DNAmapping/Snakefile

@@ -84,6 +84,10 @@ else:
     elif aligner == "bwa-mem2":
         include: os.path.join(maindir, "shared", "rules", "bwa-mem2.snakefile")
 
+if "allelic-whatshap" in mode:
+    include: os.path.join(maindir, "shared", "rules", "whatshap.snakefile")
+    include: os.path.join(maindir, "shared", "rules", "deepTools_qc_allelic.snakefile")
+
 ### conditional/optional rules #################################################
 ################################################################################
 
@@ -126,12 +130,6 @@ def run_deepTools_qc():
     file_list.append(expand("deepTools_qc/estimateReadFiltering/{sample}_filtering_estimation.txt",sample = samples))
     return (file_list)
 
-def run_Qualimap():
-    file_list = []
-    if qualimap:
-        file_list += expand("Qualimap_qc/{sample}.filtered.bamqc_report.html", sample = samples)
-        file_list += expand("Qualimap_qc/{sample}.filtered.bamqc_results.txt", sample = samples)
-    return (file_list)
 # allele specific
 def make_nmasked_genome():
     if allele_mode == 'create_and_map':
@@ -160,6 +158,21 @@ def run_allelesp_mapping():
     else:
         return([])
 
+def run_whatshap():
+    if "allelic-whatshap" in mode:
+        allele_suffix = ['allele_flagged', 'genome1', 'genome2', 'unassigned']
+        file_list = [
+        expand("allelic_bams/{sample}.{suffix}.sorted.bam", sample = samples,
+                        suffix = allele_suffix),
+        expand("allelic_bams/{sample}.{suffix}.sorted.bam.bai", sample = samples,
+                        suffix = allele_suffix),
+        expand("bamCoverage/allele_specific/{sample}.{suffix}.seq_depth_norm.bw", sample = samples,
+                        suffix = ['genome1', 'genome2'])
+        ]
+        return(file_list)
+    else:
+        return([])
+
 ### execute before workflow starts #############################################
 ################################################################################
 onstart:
@@ -196,9 +209,9 @@ rule all:
         expand("bamCoverage/{sample}.filtered.seq_depth_norm.bw", sample = samples),
         run_computeGCBias(GCBias),
         run_deepTools_qc(),
-        run_Qualimap(),
         make_nmasked_genome(),
         run_allelesp_mapping(),
+        run_whatshap(),
         "multiQC/multiqc_report.html"
 
 ### execute after workflow finished ############################################

snakePipes/workflows/DNAmapping/defaults.yaml

@@ -64,3 +64,5 @@ UMIDedupOpts:
 fragmentLength: 200
 qualimap: false
 verbose: false
+#phased vcf
+pvcf:

snakePipes/workflows/DNAmapping/internals.snakefile

@@ -27,6 +27,15 @@ if trim:
 
 ### Initialization #############################################################
 
+if "allelic-whatshap" in mode:
+    if "allelic-mapping" in mode:
+        sys.exit("Allelic-mapping and allelic-whatshap modes are not compatible. Please choose one or another.")
+    if not pvcf:
+        sys.exit("Allelic-whatshap mode was specified but no phased vcf file was provided. Please provided the path to a phased vcf file.")
+    if not os.path.isfile(pvcf):
+        sys.exit(f"File {pvcf} doesn't exist.")
+
+
 infiles = sorted(glob.glob(os.path.join(str(indir or ''), '*'+ext)))
 if infiles == []:
     sys.exit("Error! Samples extnesion in {} are not {}. "
@@ -71,3 +80,4 @@ else:
     f = open(filter_rules, "w")
     f.write(filt)
     f.close()
+

snakePipes/workflows/mRNAseq/mRNAseq.py

@@ -103,7 +103,7 @@ def parse_args(defaults={"verbose": False, "configFile": None,
                           default=defaults["LRT"])
 
 
-    optional.add_argument("--phased-vcf",
+    optional.add_argument("--phasedVcf",
                           dest="pvcf",
                           help="Phased vcf required for whatshap haplotagging. (default: '%(default)s')",
                           default=defaults["pvcf"])

tests/test_jobcounts.py

@@ -549,6 +549,25 @@ def test_allelic_2strains(self,ifs):
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
         assert parseSpOut(_p) == 143
+    def test_whatshap_allelic(self, ifs):
+        ci = [
+            "DNAmapping",
+            '-i',
+            ifs / 'PE',
+            '-o',
+            ifs / 'outdir',
+            '--snakemakeOptions',
+            SMKOPTS,
+            ifs / 'org.yaml',
+            '-m',
+            'allelic-whatshap',
+            '--phasedVcf',
+            ifs / 'allelic_input' / 'file.vcf.gz'
+        ]
+        print(' '.join([str(i) for i in ci]))
+        _p = sp.run(ci, capture_output=True, text=True)
+        assert _p.returncode == 0
+        assert parseSpOut(_p) == 180
 
 class TestChIPseq:
     def test_default(self, ifs):
@@ -1832,7 +1851,7 @@ def test_whatshap_allelic(self, ifs):
             ifs / 'sampleSheet.tsv',
             '-m',
             'allelic-whatshap,deepTools_qc',
-            '--phased-vcf',
+            '--phasedVcf',
             ifs / 'allelic_input' / 'file.vcf.gz'
         ]
         print(' '.join([str(i) for i in ci]))
@@ -1856,7 +1875,7 @@ def test_whatshap_allelic_fromBAM(self, ifs):
             ifs / 'sampleSheet.tsv',
             '-m',
             'allelic-whatshap,deepTools_qc',
-            '--phased-vcf',
+            '--phasedVcf',
             ifs / 'allelic_input' / 'file.vcf.gz'
         ]
         print(' '.join([str(i) for i in ci]))