Add build of results

cnvpipelines.py

@@ -792,7 +792,8 @@ class CnvPipeline:
 
     def run_simulation(self, nb_inds, reference, nstretches, sv_list, coverage, force_polymorphism,
                        haploid, proba_del, proba_inv, read_len, insert_len_mean, insert_len_sd, min_deletions,
-                       min_inversions, max_try, genotypes, tools, species, max_n_stretches, force_wdir=False, **kwargs):
+                       min_inversions, max_try, genotypes, tools, species, max_n_stretches,
+                       overlap_cutoff, left_precision, right_precision, force_wdir=False, **kwargs):
         """
         Run simulation
         :param nb_inds:
@@ -886,7 +887,10 @@ class CnvPipeline:
                 "no_detection_index": True,
                 "build_refbundle": "genomestrip" in tools,
                 "species": species,
-                "maxNstretches": max_n_stretches
+                "maxNstretches": max_n_stretches,
+                "overlap_cutoff": overlap_cutoff,
+                "left_precision": left_precision,
+                "right_precision": right_precision
             }
 
             if "genomestrip" in tools:
@@ -1178,6 +1182,9 @@ if __name__ == "__main__":
                                   help="[refbundle] Species name, according to the NCBI Taxonomy database")
     simul_parser.add_argument('-mn', '--max-n-stretches', type=int, required=False, default=100,
                                   help="[refbundle] Max size of nstretches to consider them as it")
+    simul_parser.add_argument('--overlap_cutoff', type=float, default=0.5, help='[results] cutoff for reciprocal overlap')
+    simul_parser.add_argument('--left-precision', type=int, default=-1, help='[results] left breakpoint precision')
+    simul_parser.add_argument('--right-precision', type=int, default=-1, help='[results] right breakpoint precision')
     add_run_and_rerun_options(simul_parser)
     simul_parser.add_argument('-w', '--wdir', type=str, required=True,
                               help="Output folder where data will be stored")

popsim @ 1903b654

-Subproject commit 3a817c8b952b74d2056a9ee41be30554c746f7df
+Subproject commit 1903b654d5395769a250cce7228c1e28aa6a86ec

snakecnv/popsim.snk

@@ -26,6 +26,9 @@ MIN_DELETIONS = config["min_deletions"]
 MIN_INVERSIONS = config["min_inversions"]
 MAX_TRY = config["max_try"]
 GENOTYPES = config["genotypes"] if "genotypes" in config else None
+OVERLAP_CUTOFF = config["overlap_cutoff"]
+LEFT_PRECISION = config["left_precision"]
+RIGHT_PRECISION = config["right_precision"]
 
 
 # Functions here
@@ -70,7 +73,9 @@ localrules: all_popsim
 
 rule all_popsim:
     input:
-        set_all_inputs
+        expand("results/{svtype}/Summarized_results_{svtype}.ipynb", svtype=config["varianttypes"]),
+        expand("results/{svtype}/Summarized_results_{svtype}.html", svtype=config["varianttypes"]),
+        #set_all_inputs
         #expand("bams/{sample}.bam", sample=["indiv%d" % i for i in range(1, NB_INDS+1)])
         #*bp_outputs
 
@@ -78,7 +83,8 @@ rule buildpop:
     input:
         **bp_inputs
     output:
-        *bp_outputs
+        *bp_outputs,
+        os.path.join("pop", "genotypes.vcf.gz")
     params:
         nb_inds = NB_INDS,
         out_dir = "pop",
@@ -134,4 +140,59 @@ rule linkdatabams:
         bai="data/bams/{sample}.bam.bai"
     run:
         os.symlink(os.path.abspath(input.bam), output.bam)
-        os.symlink(os.path.abspath(input.bai), output.bai)
\ No newline at end of file
+        os.symlink(os.path.abspath(input.bai), output.bai)
+
+rule buildinfilesresults:
+    input:
+        genotypes = expand("{batch}/genotypes/{{svtype}}/svtyper_{chrom}_{{svtype}}_genotypes.vcf.gz",
+                           chrom=chromosomes,
+                           batch=samples.keys()),
+        filtered = expand("{batch}/filtered/{{svtype}}/svtyper_{chrom}_{{svtype}}_genotypes_filtered.vcf.gz",
+                          chrom=chromosomes,
+                          batch=samples.keys())
+    output:
+        genotypes = "list_files/tools_results_files_{svtype}.list",
+        filtered = "list_files/filtered_results_files_{svtype}.list"
+    threads:
+        1
+    run:
+        with open(output.genotypes, "w") as o_gt:
+            for i_gt in input.genotypes:
+                o_gt.write("%s\n" % i_gt)
+
+rule buildresults:
+    input:
+        genotypes = "list_files/tools_results_files_{svtype}.list",
+        filtered = "list_files/filtered_results_files_{svtype}.list",
+        true_vcf = "pop/genotypes.vcf.gz"
+    output:
+        "results/{svtype}/Summarized_results_{svtype}.ipynb",
+        "results/{svtype}/Summarized_results_{svtype}.html"
+    params:
+        overlap_cutoff = OVERLAP_CUTOFF,
+        left_precision = LEFT_PRECISION,
+        right_precision = RIGHT_PRECISION,
+        svlist = SV_LIST,
+        haploid = HAPLOID
+    threads:
+        1
+    log:
+        stdout = "logs/results_{svtype}.log"
+    run:
+        command = ["build_results.py",
+                   "-v", input.genotypes,
+                   "-t", input.true_vcf,
+                   "-f", input.filtered_vcf,
+                   "-y", wildcards.svtype,
+                   "--overlap_cutoff", str(params.overlap_cutoff),
+                   "--left-precision", str(params.overlap_cutoff),
+                   "--right-precision", str(params.right_precision),
+                   "-o", os.path.join("results", wildcards.svtype),
+                   "--no-xls"]
+
+        if params.haploid:
+            command.append("--haploid")
+        if params.svlist is not None:
+            command += ["-r", params.svlist]
+
+        shell(" ".join(command) + "1> %s 2> %s" % (log.stdout, log.stderr))