relations.smake

##
## Relatetness related rules
##

#
# Detect IBD segemnts with germline
#

DEF_MIN_SEG_LEN=3.0

def min_seg_len():
    return config.get('min_seg_len', DEF_MIN_SEG_LEN)

def ersa_mask():
    return config.get('ersa_mask', ref_path('ersa-mask.tsv'))

rule germline:
    input:
        ped="{dir}/{file}.ped", 
        map="{dir}/{file}_interpolated.map"
    output:
        "{dir}_GRM/{file}.match"
    params:
        prefix="{dir}_GRM/{file}",
        min_seg=min_seg_len()
    shell:
        "germlinew -input {input.ped} {input.map} -output {params.prefix} -bits 128 -min_m {params.min_seg} -err_het 1 -err_hom 2 -g_extend -w_extend"

#
# Convert vcf to plinke format.
#

rule vcf2plink:
    input:
         "{dir}/{file}.vcf.gz"
    output:
         "{dir}/{file}.map",
         "{dir}/{file}.ped",
         "{dir}/{file}_interpolated.map"
    params:
        prefix="{dir}/{file}",
        gentic_map=ref_path('plink.chrALL.GRCh37.map.gz')
    shell:
        "gzvcf2plink -i {input} -o {params.prefix} -r {params.gentic_map}"

#
# Filter out polulation IBD segments using ERSA mask
#

rule FPI_filter_population_ibd:
    input:
        "{filename}.match.gz"
    output:
        "{filename}_FPI.match.gz"
    params:
        ersa_mask=ersa_mask()
    shell:
        "maskSegments -v -i {input} -o {output} -m {params.ersa_mask} -u cM"


#
# Estimate pairwise IBD0 and relatedness degree from IBD segments (in germline format)
#

rule estimateIBD0:
    input:
        "{filename}.match.gz"
    output:
        degrees="{filename}_IBD.csv",
        segments="{filename}_IBD-segments.match.gz"
    params:
        min_seg=min_seg_len()
    shell:
        "ibd2degree -v -i {input} -o {output.degrees} -s {output.segments} -m {params.min_seg}"

#
# Compare estimared degree relatedness vs reported (true) one.
#

rule compare_vs_true_rel:
    input:
        estRel="{filename}.csv", 
        trueRel=config['rel_true'],
        rmd=os.path.join(PWD, "notebooks/compare_vs_true_rel.Rmd")
    output:
        "{filename}_RVT.html"
    script:
        "notebooks/compare_vs_true_rel.Rmd"