<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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.
Step1:  Read in matched Gene expression data.
Step2:  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.
Step3:  Pathway and Gene Annotation Analysis
Step4:  Overexpressed pathways
Step5:  Underexpressed pathways
Step6:  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.
<ASSISTANT_TASK:>
Python Code:
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
<END_TASK>