import NotebookImport
from metaPCNA import *

"""
Explanation: Exploration of GABA Receptors
End of explanation
"""


f_win.order().tail()

switch_plot('GABRD')

"""
Explanation: GABRD is the top hit for tumor-associated, proliferation independent genes.
End of explanation
"""


gabr = [g for g in rna_df.index if g.startswith('GABR')]
f = dx_rna.ix[gabr].dropna()
f.join(f_win).sort(f_win.name)

paired_boxplot_tumor_normal(matched_rna.ix[gabr].clip(-9,10).T, 
                            order=list(f.frac.order().index))
prettify_ax(plt.gca())

paired_boxplot_tumor_normal(matched_rna.ix[gabr, ti(codes.str.startswith('KIRC'))].clip(-9,10).T, 
                            order=list(f.frac.order().index), sig=False)
prettify_ax(plt.gca())

paired_boxplot_tumor_normal(matched_rna.ix[gabr, ti(codes.str.startswith('BRCA'))].clip(-9,10).T, 
                            order=list(f.frac.order().index), sig=False)
prettify_ax(plt.gca())

"""
Explanation: GABRD is one of many GABA receptor subunits. Many of which are down-regulated.
End of explanation
"""


switch_plot('GABRA2')

switch_plot('GABRA4')

"""
Explanation: GABRA2 and GABRA4 seem to be downregulated most of the time.
End of explanation
"""


rr = df_s.corrwith(df_s.ix['GABRD'], 1).dropna()
enr = screen_feature(rr.order(), rev_kruskal, gene_sets.T, align=False)
enr.head(4)

violin_plot_pandas(gene_sets['NABA_CORE_MATRISOME'], rr.replace(1, np.nan))

"""
Explanation: Looking for gene sets concordinant with GABRD differential expression.
End of explanation
"""


RERUN = False
try:
    assert RERUN == False
    corr_profile_t = pd.read_hdf(RNA_SUBREAD_STORE, 'dfs_correlations')
    corr_profile_m = pd.read_hdf(RNA_SUBREAD_STORE, 'dx_correlations')
    corr_change = corr_profile_m - corr_profile_t
except:
    corr_profile_t = df_s.T.corr()
    corr_profile_t.to_hdf(RNA_SUBREAD_STORE, 'dfs_correlations')

    corr_profile_m = dx.T.corr()
    corr_profile_m.to_hdf(RNA_SUBREAD_STORE, 'dx_correlations')

rr = screen_feature(corr_change.ix['GABRD'], rev_kruskal, gene_sets.T, align=False)
rr.head()

%%prun
tstats = {}
for gs in rr[rr.q < 10e-10].index:
    mm = corr_change.groupby(gene_sets[gs]).mean()
    vv = corr_change.groupby(gene_sets[gs]).var()
    tt = ((mm.ix[1] - mm.ix[0]) / 
          np.sqrt(vv.div(gene_sets[gs].value_counts(), axis=0)).sum())
    tt = tt.dropna().order()
    tstats[gs] = tt
tstats = pd.DataFrame(tstats)

tstats.rank(ascending=False).ix['GABRD'].order().head(10)

tstats['BIOCARTA_NKT_PATHWAY'].order().tail()

corr_change.ix['GABRD'].ix[ti(gene_sets['BIOCARTA_NKT_PATHWAY'] > 0)].order()

pd.crosstab(gene_sets['BIOCARTA_NKT_PATHWAY'], 
            gene_sets['KEGG_HEMATOPOIETIC_CELL_LINEAGE'])

violin_plot_pandas(gene_sets['KEGG_HEMATOPOIETIC_CELL_LINEAGE'], 
                   corr_change.ix['GABRD'].dropna())

violin_plot_pandas(gene_sets['BIOCARTA_NKT_PATHWAY'], corr_change.ix['GABRD'].dropna())

gs = 'BIOCARTA_NKT_PATHWAY'
fig, ax = subplots()
series_scatter(corr_profile_m['GABRD'], corr_profile_t['GABRD'], s=10, ax=ax, ann=None,
               alpha=.1)
series_scatter(corr_profile_m['GABRD'], corr_profile_t['GABRD'].ix[ti(gene_sets[gs]>0)], 
               s=10, ax=ax, color='red', alpha=1, ann=None)

cc = combine(gene_sets['KEGG_HEMATOPOIETIC_CELL_LINEAGE'],
            gene_sets['BIOCARTA_NKT_PATHWAY'])
violin_plot_pandas(cc, corr_change.ix['GABRD'].dropna())

"""
Explanation: I'm storing intermediate results here because they can take a bit to compute. If anything gets changed, you are going to want to restore the results by setting the RERUN flag.
End of explanation
"""


