import sys
import os
from itertools import chain
from collections import defaultdict
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
pd.set_option('precision', 3)
pd.set_option('display.max_colwidth', -1)
import gseapy as gp
import goatools

obodag = goatools.obo_parser.GODag('go-basic.obo')

sys.path.insert(0, "/home/gstupp/projects/metaproteomics")
from metaproteomics import utils
#from metaproteomics.analysis import build_loci

BASE = '../out/'
grouped_loci = utils.load(os.path.join(BASE,"grouped_loci_filt_norm.pkl.gz"))

def make_go2Gene_map(grouped_loci, ontology='MF'):    

    ontology_map = {'MF': 'molecular_function', 'BP': 'biological_process', 'CC': 'cellular_component'}
    
    out = defaultdict(set)    
    for l in grouped_loci:
        if 'go' in l.annotations:
            for go in l.annotations['go']:
                if obodag[go].namespace == ontology_map[ontology]:
                    out[go].add(l.cluster_id)
                    for parent in obodag[go].get_all_parents():
                        if obodag[parent].namespace == ontology_map[ontology]:
                            out[parent].add(l.cluster_id)
                
    return dict(out)

def filter_go2gene_map(go_locus):
    
    # Remove "very broad" gene sets. Arbitrary definition: gene sets that emcompass >50% of all IDs
    all_ids = set(chain(*go_locus.values()))
    go_locus = {key: value for (key, value) in go_locus.items() if len(value) / len(all_ids) <= 0.5}

    # Remove terms with less than 5 members: changed from 10 to 5 becasue small #s of proteins compared
    # to what you would find wiht genes
    go_locus = {key: value for (key, value) in go_locus.items() if len(value) >= 5}

    # Remove child terms with identical gene sets as their parents
    to_remove = set()
    for parent in go_locus.keys():
        # If child term has exact same members as parent, remove
        child_ids = [x.id for x in obodag[parent].children if x.id in go_locus.keys()]
        for child in child_ids:
            if go_locus[child] == go_locus[parent]:
                to_remove.add(child)
    go_locus = {key: value for (key, value) in go_locus.items() if key not in to_remove}

    # Remove sibling terms with identical gene sets
    to_remove = set()
    to_keep = set()
    for brother in go_locus.keys():
        to_keep.add(brother) # make sure filtered out siblings don't filter out ones we want to keep
        for parent in obodag[brother].parents:
            siblings = set([y.id for y in parent.children])
            siblings.remove(brother)
            for sibling in siblings:
                if sibling in go_locus.keys() and go_locus[brother] == go_locus[sibling] and not sibling in to_keep:
                    to_remove.add(sibling)
    go_locus = {key: value for (key, value) in go_locus.items() if key not in to_remove}

    return go_locus

def gomap_to_csv(go2gene, out_file = 'test.tsv'):

    out = ""
    for term, loci in go2gene.items():
        out += "{}\t".format(term)
        out += "{}\t".format(obodag[term].name)
        out += '\t'.join(list(map(str,loci)))
        out += '\n'
        
    with open(out_file, 'w') as fout:
        fout.write(out)

def run_go_gsea(rank_df, g2g_map, seed, outdir='tmp'):
    """
    A ranked df and go2gene mapping returns the result dataframe for GSEA against all go-Terms
    
    loci must be grouped such that avg_ratio and p-values are correct for 1 phenotype
    see rt_unenr_grouped_loci above for example
    """
    import gseapy as gp
        
    # save the go 2 gene map, since gseapy doesn't seem to be able to use one already in memory
    gomap_to_csv(g2g_map, 'temp.gmt')
    
    res = gp.prerank(rnk=rank_df, gene_sets='temp.gmt', outdir=outdir, min_size = 5, max_size=500, 
                     permutation_n = 10000, graph_num = 20, seed=seed)
      
    def get_go_name(term):
        return obodag[term].name
    
    res['name'] = res.index.map(get_go_name)
    
    return res.sort_values('nes', ascending=False)

def plot_gsea_result(row, rank):
    return gp.plot.gsea_plot(rank, row['name'], row.hit_index, row.nes, row.pval, row.fdr, row.rank_ES, phenoPos='Tcell', phenoNeg='RAG')

mf_map = make_go2Gene_map(grouped_loci)
mf_map_f = filter_go2gene_map(mf_map)

