fix collect sites

Zilong-Li · Jun 22, 2023 · 7c2fefd · 7c2fefd
1 parent a33527a
commit 7c2fefd
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 20 deletions.
diff --git a/workflow/rules/accuracy.smk b/workflow/rules/accuracy.smk
@@ -2,11 +2,7 @@
 rule collect_truth_gts:
     """would be better to use sites in subrefs"""
     input:
-        sites=lambda wildcards: expand(
-            rules.subset_refpanel.output.sites,
-            size=config["refsize"],
-            allow_missing=True,
-        ),
+        sites=rules.concat_refpanel_sites_by_region2.output.sites,
     output:
         gt=os.path.join(OUTDIR_PANEL, "truth.gts.{chrom}.txt"),
         af=os.path.join(OUTDIR_PANEL, "af.input.panel.{chrom}.txt"),

diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
@@ -54,7 +54,7 @@ def get_all_results():
 
 def get_subset_refs():
     return expand(
-        rules.subset_refpanel.output,
+        rules.subset_refpanel_by_region2.output,
         chrom=config["chroms"],
         size=config["refsize"],
     )

diff --git a/workflow/rules/glimpse.smk b/workflow/rules/glimpse.smk
@@ -1,7 +1,7 @@
 
 rule glimpse2_prepare_panel:
     input:
-        refvcf=rules.subset_refpanel.output.vcf,
+        refvcf=rules.subset_refpanel_by_region2.output.vcf,
     output:
         os.path.join(OUTDIR_GLIMPSE2, "bin", "{chrom}.size{size}.chunks{chunkid}.bin"),
     log:
@@ -149,12 +149,12 @@ rule glimpse_prepare_glvcf:
     input:
         bams=rules.bamlist.output,
         sites=lambda wildcards: expand(
-            rules.subset_refpanel.output.sites,
+            rules.subset_refpanel_by_region2.output.sites,
             size=config["refsize"],
             allow_missing=True,
         ),
         tsv=lambda wildcards: expand(
-            rules.subset_refpanel.output.tsv,
+            rules.subset_refpanel_by_region2.output.tsv,
             size=config["refsize"],
             allow_missing=True,
         ),
@@ -184,7 +184,7 @@ rule glimpse_prepare_glvcf:
 
 rule glimpse_phase:
     input:
-        refvcf=rules.subset_refpanel.output.vcf,
+        refvcf=rules.subset_refpanel_by_region2.output.vcf,
         glvcf=rules.glimpse_prepare_glvcf.output.vcf,
     output:
         temp(

diff --git a/workflow/rules/quilt.smk b/workflow/rules/quilt.smk
@@ -1,9 +1,9 @@
 
 rule quilt_prepare_regular:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
-        # hap=rules.subset_refpanel.output.hap,
-        # leg=rules.subset_refpanel.output.leg,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
+        # hap=rules.subset_refpanel_by_region2.output.hap,
+        # leg=rules.subset_refpanel_by_region2.output.leg,
     output:
         os.path.join(
             OUTDIR_PANEL,
@@ -72,7 +72,7 @@ rule quilt_prepare_regular:
 
 rule quilt_prepare_mspbwt:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
     output:
         os.path.join(
             OUTDIR_PANEL,
@@ -144,7 +144,7 @@ rule quilt_prepare_mspbwt:
 
 rule quilt_prepare_zilong:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
     output:
         os.path.join(
             OUTDIR_PANEL,
@@ -221,7 +221,7 @@ rule quilt_prepare_zilong:
 
 rule quilt_run_regular:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
         bams=rules.bamlist.output,
         rdata=rules.quilt_prepare_regular.output,
     output:
@@ -325,7 +325,7 @@ rule quilt_ligate_regular:
 
 rule quilt_run_mspbwt:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
         bams=rules.bamlist.output,
         rdata=rules.quilt_prepare_mspbwt.output,
     output:
@@ -429,7 +429,7 @@ rule quilt_ligate_mspbwt:
 
 rule quilt_run_zilong:
     input:
-        vcf=rules.subset_refpanel.output.vcf,
+        vcf=rules.subset_refpanel_by_region2.output.vcf,
         bams=rules.bamlist.output,
         rdata=rules.quilt_prepare_zilong.output,
     output:

diff --git a/workflow/rules/refpanels.smk b/workflow/rules/refpanels.smk
@@ -44,7 +44,7 @@ rule subset_refpanel_by_chrom:
         """
 
 
-rule subset_refpanel:
+rule subset_refpanel_by_region2:
     input:
         rules.subset_sample_list.output.samples,
     output:
@@ -63,7 +63,7 @@ rule subset_refpanel:
             OUTDIR_PANEL, "panelsize{size}", "vcfs" "{chrom}.{start}.{end}.tsv.vcf.gz"
         ),
     params:
-        N="subset_refpanel",
+        N="subset_refpanel_by_region2",
         prefix=lambda wildcards, output: os.path.splitext(output[0])[0],
         vcf=lambda wildcards: REFPANEL[wildcards.chrom]["vcf"],
         start=lambda wildcards: max(1, wildcards.start - config["glimpse"]["buffer"]),
@@ -83,3 +83,28 @@ rule subset_refpanel:
             bcftools query -f'%CHROM\t%POS\t%REF,%ALT\n' {output.sites} | bgzip -c > {output.tsv} && tabix -s1 -b2 -e2 {output.tsv}
         )  &> {log}
         """
+
+
+rule concat_refpanel_sites_by_region2:
+    input:
+        rules.subset_refpanel_by_region2.output.sites,
+    output:
+        sites=os.path.join(
+            OUTDIR_PANEL,
+            "panelsize{size}",
+            "vcfs" "{chrom}.sites.vcf.gz",
+        ),
+        tsv=os.path.join(OUTDIR_PANEL, "panelsize{size}", "vcfs" "{chrom}.tsv.vcf.gz"),
+    log:
+        os.path.join(
+            OUTDIR_PANEL, "panelsize{size}", "vcfs" "{chrom}.sites.vcf.gz.llog"
+        ),
+    conda:
+        "../envs/pandas.yaml"
+    shell:
+        """
+        ( \
+            bcftools concat -D --threads 4 -Oz -o {output.sites} {input} && bcftools index -f {output.sites} \
+            bcftools query -f'%CHROM\t%POS\t%REF,%ALT\n' {output.sites} | bgzip -c > {output.tsv} && tabix -s1 -b2 -e2 {output.tsv}
+        ) & > {log}
+        """