Created
August 13, 2018 21:16
-
-
Save cplaisier/3654cdb0de6da6fa9a2bd7fca4aadff4 to your computer and use it in GitHub Desktop.
Look at lines 29, 30, 263-280, 342-353, 371-382
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ########################################################## | |
| ## Consistilator: plotNetworkMotifs.py ## | |
| ## ______ ______ __ __ ## | |
| ## /\ __ \ /\ ___\ /\ \/\ \ ## | |
| ## \ \ __ \ \ \___ \ \ \ \_\ \ ## | |
| ## \ \_\ \_\ \/\_____\ \ \_____\ ## | |
| ## \/_/\/_/ \/_____/ \/_____/ ## | |
| ## @Developed by: Plaisier Lab ## | |
| ## (https://plaisierlab.engineering.asu.edu/) ## | |
| ## Arizona State University ## | |
| ## 242 ISTB1, 550 E Orange St ## | |
| ## Tempe, AZ 85281 ## | |
| ## @Author: Chris Plaisier ## | |
| ## @License: GNU GPLv3 ## | |
| ## ## | |
| ## If this program is used in your analysis please ## | |
| ## mention who built it. Thanks. :-) ## | |
| ########################################################## | |
| from subprocess import * | |
| import numpy as np | |
| import pandas as pd | |
| import networkx as nx | |
| #from networkx.algorithms import community | |
| from networkx.algorithms import clique | |
| import matplotlib.pyplot as plt | |
| from scipy.stats import pearsonr | |
| from lifelines import KaplanMeierFitter | |
| from lifelines.statistics import pairwise_logrank_test | |
| from matplotlib.backends.backend_pdf import PdfPages | |
| import matplotlib.gridspec as gridspec | |
| import boolean2 | |
| from boolean2 import util, state, network | |
| import palettable as pal | |
| #import copy | |
| #from multiprocessing import Pool, cpu_count, Manager | |
| from matplotlib.patches import FancyArrowPatch, Circle | |
| def simulation(model, trans): | |
| "One simulation step will update the transition graph" | |
| # generates all states, set limit to a value to keep only the first that many states | |
| # when limit is a number it will take the first that many initial states | |
| initializer = state.all_initial_states( model.nodes, limit=None ) | |
| # the data is the inital data, the func is the initializer | |
| for data, initfunc in initializer: | |
| model.initialize(missing=initfunc) | |
| model.iterate(100) | |
| trans.add( model.states, times=range(100) ) | |
| return trans | |
| ## Plot network information | |
| ## _____________________________________ | |
| ## | | | | | | |
| ## | NetMot | R^2 | And | Or | | |
| ## | | | Attr. | Attr. | | |
| ## _____________________________________ | |
| ## | | | Surv. | Surv. | [ TCGA LGG_GBM | |
| ## | Dist | Surv. | And | Or | X3 [ REMBRANDT | |
| ## | States | States | Attr. | Attr. | [ French, et al. | |
| ## _____________________________________ | |
| ## | | | Surv. | Surv. | [ TCGA LGG_GBM | |
| ## | Dist | Surv. | And | Or | X3 [ REMBRANDT | |
| ## | States | States | Attr. | Attr. | [ French, et al. | |
| ## _____________________________________ | |
| #def plotNetworks(geneSets, ids, consistencies, networks, rsq, exp, binExp, pheno, symbol2entrez): | |
| #def plotNetworks(geneSets, ids, consistencies, networks, rules_and, rules_or, rsq, exp, binExp, pheno, exp_lgg, binExp_lgg, pheno_lgg, symbol2entrez): | |
| def plotNetworks(geneSets, ids, consistencies, networks, rules_and, rules_or, rsq, binExp_lgg, pheno_lgg, binExp_sc, samples, symbol2entrez): | |
| nodeColors = ['k','r','g'] | |
| #pp = PdfPages('gbmNetMotifs_3node_scRNA_seq_grade.pdf') | |
| for set1 in range(len(geneSets)): | |
| fig = plt.figure(figsize=(20,17)) | |
| plt.rcParams['legend.fontsize'] = 8 | |
| grid = plt.GridSpec(7,5, wspace=0.25, hspace=0.25, left=0.075, right=0.95, bottom=0.05, top=0.95) | |
| # Make boolean networks first to ensure order is conserved across all aspects | |
| # And | |
| model_and = boolean2.Model( text=rules_and[set1], mode='sync') | |
| trans_and = network.TransGraph( logfile='threenodes.log', verbose=True ) | |
| simulation( model_and, trans_and) | |
| att_and_graph = nx.DiGraph() | |
| for state_edge in trans_and.graph.adjacency(): | |
| att_and_graph.add_edge(state_edge[0],state_edge[1].keys()[0]) | |
| print state_edge[0],state_edge[1].keys()[0] | |
| attractors_and = sorted(nx.connected_components(att_and_graph.to_undirected()), key=len, reverse=True) | |
| # Or | |
| model_or = boolean2.Model( text=rules_or[set1], mode='sync') | |
| trans_or = network.TransGraph( logfile='threenodes.log', verbose=True ) | |
| simulation (model_or, trans_or) | |
| att_or_graph = nx.DiGraph() | |
| for state_edge in trans_or.graph.adjacency(): | |
| att_or_graph.add_edge(state_edge[0],state_edge[1].keys()[0]) | |
| print state_edge[0],state_edge[1].keys()[0] | |
| attractors_or = sorted(nx.connected_components(att_or_graph.to_undirected()), key=len, reverse=True) | |
| if not model_and.states[0].keys()==model_or.states[0].keys(): | |
| print geneSets[set1], model_and.states[0].keys(), model_or.states[0].keys() | |
| break | |
| nodes = dict(zip(model_and.states[0].keys(),nodeColors)) | |
| # Network plot [0,0] | |
| plt.subplot(grid[0,0]) | |
| G = networks[set1] | |
| node_colors = [nodes[i] for i in list(G.nodes)] | |
| edges = G.edges() | |
| edge_colors = [G[u][v]['color'] for u,v in edges] | |
| nx.draw(G, pos=nx.circular_layout(G),label_pos=3,with_labels=False,node_size=500,node_color=node_colors,edge_color=edge_colors,width=3,arrowsize=25,font_color='w') | |
| #createing label offset on network figure | |
| label_ratio = 0.27 | |
| pos_labels = {} | |
| #For each node in the Graph | |
| pos = nx.circular_layout(G) | |
| p=0 | |
| for aNode in G.nodes(): | |
| #Get the node's position from the layout | |
| x,y = pos[aNode] | |
| #Set Offset | |
| if p==0: | |
| pos_labels[aNode] = (x-label_ratio, y) | |
| else: | |
| pos_labels[aNode] = (x+label_ratio, y) | |
| p+=1 | |
| nx.draw_networkx_labels(G,pos=pos_labels,fontsize=3) | |
| plt.title('Network Motif',fontdict={'fontsize':8}) | |
| # Plot R squared values [0,1] | |
| ax = plt.subplot(grid[0,1]) | |
| index = np.arange(len(model_and.states[0].keys())) | |
| print index | |
| plt.xticks(index, model_and.states[0].keys()) | |
| g1, g2, g3 = plt.bar(index, [consistencies[set1]['all'][j] for j in model_and.states[0].keys()]) | |
| print consistencies[set1]['all'] | |
| plt.ylim((0,1)) | |
| plt.ylabel('$R^2$') | |
| #plt.xlabel('Gene') | |
| ax.axhline(rsq,color='k',linestyle='--',alpha=0.5) | |
| g1.set_facecolor(nodes[list(model_and.states[0].keys())[0]]) | |
| g2.set_facecolor(nodes[list(model_and.states[0].keys())[1]]) | |
| g3.set_facecolor(nodes[list(model_and.states[0].keys())[2]]) | |
| plt.title('$R^2$ for Inputs',fontdict={'fontsize':8}) | |
| plt.tight_layout() | |
| # And attractors [0,2] | |
| print attractors_and | |
| plt.subplot(grid[0,2]) | |
| plt.xticks([]) | |
| plt.yticks([]) | |
| positions={'000':[-1,-1],'001':[-1,0],'010':[-1,1],'100':[0,-1],'110':[0,0],'011':[0,1],'101':[1,-1],'111':[1,0]} | |
| pos = nx.spring_layout(att_and_graph,k=0.8,pos=positions,iterations=10) | |
| #pos = nx.nx_agraph.graphviz_layout(att_and_graph, prog='dot') ### TODO ### See if we can get this to work. | |
| selfies = list(att_and_graph.nodes_with_selfloops()) | |
| nx.draw_networkx_nodes(att_and_graph,pos,nodelist=[i for i in att_and_graph.nodes if not i in selfies],node_color='w',node_size=100) | |
| nx.draw_networkx_nodes(att_and_graph,pos,nodelist=selfies,node_color='#fa9fb5',node_size=100) | |
| nx.draw_networkx_edges(att_and_graph,pos,fontsize=3, alpha=0.5) | |
| nx.draw_networkx_labels(att_and_graph,pos,fontsize=3,font_color='k') | |
| plt.title('AND Attractors',fontdict={'fontsize':8}) | |
| # Or attractors [0,3] | |
| print attractors_or | |
| plt.subplot(grid[0,3]) | |
| plt.xticks([]) | |
| plt.yticks([]) | |
| positions={'000':[-1,-1],'001':[-1,0],'010':[-1,1],'100':[0,-1],'110':[0,0],'011':[0,1],'101':[1,-1],'111':[1,0]} | |
| pos = nx.spring_layout(att_or_graph,k=0.8,pos=positions,iterations=10) | |
| selfies = list(att_or_graph.nodes_with_selfloops()) | |
| nx.draw_networkx_nodes(att_or_graph,pos,nodelist=[i for i in att_or_graph.nodes if not i in selfies],node_color='w',node_size=100) | |
| nx.draw_networkx_nodes(att_or_graph,pos,nodelist=selfies,node_color='#fa9fb5',node_size=100) | |
| nx.draw_networkx_edges(att_or_graph,pos,fontsize=3, alpha=0.5) | |
| nx.draw_networkx_labels(att_or_graph,pos,fontsize=3,font_color='k') | |
| plt.title('OR Attractors',fontdict={'fontsize':8}) | |
| # Make binary data for each | |
| states = ['000', '100','010','001','110','011','101','111'] | |
| states_colors = dict(zip(states,['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown', 'tab:pink', 'tab:gray'])) | |
| plt.tight_layout() | |
| ### LGG+GBM samples | |
| cohorts = ['TCGA','REMBRANDT','French'] | |
| for cohort in [0,1,2]: | |
| print cohort, cohorts[cohort] | |
| # Bin for LGG+GBM | |
| if len([i for i in binExp_lgg[cohorts[cohort]].index if i in model_and.states[0].keys()])==3: | |
| tmp = binExp_lgg[cohorts[cohort]].loc[model_and.states[0].keys()].transpose() | |
| tmp1 = dict(zip(tmp.index,[''.join([str(int(k)) for k in tmp.loc[j].values]) for j in tmp.index])) | |
| tmp2 = [] | |
| for j in pheno_lgg[cohorts[cohort]].index: | |
| if not j in tmp1: | |
| tmp2.append(np.nan) | |
| else: | |
| tmp2.append(tmp1[j]) | |
| pheno2 = pheno_lgg[cohorts[cohort]].assign(binSet=tmp2) | |
| # Distribution of LGG+GBM states [1,0][0] as heatmap | |
| grades = ['All','IV','III','II'] | |
| #ax = plt.subplot(inner[0]) | |
| groups = pheno2['binSet'] | |
| dist1 = groups.value_counts() | |
| all1 = [] | |
| total = float(sum(list(dist1))) | |
| for i in states: | |
| if i in dist1: | |
| all1.append(float(dist1[i])/total) | |
| else: | |
| all1.append(0) | |
| # gIV | |
| subset1 = pheno2.loc[pheno2['grade']==4].index | |
| groups = pheno2['binSet'].loc[subset1] | |
| dist1 = groups.value_counts() | |
| gIV = [] | |
| total = float(sum(list(dist1))) | |
| for i in states: | |
| if i in dist1: | |
| gIV.append(float(dist1[i])/total) | |
| else: | |
| gIV.append(0) | |
| # gIII | |
| subset1 = pheno2.loc[pheno2['grade']==3].index | |
| groups = pheno2['binSet'].loc[subset1] | |
| dist1 = groups.value_counts() | |
| gIII = [] | |
| total = float(sum(list(dist1))) | |
| for i in states: | |
| if i in dist1: | |
| gIII.append(float(dist1[i])/total) | |
| else: | |
| gIII.append(0) | |
| # gII | |
| subset1 = pheno2.loc[pheno2['grade']==2].index | |
| groups = pheno2['binSet'].loc[subset1] | |
| dist1 = groups.value_counts() | |
| gII = [] | |
| total = float(sum(list(dist1))) | |
| for i in states: | |
| if i in dist1: | |
| gII.append(float(dist1[i])/total) | |
| else: | |
| gII.append(0) | |
| # Heatmap | |
| props1 = np.array([all1, gIV, gIII, gII]) | |
| ax = plt.subplot(grid[cohort+1,0]) | |
| im = ax.imshow(props1, cmap=plt.get_cmap('Reds'), vmin=0, vmax=0.5, aspect='auto') | |
| ax.set_xticks(np.arange(len(states))) | |
| ax.set_yticks(np.arange(len(grades))) | |
| ax.set_xticklabels(states) | |
| ax.set_yticklabels(grades) | |
| plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") | |
| for i in range(len(grades)): | |
| for j in range(len(states)): | |
| if props1[i,j]>0.25: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.2g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| ax.set_title(cohorts[cohort],fontdict={'fontsize':8}) | |
| # LGG+GBM Survival [1,1] | |
| ax = plt.subplot(grid[cohort+1,1]) | |
| kmf = KaplanMeierFitter() | |
| pheno3 = pheno2[['SURVIVAL','DEAD','binSet']].dropna(axis='rows') | |
| T = pheno3['SURVIVAL'] | |
| E = pheno3['DEAD']=='Dead' | |
| groups = pheno3['binSet'] | |
| dist1 = groups.value_counts() | |
| for k in states: | |
| if k in dist1: | |
| ix = (groups==k) | |
| if len(T[ix]) > 0: | |
| kmf.fit(T[ix], E[ix], label=k) | |
| kmf.plot(ax=ax, ci_show=False, color=states_colors[k]) | |
| plt.title('Survival Bin. States ('+cohorts[cohort]+')',fontdict={'fontsize':8}) | |
| print pairwise_logrank_test(T, groups, E) | |
| tmp = pairwise_logrank_test(T, groups, E) | |
| import pdb; pdb.set_trace() | |
| # Plot both attractors | |
| inner = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=grid[cohort+1,2], wspace=0.3, hspace=0.7) | |
| ## Attractor distributions | |
| # And | |
| attrSets = [', '.join(list(i)) for i in attractors_and] | |
| attrDict = dict(zip(attrSets,attractors_and)) | |
| tmp1 = [] | |
| for i in pheno2['binSet']: | |
| for j in attrDict: | |
| if i in attrDict[j]: | |
| tmp1.append(j) | |
| attrAnd = [float(i)/float(len(tmp1)) for i in pd.Series(tmp1).value_counts()] | |
| # Heatmap | |
| ax = plt.subplot(inner[0]) | |
| props1 = np.array([attrAnd]) | |
| print props1 | |
| im = ax.imshow(props1, cmap=plt.get_cmap('Reds'), vmin=0, vmax=1, aspect='auto') | |
| ax.set_xticks(np.arange(len(attrSets))) | |
| ax.set_yticks(np.arange(1)) | |
| ax.set_xticklabels(attrSets) | |
| ax.set_yticklabels(['And']) | |
| plt.setp(ax.get_xticklabels(), rotation=15, ha="center", rotation_mode="anchor") | |
| for i in range(1): | |
| for j in range(len(attrSets)): | |
| if props1[i,j]>0.7: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.2g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| ax.set_title(cohorts[cohort],fontdict={'fontsize':8}) | |
| # Or | |
| attrSets = [', '.join(list(i)) for i in attractors_or] | |
| attrDict = dict(zip(attrSets,attractors_or)) | |
| tmp1 = [] | |
| for i in pheno2['binSet']: | |
| for j in attrDict: | |
| if i in attrDict[j]: | |
| tmp1.append(j) | |
| attrOr = [float(i)/float(len(tmp1)) for i in pd.Series(tmp1).value_counts()] | |
| # Heatmap | |
| ax = plt.subplot(inner[1]) | |
| props1 = np.array([attrOr]) | |
| print props1 | |
| im = ax.imshow(props1, cmap=plt.get_cmap('Reds'), vmin=0, vmax=1, aspect='auto') | |
| ax.set_xticks(np.arange(len(attrSets))) | |
| ax.set_yticks(np.arange(1)) | |
| ax.set_xticklabels(attrSets) | |
| ax.set_yticklabels(['Or']) | |
| plt.setp(ax.get_xticklabels(), rotation=15, ha="center", rotation_mode="anchor") | |
| for i in range(1): | |
| for j in range(len(attrSets)): | |
| if props1[i,j]>0.7: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.2g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| #ax.set_title(cohorts[cohort],fontdict={'fontsize':8}) | |
| # LGG+GBM Survival grouped by AND attractor states [1,1] | |
| ax = plt.subplot(grid[cohort+1,3]) | |
| kmf = KaplanMeierFitter() | |
| pheno3 = pheno2[['SURVIVAL','DEAD','binSet']].dropna(axis='rows') | |
| T = pheno3['SURVIVAL'] | |
| E = pheno3['DEAD']=='Dead' | |
| groups = pheno3['binSet'] | |
| dist1 = groups.value_counts() | |
| attrSets = [', '.join(list(i)) for i in attractors_and] | |
| attrDict = dict(zip(attrSets,attractors_and)) | |
| attrColors = dict(zip(attrSets,pal.tableau.Tableau_10.hex_colors[0:len(attractors_and)])) | |
| print 'andAttractors', pairwise_logrank_test(T, groups, E) | |
| for k in attrDict: | |
| ix = [] | |
| for i in groups: | |
| if not i in attrDict[k]: | |
| ix.append(False) | |
| else: | |
| ix.append(True) | |
| #print ix | |
| #print k, ix, T[ix], E[ix] | |
| if len(T[ix])>0: | |
| kmf.fit(T[ix], E[ix], label=k) | |
| kmf.plot(ax=ax, ci_show=False, color=attrColors[k]) | |
| plt.title('Survival AND Attractor States ('+cohorts[cohort]+')',fontdict={'fontsize':8}) | |
| #print dir(ax) | |
| #print ax._get_lines().get_next_color() | |
| # LGG+GBM Survival grouped by OR attractor states [1,1] | |
| ax = plt.subplot(grid[cohort+1,4]) | |
| kmf = KaplanMeierFitter() | |
| pheno3 = pheno2[['SURVIVAL','DEAD','binSet']].dropna(axis='rows') | |
| T = pheno3['SURVIVAL'] | |
| E = pheno3['DEAD']=='Dead' | |
| groups = pheno3['binSet'] | |
| dist1 = groups.value_counts() | |
| attrSets = [', '.join(list(i)) for i in attractors_or] | |
| attrDict = dict(zip(attrSets,attractors_or)) | |
| attrColors = dict(zip(attrSets,pal.tableau.Tableau_10.hex_colors[0:len(attractors_or)])) | |
| print 'orAttractors', pairwise_logrank_test(T, groups, E) | |
| for k in attrDict: | |
| ix = [] | |
| for i in groups: | |
| if not i in attrDict[k]: | |
| ix.append(False) | |
| else: | |
| ix.append(True) | |
| #print ix | |
| #print k, ix, T[ix], E[ix] | |
| if len(T[ix])>0: | |
| kmf.fit(T[ix], E[ix], label=k) | |
| kmf.plot(ax=ax, ci_show=False, color=attrColors[k]) | |
| plt.title('Survival OR Attractor States ('+cohorts[cohort]+')',fontdict={'fontsize':8}) | |
| ## Plot single cell data | |
| #scData = ['GSE57872','GSE84465','GSE89567','GSE102130','GSE70630'] | |
| scData = ['IV','III','II','H3K27M'] | |
| for dataset in range(len(scData)): | |
| sources = sorted(list(set(samples[scData[dataset]]))) | |
| print len(sources), sources | |
| # Concatenate states | |
| tmp = binExp_sc[scData[dataset]].loc[model_and.states[0].keys()].transpose().dropna() | |
| inc_samps = [i for i in range(len(binExp_sc[scData[dataset]].columns)) if binExp_sc[scData[dataset]].columns[i] in tmp.index] | |
| #print scData[dataset], inc_samps | |
| tmp_samples = [samples[scData[dataset]][i] for i in inc_samps] | |
| if tmp.shape[0]>0: | |
| tmp1 = tmp.assign(binSet=[''.join([str(int(k)) for k in tmp.loc[j].values]) for j in tmp.index]) | |
| # Distribution of LGG+GBM states [1,0][0] | |
| allLines = [] | |
| allLines_and = [] | |
| allLines_or = [] | |
| for s1 in range(len(sources)): | |
| cur = [i for i in range(len(tmp_samples)) if tmp_samples[i]==sources[s1]] | |
| if len(cur)>0: | |
| #print cur, tmp1['binSet'].shape | |
| groups = tmp1['binSet'].iloc[cur] | |
| dist1 = groups.value_counts() | |
| dist2 = [] | |
| total = float(sum(list(dist1))) | |
| for i in states: | |
| if i in dist1: | |
| dist2.append(float(dist1[i])/total) | |
| else: | |
| dist2.append(0) | |
| allLines.append(dist2) | |
| # And | |
| attrSets_and = [', '.join(list(i)) for i in attractors_and] | |
| attrDict = dict(zip(attrSets_and,attractors_and)) | |
| tmp2 = [] | |
| for i in tmp1['binSet'].iloc[cur]: | |
| for j in attrDict: | |
| if i in attrDict[j]: | |
| tmp2.append(j) | |
| attrAnd = [] | |
| dist1 = pd.Series(tmp2).value_counts() | |
| total = float(sum(list(dist1))) | |
| for i in attrSets_and: | |
| if i in dist1: | |
| attrAnd.append(float(dist1[i])/total) | |
| else: | |
| attrAnd.append(0) | |
| allLines_and.append(attrAnd) | |
| # Or | |
| attrSets_or = [', '.join(list(i)) for i in attractors_or] | |
| attrDict = dict(zip(attrSets_or,attractors_or)) | |
| tmp2 = [] | |
| for i in tmp1['binSet'].iloc[cur]: | |
| for j in attrDict: | |
| if i in attrDict[j]: | |
| tmp2.append(j) | |
| attrOr = [] | |
| dist1 = pd.Series(tmp2).value_counts() | |
| total = float(sum(list(dist1))) | |
| for i in attrSets_or: | |
| if i in dist1: | |
| attrOr.append(float(dist1[i])/total) | |
| else: | |
| attrOr.append(0) | |
| allLines_or.append(attrOr) | |
| else: | |
| allLines.append([np.NaN,np.NaN,np.NaN,np.NaN,np.NaN,np.NaN,np.NaN,np.NaN]) | |
| allLines_and.append([np.NaN for i in range(len(attractors_and))]) | |
| allLines_or.append([np.NaN for i in range(len(attractors_or))]) | |
| # Heatmap | |
| props1 = np.array(allLines) | |
| print(props1) | |
| ax = plt.subplot(grid[4,dataset]) | |
| masked_array = np.ma.array(allLines, mask=np.isnan(allLines)) | |
| print masked_array | |
| cmap = plt.get_cmap('Reds') | |
| cmap.set_bad('white',1.) | |
| im = ax.imshow(masked_array, interpolation='nearest', cmap=cmap, vmin=0, vmax=0.5, aspect='auto') | |
| #im = ax.imshow(props1, cmap=plt.get_cmap('Reds')) | |
| ax.set_xticks(np.arange(len(states))) | |
| ax.set_yticks(np.arange(len(sources))) | |
| ax.set_xticklabels(states) | |
| ax.set_yticklabels(sources) | |
| plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") | |
| for i in range(len(sources)): | |
| for j in range(len(states)): | |
| if props1[i,j]>0.25: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.1g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| ax.set_title(scData[dataset],fontdict={'fontsize':8}) | |
| ## Attractor distributions | |
| # Heatmap | |
| print allLines | |
| print allLines_and | |
| props1 = np.array(allLines_and) | |
| #print props1 | |
| ax = plt.subplot(grid[5,dataset]) | |
| #print allLines_and | |
| masked_array = np.ma.array (allLines_and, mask=np.isnan(allLines_and)) | |
| cmap = plt.get_cmap('Reds') | |
| cmap.set_bad('white',1.) | |
| im = ax.imshow(masked_array, interpolation='nearest', cmap=cmap, vmin=0, vmax=0.5, aspect='auto') | |
| #im = ax.imshow(props1, cmap=plt.get_cmap('Reds'), vmin=0, vmax=1, aspect='auto') | |
| ax.set_xticks(np.arange(len(attrSets_and))) | |
| ax.set_yticks(np.arange(len(sources))) | |
| ax.set_xticklabels(attrSets_and) | |
| ax.set_yticklabels(sources) | |
| plt.setp(ax.get_xticklabels(), rotation=15, ha="center", rotation_mode="anchor") | |
| for i in range(len(sources)): | |
| for j in range(len(attrSets_and)): | |
| if props1[i,j]>0.7: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.2g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| ax.set_title(scData[dataset]+' AND',fontdict={'fontsize':8}) | |
| # Heatmap | |
| #props1 = np.array(allLines_or) | |
| #print props1 | |
| ax = plt.subplot(grid[6,dataset]) | |
| masked_array = np.ma.array (allLines_or, mask=np.isnan(allLines_or)) | |
| cmap = plt.get_cmap('Reds') | |
| cmap.set_bad('white',1.) | |
| im = ax.imshow(masked_array, interpolation='nearest', cmap=cmap, vmin=0, vmax=0.5, aspect='auto') | |
| #im = ax.imshow(props1, cmap=plt.get_cmap('Reds'), vmin=0, vmax=1, aspect='auto') | |
| ax.set_xticks(np.arange(len(attrSets_or))) | |
| ax.set_yticks(np.arange(len(sources))) | |
| ax.set_xticklabels(attrSets_or) | |
| ax.set_yticklabels(sources) | |
| plt.setp(ax.get_xticklabels(), rotation=15, ha="center", rotation_mode="anchor") | |
| for i in range(len(sources)): | |
| for j in range(len(attrSets_or)): | |
| if props1[i,j]>0.7: | |
| col1 = 'w' | |
| else: | |
| col1 = 'k' | |
| text = ax.text(j, i, '%s' % float('%.2g' % props1[i, j]), ha="center", va="center", color=col1, fontdict={'fontsize':8}) | |
| ax.set_title(scData[dataset]+' OR',fontdict={'fontsize':8}) | |
| plt.show() | |
| #pp.savefig(fig) | |
| break | |
| #pp.close() | |
| # Read in Biotapesty file | |
| inEdges = {} | |
| inFile = open('biotapestry_GBM.csv','r') | |
| while 1: | |
| line = inFile.readline() | |
| if not line: | |
| break | |
| splitUp = line.strip().split(',') | |
| if splitUp[0]=='"# Standard Interactions"': | |
| inFile.readline() # Remove header | |
| break | |
| while 1: | |
| line = inFile.readline() | |
| if not line: | |
| break | |
| splitUp = line.strip().split(',') | |
| node1 = splitUp[3].strip('"') | |
| node2 = splitUp[5].strip('"') | |
| for i in node2.split(';'): | |
| if not i in inEdges: | |
| inEdges[i] = {} | |
| if not node1 in inEdges[i]: | |
| inEdges[i][node1] = splitUp[6].strip('"') | |
| inFile.close() | |
| # Load up genesets and netMatrices | |
| data = {} | |
| consistent = {} | |
| gene2probe = {} | |
| subsets = ['all'] | |
| rsq = 0.8 | |
| for subset in subsets: | |
| inFile = open('results_'+subset+'.csv','r') | |
| inFile.readline() # Get rid of header | |
| while 1: | |
| line = inFile.readline() | |
| if not line: | |
| break | |
| splitUp = line.strip().split(',') | |
| if not splitUp[1] in data: | |
| data[splitUp[1]] = {} | |
| consistent[splitUp[1]] = {} | |
| if not splitUp[3] in data[splitUp[1]]: | |
| data[splitUp[1]][splitUp[3]] = {} | |
| consistent[splitUp[1]][splitUp[3]] = {} | |
| data[splitUp[1]][splitUp[3]][subset] = splitUp | |
| if not (splitUp[6]=='Inconsistent' or splitUp[9]=='Inconsistent' or splitUp[12]=='Inconsistent') and (splitUp[6]=='NA' or float(splitUp[6])>=rsq) and (splitUp[9]=='NA' or float(splitUp[9])>=rsq) and (splitUp[12]=='NA' or float(splitUp[12])>=rsq): | |
| consistent[splitUp[1]][splitUp[3]][subset] = 'Yes' | |
| else: | |
| consistent[splitUp[1]][splitUp[3]][subset] = 'No' | |
| inFile.close() | |
| # To translate gene ids later | |
| symbol2entrez = {} | |
| entrez2symbol = {} | |
| with open('gene2entrezId.csv','r') as inFile: | |
| while 1: | |
| inLine = inFile.readline() | |
| if not inLine: | |
| break | |
| split = inLine.strip().split(',') | |
| entrez2symbol[split[1]] = split[0] | |
| symbol2entrez[split[0]] = split[1] | |
| # Read in binarized data | |
| binExp_lgg = {} | |
| binExp_lgg['TCGA'] = pd.read_csv('binGeneExp_Vgt0_Median_lgg.csv', header=0, index_col=0)#.transpose() | |
| binExp_lgg['TCGA'].columns = [i.replace('.','-') for i in binExp_lgg['TCGA'].columns.values] | |
| binExp_lgg['REMBRANDT'] = pd.read_csv('gbm_lgg/binGeneExp_REMBRANDT_VMedian.csv', header=0, index_col=0)#.transpose() | |
| binExp_lgg['REMBRANDT'] = binExp_lgg['REMBRANDT'].dropna() | |
| binExp_lgg['French'] = pd.read_csv('gbm_lgg/binGeneExp_French_VMedian.csv', header=0, index_col=0)#.transpose() | |
| binExp_lgg['French'] = binExp_lgg['French'].dropna() | |
| # Phenotypes | |
| pheno_lgg = {} | |
| pheno_lgg['TCGA'] = pd.read_csv('phenotypes_lgg.csv', header=0, index_col=0) | |
| pheno_lgg['REMBRANDT'] = pd.read_csv('gbm_lgg/phenotypes_REMBRANDT.csv', header=0, index_col=0) | |
| pheno_lgg['French'] = pd.read_csv('gbm_lgg/phenotypes_French.csv', header=0, index_col=0) | |
| # Read in single cell RNA-seq data | |
| binExp_sc = {} | |
| samples = {} | |
| tmp = pd.read_csv('singleCell/binGeneExp_scGBM_Vgt0_GSE57872.csv', header=0, index_col=0)#.transpose() | |
| binExp_sc['IV'] = tmp | |
| samples['IV'] = [i.split('_')[0] for i in tmp.columns] | |
| tmp = pd.read_csv('singleCell/binGeneExp_scGBM_Vgt0_GSE84465.csv', header=0, index_col=0)#.transpose() | |
| binExp_sc['IV'] = pd.concat([binExp_sc['IV'],tmp], axis=1) | |
| samples['IV'] = samples['IV']+['BT_S2']*1169+['BT_S1']*489+['BT_S4']*1542+['BT_S6']*389 | |
| tmp = pd.read_csv('singleCell/binGeneExp_scGBM_Vgt0_GSE89567.csv', header=0, index_col=0)#.transpose() | |
| tmp2 = [j.split('.')[0] for j in [i.split('_')[0] for i in tmp.columns]] | |
| for i in range(len(tmp2)): | |
| if tmp2[i]=='X57': | |
| tmp2[i] = 'MGH57' | |
| if tmp2[i]=='mgh103': | |
| tmp2[i] = 'MGH103' | |
| gbmIV = [i for i in range(len(tmp2)) if tmp2[i] in ['MGH45','MGH57']] | |
| binExp_sc['IV'] = pd.concat([binExp_sc['IV'],tmp[tmp.columns[gbmIV]]], axis=1) | |
| samples['IV'] = samples['IV']+[tmp2[i] for i in gbmIV] | |
| gliomaIII = [i for i in range(len(tmp2)) if tmp2[i] in ['MGH42','MGH43','MGH44','MGH56','MGH61','MGH64','MGH103']] | |
| binExp_sc['III'] = tmp[tmp.columns[gliomaIII]] | |
| samples['III'] = [tmp2[i] for i in gliomaIII] | |
| gliomaII = [i for i in range(len(tmp2)) if tmp2[i] in ['MGH107neg','MGH107pos']] | |
| binExp_sc['II'] = tmp[tmp.columns[gliomaII]] | |
| samples['II'] = [tmp2[i] for i in gliomaII] | |
| binExp_sc['H3K27M'] = pd.read_csv('singleCell/binGeneExp_scGBM_K27M_Vgt0_GSE102130.csv', header=0, index_col=0)#.transpose() | |
| samples['H3K27M'] = [i.split('_')[0].split('.')[0] for i in binExp_sc['H3K27M'].columns] | |
| tmp = pd.read_csv('singleCell/binGeneExp_scOligo_Vgt0_GSE70630.csv', header=0, index_col=0)#.transpose() | |
| tmp2 = [i.split('_')[0] for i in tmp.columns] | |
| for i in range(len(tmp2)): | |
| if tmp2[i]=='X93': | |
| tmp2[i] = 'MGH93' | |
| if tmp2[i]=='X97': | |
| tmp2[i] = 'MGH97' | |
| gliomaII = [i for i in range(len(tmp2)) if tmp2[i] in ['MGH36','MGH53','MGH54','MGH60','MGH93','MGH97']] | |
| binExp_sc['II'] = pd.concat([binExp_sc['II'],tmp[tmp.columns[gliomaII]]], axis=1) | |
| samples['II'] = samples['II']+[tmp2[i] for i in gliomaII] | |
| # Gather data | |
| geneSets = [] | |
| ids = [] | |
| netMatrices = [] | |
| consistencies = [] | |
| networks = [] | |
| rules_and = [] | |
| rules_or = [] | |
| for netMotif in data: | |
| for instance in data[netMotif]: | |
| # Only plot if significant amount of variance explained | |
| if len([i for i in subsets if consistent[netMotif][instance][i]=='Yes']) > 0: | |
| genes = data[netMotif][instance]['all'][3].split(';') | |
| if len(set(genes))==len(genes): | |
| geneSets.append(genes) | |
| ids.append('id'+data[netMotif][instance]['all'][1]+' '+data[netMotif][instance]['all'][2]) | |
| tmp = {'all':dict(zip(genes, [data[netMotif][instance]['all'][i] for i in [6,9,12]]))} | |
| for i in tmp: | |
| for j in tmp[i]: | |
| if tmp[i][j]=='Inconsistent': | |
| tmp[i][j] = -0.25 | |
| elif not tmp[i][j]=='NA': | |
| tmp[i][j] = float(tmp[i][j]) | |
| else: | |
| tmp[i][j] = 0 | |
| consistencies.append(tmp) | |
| # Make networks | |
| G = nx.DiGraph() | |
| tmp_and = "" | |
| tmp_or = "" | |
| for gene1 in genes: | |
| tmp_pos = [] | |
| tmp_neg = [] | |
| if gene1 in inEdges and len([i for i in inEdges[gene1] if (i in genes and not i==gene1)])>0: | |
| for gene2 in inEdges[gene1]: | |
| if gene2 in genes and not gene2==gene1: | |
| if inEdges[gene1][gene2]=='positive': | |
| G.add_edge(gene2,gene1,color='g') | |
| tmp_pos.append(str(gene2)) | |
| elif inEdges[gene1][gene2]=='negative': | |
| G.add_edge(gene2,gene1,color='r') | |
| tmp_neg.append(str(gene2)) | |
| else: | |
| print gene1, genes | |
| tmp_pos.append(gene1) | |
| if not len(tmp)==0: | |
| tmp_neg2_and = [] | |
| tmp_neg2_or = [] | |
| if len(tmp_neg)>0: | |
| tmp_neg2_and = ['((not ('+' and '.join(tmp_neg)+')) and '+gene1+')'] | |
| tmp_neg2_or = ['((not ('+' or '.join(tmp_neg)+')) and '+gene1+')'] | |
| tmp_and += '\n'+str(gene1)+'* = '+' and '.join(tmp_pos+tmp_neg2_and) | |
| tmp_or += '\n'+str(gene1)+'* = '+' or '.join(tmp_pos+tmp_neg2_or) | |
| rules_and.append(tmp_and) | |
| print tmp_and | |
| rules_or.append(tmp_or) | |
| print tmp_or | |
| networks.append(G) | |
| # Plot them | |
| plotNetworks(geneSets, ids, consistencies, networks, rules_and, rules_or, rsq, binExp_lgg, pheno_lgg, binExp_sc, samples, symbol2entrez) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment