Snakefile_whatshap_WG

import glob
import peds 
import uuid

### configure
output_dir="output"
wgs_interval= "ref/hg38.wgs_interval.bed"
hg38_ref="ref/Homo_sapiens_assembly38.fasta"
ref_dict=os.path.splitext(hg38_ref)[0]+'.dict'

### Define trios
ped_dir = "ped_files"
ped_files = glob.glob(ped_dir + "/*.ped")

families = {}
for fn in ped_files:
    f=peds.open_ped(fn)[0]
    families[f.id]=f

fam_ids = [line.strip() for line in open("8trios.lst", "r")]

child_ids = [[person.id for person in families[fid] if families[fid].get_father(person) ][0] for fid in fam_ids]

CHILD_DICT=dict(zip(fam_ids,child_ids))
# print (CHILD_DICT)
# exit(0)

rule all:
    input:
        expand(output_dir +"/phase_child/{fam}_child.phased.vcf.gz", fam=fam_ids),
        expand(output_dir +"/phase_trios/{fam}.phased.vcf.gz", fam=fam_ids)

rule phase_child:
    input: 
        bams=lambda w: expand(output_dir +"/cram/{id}.cram", id=[CHILD_DICT[w.fam]]),
        bais=lambda w: expand(output_dir +"/cram/{id}.cram.crai", id=[CHILD_DICT[w.fam]]),
        ref=hg38_ref,
        vcf=output_dir +"/glnexus/{fam}.dv_combined.vcf.gz"
    output: 
        vcf= output_dir +"/phase_child/{fam}_child.phased.vcf.gz"
    benchmark:
        output_dir +"/benchmark/phase_child/{fam}.tsv"
    resources: 
        threads=2,
        mem_mb = 40*1000,
        runtime= "7d"
    shell: """
        whatshap phase -o {output.vcf} --tag=PS --indels --reference={input.ref} {input.vcf} {input.bams}
    """

rule phase_trios:
    input: 
        bams=lambda w: expand(output_dir +"/cram/{id}.cram", id=[person.id for person in families[w.fam]]),
        bais=lambda w: expand(output_dir +"/cram/{id}.cram.crai", id=[person.id for person in families[w.fam]]),
        ped=ped_dir + "/{fam}.ped",
        ref=hg38_ref,
        vcf=output_dir +"/glnexus/{fam}.dv_combined.vcf.gz"
    output: 
        vcf= output_dir +"/phase_trios/{fam}.phased.vcf.gz"
    benchmark:
        output_dir +"/benchmark/phase_trios/{fam}.tsv"
    resources: 
        threads=2,
        mem_mb = 160*1000,
        runtime= "7d"
    shell: """
        whatshap phase -o {output.vcf} --tag=PS --indels --ped {input.ped} --reference={input.ref} {input.vcf} {input.bams} 
    """

### add whatshap stats