notebook.community

Edit and run



In [21]:

    
import theano
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import MDBN



In [22]:

    
import AML
me_DBN, ge_DBN, dm_DBN, top_DBN = AML.load_network('Exp_2016-12-11_1841_run_0.npz','../MDBN_run')









    



Adding a layer with 559 input and 40 outputs
Adding a layer with 19940 input and 400 outputs
Adding a layer with 400 input and 40 outputs
Adding a layer with 80 input and 24 outputs
Adding a layer with 24 input and 3 outputs



In [23]:

    
datafiles = AML.prepare_AML_TCGA_datafiles('../data')



In [24]:

    
ME_output, _ = me_DBN.MLP_output_from_datafile(datafiles['ME'], datadir='../data')
GE_output, _ = ge_DBN.MLP_output_from_datafile(datafiles['GE'], datadir='../data')



In [6]:

    
joint_layer = np.concatenate([ME_output, GE_output],axis=1)
plt.imshow(joint_layer, cmap='gray',interpolation='none')
plt.axis('tight')









    Out[6]:





(-0.5, 79.5, 172.5, -0.5)



In [7]:

    
top_output = top_DBN.get_output(theano.shared(joint_layer,borrow=True))
plt.imshow((top_output>0.8)*np.ones_like(top_output)-(top_output<0.2)*np.ones_like(top_output),interpolation='none',extent=[0,3,385,0])
plt.colorbar()
plt.axis('tight')
plt.xticks(np.arange(0.5,3.5,1),('0','1','2'))









    Out[7]:





([<matplotlib.axis.XTick at 0x12de4e310>,
  <matplotlib.axis.XTick at 0x12de6bc10>,
  <matplotlib.axis.XTick at 0x12baaae50>],
 <a list of 3 Text xticklabel objects>)



In [8]:

    
plt.imshow(top_output, interpolation='none',extent=[0,3,385,0])
plt.axis('tight')
plt.colorbar()
plt.xticks(np.arange(0.5,3.5,1),('0','1','2'))









    Out[8]:





([<matplotlib.axis.XTick at 0x12bb77950>,
  <matplotlib.axis.XTick at 0x12de9c090>,
  <matplotlib.axis.XTick at 0x12bb37cd0>],
 <a list of 3 Text xticklabel objects>)



In [9]:

    
plt.hist(top_output)









    Out[9]:





([array([ 127.,    3.,    0.,    2.,    1.,    0.,    1.,    2.,    3.,   34.]),
  array([  53.,    2.,    0.,    3.,    3.,    3.,    1.,    5.,    2.,  101.]),
  array([ 69.,   1.,   1.,   0.,   0.,   0.,   0.,   4.,   0.,  98.])],
 array([  1.78998524e-10,   1.00000000e-01,   2.00000000e-01,
          3.00000000e-01,   4.00000000e-01,   5.00000000e-01,
          6.00000000e-01,   7.00000000e-01,   8.00000000e-01,
          9.00000000e-01,   1.00000000e+00]),
 <a list of 3 Lists of Patches objects>)



In [10]:

    
plt.hist(top_output[:,0])









    Out[10]:





(array([ 127.,    3.,    0.,    2.,    1.,    0.,    1.,    2.,    3.,   34.]),
 array([  1.78998524e-10,   9.99999406e-02,   1.99999881e-01,
          2.99999821e-01,   3.99999762e-01,   4.99999702e-01,
          5.99999642e-01,   6.99999583e-01,   7.99999523e-01,
          8.99999464e-01,   9.99999404e-01]),
 <a list of 10 Patch objects>)



In [11]:

    
code = (top_output[:,0:2] > 0.5) * np.ones_like(top_output[:,0:2])



In [12]:

    
from utils import find_unique_classes
U = find_unique_classes(code)
cl = U[0]



In [13]:

    
cl









    Out[13]:





