Snakefile

###############################################################################
#                                  Parameters                                 #
###############################################################################

# NOTE: Setting this to "True" will recalibrate CAPS, which is time-
# and resource-intensive. Re-use the files in "model/" if possible.
recalibrate_CAPS = False

colors2 = ["#999899", "#B11116"]
colors3 = ["#999899", "#B11116", "#F04B5D"]

# Pairs of contexts (defined by their mutability values) to compare
# MAPS and CAPS on, based on maximum differences in the residuals from
# the MAPS model
mutability_pairs = [
    # Non-CpGs
    (3.1245e-09, 7.7544e-09),
    (6.7244e-09, 9.9720e-09),
    (1.0799e-09, 1.0244e-09),
    (9.4166e-10, 2.9242e-09),
    (2.9409e-09, 8.8520e-09),
    (2.4639e-09, 4.9426e-09),
    # CpGs
    (1.4169e-07, 1.0005e-07),
    (1.0750e-07, 7.2389e-08),
]
mutability_pairs1, mutability_pairs2 = zip(*mutability_pairs)

###############################################################################
#                                   Scripts                                   #
###############################################################################

caps_model = "../model/greta_model.R"
calculate_caps = "calculate_caps.R"
calculate_maps = "calculate_maps.R"
compare_scores = "compare_scores.R"
viz_scores = "viz_scores.R"
viz_model = "viz_model.R"
counts = "counts.R"
contexts = "contexts.R"
filtering = "filtering.R"
join_groups = "join_groups.R"
viz_scores_multi = "viz_scores_multi.R"
annotate_with_mutability_bands = "annotate_with_mutability_bands.R"
Whiffin2020_uORFs_preprocessing = "Whiffin2020_uORFs_preprocessing.R"
filter_alphamissense = "filter_alphamissense.R"
calculate_avg_alphamissense = "calculate_avg_alphamissense.R"

###############################################################################
#                                 Input files                                 #
###############################################################################

# General #####################################################################
if config["gcp"] == False:
    # Synonymous variants grouped by LOEUF bin
    syn_constrained = "../files/syn_constrained.tsv"
    # Mutability table
    mutation_ht = "../files/mutation_ht.tsv"
    # All variants grouped by their most severe consequence
    by_csq_unprocessed = "../files/by_csq.tsv"
    # All variants grouped by their most severe consequence (genomes)
    by_csq_genomes_unprocessed = "../files/by_csq_genomes.tsv"
    # Near-splice variants
    near_splice_vars_unprocessed = "../files/near_splice_vars.tsv"
    near_splice_vars_genomes_unprocessed = "../files/near_splice_vars_genomes.tsv"
    # Missense variants annotated with some AlphaMissense information
    missense_vars_with_alphamissense_unprocessed = (
        "../files/missense_vars_with_alphamissense.tsv"
    )
    # Whiffin2020
    syn_by_context_Whiffin2020 = "../files/syn_by_context_Whiffin2020.tsv"
    intronic_by_context_Whiffin2020 = "../files/intronic_by_context_Whiffin2020.tsv"
else:
    from snakemake.remote.GS import RemoteProvider

    GS = RemoteProvider()
    # Synonymous variants grouped by LOEUF bin
    syn_constrained = GS.remote(config["gcp_rootdir"] + "syn_constrained.tsv")
    # Mutability table
    mutation_ht = GS.remote(config["gcp_rootdir"] + "mutation_ht.tsv")
    # All variants grouped by their most severe consequence
    by_csq_unprocessed = GS.remote(config["gcp_rootdir"] + "by_csq.tsv")
    # All variants grouped by their most severe consequence (genomes)
    by_csq_genomes_unprocessed = GS.remote(config["gcp_rootdir"] + "by_csq_genomes.tsv")
    # Near-splice variants
    near_splice_vars_unprocessed = GS.remote(
        config["gcp_rootdir"] + "near_splice_vars.tsv"
    )
    near_splice_vars_genomes_unprocessed = GS.remote(
        config["gcp_rootdir"] + "near_splice_vars_genomes.tsv"
    )
    # Missense variants annotated with some AlphaMissense information
    missense_vars_with_alphamissense_unprocessed = GS.remote(
        config["gcp_rootdir"] + "missense_vars_with_alphamissense.tsv"
    )
    # Whiffin2020
    syn_by_context_Whiffin2020 = GS.remote(
        config["gcp_rootdir"] + "syn_by_context_Whiffin2020.tsv"
    )
    intronic_by_context_Whiffin2020 = GS.remote(
        config["gcp_rootdir"] + "intronic_by_context_Whiffin2020.tsv"
    )

by_csq = "by_csq.tsv"
by_csq_genomes = "by_csq_genomes.tsv"

near_acceptor_vars = "near_acceptor_vars.tsv"
near_acceptor_genomes_vars = "near_acceptor_genomes_vars.tsv"
near_donor_vars = "near_donor_vars.tsv"
near_donor_genomes_vars = "near_donor_genomes_vars.tsv"

missense_vars_with_alphamissense = "missense_vars_with_alphamissense.tsv"

Whiffin2020_stop_unprocessed = (
    "uORFs/data_files/stop-removing_all_possible_annotated.txt"
)
Whiffin2020_uAUG_unprocessed = (
    "uORFs/data_files/uAUG-creating_all_possible_annotated.txt"
)

# Model #######################################################################
# Synonymous variants with 4 missing contexts added
syn_by_context_p4 = "../model/phat.tsv"
syn_by_context_ppd = "../model/phat_sim.tsv"
# Synonymous variants with 4 missing contexts added (genomes)
syn_by_context_p4_genomes = "../model/phat_genomes.tsv"
syn_by_context_ppd_genomes = "../model/phat_sim_genomes.tsv"
# Synonymous variants with 4 missing contexts added (Whiffin2020 QC)
syn_by_context_p4_Whiffin2020 = "../model/phat_Whiffin2020.tsv"
syn_by_context_ppd_Whiffin2020 = "../model/phat_sim_Whiffin2020.tsv"

###############################################################################
#                               Auxiliary files                               #
###############################################################################

