import NotebookImport
from Imports import *

import seaborn as sns
sns.set_context('paper',font_scale=1.5)
sns.set_style('white')

"""
Explanation: Differential Analysis with both GEO and RNA-Seq
Import everything from the imports notebook. This reads in all of the expression data as well as the functions needed to analyse differential expression data.
End of explanation
"""


matched_rna = pd.read_hdf('/data_ssd/RNASeq_2014_07_15.h5', 'matched_tn')
rna_microarray = pd.read_hdf('/data_ssd/GEO_microarray_dx.h5', 'data')
matched_rna = rna_microarray.join(matched_rna)

"""
Explanation: Read in matched Gene expression data.
End of explanation
"""


dx_rna = binomial_test_screen(matched_rna, fc=1.)
dx_rna = dx_rna[dx_rna.num_dx > 300]

dx_rna.frac.hist(bins=30)

dx_rna.ix[['ADH1A','ADH1B','ADH1C']]

dx_rna.shape

dx_rna.p.rank().ix[['ADH1A','ADH1B','ADH1C']]

dx_rna.sort('p').head(10)

paired_bp_tn_split(matched_rna.ix['ADH1B'], codes, data_type='mRNA')

"""
Explanation: Run a simple screen for DX probes

Here we take the matched data and run a basic screen
fc = 1 means that we have no foldchange buffer for a gene to be considered over or underexpressed in a patient 
If there are ties or missing data, I omit these from the test.  This can cause underpowered tests which have extreme test statistics but weak p-values.  For this reason I filter all gene/probes/markers with a sample size of less than 300 patients.
End of explanation
"""


gs2 = gene_sets.ix[dx_rna.index].fillna(0)

rr = screen_feature(dx_rna.frac, rev_kruskal, gs2.T, 
                    align=False)
fp = (1.*gene_sets.T * dx_rna.frac).T.dropna().replace(0, np.nan).mean().order()
fp.name = 'mean frac'

"""
Explanation: Pathway and Gene Annotation Analysis
End of explanation
"""


rr.ix[ti(fp > .5)].join(fp).sort('p').head()

"""
Explanation: Overexpressed pathways
End of explanation
"""


rr.ix[ti(fp < .5)].join(fp).sort('p').head()

"""
Explanation: Underexpressed pathways
End of explanation
"""


def fig_1f(ax):
    v = pd.concat([dx_rna.frac, 
                   dx_rna.frac.ix[ti(gs2['REACTOME_CELL_CYCLE_MITOTIC']>0)],
                   dx_rna.frac.ix[ti(gs2['KEGG_FATTY_ACID_METABOLISM']>0)]])

    v1 = pd.concat([pd.Series('All Genes', dx_rna.frac.index), 
                    pd.Series('Cell Cycle\nMitotic',
                              ti(gs2['REACTOME_CELL_CYCLE_MITOTIC']>0)),
                    pd.Series('Fatty Acid\nMetabolism', 
                              ti(gs2['KEGG_FATTY_ACID_METABOLISM']>0))])
    v1.name = ''
    v.name = 'Fraction Overexpressed'
    violin_plot_pandas(v1, v, ann=None, ax=ax)
    prettify_ax(ax)
    return ax

#Do not import
fig, ax = subplots(1,1, figsize=(5,3))
fig_1f(ax);

"""
Explanation: I am folling up on Fatty Acid Metabolism as opposed to biological oxidations, because it has a larger effect size, although the smaller gene-set size gives it a less extreme p-value.
End of explanation
"""


