clair3 defaults

EichlerLab · Nov 17, 2023 · be97c44 · be97c44
1 parent 800aac7
commit be97c44
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 4 deletions.
diff --git a/202201-6445530 b/202201-6445530
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
@@ -34,6 +34,17 @@ def find_clair_chrs(wildcards):
         chrom=chroms,
     )
 
+def find_clair_gvcf(wildcards):
+    if config.get("CHRS") == None:
+        with open(f"{REF}.fai", "r") as infile:
+            chroms = [line.split("\t")[0] for line in infile]
+    else:
+        chroms = config.get("CHRS")
+    return expand(
+        "tmp/variants/{{sample}}/{chrom}/merge_output.gvcf.gz",
+        chrom=chroms,
+    )
+
 
 def concatenate_fastq(wildcards):
     FOFN = manifest_df.at[wildcards.sample, "FOFN"]
@@ -64,4 +75,4 @@ def find_aln_bai(wildcards):
         return rules.index_aln.output.merged_bai
     else:
         print("phase state must be either minimap2 (unphased) or longphase (phased)")
-        exit(1)
+        exit(1)
diff --git a/workflow/rules/variant_calling.smk b/workflow/rules/variant_calling.smk
@@ -6,31 +6,39 @@ rule clair_chr:
         ref=REF,
     output:
         vcf="tmp/variants/{sample}/{chrom}/phased_merge_output.vcf.gz",
+        gvcf="tmp/variants/{sample}/{chrom}/merge_output.gvcf.gz",
     log:
         "log/{sample}_{chrom}.clair.log",
     conda:
         "../envs/clair3.yaml"
     threads: 8
     resources:
-        mem=2,
+        mem=32,
         hrs=24,
         disk_free=1,
     shell:
         """
-        run_clair3.sh --bam_fn={input.merged_bam} --sample_name={wildcards.sample} --ref_fn={input.ref} --threads={threads} --platform=ont --model_path=$(dirname $( which run_clair3.sh  ) )/models/ont_guppy5 --output=$(dirname {output.vcf}) --ctg_name={wildcards.chrom} --enable_phasing
+        run_clair3.sh --bam_fn={input.merged_bam} --sample_name={wildcards.sample} --ref_fn={input.ref} --threads={threads} --platform=ont --model_path=$(dirname $( which run_clair3.sh  ) )/models/ont_guppy5 --output=$(dirname {output.vcf}) --ctg_name={wildcards.chrom} --enable_phasing --gvcf
         if [[ ( -f $( echo {output.vcf} | sed 's/phased_//' ) ) && ( ! -f {output.vcf} ) ]]; then
             cp $( echo {output.vcf} | sed 's/phased_//' ) {output.vcf}
         elif [[ ( ! -f $( echo {output.vcf} | sed 's/phased_//' ) ) && ( ! -f {output.vcf} ) ]]; then
             zcat $(dirname {output.vcf})/pileup.vcf.gz | grep ^# | bgzip -c > {output.vcf}
         fi
+        if [[ ( -f $( echo {output.vcf} | sed 's/phased_//' ) ) && ( ! -f {output.gvcf} ) ]]; then
+            cp $( echo {output.vcf} | sed 's/phased_//' ) {output.gvcf}
+        elif [[ ( ! -f $( echo {output.vcf} | sed 's/phased_//' ) ) && ( ! -f {output.gvcf} ) ]]; then
+            zcat $(dirname {output.vcf})/pileup.vcf.gz | grep ^# | bgzip -c > {output.gvcf}
+        fi
         """
 
 
 rule concat_clair:
     input:
         vcf=find_clair_chrs,
+        gvcf=find_clair_gvcf,
     output:
         vcf=temp("variants/{sample}/{sample}.clair3.vcf"),
+        gvcf="variants/{sample}/{sample}.clair3.gvcf.gz",
     envmodules:
         "modules",
         "modules-init",
@@ -48,6 +56,7 @@ rule concat_clair:
     shell:
         """
         bcftools concat -O v -o {output.vcf} {input.vcf}
+        bcftools concat -O v {input.gvcf} | bcftools sort -O v | bgzip -c > {output.gvcf} 
         """
 
 
@@ -69,7 +78,7 @@ rule sniffles:
     conda:
         "../envs/sniffles.yaml"
     resources:
-        mem=10,
+        mem=24,
         hrs=24,
         disk_free=1,
     threads: 1