# General #####################################################################
# Synonymous variants grouped by context/mutability (exomes)
syn_by_context_and_coverage = "syn_by_context_and_coverage.tsv"
# Intronic variants grouped by context/mutability (exomes)
intronic_by_context_and_coverage = "intronic_by_context_and_coverage.tsv"
# Synonymous variants grouped by context/mutability (genomes)
syn_by_context_and_coverage_genomes = "syn_by_context_and_coverage_genomes.tsv"
# Intronic variants grouped by context/mutability (genomes)
intronic_by_context_and_coverage_genomes = (
    "intronic_by_context_and_coverage_genomes.tsv"
)
# CAPS scores for near-splice variants
caps_near_acceptor_scores = "near_acceptor_caps.tsv"
caps_near_acceptor_genomes_scores = "near_acceptor_genomes_caps.tsv"
caps_near_donor_scores = "near_donor_caps.tsv"
caps_near_donor_genomes_scores = "near_donor_genomes_caps.tsv"
# CAPS scores for all major functional classes of SNVs
caps_by_csq_scores = "all_csq_caps.tsv"
caps_by_csq_genomes_scores = "all_csq_genomes_caps.tsv"
# CAPS-PPD scores for all major functional classes of SNVs
caps_ppd_by_csq_scores = "all_csq_caps_ppd.tsv"
caps_ppd_by_csq_genomes_scores = "all_csq_genomes_caps_ppd.tsv"
# Synonymous variants annotated with mutability bands
syn_w_mutability_bands = "syn_w_mutability_bands.tsv"
# All variants annotated with mutability bands
all_w_mutability_bands = "all_w_mutability_bands.tsv"
# Synonymous variants by mutability band (CAPS score)
syn_vars_caps_by_mutability_band_scores = "syn_mutband_caps.tsv"
# Synonymous variants by mutability band (MAPS score)
syn_vars_maps_by_mutability_band_scores = "syn_mutband_maps.tsv"
# All variant by mutability band (MAPS score)
all_vars_maps_by_mutability_band_scores = "all_mutband_maps.tsv"
# All variant by mutability band (MAPS score, square-root transformed mutability model)
all_vars_sqrt_by_mutability_band_scores = "all_mutband_sqrt.tsv"
# All variant by mutability band (CAPS score)
all_vars_caps_by_mutability_band_scores = "all_mutband_caps.tsv"
# Near-splice variants (MAPS score)
maps_near_acceptor_scores = "near_acceptor_maps.tsv"
maps_near_acceptor_genomes_scores = "near_acceptor_genomes_maps.tsv"
maps_near_donor_scores = "near_donor_maps.tsv"
maps_near_donor_genomes_scores = "near_donor_genomes_maps.tsv"
# All variants by most severe consequences (MAPS score)
maps_by_csq_scores = "all_csq_maps.tsv"
maps_by_csq_genomes_scores = "all_csq_genomes_maps.tsv"
# Synonymous variants by Top % LOEUF group
syn_constrained_joined_groups = "syn_constrained_joined_groups.tsv"
syn_constrained_tops_caps_scores = "syn_constrained_caps.tsv"
syn_constrained_tops_maps_scores = "syn_constrained_maps.tsv"
# Model #######################################################################
theta_sample = "../model/theta_sample.tsv"
theta_sample_genomes = "../model/theta_sample_genomes.tsv"
theta_sample_Whiffin2020 = "../model/theta_sample_Whiffin2020.tsv"
# Whiffin2020 #################################################################
Whiffin2020_stop = "stop-removing_all_possible_annotated.tsv"
Whiffin2020_uAUG = "uAUG-creating_all_possible_annotated.tsv"
Whiffin2020_uORFs_effects_stop_maps_scores = (
    "Whiffin2020_uORFs_effects_stop_maps_scores.tsv"
)
Whiffin2020_uORFs_effects_stop_caps_scores = (
    "Whiffin2020_uORFs_effects_stop_caps_scores.tsv"
)
Whiffin2020_uORFs_effects_uAUG_maps_scores = (
    "Whiffin2020_uORFs_effects_uAUG_maps_scores.tsv"
)
Whiffin2020_uORFs_effects_uAUG_caps_scores = (
    "Whiffin2020_uORFs_effects_uAUG_caps_scores.tsv"
)
Whiffin2020_uORFs_kozak_stop_maps_scores = (
    "Whiffin2020_uORFs_kozak_stop_maps_scores.tsv"
)
Whiffin2020_uORFs_kozak_stop_caps_scores = (
    "Whiffin2020_uORFs_kozak_stop_caps_scores.tsv"
)
Whiffin2020_uORFs_kozak_uAUG_maps_scores = (
    "Whiffin2020_uORFs_kozak_uAUG_maps_scores.tsv"
)
Whiffin2020_uORFs_kozak_uAUG_caps_scores = (
    "Whiffin2020_uORFs_kozak_uAUG_caps_scores.tsv"
)
Whiffin2020_uORFs_distanceCDS_uAUG_maps_scores = (
    "Whiffin2020_uORFs_distanceCDS_uAUG_maps_scores.tsv"
)
Whiffin2020_uORFs_distanceCDS_uAUG_caps_scores = (
    "Whiffin2020_uORFs_distanceCDS_uAUG_caps_scores.tsv"
)

###############################################################################
#                                 Output files                                #
###############################################################################

# Counts by variant class
counts_by_csq = "all_csq_counts.tex"
# Counts by variant class (genomes)
counts_by_csq_genomes = "all_csq_counts_genomes.tex"
# Contexts by variant class
contexts_by_csq = "all_csq_contexts.tex"
# Contexts by variant class (genomes)
contexts_by_csq_genomes = "all_csq_contexts_genomes.tex"
# CAPS and MAPS scores for near-splice variants
scores2_near_acceptor_plot = "near_acceptor_scores2.pdf"
scores2_near_donor_plot = "near_donor_scores2.pdf"
scores2_near_acceptor_genomes_plot = "near_acceptor_genomes_scores2.pdf"
scores2_near_donor_genomes_plot = "near_donor_genomes_scores2.pdf"
# CAPS and MAPS scores for all major functional classes of SNVs
scores2_by_csq_plot = "all_csq_scores2.pdf"
scores2_by_csq_genomes_plot = "all_csq_scores2_genomes.pdf"
scores3_by_csq_plot = "all_csq_scores3.pdf"
scores3_by_csq_genomes_plot = "all_csq_scores3_genomes.pdf"
# Plot of CAPS scores for all major functional classes of SNVs
caps_by_csq_plot = "all_csq_caps.pdf"
# CAPS-PPD scores for all major functional classes of SNVs
caps_ppd_by_csq_plot = "all_csq_caps_ppd.pdf"
# Synonymous variants by mutability band (CAPS score)
syn_vars_caps_by_mutability_band_plot = "syn_mutband_caps.pdf"
# Synonymous variants by mutability band (MAPS score)
syn_vars_maps_by_mutability_band_plot = "syn_mutband_maps.pdf"
# All variant by mutability band (MAPS score)
all_vars_maps_by_mutability_band_plot = "all_mutband_maps.pdf"
# All variant by mutability band (MAPS score, square-root transformed mutability model)
all_vars_sqrt_by_mutability_band_plot = "all_mutband_sqrt.pdf"
# All variant by mutability band (CAPS score)
all_vars_caps_by_mutability_band_plot = "all_mutband_caps.pdf"
# CAPS and MAPS scores by mutability band
all_vars_scores2_by_mutability_band_plot = "all_mutband_scores2.pdf"
# All variants by most severe consequences (MAPS score)
maps_by_csq_plot = "all_csq_maps.pdf"
# Table summarising the difference in scores of the 3 models
scores3_by_csq_scores_comparison = "scores3_by_csq_scores_comparison.tex"
# Synonymous variants by Top % LOEUF group
syn_constrained_caps_plot = "syn_constrained_caps.pdf"
syn_constrained_scores2_plot = "syn_constrained_scores2.pdf"
# Fit of the MAPS model calibrated on synonymous variants
MAPS_model_fit_on_syn = "MAPS_model_fit_on_syn.pdf"
MAPS_model_fit_on_syn_genomes = "MAPS_model_fit_on_syn_genomes.pdf"
# Whiffin2020 #################################################################
Whiffin2020_uORFs_effects_uAUG_plot = "Whiffin2020_uORFs_effects_uAUG.pdf"
Whiffin2020_uORFs_effects_stop_plot = "Whiffin2020_uORFs_effects_stop.pdf"
Whiffin2020_uORFs_kozak_uAUG_plot = "Whiffin2020_uORFs_kozak_uAUG.pdf"
Whiffin2020_uORFs_distanceCDS_uAUG_plot = "Whiffin2020_uORFs_distanceCDS_uAUG.pdf"
Whiffin2020_uORFs_kozak_stop_plot = "Whiffin2020_uORFs_kozak_stop.pdf"


