import logging
import os
import time

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pandas.tools.plotting import scatter_matrix
from matplotlib_venn import venn2
import seaborn as sns

import pybel
import pybel_tools as pbt
from pybel_tools.analysis.concordance import *
from pybel.constants import *
from pybel.canonicalize import calculate_canonical_name
from pybel_tools.visualization import to_jupyter

%config InlineBackend.figure_format = 'svg'
%matplotlib inline

"""
Explanation: Abstract
Author: Charles Tapley Hoyt
Estimated Run Time: 4 minutes
This notebook calculates the concordance between differential gene expression and the causal relationships between their protein products. This analysis assumes the forward hypothesis, under which gene expression is believed to be correlated with protein activity. While there are many faults to this assumption, it directly enables a very simple analysis of a knowledge assembly. 
Notebook Imports
End of explanation
"""


time.asctime()

# seed the random number generator
import random
random.seed(127)

pybel.__version__

pbt.__version__

"""
Explanation: Notebook Provenance
The time of execution, random number generator seed, and the versions of the software packages used are displayed explicitly.
End of explanation
"""


bms_base = os.environ['BMS_BASE']

"""
Explanation: Local Path Definitions
To make this notebook interoperable across many machines, locations to the repositories that contain the data used in this notebook are referenced from the environment, set in ~/.bashrc to point to the place where the repositories have been cloned. Assuming the repositories have been git clone'd into the ~/dev folder, the entries in ~/.bashrc should look like:
bash
...
export BMS_BASE=~/dev/bms
...
BMS
The biological model store (BMS) is the internal Fraunhofer SCAI repository for keeping BEL models under version control. It can be downloaded from https://tor-2.scai.fraunhofer.de/gf/project/bms/
End of explanation
"""


owncloud_base = os.environ['OWNCLOUD_BASE']

"""
Explanation: OwnCloud
The differential gene expression data used in this notebook is currently not published, and is obfuscated with reference through our team's internal data storage system with OwnCloud.
End of explanation
"""


pickle_path = os.path.join(bms_base, 'aetionomy', 'alzheimers', 'alzheimers.gpickle')

graph = pybel.from_pickle(pickle_path)

graph.version

"""
Explanation: Data
Alzheimer's Disease Knowledge Assembly
The Alzheimer's Disease Knowledge Assembly has been precompiled with the following command line script, and will be loaded from this format for improved performance. In general, derived data, such as the gpickle representation of a BEL script, are not saved under version control to ensure that the most up-to-date data is always used.
sh
pybel convert --path "$BMS_BASE/aetionomy/alzheimers.bel" --pickle "$BMS_BASE/aetionomy/alzheimers.gpickle"
The BEL script can also be compiled from inside this notebook with the following python code:
```python



import os
import pybel
Input from BEL script
bel_path = os.path.join(bms_base, 'aetionomy', 'alzheimers.bel')
graph = pybel.from_path(bel_path)
Output to gpickle for fast loading later
pickle_path = os.path.join(bms_base, 'aetionomy', 'alzheimers.gpickle')
pybel.to_pickle(graph, pickle_path)
```
End of explanation
"""


pbt.filters.remove_nodes_by_namespace(graph, {'MGI', 'RGD'})

"""
Explanation: All orthologies are discarded before analysis.
End of explanation
"""


pbt.mutation.collapse_by_central_dogma_to_genes(graph)

"""
Explanation: Following the forward hypothesis, the knowledge graph is collapsed to genes using  pbt.mutation.collapse_by_central_dogma_to_genes.
End of explanation
"""


pbt.mutation.rewire_variants_to_genes(graph)

pbt.summary.print_summary(graph)

pbt.summary.plot_summary(graph, plt, figsize=(10, 4))
plt.show()

"""
Explanation: An additional assumption is made about the activities of variants of proteins and genes. All variants are collapsed to the reference gene.
End of explanation
"""


data_path = os.path.join(owncloud_base, 'alzheimers', 'SevAD.csv')
target_columns = ['Gene.symbol', 'logFC']

df = pd.read_csv(data_path)
df = df.loc[df['Gene.symbol'].notnull(), target_columns]
df.head()