array([ 0.,  1.,  1.,  3.,  1.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,
        1.,  1.,  3.,  3.,  1.,  0.,  0.,  1.,  1.,  1.,  2.,  0.,  1.,
        1.,  0.,  1.,  1.,  3.,  1.,  0.,  0.,  1.,  3.,  1.,  1.,  1.,
        2.,  1.,  3.,  1.,  1.,  1.,  2.,  0.,  3.,  1.,  1.,  1.,  0.,
        0.,  1.,  3.,  0.,  1.,  3.,  1.,  1.,  1.,  0.,  1.,  1.,  1.,
        0.,  0.,  1.,  1.,  3.,  1.,  3.,  1.,  2.,  1.,  1.,  3.,  1.,
        1.,  1.,  0.,  3.,  3.,  1.,  3.,  1.,  3.,  2.,  2.,  0.,  1.,
        1.,  1.,  1.,  2.,  1.,  1.,  0.,  2.,  0.,  1.,  1.,  0.,  2.,
        2.,  3.,  0.,  1.,  2.,  1.,  1.,  2.,  3.,  0.,  3.,  2.,  0.,
        2.,  1.,  3.,  1.,  1.,  1.,  1.,  1.,  0.,  1.,  3.,  0.,  1.,
        1.,  1.,  1.,  3.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  1.,  0.,
        0.,  0.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,
        1.,  1.,  0.,  3.,  3.,  0.,  0.,  0.,  0.,  2.,  0.,  1.,  0.,
        0.,  3.,  0.,  0.])



In [15]:

    
plt.hist(cl,bins=4)









    Out[15]:





(array([ 46.,  87.,  15.,  25.]),
 array([ 0.  ,  0.75,  1.5 ,  2.25,  3.  ]),
 <a list of 4 Patch objects>)

Check Survival curves for the different classes



In [16]:

    
import csv
id=[]
with open('../data/'+datafiles['ME']) as f:
    my_csv = csv.reader(f,delimiter='\t')
    id = my_csv.next()



In [17]:

    
stat={}
with open('../data/AML/AML_clinical_data2.csv') as f:
    reader = csv.reader(f, delimiter=',')
    for row in reader:
        patient_id=row[0]
        stat[patient_id]=(row[4],row[7],row[6])



In [18]:

    
import re
time_list = []
event_list = []
group_list = []
print('The following case IDs were  not found in clinical data')
for index, key in enumerate(id[1:]):
    m = re.match('TCGA-\w+-\d+', key)
    patient_id = m.group(0)
    if patient_id in stat:
        patient_stat = stat[patient_id]
        add_group = True
        try:
            time_list.append(float(patient_stat[2]))
            event_list.append(1)
        except ValueError:
            try:
                time_list.append(float(patient_stat[1]))
                event_list.append(0)
            except ValueError:
                print('No data for %s' % patient_id)
                add_group = False
        if add_group:
            group_list.append(cl[index])
    else:
        print(patient_id)









    



The following case IDs were  not found in clinical data
No data for TCGA-AB-2887
No data for TCGA-AB-2891
No data for TCGA-AB-2918
No data for TCGA-AB-2921
No data for TCGA-AB-2930
No data for TCGA-AB-2940
No data for TCGA-AB-2943
No data for TCGA-AB-2944
No data for TCGA-AB-2946
No data for TCGA-AB-2975



In [19]:

    
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(time_list,event_observed=event_list)
kmf.plot()









    Out[19]:





<matplotlib.axes._subplots.AxesSubplot at 0x12d90de50>



In [20]:

    
T=np.array(time_list)
E=np.array(event_list)
ix = (np.array(group_list) == 0)
kmf.fit(T[ix], E[ix], label='group 0')
ax=kmf.plot()
for i in range(1,4):
    ix=(np.array(group_list)==i)
    kmf.fit(T[ix], E[ix], label='group %d' % i)
    kmf.plot(ax=ax)



In [ ]: