Snakefile

import itertools
import matplotlib
import os
matplotlib.use('Agg')
import pandas as pd
import seaborn as sns

import megago.megago as mg

SAMPLES, = glob_wildcards("../megago_comparison/go-per-sample/go_terms_{sample}.csv")


rule all:
    input:
        molecular_function="figures/similarity_mf.svg",
        biological_process="figures/similarity_bp.svg",
        cellular_component="figures/similarity_cc.svg"

rule visualization:
    input:
        molecular_function="aggregated_sim/similarity_mf.csv",
        biological_process="aggregated_sim/similarity_bp.csv",
        cellular_component="aggregated_sim/similarity_cc.csv"
    output:
        molecular_function="figures/similarity_mf.svg",
        biological_process="figures/similarity_bp.svg",
        cellular_component="figures/similarity_cc.svg"
    run:
        df_dict = {namespace: pd.read_csv(f, index_col=0) for namespace, f in input.items()}
        index_sorted = sorted(df_dict["molecular_function"].index, key=lambda s: s[-2]+s[:3])
        width = 6
        for i, (namespace, df) in enumerate(df_dict.items()):
            cluster_grid = sns.clustermap(df, xticklabels=True, yticklabels=True, method="average")
            cluster_grid.savefig(output[namespace])

rule aggregate_similaries:
    input:
        expand("similarity/{sample_comb}.csv",
               sample_comb=[f"{c[0]}-vs-{c[1]}" for c in itertools.product(SAMPLES, repeat=2)])
    output:
        molecular_function="aggregated_sim/similarity_mf.csv",
        biological_process="aggregated_sim/similarity_bp.csv",
        cellular_component="aggregated_sim/similarity_cc.csv"
    run:
        namespaces = ["molecular_function", "biological_process", "cellular_component"]
        df_dict = {namespace: pd.DataFrame(index=sorted(SAMPLES), columns=sorted(SAMPLES)) for namespace in namespaces}
        for f_in in input:
            filename_wo_ext = f_in.split("/")[-1].split(".")[0]
            s1, s2 = filename_wo_ext.split("-vs-")
            df_s = pd.read_csv(f_in, index_col=0)
            for namespace, df_agg in df_dict.items():
                sim = df_s.loc[namespace, "SIMILARITY"]
                df_agg.loc[s1, s2] = sim
        for namespace, df_agg in df_dict.items():
            f_out = output[namespace]
            df_agg.to_csv(f_out)

rule calc_pairwise_sample_similarity:
    input:
        "../megago_comparison/go-per-sample/go_terms_{sample_a}.csv",
        "../megago_comparison/go-per-sample/go_terms_{sample_b}.csv"
    output:
        "similarity/{sample_a}-vs-{sample_b}.csv",
        "similarity/{sample_b}-vs-{sample_a}.csv"
    run:
        with open(input[0], "r") as f_a:
            go_a_raw = set(mg.read_input(f_a))
        with open(input[1], "r") as f_b:
            go_b_raw = set(mg.read_input(f_b))
        go_dag = mg.GODag(mg.GO_DAG_FILE_PATH, prt=open(os.devnull, 'w'))
        go_a_bp, go_a_cc, go_a_mf = mg.split_per_domain(go_a_raw, go_dag)
        go_b_bp, go_b_cc, go_b_mf = mg.split_per_domain(go_b_raw, go_dag)

        with open(output[0], "w") as ab:
            with open(output[1], "w") as ba:
                ab.write("DOMAIN,SIMILARITY\n")
                ba.write("DOMAIN,SIMILARITY\n")
                for cat, go_a, go_b in zip(["biological_process", "cellular_component", "molecular_function"],
                                           [go_a_bp, go_a_cc, go_a_mf],
                                           [go_b_bp, go_b_cc, go_b_mf]
                                            ):
                    sim_ab = len(set(go_a).intersection(set(go_b))) / len(set(go_a))
                    sim_ba = len(set(go_a).intersection(set(go_b))) / len(set(go_b))
                    ab.write(f"{cat},{sim_ab}\n")
                    ba.write(f"{cat},{sim_ba}\n")