shell:
    "rm -r input; rm -r plink/*; rm -r data; rm *txt; rm accessory/*"

run:
    shell(f"{CMD_PREFIX} cat plink/{BASE}.genome.* | gzip > {BASE}.genome.gz; mv {BASE}.genome.gz plink")

run:
    shell(f"zcat {{input[1]}} | awk {{params.awk_string}}; " #Extract related individuals
    f"{CMD_PREFIX} Rscript scripts/filter_relateds.R -r accessory/related_IBS_ctrls.txt -o accessory/excluded_related_ctrls.txt; " #This script cycles through related pairs and keeps one 
    f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDprune_ctrls --remove accessory/excluded_related_ctrls.txt --keep-allele-order --make-bed --out plink/{BASE}_LDp_ctrls_IBDflt {plnk_rsrc(rule, plink_run_specs)}")

run:
    shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_ctrls_IBDflt --keep plink/DS/{BASE}_LDprune_ctrls_DS{{wildcards.round}}.txt --maf {config['sHWE']['maf']} "
          f"--geno 0.0 --keep-allele-order --make-bed --out plink/DS/{BASE}_LDprune_ctrls_DS{{wildcards.round}} {plnk_rsrc(rule, plink_run_specs)}")

run:
    with open('accessory/sHWE_downsampled_snps.txt', 'r') as tested: #count up the SNPs that were combined in previous rule.
        for num_tested, line in enumerate(tested): pass
    if num_tested < int(config['sHWE']['test_threshold']): #If not enough SNPs will be tested, print error message and DO NOT create output file.  The documentatino of the pipeline encourages the user to adjust the parameters of the downsampling (more rounds, fewer individuals, lower threshold, etc.)
        print(f"Failed Checkpoint\nReason: Fewer than {config['sHWE']['test_threshold']} unique SNPs downsampled for sHWE (number downsampled = {num_tested})\n"
              f"Increase the value set for \'downsample_rounds\' entry in the \'sHWE\' section of the config file or reduce the sHWE 'test_threshold'")
        shell("cp accessory/sHWE_downsampled_snps.txt accessory/.sHWE_downsampled_snps.txt; rm accessory/sHWE_downsampled_snps.txt") #This will remove output of prior rule, which will prompt the re-creation of accessory/sHWE_downsampled_snps.txt.  This is necessary to prompt the downsampling process when the user changes parameters of downsampling in config file.
    else:
        shell('touch accessory/.sHWE_downsample_pass.txt') #If pass, then create hidden file required by next rule

run:
    max_entropy = -9
    best_popnum = -9
    with open("sHWE/entropy.txt", 'w') as ent_file:
        for file in input: #Loop over input files and calculate entropy of each
            if not os.stat(file).st_size == 0: #Make sure file is not empty
                popnum = file.split("_")[-1] #Extract population number from filename
                dat = pd.read_table(file, header=None) # read in vector of p-values from sHWE
                binned = pd.cut(dat.iloc[:,0], params.bins).value_counts(sort=False) # bin and count vector values
                entropy = scipy.stats.entropy(binned.iloc[1:,]) # Calculate entropy of bin counts, excluding first bin that corrresponds to bin containing true positives (i.e. sites truly out of HWE)
                ent_file.write(f"{popnum}\t{entropy}\n")
                if entropy > max_entropy: #Store current value if better than prior best result.
                    best_popnum = popnum
                    max_entropy = entropy
            else:
                print(f"sHWE calculation failed for {popnum}; debug or remove this value from sHWE -> pop_nums in the config.yml file")
    assert(max_entropy != -9 and best_popnum != -9) #Make sure that we have successfully identified an optimum population number
    with open("sHWE/.optimum_popnum.txt", 'w') as outfile: outfile.write(str(best_popnum)) #Write the optimum population number to hidden file.