rule all:
    input:
        # Misc files ##########################################################
        counts_by_csq,
        counts_by_csq_genomes,
        contexts_by_csq,
        contexts_by_csq_genomes,
        # Plots ###############################################################
        MAPS_model_fit_on_syn,
        MAPS_model_fit_on_syn_genomes,
        caps_by_csq_plot,
        all_vars_scores2_by_mutability_band_plot,
        all_vars_maps_by_mutability_band_plot,
        all_vars_sqrt_by_mutability_band_plot,
        all_vars_caps_by_mutability_band_plot,
        scores2_near_acceptor_plot,
        scores2_near_donor_plot,
        scores2_near_acceptor_genomes_plot,
        scores2_near_donor_genomes_plot,
        scores2_by_csq_plot,
        scores3_by_csq_plot,
        scores3_by_csq_genomes_plot,
        scores2_by_csq_genomes_plot,
        expand(
            "CAPSvsMAPS_for_contexts_{mu1}_{mu2}.pdf",
            zip,
            mu1=mutability_pairs1,
            mu2=mutability_pairs2,
        ),
        syn_constrained_caps_plot,
        syn_constrained_scores2_plot,
        Whiffin2020_uORFs_kozak_uAUG_plot,
        Whiffin2020_uORFs_distanceCDS_uAUG_plot,
        Whiffin2020_uORFs_kozak_stop_plot,
        Whiffin2020_uORFs_effects_uAUG_plot,
        Whiffin2020_uORFs_effects_stop_plot,
        # Scores ##############################################################
        scores3_by_csq_scores_comparison,


# QC ##########################################################################


rule process_by_csq_variants:
    input:
        In=by_csq_unprocessed,
    output:
        Out=by_csq,
    params:
        min_coverage=30,
        add_essential_splice_cat=True,
    script:
        filtering


rule process_by_csq_variants_genomes:
    input:
        In=by_csq_genomes_unprocessed,
    output:
        Out=by_csq_genomes,
    params:
        min_coverage=30,
        add_essential_splice_cat=True,
    script:
        filtering


rule get_near_acceptor_variants:
    input:
        In=near_splice_vars_unprocessed,
    output:
        Out=near_acceptor_vars,
    params:
        min_coverage=30,
        variable="splice_type",
        variable_value="splice_acceptor_variant",
    script:
        filtering


rule get_near_acceptor_variants_genomes:
    input:
        In=near_splice_vars_genomes_unprocessed,
    output:
        Out=near_acceptor_genomes_vars,
    params:
        min_coverage=30,
        variable="splice_type",
        variable_value="splice_acceptor_variant",
    script:
        filtering


rule get_near_donor_variants:
    input:
        In=near_splice_vars_unprocessed,
    output:
        Out=near_donor_vars,
    params:
        min_coverage=30,
        variable="splice_type",
        variable_value="splice_donor_variant",
    script:
        filtering


rule get_near_donor_variants_genomes:
    input:
        In=near_splice_vars_genomes_unprocessed,
    output:
        Out=near_donor_genomes_vars,
    params:
        min_coverage=30,
        variable="splice_type",
        variable_value="splice_donor_variant",
    script:
        filtering

rule process_missense_vars_with_alphamissense:
    input:
        In=missense_vars_with_alphamissense_unprocessed,
    output:
        Out=missense_vars_with_alphamissense,
    params:
        min_coverage=30,
    script:
        filtering


###############################################################################
#                                  MAPS model                                 #
###############################################################################


rule MAPS_model_fit_on_syn:
    input:
        mutation_ht=mutation_ht,
        calibrate_on=syn_by_context_and_coverage,
    params:
        transformation="+",
        legend_title="",
        colors=colors3,
        text_size=22,
    output:
        plot=MAPS_model_fit_on_syn,
    script:
        viz_model


rule MAPS_model_fit_on_syn_genomes:
    input:
        mutation_ht=mutation_ht,
        calibrate_on=syn_by_context_and_coverage_genomes,
    params:
        transformation="+",
        legend_title="",
        colors=colors3,
        text_size=22,
    output:
        plot=MAPS_model_fit_on_syn_genomes,
    script:
        viz_model


###############################################################################
#                                  CAPS model                                 #
###############################################################################

if recalibrate_CAPS == True:

    rule caps_model:
        input:
            syn_vars=syn_by_context_and_coverage,
            intron_vars=intronic_by_context_and_coverage,
        output:
            phat=syn_by_context_p4,
            phat_sim=syn_by_context_ppd,
            theta_sample=theta_sample,
        params:
            group_by_context=True,
        resources:
            mem_mb=3072,
        script:
            caps_model

    rule caps_model_genomes:
        input:
            syn_vars=syn_by_context_and_coverage_genomes,
            intron_vars=intronic_by_context_and_coverage_genomes,
        output:
            phat=syn_by_context_p4_genomes,
            phat_sim=syn_by_context_ppd_genomes,
            theta_sample=theta_sample_genomes,
        params:
            group_by_context=True,
        resources:
            mem_mb=3072,
        script:
            caps_model

    rule caps_model_Whiffin2020:
        input:
            syn_vars=syn_by_context_Whiffin2020,
            intron_vars=intronic_by_context_Whiffin2020,
        output:
            phat=syn_by_context_p4_Whiffin2020,
            phat_sim=syn_by_context_ppd_Whiffin2020,
            theta_sample=theta_sample_Whiffin2020,
        params:
            group_by_context=False,
        resources:
            mem_mb=3072,
        script:
            caps_model


###############################################################################


rule syn_by_context_and_coverage:
    input:
        In=by_csq,
    output:
        Out=syn_by_context_and_coverage,
    params:
        variable="worst_csq",
        variable_value="synonymous_variant",
    script:
        filtering


rule syn_by_context_and_coverage_genomes:
    input:
        In=by_csq_genomes,
    output:
        Out=syn_by_context_and_coverage_genomes,
    params:
        variable="worst_csq",
        variable_value="synonymous_variant",
    script:
        filtering