print('Unfiltered: {}\tFiltered: {}'.format(len(mf_map), len(mf_map_f)))

bp_map = make_go2Gene_map(grouped_loci, 'BP')
bp_map_f = filter_go2gene_map(bp_map)

print('Unfiltered: {}\tFiltered: {}'.format(len(bp_map), len(bp_map_f)))

cc_map = make_go2Gene_map(grouped_loci, 'CC')
cc_map_f = filter_go2gene_map(cc_map)

print('Unfiltered: {}\tFiltered: {}'.format(len(cc_map), len(cc_map_f)))

"""
Explanation: gseapy analysis of CM7_CM1E2d56col_unenr123_rawextract_2017
gsea algorithm: http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html
gseapy: https://media.readthedocs.org/pdf/gseapy/latest/gseapy.pdf
https://github.com/stuppie/CM7_CM1E2d56col_unenr123_rawextract_2017/blob/master/scripts/gseapy.ipynb
End of explanation
"""


out_dir = "RT_control_hm_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_control_results_named_annot.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)
mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
#cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_control = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_control = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
#cc_rt = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
#cc_control = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_control

bp_rt

bp_control

"""
Explanation: RT vs Control Human/Mouse
End of explanation
"""


out_dir = "RT_control_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_control_results_named_annot.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[~df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)

mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_control = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_control = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
cc_rt = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
cc_control = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_control

bp_rt

bp_control

cc_rt

cc_control

"""
Explanation: RT vs Control Non-Human/Mouse
End of explanation
"""


out_dir = "Rag_WT_hm_gsea"
df = pd.read_csv(os.path.join(BASE,"Rag_WT_results_named_annot.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)
mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
#cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rag = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_wt = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rag = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_wt = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
#cc_rt = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
#cc_control = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rag

mf_wt

bp_rag

bp_wt

"""
Explanation: Rag vs WT Human/Mouse
End of explanation
"""


out_dir = "Rag_WT_gsea"
df = pd.read_csv(os.path.join(BASE,"Rag_WT_results_named_annot.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[~df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)

mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rag = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_wt = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rag = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_wt = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
cc_rag = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
cc_wt = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rag

mf_wt

bp_rag

bp_wt

cc_rag

cc_wt

"""
Explanation: Rag vs WT Non-Human/Mouse
End of explanation
"""


out_dir = "RT_WT_hm_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_WT_deseq_results.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)
mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_wt = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_wt = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_wt

bp_rt

bp_wt

"""
Explanation: RT vs WT Human/Mouse
End of explanation
"""


out_dir = "RT_WT_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_WT_deseq_results.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[~df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)

mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_wt = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_wt = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
cc_rt = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
cc_wt = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_wt

bp_rt

bp_wt

cc_rt

cc_wt

"""
Explanation: RT vs WT Non-Human/Mouse
End of explanation
"""


out_dir = "RT_Rag_hm_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_Rag_deseq_results.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)
mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_rag = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_rag = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_rag

bp_rt

bp_rag

"""
Explanation: RT vs Rag Human/Mouse
End of explanation
"""


out_dir = "RT_Rag_gsea"
df = pd.read_csv(os.path.join(BASE,"RT_Rag_deseq_results.csv"))
df = df[(df.padj.abs()<=0.2)]
df = df[~df.human_mouse]
df['log2FoldChange'] = -1 * df['log2FoldChange']
rank_df = df[['Unnamed: 0', 'log2FoldChange']].rename(columns={'Unnamed: 0': 'gene_name', 'log2FoldChange': 'rank'})
rank_df = rank_df.sort_values('rank').reset_index(drop=True)

mf_res = run_go_gsea(rank_df, mf_map_f, seed=1111, outdir=out_dir)
bp_res = run_go_gsea(rank_df, bp_map_f, seed=1111, outdir=out_dir)
cc_res = run_go_gsea(rank_df, cc_map_f, seed=1111, outdir=out_dir)
mf_rt = mf_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
mf_rag = mf_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
bp_rt = bp_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
bp_rag = bp_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)
cc_rt = cc_res.query('nes > 0 and pval < 0.05').sort_values('nes', ascending=False)
cc_rag = cc_res.query('nes < 0 and pval < 0.05').sort_values('nes', ascending=True)

mf_rt

mf_rag

bp_rt

bp_rag

cc_rt

cc_rag

"""
Explanation: RT vs Rag Non-Human/Mouse
End of explanation
"""