run:
    popfile = open(f"sHWE/.optimum_popnum.txt", 'r') #Grab the estimated population number from the hidden file
    popnum = popfile.readline()
    popfile.close()
    shell(f"{CMD_PREFIX} Rscript {CODE}/run_sHWE.R -p plink/DS/{BASE}_LDprune_ctrls_DS{{wildcards.f}} -d {popnum} -t {{params.threshold}} -o sHWE/{BASE}_sHWE-DS{{wildcards.f}} || touch sHWE/{BASE}_sHWE-DS{{wildcards.f}}")

run:
    file_list = [x for x in glob.glob(f"plink/DS/{BASE}_LDprune_ctrls_DS*.bim")
                 if os.path.exists(f"sHWE/{BASE}_sHWE-DS{int(x.strip('.bim').replace('DS', '').split('_')[-1])}.rmv.bim")] #Create file list of all successful runs of sHWE (the .rmv.bim file will be present if sHWE completed successfully)
    with open('accessory/.tested_file_list.txt', 'w') as tested: #Write file list to hidden file
        for file in file_list: tested.write(f"{file}\n")
    if file_list: shell(f"find plink/DS/ -name \'{BASE}_LDprune_ctrls_DS*.bim\' | grep -v -w accessory/.tested_file_list.txt | xargs cat | cut -d \" \" -f 2 | sort -T {os.getcwd()} | uniq > accessory/sHWE_tested_snps.txt") #This command prints markers from the input .bim files IFF they were successfully tested (in our file list)

run:
    with open('accessory/sHWE_tested_snps.txt', 'r') as tested: #Count up tested SNPs in file
        for num_tested, line in enumerate(tested): pass
    if num_tested < int(config['sHWE']['test_threshold']):
        print(f"Failed Checkpoint\nReason: Only {num_tested} markers were tested for sHWE, which is fewer than requested threshold ({config['sHWE']['test_threshold']})\n"
              f"Increase the value set for \'downsample_rounds\' entry in the \'sHWE\' section of the config file")
        shell("cp accessory/sHWE_tested_snps.txt accessory/.sHWE_tested_snps.txt; rm accessory/sHWE_tested_snps.txt; rm accessory/sHWE_downsampled_snps.txt")
    else:
        shell('touch accessory/.sHWE_pass.txt')

run:
    if config['IBD']['filter_relateds'] == 'true':
        shell(f"zcat {{input[1]}} | awk {{params.awk_string}}; {CMD_PREFIX} Rscript {CODE}/filter_relateds.R -r accessory/related_IBS_results.txt -o accessory/excluded_related_samples.txt; "
              f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_sHWE --remove accessory/excluded_related_samples.txt --keep-allele-order --make-bed --out plink/{BASE}_LDp_sHWE_IBDflt {plnk_rsrc(rule, plink_run_specs)}")
    else:
        shell(f"zcat {{input[1]}} | awk {{params.awk_string}}; {CMD_PREFIX} Rscript {CODE}/filter_relateds.R -r accessory/related_IBS_results.txt -o accessory/excluded_related_samples.txt; "
              f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_sHWE --keep-allele-order --make-bed --out plink/{BASE}_LDp_sHWE_IBDflt {plnk_rsrc(rule, plink_run_specs)}")

run:
    if config['match_controls'] == 'true': #Run control-matching and parse results into a format that PLINK can use to extract matched samples.
        shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_sHWE_IBDflt --allow-no-sex --read-genome {{input[1]}} --cluster cc --mcc 1 {config['control_number']} --out plink/{BASE} {plnk_rsrc(rule, plink_run_specs)}; "
              f"{CMD_PREFIX} " + "awk '{{if(NF!=2) for (i = 2; i <= NF; i++) print $i}}'" + f" plink/{BASE}.cluster1 " + f"| cut -d \"(\" -f 1 | sed 's/_/\\t/' > accessory/samples2.txt; "
              "awk '{{split($2,a,\"_\"); printf(\"%s\\n%s\\n\",$2,a[1])}}' accessory/samples2.txt | sed 's/#//g' > accessory/samples.txt") #The sed command is a holdover from dealing with pound signs in sample names in StJude data
        shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_sub_mflt --keep accessory/samples2.txt --make-bed --out input/{BASE}_CtlMat {plnk_rsrc(rule, plink_run_specs)}; "
              f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_sHWE_IBDflt --keep accessory/samples2.txt --make-bed --out input/{BASE}_PCA {plnk_rsrc(rule, plink_run_specs)}")
    else:
        shell("{CMD_PREFIX} awk '{{split($2,a,\"_\"); printf(\"%s\\n%s\\n\",$2,a[1])}}' plink/" + BASE + "_LDp_sHWE_IBDflt.fam | sed 's/#//g' > accessory/samples.txt")
        shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_sub_mflt --make-bed --out input/{BASE}_CtlMat {plnk_rsrc(rule, plink_run_specs)}; "
              f"{CMD_PREFIX} plink --bfile plink/{BASE}_LDp_sHWE_IBDflt --make-bed --out input/{BASE}_PCA {plnk_rsrc(rule, plink_run_specs)}")