rule intronic_by_context_and_coverage:
    input:
        In=by_csq,
    output:
        Out=intronic_by_context_and_coverage,
    params:
        variable="worst_csq",
        variable_value="intron_variant",
    script:
        filtering


rule intronic_by_context_and_coverage_genomes:
    input:
        In=by_csq_genomes,
    output:
        Out=intronic_by_context_and_coverage_genomes,
    params:
        variable="worst_csq",
        variable_value="intron_variant",
    script:
        filtering


###############################################################################


rule all_vars_scores2_by_mutability_band_plot:
    input:
        scores=[
            all_vars_caps_by_mutability_band_scores,
            all_vars_maps_by_mutability_band_scores,
        ],
    params:
        xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        new_xlab_labels=["Highest", "High", "Low", "Lowest"],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        aspect_ratio=0.7,
        text_size=20,
        xlab_size=20,
        point_size=1,
        colors=colors2,
        legend_title="",
        xlab="Mutability band",
        ylab="Selection score",
        xlab_angle=0,
        xlab_hjust=0.5,
        xlab_vjust=0,
    output:
        plot=all_vars_scores2_by_mutability_band_plot,
    script:
        viz_scores_multi


rule scores2_by_csq_plot:
    input:
        scores=[caps_by_csq_scores, maps_by_csq_scores],
    params:
        xlab_labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        aspect_ratio=0.5,
        text_size=16,
        xlab_size=16,
        point_size=1,
        xlab="",
        ylab="Selection score",
        xlab_angle=20,
        xlab_hjust=1,
        xlab_vjust=1,
        legend_title="",
    output:
        plot=scores2_by_csq_plot,
    script:
        viz_scores_multi


rule scores2_by_csq_genomes_plot:
    input:
        scores=[caps_by_csq_genomes_scores, maps_by_csq_genomes_scores],
    params:
        xlab_labels_set=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intergenic",
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        xlab="",
        ylab="Selection score",
        aspect_ratio=0.9,  #aspect_ratio=0.2,
        text_size=22,  #text_size=15,
        xlab_size=20,  #xlab_size=15,
        point_size=1,  #point_size=0.5
        xlab_angle=30,
        xlab_hjust=1,
        xlab_vjust=1,
        legend_title="",
    output:
        plot=scores2_by_csq_genomes_plot,
    script:
        viz_scores_multi


rule scores3_by_csq_plot:
    input:
        scores=[caps_ppd_by_csq_scores, caps_by_csq_scores, maps_by_csq_scores],
    params:
        xlab_labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        model_labels=["maps", "caps", "caps_ppd"],
        new_model_labels=["MAPS", "CAPS", "CAPS-PPD"],
        text_size=15,
        xlab_size=15,
        point_size=1,
        aspect_ratio=0.6,
        dodge_width=0.8,
        colors=colors3,
        xlab="",
        ylab="Selection score",
        xlab_angle=20,
        xlab_hjust=1,
        xlab_vjust=1,
        legend_title="",
    output:
        plot=scores3_by_csq_plot,
    script:
        viz_scores_multi


rule scores3_by_csq_genomes_plot:
    input:
        scores=[
            caps_ppd_by_csq_genomes_scores,
            caps_by_csq_genomes_scores,
            maps_by_csq_genomes_scores,
        ],
    params:
        xlab_labels_set=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intergenic",
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        model_labels=["maps", "caps", "caps_ppd"],
        new_model_labels=["MAPS", "CAPS", "CAPS-PPD"],
        dodge_width=0.95,
        colors=colors3,
        xlab="",
        ylab="Selection score",
        xlab_angle=45,
        xlab_hjust=1,
        xlab_vjust=1,
        legend_title="",
    output:
        plot=scores3_by_csq_genomes_plot,
    script:
        viz_scores_multi


rule scores2_near_acceptor_plot:
    input:
        scores=[caps_near_acceptor_scores, maps_near_acceptor_scores],
    params:
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        ylab="Selection score",
        NA_omit=True,
        point_size=0.6,
        point_alpha=0.7,
        dodge_width=1,
        ylim_min=-0.08,
        ylim_max=0.18,
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=0.5,
        legend_title="",
        xlab_size=14,
        aspect_ratio=0.3,
        xlab="Position relative to exon boundary",
        xlab_labels=[
            -24,
            -23,
            -22,
            -21,
            -20,
            -19,
            -18,
            -17,
            -16,
            -15,
            -14,
            -13,
            -12,
            -11,
            -10,
            -9,
            -8,
            -7,
            -6,
            -5,
            -4,
            -3,
            -2,
            -1,
            0,
            1,
            2,
            3,
            4,
            5,
            6,
            7,
            8,
            9,
            10,
            11,
        ],
        new_xlab_labels=[
            "-25",
            "-24",
            "-23",
            "-22",
            "-21",
            "-20",
            "-19",
            "-18",
            "-17",
            "-16",
            "-15",
            "-14",
            "-13",
            "-12",
            "-11",
            "-10",
            "-9",
            "-8",
            "-7",
            "-6",
            "-5",
            "-4",
            "-3",
            "-2",
            "-1",
            "0",
            "+1",
            "+2",
            "+3",
            "+4",
            "+5",
            "+6",
            "+7",
            "+8",
            "+9",
            "+10",
        ],
    output:
        plot=scores2_near_acceptor_plot,
    script:
        viz_scores_multi


rule scores2_near_acceptor_genomes_plot:
    input:
        scores=[caps_near_acceptor_genomes_scores, maps_near_acceptor_genomes_scores],
    params:
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        ylab="Selection score",
        NA_omit=True,
        point_size=0.6,
        point_alpha=0.7,
        dodge_width=0.2,
        ylim_min=-0.2,
        ylim_max=0.3,
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=0.5,
        legend_title="",
        xlab_size=14,
        aspect_ratio=0.3,
        xlab="Position relative to exon boundary",
        xlab_labels=[
            -24,
            -23,
            -22,
            -21,
            -20,
            -19,
            -18,
            -17,
            -16,
            -15,
            -14,
            -13,
            -12,
            -11,
            -10,
            -9,
            -8,
            -7,
            -6,
            -5,
            -4,
            -3,
            -2,
            -1,
            0,
            1,
            2,
            3,
            4,
            5,
            6,
            7,
            8,
            9,
            10,
            11,
        ],
        new_xlab_labels=[
            "-25",
            "-24",
            "-23",
            "-22",
            "-21",
            "-20",
            "-19",
            "-18",
            "-17",
            "-16",
            "-15",
            "-14",
            "-13",
            "-12",
            "-11",
            "-10",
            "-9",
            "-8",
            "-7",
            "-6",
            "-5",
            "-4",
            "-3",
            "-2",
            "-1",
            "0",
            "+1",
            "+2",
            "+3",
            "+4",
            "+5",
            "+6",
            "+7",
            "+8",
            "+9",
            "+10",
        ],
    output:
        plot=scores2_near_acceptor_genomes_plot,
    script:
        viz_scores_multi