"""
Explanation: Differential Gene Expression
Differential gene expression data can be obtained from many sources, including ADNI and other large clinical studies. This analysis is concerned with the log-fold-changes on each gene, and not necessarily the p-value. This is better as a data-driven process becuase it does not require a model or multiple hypothesis testing on raw data.
End of explanation
"""


data = {k: v for _, k, v in df.itertuples()}

hgnc_names = pbt.summary.get_names_by_namespace(graph, 'HGNC')
df_names = set(df['Gene.symbol'])
overlapping_hgnc_names = hgnc_names & df_names

fix, (rax, lax) = plt.subplots(1, 2, figsize=(10, 3))

lax.set_title('Distribution of Log-Fold-Changes in E-GEOD-5281')
lax.set_xlabel('Log-Fold-Change')
lax.set_ylabel('Frequency')
lax.hist(list(data.values()))

rax.set_title('Gene Overlap')
venn2([hgnc_names, df_names], set_labels=["Alzheimer's Disease KA", 'E-GEOD-5281'], ax=rax)
plt.tight_layout()
plt.show()

"""
Explanation: A histogram of the log-fold-changes shows that the data are normally distributed, as expecOn the left, the number of shared elements in the knowledge assembly and differential gene data set are counted.
On the right, a histogram of the log-fold-changes shows that the data are normally distributed, as expected for differential gene expression data.ted for differential gene expression data.
End of explanation
"""


key = 'weight'
cutoff = 0.3

pbt.integration.overlay_type_data(graph, data, key, GENE, 'HGNC', overwrite=False, impute=0)

graph.edge[graph.nodes()[25]].keys()

"""
Explanation: Analysis
Data Integration
Finally, the differential gene expression data are ovelayed on the BEL graph with pbt.integration.overlay_type_data
End of explanation
"""


cutoffs = np.linspace(0, 3, 50)

plt.plot(cutoffs, [
    calculate_concordance(graph, key, cutoff=cutoff)
    for cutoff in cutoffs
])

plt.title('Effect of Cutoff on Concordance \nfor {}'.format(graph))
plt.ylabel('Concordance')
plt.xlabel('Cutoff')
plt.show()

"""
Explanation: Concordance
The effect of the cutoff is explored when calculating the concordance of the full knowledge assembly.
End of explanation
"""


concordance_df = pd.DataFrame.from_dict(
    {
        value.replace(' subgraph', ''): calculate_concordance_helper(subgraph, key, cutoff)
        for value, subgraph in pbt.selection.get_subgraphs_by_annotation(graph, 'Subgraph').items()
    }, 
    orient='index'
)

concordance_df.columns = 'Correct', 'Incorrect', 'Ambiguous', 'Unassigned'

"""
Explanation: Stratify by Subgraph
The distribution of the concordance values across the stratified subgraphs is displayed below.
End of explanation
"""


has_correct = concordance_df['Correct'] > 0
has_incorrect = concordance_df['Incorrect'] > 0

concordance_df[has_correct | has_incorrect].sort_values('Correct', ascending=False).head(20)

"""
Explanation: subgraphs that contain at least one correct or one incorrect relation are shown.
End of explanation
"""


concordance_df['Concordance'] = concordance_df['Correct'] / (concordance_df['Correct'] + concordance_df['Incorrect'])

concordance_df.sort_values('Concordance', ascending=False).head(20)

"""
Explanation: The 20 highest concording subgraphs are shown. Even without statistical analysis, these data suggest that high concorance values are most likely due to random chance.
End of explanation
"""


%%time
results = calculate_concordance_by_annotation(graph, 'Subgraph', key, cutoff)

"""
Explanation: While these calculations were presented in a dataframe, they can be directly acquired with calculate_concordance_by_annotation.
End of explanation
"""


sns.distplot([
    result
    for result in results.values() 
    if result != -1
])

plt.title('Distribution of concordance values at $C={}$'.format(cutoff))
plt.xlabel('Concordance')
plt.ylabel('Frequency')
plt.xlim([0, 1])
plt.show()

"""
Explanation: The distribution of concordance values over all subgraphs is shown.
End of explanation
"""