run:
    shell(f"head -n 1 {BASE}_chr1.admixmap.txt > header.txt; tail -q -n +2  {BASE}_chr*.admixmap.txt > body.txt; cat header.txt body.txt > {BASE}.admixmap.txt; rm header.txt; rm body.txt")

run:
    cmd = f"{CMD_PREFIX} Rscript {CODE}/genesis_prep.R -p plink/{BASE}_sub_mflt -k {{input[1]}} -n {config['gwas']['pc_num']} -o input/{BASE} -t {config['kinship_threshold']}"
    if os.path.exists(config['covariate_file']): #Not tested
            cmd += f" -C {config['covariate_file']}"
    shell(cmd)

run:
    cmd = f"{CMD_PREFIX} Rscript {CODE}/genesis_gwas.R -g {{input[0]}} -p {{input[1]}} -k {{input[2]}} -n {config['gwas']['pc_num']} -C {config['gwas']['other_predictors']} -O {os.getcwd()} -o {BASE}.genesis.txt -c {{threads}}"
    shell(cmd)

run:
    with open("scripts/gather_report_data.sh", 'w') as report_cmds:
        report_line = f"echo \'rmarkdown::render(\"scripts/admixMap-report.Rmd\", output_file=\"{BASE}-Mapping-report.pdf\", " \
                      f"params=list(genesis_rslt=\"{BASE}.genesis.txt\", " \
                      f"fam_file=\"plink/{BASE}_sub_mflt.fam\", " \
                      f"input_directory=\"input/\", " \
                      f"blink_rslt=\"-9\", " \
                      f"gen_since_admix=\"{config['admixMapping']['gen_since_admixture']}\", " \
                      f"gwas_threshold=\"{config['gwas']['sig_threshold']}\", " \
                      f"ancestry_comps=\"{config['admixMapping']['ancestry_predictors']}\", " \ 
                      f"pruned_bim=\"plink/{BASE}_LDp_sHWE.bim\", " \
                      f"sHWE_entropy=\"sHWE/entropy.txt\", " \
                      f"sHWE_markers=\"accessory/sHWE_tested_snps.txt\", " \
                      f"ibd_ctrls=\"plink/{BASE}_LDp_ctrls_IBDflt.fam\", " \
                      f"eigenval=\"input/{BASE}.pcair.varprop\", " \
                      f"phenodat=\"input/{BASE}.genesis.pdat\", " \
                      f"kinplot=\"input/{BASE}.kinplot.png\", " \
                      f"mflt=\"accessory/missing_filter.txt\", " \
                      f"rulegraph_file=\"{BASE}-rulegraph.png\", "
        if f"{BASE}.admixmap.txt" in input:
            report_line += f"admixMap_rslt=\"{BASE}.admixmap.txt\", global_ancestry=\"{BASE}.globalancestry.txt\","
        else:
            report_line += f"admixMap_rslt=\"-9\", global_ancestry=\"-9\","
            shell(f"touch {BASE}.admixmap.sig.txt")
        if f"{BASE}.dos.genesis.txt" in input: report_line += f"dosage_rslt=\"{BASE}.dos.genesis.txt\","
        else:
            report_line += f"dosage_rslt=\"-9\","
            shell(f"touch {BASE}.dos.genesis.sig.txt")
        report_line += f"config_file=\"workflow/config.yml\"))\' | R --vanilla"
        report_cmds.write(report_line)
    shell(f"{CMD_PREFIX} sh scripts/gather_report_data.sh; mv scripts/{BASE}-Mapping-report.pdf {BASE}-Mapping-report.pdf; mv input/*kinplot.png figures")