rule scores2_near_donor_plot:
    input:
        scores=[caps_near_donor_scores, maps_near_donor_scores],
    params:
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        ylab="Selection score",
        NA_omit=True,
        point_size=0.6,
        point_alpha=0.7,
        dodge_width=0.75,
        ylim_min=-0.08,
        ylim_max=0.18,
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=0.5,
        legend_title="",
        xlab_labels=[
            -11,
            -10,
            -9,
            -8,
            -7,
            -6,
            -5,
            -4,
            -3,
            -2,
            -1,
            0,
            1,
            2,
            3,
            4,
            5,
            6,
            7,
            8,
            9,
        ],
        new_xlab_labels=[
            "-10",
            "-9",
            "-8",
            "-7",
            "-6",
            "-5",
            "-4",
            "-3",
            "-2",
            "-1",
            "0",
            "+1",
            "+2",
            "+3",
            "+4",
            "+5",
            "+6",
            "+7",
            "+8",
            "+9",
            "+10",
        ],
        xlab_size=14,
        aspect_ratio=0.3,
        xlab="Position relative to exon boundary",
    output:
        plot=scores2_near_donor_plot,
    script:
        viz_scores_multi


rule scores2_near_donor_genomes_plot:
    input:
        scores=[caps_near_donor_genomes_scores, maps_near_donor_genomes_scores],
    params:
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        ylab="Selection score",
        NA_omit=True,
        point_size=0.6,
        point_alpha=0.7,
        dodge_width=0.2,
        ylim_min=-0.2,
        ylim_max=0.3,
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=0.5,
        legend_title="",
        xlab_labels=[
            -11,
            -10,
            -9,
            -8,
            -7,
            -6,
            -5,
            -4,
            -3,
            -2,
            -1,
            0,
            1,
            2,
            3,
            4,
            5,
            6,
            7,
            8,
            9,
        ],
        new_xlab_labels=[
            "-10",
            "-9",
            "-8",
            "-7",
            "-6",
            "-5",
            "-4",
            "-3",
            "-2",
            "-1",
            "0",
            "+1",
            "+2",
            "+3",
            "+4",
            "+5",
            "+6",
            "+7",
            "+8",
            "+9",
            "+10",
        ],
        xlab_size=14,
        aspect_ratio=0.3,
        xlab="Position relative to exon boundary",
    output:
        plot=scores2_near_donor_genomes_plot,
    script:
        viz_scores_multi


rule scores3_by_csq_scores_comparison:
    input:
        scores=[caps_ppd_by_csq_scores, caps_by_csq_scores, maps_by_csq_scores],
    params:
        labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        model_labels=["maps", "caps", "caps_ppd"],
        new_model_labels=["MAPS", "CAPS", "CAPS-PPD"],
        caption="Increase in the size of confidence intervals in CAPS and CAPS-PPD compared to MAPS by variant class.",
        label="table:CIs_by_csq",
        variable="worst_csq",
    output:
        table=scores3_by_csq_scores_comparison,
    script:
        compare_scores


rule syn_vars_maps_by_mutability_band_plot:
    input:
        scores=syn_vars_maps_by_mutability_band_scores,
    output:
        plot=syn_vars_maps_by_mutability_band_plot,
    params:
        xlab_angle=45,
        ylim_min=-0.02,
        ylim_max=0.04,
        xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        new_xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        score_name="maps",
    script:
        viz_scores


rule syn_vars_maps_by_mutability_band:
    input:
        variants=syn_w_mutability_bands,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="mutability_band",
        transformation="+",
    output:
        scores=syn_vars_maps_by_mutability_band_scores,
    script:
        calculate_maps


rule syn_vars_caps_by_mutability_band:
    input:
        variants=syn_w_mutability_bands,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="mutability_band",
    output:
        scores=syn_vars_caps_by_mutability_band_scores,
    script:
        calculate_caps


rule maps_near_acceptor_scores:
    input:
        variants=near_acceptor_vars,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="splice_distance",
        transformation="+",
    output:
        scores=maps_near_acceptor_scores,
    script:
        calculate_maps


rule maps_near_acceptor_genomes_scores:
    input:
        variants=near_acceptor_genomes_vars,
        calibrate_on=syn_by_context_and_coverage_genomes,
    params:
        extra="splice_distance",
        transformation="+",
    output:
        scores=maps_near_acceptor_genomes_scores,
    script:
        calculate_maps


rule maps_near_donor_scores:
    input:
        variants=near_donor_vars,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="splice_distance",
        transformation="+",
    output:
        scores=maps_near_donor_scores,
    script:
        calculate_maps


rule maps_near_donor_genomes_scores:
    input:
        variants=near_donor_genomes_vars,
        calibrate_on=syn_by_context_and_coverage_genomes,
    params:
        extra="splice_distance",
        transformation="+",
    output:
        scores=maps_near_donor_genomes_scores,
    script:
        calculate_maps


rule maps_by_csq_scores:
    input:
        variants=by_csq,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="worst_csq",
        transformation="+",
    output:
        scores=maps_by_csq_scores,
    script:
        calculate_maps


rule maps_by_csq_genomes_scores:
    input:
        variants=by_csq_genomes,
        calibrate_on=syn_by_context_and_coverage_genomes,
    params:
        extra="worst_csq",
        transformation="+",
    output:
        scores=maps_by_csq_genomes_scores,
    script:
        calculate_maps


rule maps_by_csq_plot:
    input:
        scores=maps_by_csq_scores,
    output:
        plot=maps_by_csq_plot,
    params:
        ylim_max=0.16,
        ylim_min=-0.04,
        xlab_angle=45,
        xlab_labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        score_name="maps",
    script:
        viz_scores


rule annotate_all_vars_with_mutability_bands:
    input:
        variants=by_csq,
    output:
        annotated_variants=all_w_mutability_bands,
    script:
        annotate_with_mutability_bands


rule all_vars_caps_by_mutability_band_plot:
    input:
        scores=all_vars_caps_by_mutability_band_scores,
    output:
        plot=all_vars_caps_by_mutability_band_plot,
    params:
        xlab_angle=45,
        ylim_min=0,
        ylim_max=0.06,
        xlab="Mutability band",
        xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        new_xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        score_name="caps",
    script:
        viz_scores


rule all_vars_maps_by_mutability_band_plot:
    input:
        scores=all_vars_maps_by_mutability_band_scores,
    output:
        plot=all_vars_maps_by_mutability_band_plot,
    params:
        xlab_angle=45,
        ylim_min=0,
        ylim_max=0.06,
        xlab="Mutability band",
        xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        new_xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        score_name="maps",
    script:
        viz_scores


rule all_vars_sqrt_by_mutability_band_plot:
    input:
        scores=all_vars_sqrt_by_mutability_band_scores,
    output:
        plot=all_vars_sqrt_by_mutability_band_plot,
    params:
        xlab_angle=45,
        ylim_min=0,
        ylim_max=0.06,
        xlab="Mutability band",
        xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        new_xlab_labels=["Highest", "Higher", "Lower", "Lowest"],
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        score_name="maps",
    script:
        viz_scores


rule all_vars_caps_by_mutability_band:
    input:
        variants=all_w_mutability_bands,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="mutability_band",
    output:
        scores=all_vars_caps_by_mutability_band_scores,
    script:
        calculate_caps


rule all_vars_maps_by_mutability_band:
    input:
        variants=all_w_mutability_bands,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="mutability_band",
        transformation="+",
    output:
        scores=all_vars_maps_by_mutability_band_scores,
    script:
        calculate_maps


rule all_vars_sqrt_by_mutability_band:
    input:
        variants=all_w_mutability_bands,
        calibrate_on=syn_by_context_and_coverage,
    params:
        extra="mutability_band",
        transformation="sqrt",
    output:
        scores=all_vars_sqrt_by_mutability_band_scores,
    script:
        calculate_maps


###############################################################################


rule caps_near_acceptor_scores:
    input:
        variants=near_acceptor_vars,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="splice_distance",
    output:
        scores=caps_near_acceptor_scores,
    script:
        calculate_caps


rule caps_near_acceptor_genomes_scores:
    input:
        variants=near_acceptor_genomes_vars,
        exp_variants=syn_by_context_p4_genomes,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="splice_distance",
    output:
        scores=caps_near_acceptor_genomes_scores,
    script:
        calculate_caps


rule caps_near_donor_scores:
    input:
        variants=near_donor_vars,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="splice_distance",
    output:
        scores=caps_near_donor_scores,
    script:
        calculate_caps


rule caps_near_donor_genomes_scores:
    input:
        variants=near_donor_genomes_vars,
        exp_variants=syn_by_context_p4_genomes,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="splice_distance",
    output:
        scores=caps_near_donor_genomes_scores,
    script:
        calculate_caps


rule caps_by_csq_scores:
    input:
        variants=by_csq,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="worst_csq",
    output:
        scores=caps_by_csq_scores,
    script:
        calculate_caps


rule caps_by_csq_genomes_scores:
    input:
        variants=by_csq_genomes,
        exp_variants=syn_by_context_p4_genomes,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="worst_csq",
    output:
        scores=caps_by_csq_genomes_scores,
    script:
        calculate_caps


rule caps_ppd_by_csq_scores:
    input:
        variants=by_csq,
        exp_variants=syn_by_context_ppd,
    params:
        phat_method="PPD",
        variable="worst_csq",
    output:
        scores=caps_ppd_by_csq_scores,
    script:
        calculate_caps


rule caps_ppd_by_csq_genomes_scores:
    input:
        variants=by_csq_genomes,
        exp_variants=syn_by_context_ppd_genomes,
    params:
        phat_method="PPD",
        variable="worst_csq",
    output:
        scores=caps_ppd_by_csq_genomes_scores,
    script:
        calculate_caps


rule caps_by_csq_plot:
    input:
        scores=caps_by_csq_scores,
    output:
        plot=caps_by_csq_plot,
    params:
        ylim_max=0.16,
        aspect_ratio=0.6,
        ylim_min=-0.04,
        xlab_angle=30,
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        xlab_labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        score_name="caps",
    script:
        viz_scores


rule caps_ppd_by_csq_plot:
    input:
        scores=caps_ppd_by_csq_scores,
    output:
        plot=caps_ppd_by_csq_plot,
    params:
        ylim_max=0.16,
        ylim_min=-0.04,
        xlab_angle=45,
        xlab_vjust=1,
        xlab_hjust=1,
        reorder_xlab_by_score=False,
        xlab_labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        xlab_labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_xlab_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
        score_name="caps_ppd",
    script:
        viz_scores


######################################################################


rule syn_loeuf_bins_to_tops:
    input:
        groups=syn_constrained,
    output:
        joined_groups=syn_constrained_joined_groups,
    params:
        variable="loeuf",
        desc=True,
    script:
        join_groups


rule calculate_caps_by_loeuf_group_syn:
    input:
        variants=syn_constrained_joined_groups,
        exp_variants=syn_by_context_p4,
    params:
        phat_method="Var",
        variable="loeuf",
    output:
        scores=syn_constrained_tops_caps_scores,
    script:
        calculate_caps


rule calculate_maps_by_loeuf_group_syn:
    input:
        calibrate_on=syn_by_context_and_coverage,
        variants=syn_constrained_joined_groups,
    params:
        extra="loeuf",
        transformation="+",
    output:
        scores=syn_constrained_tops_maps_scores,
    script:
        calculate_maps


rule syn_loeuf_scores2_groups_plot:
    input:
        scores=[
            syn_constrained_tops_caps_scores,
            syn_constrained_tops_maps_scores,
        ],
    params:
        ylab="Selection score",
        xlab="LOEUF top group",
        reorder_xlab_by_score=False,
        xlab_labels=["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"],
        new_xlab_labels=[
            "10%",
            "20%",
            "30%",
            "40%",
            "50%",
            "60%",
            "70%",
            "80%",
            "90%",
            "100%",
        ],
        aspect_ratio=0.4,
        text_size=15,
        xlab_size=15,
        point_size=0.6,
        dodge_width=0.5,
        point_alpha=0.7,
        xlab_angle=20,
        xlab_vjust=1,
        xlab_hjust=1,
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
    output:
        plot=syn_constrained_scores2_plot,
    script:
        viz_scores_multi


rule syn_loeuf_groups_caps_plot:
    input:
        scores=syn_constrained_tops_caps_scores,
    output:
        plot=syn_constrained_caps_plot,
    params:
        xlab="LOEUF top group",
        reorder_xlab_by_score=False,
        aspect_ratio=0.5,
        xlab_labels=["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"],
        new_xlab_labels=[
            "10%",
            "20%",
            "30%",
            "40%",
            "50%",
            "60%",
            "70%",
            "80%",
            "90%",
            "100%",
        ],
        xlab_angle=45,
        xlab_vjust=1,
        xlab_hjust=1,
        score_name="caps",
    script:
        viz_scores


###############################################################################