run:
    shell(f"cut -f 1 {BASE}.genesis.sig.txt > {BASE}.genesis.sigID.txt")
    shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_sub_mflt --extract {BASE}.genesis.sigID.txt --keep-allele-order --freq case-control --allow-no-sex --recode vcf --out {BASE}.genesis.sig {plnk_rsrc(rule, plink_run_specs)}")
    cmd = f"{CMD_PREFIX} /usr/lib/jvm/java-8-openjdk-amd64/bin/java -Xmx8G -jar /snpEff/snpEff.jar ann -t GRCh37.75 {BASE}.genesis.sig.vcf | " \
          "awk '{{split($8,a,\"|\"); split(a[1],b,\"=\"); print $1,$2,$4,$5,a[4],a[3],a[2],b[2]}}' | grep -v \"#\" > {BASE}.genesis.sig.preannt.txt"
    shell(cmd)
    for soft in ['genesis']:
        cmd = f"{CMD_PREFIX} Rscript {CODE}/merge_annt.R -i {BASE}.{soft}.txt -s {BASE}.{soft}.sig.preannt.txt -r {config['annotation']['rsIDs']} -f {BASE}.{soft}.sig.frq.cc -o {BASE}.{soft}.sig.annt.txt"
        if os.path.exists(config['annotation']['typed_key']): cmd += f" -t {config['annotation']['typed_key']}"
        shell(cmd)

run:
    if os.path.exists(config['conditional_analysis']['snp_list']):
        if not os.path.exists('conditional-analysis'): os.mkdir('conditional-analysis')
        chrom, pos = f"{wildcards.marker}".split("-")
        shell(f"{CMD_PREFIX} plink --bfile plink/{BASE}_sub_mflt --snp {wildcards.marker.replace('-',':')} --recodeAD --out conditional-analysis/{wildcards.marker} {plnk_rsrc(rule, plink_run_specs)} || true")
        if os.path.exists(f"conditional-analysis/{wildcards.marker}.raw"):
            shell(f"{CMD_PREFIX} Rscript {CODE}/genesis_conditional_analysis.R -a {chrom} -b {pos} -d {config['conditional_analysis']['distance']} "
                  f"-m conditional-analysis/{wildcards.marker}.raw "
                  f"-g {{input[0]}} -p {{input[1]}} -k {{input[2]}} -n {config['gwas']['pc_num']} -C {config['gwas']['other_predictors']} "
                  f"-o conditional-analysis/{BASE}-{wildcards.marker}.genesis.txt")
        else:
            print(f"Did not find {wildcards.marker} in plink/{BASE}_sub_mflt")
            shell(f"touch conditional-analysis/{BASE}-{wildcards.marker}.genesis.txt")

run:
    if not os.path.exists('fine-mapping'): os.mkdir('fine-mapping')
    chrom, pos = f"{wildcards.marker}".split("-")
    shell(f"{CMD_PREFIX} plink2 --pfile input/{BASE}_sub --snp {wildcards.marker.replace('-',':')} --window {int(config['fine_mapping']['distance'])} "
          f"--export A --out fine-mapping/{wildcards.marker} {plnk_rsrc(rule, plink_run_specs)} || true")
    #NOTE: phenotype data is not specified in dosage input.  Need to grab and merge from data in plink folder.