rule context_by_csq_genomes:
    input:
        variants=by_csq_genomes,
    params:
        caption="Proportion of transversion contexts by variant class for genomes.",
        variable="worst_csq",
        variable_name="Class",
        labels_set=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        labels=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_labels=[
            "Intergenic",
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
    output:
        Out=contexts_by_csq_genomes,
    script:
        contexts


rule context_by_csq:
    input:
        variants=by_csq,
    params:
        caption="Proportion of transversion contexts by variant class for exomes.",
        variable="worst_csq",
        variable_name="Class",
        labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
    output:
        Out=contexts_by_csq,
    script:
        contexts


rule counts_by_csq:
    input:
        variants=by_csq,
    params:
        coverage_threshold=30,
        mode="tex",
        caption="Total number of variants observed in each functional class for exomes.",
        variable="worst_csq",
        labels_set=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        labels=[
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_labels=[
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
    output:
        Out=counts_by_csq,
    script:
        counts


rule counts_by_csq_genomes:
    input:
        variants=by_csq_genomes,
    params:
        coverage_threshold=30,
        mode="tex",
        caption="Total number of variants observed in each functional class for genomes.",
        variable="worst_csq",
        labels_set=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        labels=[
            "intergenic_variant",
            "intron_variant",
            "5_prime_UTR_variant",
            "3_prime_UTR_variant",
            "synonymous_variant",
            "missense_variant",
            "essential_splice",
            "stop_gained",
        ],
        new_labels=[
            "Intergenic",
            "Intron",
            "5'UTR",
            "3'UTR",
            "Synonymous",
            "Missense",
            "Essential splice",
            "Nonsense",
        ],
    output:
        Out=counts_by_csq_genomes,
    script:
        counts


###############################################################################
#                          Whiffin2020 uORFs analysis                         #
###############################################################################


rule Whiffin2020_get_files1:
    output:
        Whiffin2020_stop_unprocessed,
    shell:
        """
        wget https://github.com/ImperialCardioGenetics/uORFs/raw/master/data_files/stop-removing_all_possible_annotated.txt.gz
        gunzip -c stop-removing_all_possible_annotated.txt.gz > {output[0]}
        """


rule Whiffin2020_get_files2:
    output:
        Whiffin2020_uAUG_unprocessed,
    shell:
        """
        wget https://github.com/ImperialCardioGenetics/uORFs/raw/master/data_files/uAUG-creating_all_possible_annotated.txt.gz
        gunzip -c uAUG-creating_all_possible_annotated.txt.gz > {output[0]}
        """


rule Whiffin2020_uORFs_preprocessing_stop:
    input:
        inFile=Whiffin2020_stop_unprocessed,
        mu=mutation_ht,
    output:
        outFile=Whiffin2020_stop,
    script:
        Whiffin2020_uORFs_preprocessing


rule Whiffin2020_uORFs_preprocessing_uAUG:
    input:
        inFile=Whiffin2020_uAUG_unprocessed,
        mu=mutation_ht,
    output:
        outFile=Whiffin2020_uAUG,
    script:
        Whiffin2020_uORFs_preprocessing


# Effects #############################################################


rule Whiffin2020_uORFs_effects_maps_stop:
    input:
        calibrate_on=syn_by_context_Whiffin2020,
        variants=Whiffin2020_stop,
    params:
        extra="effect",
        transformation="+",
    output:
        scores=Whiffin2020_uORFs_effects_stop_maps_scores,
    script:
        calculate_maps


rule Whiffin2020_uORFs_effects_caps_stop:
    input:
        variants=Whiffin2020_stop,
        exp_variants=syn_by_context_p4_Whiffin2020,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="effect",
    output:
        scores=Whiffin2020_uORFs_effects_stop_caps_scores,
    script:
        calculate_caps


rule Whiffin2020_uORFs_effects_plot_stop:
    input:
        scores=[
            Whiffin2020_uORFs_effects_stop_maps_scores,
            Whiffin2020_uORFs_effects_stop_caps_scores,
        ],
    params:
        xlab_labels_set=["uORF_elongated", "out-of-frame_oORF"],
        xlab_labels=["uORF_elongated", "out-of-frame_oORF"],
        new_xlab_labels=["uORF elongated", "oORF formed"],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
        xlab="",
        ylab="Selection score",
        aspect_ratio=0.3,
        text_size=23,
        xlab_size=23,
        xlab_angle=10,
        point_size=1,
        xlab_hjust=1,
        xlab_vjust=1,
        ylim_min=-0.07,
        ylim_max=0.15,
        # Value taken from "gs://gcp-public-data--gnomad/papers/2019-flagship-lof/v1.0/summary_results/maps_plain_genomes.txt.bgz"
        gnomAD_missense_level_genomes=0.06182,
    output:
        plot=Whiffin2020_uORFs_effects_stop_plot,
    script:
        viz_scores_multi


rule Whiffin2020_uORFs_effects_maps_uAUG:
    input:
        calibrate_on=syn_by_context_Whiffin2020,
        variants=Whiffin2020_uAUG,
    params:
        extra="effect",
        transformation="+",
    output:
        scores=Whiffin2020_uORFs_effects_uAUG_maps_scores,
    script:
        calculate_maps


rule Whiffin2020_uORFs_effects_caps_uAUG:
    input:
        variants=Whiffin2020_uAUG,
        exp_variants=syn_by_context_p4_Whiffin2020,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="effect",
    output:
        scores=Whiffin2020_uORFs_effects_uAUG_caps_scores,
    script:
        calculate_caps


rule Whiffin2020_uORFs_effects_plot_uAUG:
    input:
        scores=[
            Whiffin2020_uORFs_effects_uAUG_maps_scores,
            Whiffin2020_uORFs_effects_uAUG_caps_scores,
        ],
    params:
        xlab_labels_set=["uORF_created", "out-of-frame_oORF", "CDS_elongated"],
        xlab_labels=["uORF_created", "out-of-frame_oORF", "CDS_elongated"],
        new_xlab_labels=["uORF created", "Out-of-frame oORF", "CDS elongated"],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
        xlab="",
        ylab="Selection score",
        text_size=23,
        xlab_size=23,
        xlab_angle=10,
        aspect_ratio=0.3,
        point_size=1,
        xlab_hjust=1,
        xlab_vjust=1,
        ylim_min=-0.07,
        ylim_max=0.15,
        # Value taken from "gs://gcp-public-data--gnomad/papers/2019-flagship-lof/v1.0/summary_results/maps_plain_genomes.txt.bgz"
        gnomAD_missense_level_genomes=0.06182,
    output:
        plot=Whiffin2020_uORFs_effects_uAUG_plot,
    script:
        viz_scores_multi


# Kozak strength ######################################################


rule Whiffin2020_uORFs_kozak_maps_stop:
    input:
        calibrate_on=syn_by_context_Whiffin2020,
        variants=Whiffin2020_stop,
    params:
        extra="Kozak_strength",
        transformation="+",
    output:
        scores=Whiffin2020_uORFs_kozak_stop_maps_scores,
    script:
        calculate_maps


rule Whiffin2020_uORFs_kozak_caps_stop:
    input:
        variants=Whiffin2020_stop,
        exp_variants=syn_by_context_p4_Whiffin2020,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="Kozak_strength",
    output:
        scores=Whiffin2020_uORFs_kozak_stop_caps_scores,
    script:
        calculate_caps


rule Whiffin2020_uORFs_kozak_plot_stop:
    input:
        scores=[
            Whiffin2020_uORFs_kozak_stop_maps_scores,
            Whiffin2020_uORFs_kozak_stop_caps_scores,
        ],
    params:
        xlab_labels_set=["weak", "moderate", "strong"],
        xlab_labels=["weak", "moderate", "strong"],
        new_xlab_labels=["Weak", "Moderate", "Strong"],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
        xlab="",
        ylab="Selection score",
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=1,
        ylim_min=-0.1,
        ylim_max=0.2,
        # Value taken from "gs://gcp-public-data--gnomad/papers/2019-flagship-lof/v1.0/summary_results/maps_plain_genomes.txt.bgz"
        gnomAD_missense_level_genomes=0.06182,
    output:
        plot=Whiffin2020_uORFs_kozak_stop_plot,
    script:
        viz_scores_multi


rule Whiffin2020_uORFs_kozak_maps_uAUG:
    input:
        calibrate_on=syn_by_context_Whiffin2020,
        variants=Whiffin2020_uAUG,
    params:
        extra="Kozak_strength",
        transformation="+",
    output:
        scores=Whiffin2020_uORFs_kozak_uAUG_maps_scores,
    script:
        calculate_maps


rule Whiffin2020_uORFs_kozak_caps_uAUG:
    input:
        variants=Whiffin2020_uAUG,
        exp_variants=syn_by_context_p4_Whiffin2020,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="Kozak_strength",
    output:
        scores=Whiffin2020_uORFs_kozak_uAUG_caps_scores,
    script:
        calculate_caps


rule Whiffin2020_uORFs_kozak_plot_uAUG:
    input:
        scores=[
            Whiffin2020_uORFs_kozak_uAUG_maps_scores,
            Whiffin2020_uORFs_kozak_uAUG_caps_scores,
        ],
    params:
        xlab_labels_set=["weak", "moderate", "strong"],
        xlab_labels=["weak", "moderate", "strong"],
        new_xlab_labels=["Weak", "Moderate", "Strong"],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
        xlab="",
        ylab="Selection score",
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=1,
        ylim_min=-0.1,
        ylim_max=0.2,
        # Value taken from "gs://gcp-public-data--gnomad/papers/2019-flagship-lof/v1.0/summary_results/maps_plain_genomes.txt.bgz"
        gnomAD_missense_level_genomes=0.06182,
    output:
        plot=Whiffin2020_uORFs_kozak_uAUG_plot,
    script:
        viz_scores_multi


# Distance to CDS #####################################################


rule Whiffin2020_uORFs_distanceCDS_maps_uAUG:
    input:
        calibrate_on=syn_by_context_Whiffin2020,
        variants=Whiffin2020_uAUG,
    params:
        extra="distanceToCDS",
        transformation="+",
    output:
        scores=Whiffin2020_uORFs_distanceCDS_uAUG_maps_scores,
    script:
        calculate_maps


rule Whiffin2020_uORFs_distanceCDS_caps_uAUG:
    input:
        variants=Whiffin2020_uAUG,
        exp_variants=syn_by_context_p4_Whiffin2020,
    params:
        phat_method="Var",
        confint_method="CAPS",
        variable="distanceToCDS",
    output:
        scores=Whiffin2020_uORFs_distanceCDS_uAUG_caps_scores,
    script:
        calculate_caps


rule Whiffin2020_uORFs_distanceCDS_plot_uAUG:
    input:
        scores=[
            Whiffin2020_uORFs_distanceCDS_uAUG_maps_scores,
            Whiffin2020_uORFs_distanceCDS_uAUG_caps_scores,
        ],
    params:
        # xlab_labels_set=["", "", ""],
        # xlab_labels=["", "", ""],
        # new_xlab_labels=["", "", ""],
        model_labels=["maps", "caps"],
        new_model_labels=["MAPS", "CAPS"],
        colors=colors2,
        legend_title="",
        xlab="",
        ylab="Selection score",
        xlab_angle=90,
        xlab_hjust=1,
        xlab_vjust=1,
        ylim_min=-0.1,
        ylim_max=0.2,
        # Value taken from "gs://gcp-public-data--gnomad/papers/2019-flagship-lof/v1.0/summary_results/maps_plain_genomes.txt.bgz"
        gnomAD_missense_level_genomes=0.06182,
    output:
        plot=Whiffin2020_uORFs_distanceCDS_uAUG_plot,
    script:
        viz_scores_multi


###############################################################################
#                        CAPS vs MAPS on AlphaMissense                        #
###############################################################################

for mu1, mu2 in mutability_pairs:

    rule filter_missense_vars_with_alphamissense:
        name:
            "filter_missense_vars_with_alphamissense_{mu1}_{mu2}"
        input:
            mutation_ht=mutation_ht,
            variants=missense_vars_with_alphamissense,
            syn_by_context_and_coverage=syn_by_context_and_coverage,
        output:
            Out="CAPSvsMAPS_for_contexts_{mu1}_{mu2}.tsv",
        params:
            mu1="{mu1}",
            mu2="{mu2}",
            variant_level=True,
            only_deleterious=True,
        script:
            filter_alphamissense

    rule missense_vars_with_alphamissense_contexts_caps:
        name:
            "missense_vars_with_alphamissense_contexts_caps_{mu1}_{mu2}"
        input:
            variants="CAPSvsMAPS_for_contexts_{mu1}_{mu2}.tsv",
            exp_variants=syn_by_context_p4,
        params:
            phat_method="Var",
            confint_method="CAPS",
            variable="full_context",
        output:
            scores="missense_vars_with_alphamissense_caps_contexts_{mu1}_{mu2}.tsv",
        script:
            calculate_caps

    rule missense_vars_with_alphamissense_contexts_maps:
        name:
            "missense_vars_with_alphamissense_contexts_maps_{mu1}_{mu2}"
        input:
            variants="CAPSvsMAPS_for_contexts_{mu1}_{mu2}.tsv",
            calibrate_on=syn_by_context_and_coverage,
        params:
            extra="full_context",
            transformation="+",
        output:
            scores="missense_vars_with_alphamissense_maps_contexts_{mu1}_{mu2}.tsv",
        script:
            calculate_maps

    rule missense_vars_with_alphamissense_contexts_avg_alphamissense:
        name:
            "missense_vars_with_alphamissense_contexts_avg_alphamissense_{mu1}_{mu2}"
        input:
            variants="CAPSvsMAPS_for_contexts_{mu1}_{mu2}.tsv",
        output:
            scores="missense_vars_with_alphamissense_contexts_{mu1}_{mu2}.tsv",
        script:
            calculate_avg_alphamissense

    rule missense_vars_with_alphamissense_contexts_scores2_plot:
        name:
            "missense_vars_with_alphamissense_contexts_scores2_plot_{mu1}_{mu2}"
        input:
            scores=[
                "missense_vars_with_alphamissense_caps_contexts_{mu1}_{mu2}.tsv",
                "missense_vars_with_alphamissense_maps_contexts_{mu1}_{mu2}.tsv",
                "missense_vars_with_alphamissense_contexts_{mu1}_{mu2}.tsv",
            ],
        params:
            model_labels=["maps", "caps", "alphamissense"],
            new_model_labels=["MAPS", "CAPS", "AlphaMissense"],
            with_AlphaMissense = True,
            colors=colors2 + ["black"],
            aspect_ratio=0.5,
            text_size=20,
            xlab_size=20,
            point_size=1,
            legend_title="",
            xlab="Context",
            ylab="Score",
            xlab_angle=20,
            xlab_hjust=0.5,
            xlab_vjust=0.5,
        output:
            plot="CAPSvsMAPS_for_contexts_{mu1}_{mu2}.pdf",
        script:
            viz_scores_multi