Created
October 10, 2019 00:08
-
-
Save cplaisier/0e171d5521b796a934e6c5387789e6da to your computer and use it in GitHub Desktop.
Stem Cell Changes
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
| from __future__ import division | |
| ########################################################## | |
| ## Consistilator: plotNetworkMotifs_SCD.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 pyBinarize as pyBin | |
| #import copy | |
| #from multiprocessing import Pool, cpu_count, Manager | |
| from matplotlib.patches import FancyArrowPatch, Circle | |
| #import pdb | |
| import json | |
| import seaborn as sns | |
| from matplotlib.colors import ListedColormap | |
| 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 | |
| # Read in Biotapesty file | |
| inEdges = {} | |
| inFile = open('biotapestry_StemCellDiff.csv','r') | |
| while 1: | |
| line = inFile.readline() | |
| if not line: | |
| break | |
| splitUp = line.strip().split(',') | |
| if splitUp['Gene Expression']=='"# 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() | |
| # 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] | |
| # Load up genesets and netMatrices | |
| data = {} | |
| consistent = {} | |
| gene2probe = {} | |
| subsets = {'GSE20573':['all','nsb','nsbs'],'GSE26867':['all','lsb','lsb3i'],'GSE32658':['all','l','lsf','lsfc'],'GSE45223':['all','be','lsb','nc'],'GSE51533':['all','nsb','pip']} | |
| feedLoop = ['020200210'] #['000102220','020200210'] | |
| rsq = 0.90 | |
| #rList = [] | |
| #rDF = {} | |
| rDict = {} | |
| geneCount = 0 | |
| for dataset1 in subsets: | |
| # rDF[dataset1] = {} | |
| for subset1 in subsets[dataset1]: | |
| # rDF[dataset1][subset1] = {} | |
| # Gene Expression Series,Data Subset,Netowrk Motif,Adjacency Matrix,Instance,Node1,Node1.Inputs,Node1.Consistency,Node2,Node2.Inputs,Node2.Consistency,Node3,Node3.Inputs,Node3.Consistency | |
| inFile = open('results_'+dataset1+'_'+subset1+'.csv','r') | |
| header = inFile.readline().strip().split(',') # Get rid of header | |
| while 1: | |
| line = inFile.readline() | |
| if not line: | |
| break | |
| splitUp = dict(zip(header, line.strip().split(','))) | |
| if not (splitUp['Network Motif'] in []): # I am pretty sure this could be used as limiter like line 116 | |
| #if splitUp[3] in feedLoop: | |
| if not splitUp['Gene Expression Series'] in data: | |
| data[splitUp['Gene Expression Series']] = {} | |
| consistent[splitUp['Gene Expression Series']] = {} | |
| rDict[splitUp['Gene Expression Series']] = {} | |
| if not splitUp['Network Motif'] in data[splitUp['Gene Expression Series']]: | |
| data[splitUp['Gene Expression Series']][splitUp['Network Motif']] = {} | |
| consistent[splitUp['Gene Expression Series']][splitUp['Network Motif']] = {} | |
| if not splitUp['Instance'] in data[splitUp['Gene Expression Series']][splitUp['Network Motif']]: | |
| data[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']] = {} | |
| consistent[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']] = {} | |
| if not splitUp['Data Subset'] in data[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']]: | |
| data[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']][splitUp['Data Subset']] = {} | |
| consistent[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']][splitUp['Data Subset']] = {} | |
| # rDF[splitUp['Gene Expression Series']][splitUp['Data Subset']] = {} | |
| if not splitUp['Data Subset'] in rDict[splitUp['Gene Expression Series']]: | |
| rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']] = {} | |
| data[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']][splitUp['Data Subset']] = splitUp.values() | |
| if not splitUp['Instance'] in rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']]: | |
| rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']][splitUp['Instance']] = {} | |
| # if not (splitUp[7]=='Inconsistent' or splitUp[10]=='Inconsistent' or splitUp[13]=='Inconsistent') and (splitUp[7]=='NA' or float(splitUp[7])>=rsq) and (splitUp[10]=='NA' or float(splitUp[10])>=rsq) and (splitUp[13]=='NA' or float(splitUp[13])>=rsq): | |
| rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']][splitUp['Instance']][splitUp['Node1']] = [splitUp['Node1.Consistency']] | |
| rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']][splitUp['Instance']][splitUp['Node2']] = [splitUp['Node2.Consistency']] | |
| rDict[splitUp['Gene Expression Series']][splitUp['Data Subset']][splitUp['Instance']][splitUp['Node3']] = [splitUp['Node3.Consistency']] | |
| if not (splitUp['Node1.Consistency']=='Inconsistent' or splitUp['Node2.Consistency']=='Inconsistent' or splitUp['Node3.Consistency']=='Inconsistent') and (splitUp['Node1.Consistency']=='NA' or float(splitUp['Node1.Consistency'])>=rsq) and (splitUp['Node2.Consistency']=='NA' or float(splitUp['Node2.Consistency'])>=rsq) and (splitUp['Node3.Consistency']=='NA' or float(splitUp['Node3.Consistency'])>=rsq): | |
| consistent[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']][splitUp['Data Subset']] = 'Yes' | |
| else: | |
| consistent[splitUp['Gene Expression Series']][splitUp['Network Motif']][splitUp['Instance']][splitUp['Data Subset']] = 'No' | |
| # rDF[dataset1][subset1][splitUp[4]] = {} | |
| # rDF[dataset1][subset1][splitUp[4]][splitUp[5]] = splitUp[7] | |
| # rDF[dataset1][subset1][splitUp[4]][splitUp[8]] = splitUp[10] | |
| # rDF[dataset1][subset1][splitUp[4]][splitUp[11]] = splitUp[13] | |
| # rList.append(rDF) | |
| inFile.close() | |
| tfs1 = [int(symbol2entrez[j]) for j in set([i for gse1 in data for mot1 in data[gse1] for motTfs1 in data[gse1][mot1] for i in motTfs1.split(';')])] | |
| # Read in binarized data | |
| geneExp = {} | |
| binExp = {} | |
| for dataset1 in subsets: | |
| print('Loading '+dataset1+'...') | |
| geneExp[dataset1] = pd.DataFrame(pd.read_csv('../../'+dataset1+'/'+dataset1+'_genesExpMatrix_all.csv', header=0, index_col=0).loc[tfs1]) | |
| binExp[dataset1] = pyBin.binarize_kMeans_matrix(geneExp[dataset1]) | |
| # Phenotypes | |
| pheno = pd.read_csv('../metaData_StuderExperiments.csv', header=0, index_col=0) | |
| # Import series information | |
| with open('../series.json','r') as inFile: | |
| series = json.load(inFile) | |
| # Gather data | |
| geneSets = [] | |
| rgeneSets = [] | |
| ids = [] | |
| netMatrices = [] | |
| consistencies = [] | |
| networks = [] | |
| rules_and = [] | |
| rules_or = [] | |
| rGeneList = [] | |
| # CONSISTENCIES AND GENESETS, GENES ONLY -- NO R2 VALUE FOR EACH TREATMENT | |
| for gse1 in ['GSE20573','GSE26867', 'GSE45223', 'GSE51533', 'GSE20573','GSE32658']: | |
| print gse1 | |
| for netMotif in data[gse1]: | |
| print netMotif | |
| for instance in data[gse1][netMotif]: | |
| # Only plot if significant amount of variance explained | |
| if len([i for i in subsets[gse1] if consistent[gse1][netMotif][instance][i]=='Yes']) > 0: | |
| genes = data[gse1][netMotif][instance]['all'][4].split(';') # THis is OK because 'all' is in the others, too | |
| rgenes = data[gse1][netMotif][instance]['all'][4].split(';') | |
| rgene = data[gse1][netMotif][instance]['all'][4] | |
| if len(set(genes))==len(genes): | |
| #if len([i for i in genes if i in ['MYB','SMAD9','SPDEF']])==3: | |
| geneSets.append(genes) | |
| rgeneSets.append(rgenes) | |
| rGeneList.append(rgene) | |
| ids.append('id'+data[gse1][netMotif][instance]['all'][2]+' '+data[gse1][netMotif][instance]['all'][3]) | |
| tmp = {'all': dict(zip(genes, [data[gse1][netMotif][instance]['all'][r] for r in [7,10,13]]))} # consistencies - not needed because you built the dictionary which is already excluded? May need to revisit because it doesn't match with consistent | |
| 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, pheno, symbol2entrez) #, samples | |
| ## Plot network information | |
| ## _____________________________________ | |
| ## | | | | | | |
| ## | NetMot | And | Or | | | |
| ## | | Attr. | Attr. | | | |
| ## ________________________________________________ | |
| ## | | | | | | |
| ## R^2 | GSE | All | Treat | Treat | x 5 for each GSE | |
| ## | Heatmp | NI | 1 | ...n | | |
| ## ________________________________________________ | |
| ## | | | | | | |
| ## | Htmp | Bar% | Htmp | Bar% | subset of the GSE | |
| ## | And | And | Or | Or | | |
| ## _____________________________________ | |
| #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, pheno, symbol2entrez): | |
| #def plotNetworks(geneSets, ids, consistencies, networks, rules_and, rules_or, rsq, binExp, pheno, binExp_sc, samples, symbol2entrez): | |
| #def plotNetworks(geneSets, ids, consistencies, networks, rules_and, rules_or, rsq, binExp, pheno, symbol2entrez): # ,samples | |
| nodeColors = ['k','r','g'] | |
| pp = PdfPages('gbmNetMotifs_3node_scRNA_seq.pdf') | |
| outFile = open('pairwise_comparisons.csv','w') | |
| states = ['000', '100','010','001','110','011','101','111'] | |
| for set1 in range(len(geneSets)): | |
| #for set1 in range(len(rGeneList)): | |
| fig = plt.figure(figsize=(40,34)) | |
| plt.rcParams['legend.fontsize'] = 8 | |
| grid = plt.GridSpec(18,16, wspace=0.35, hspace=0.60, 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 ) | |
| #trans_and = network.TransGraph( 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 ) | |
| #trans_or = network.TransGraph( 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}) | |
| # And attractors [0,2] | |
| print attractors_and | |
| plt.subplot(grid[0,1]) | |
| 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) #USE MEEEEEE | |
| 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,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_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}) | |
| # Positioning for the juicy stuff | |
| counts = 1 # Initiate count for location PDF | |
| count = 0 | |
| # Filters out GSEs if they pass the rsq value | |
| rFilter = {} | |
| for gse1 in series: | |
| rFilter[gse1] = [] | |
| for netMotif in data[gse1]: | |
| for instance in data[gse1][netMotif]: | |
| if len([i for i in subsets[gse1] if consistent[gse1][netMotif][instance][i]=='Yes']) > 0: | |
| rFilter[gse1] += [instance.split(';')] | |
| # The Juicy Stuff | |
| for gse1 in series: | |
| # Verify that all the data is present for the treatment line plots | |
| if binExp[gse1].loc[int(symbol2entrez[geneSets[set1][0]])].isna().sum()==0 and binExp[gse1].loc[int(symbol2entrez[geneSets[set1][1]])].isna().sum()==0 and binExp[gse1].loc[int(symbol2entrez[geneSets[set1][2]])].isna().sum()==0: | |
| if geneSets[0] in rFilter[gse1]: | |
| rDFrame = pd.DataFrame(columns = ['treatments','genes','r2value']) | |
| rTryDict = {} | |
| rTreat = [] | |
| rGene = [] | |
| rRvalue = [] | |
| for treatments in rDict[gse1]: | |
| rTryDict[treatments] = {} | |
| for geneTriad in rDict[gse1][treatments]: | |
| if geneTriad == rGeneList[set1]: | |
| for gene in rDict[gse1][treatments][geneTriad]: | |
| rTryDict[treatments][gene] = rDict[gse1][treatments][geneTriad][gene] | |
| rGene.append(gene) | |
| rTreat.append(treatments) | |
| for rvalue in rDict[gse1][treatments][geneTriad][gene]: | |
| if rvalue =='Inconsistent': | |
| rDict[gse1][treatments][geneTriad][gene] = -0.25 | |
| rTryDict[treatments][gene] = -0.25 | |
| elif not rvalue == 'NA': | |
| rDict[gse1][treatments][geneTriad][gene] = float(rvalue) | |
| rTryDict[treatments][gene] = float(rvalue) | |
| else: | |
| rDict[gse1][treatments][geneTriad][gene] = 0 | |
| rTryDict[treatments][gene] = 0 | |
| rRvalue.append(rTryDict[treatments][gene]) | |
| rDFrame = pd.DataFrame({'treatments': rTreat, 'genes': rGene, 'r2value': rRvalue}) | |
| # Plot R squared values [0,1] | |
| plt.subplot(grid[counts, count]) | |
| plt.title(gse1 + ' r2 Values', ha = 'center') | |
| bp = sns.barplot(data = rDFrame, x = 'treatments', y = 'r2value', hue = 'genes', palette = dict(zip(list(G.nodes()),node_colors))) | |
| bp.legend_.remove() | |
| count += 1 | |
| counts += 0 | |
| #Seaborn plots for treatment over days | |
| #plt.subplot(grid[1,1]) | |
| for treatment in series[gse1]: | |
| if not treatment == 'all': | |
| dTmp = pd.DataFrame(columns=['day','TF','gexp']) | |
| for gene1 in geneSets[set1]: | |
| controls = pheno['Control'].loc[series[gse1][treatment]] | |
| pheno1 = pheno['Day'].loc[series[gse1][treatment]] | |
| gexp1 = geneExp[gse1][series[gse1][treatment]].loc[int(symbol2entrez[gene1])] | |
| avgZero = np.mean(gexp1[controls==1]) | |
| stdGexp1 = [i - avgZero for i in gexp1] | |
| dTmp_2 = pd.DataFrame({'gexp': stdGexp1, | |
| 'day': pheno1, | |
| 'TF': [gene1]*len(series[gse1][treatment])}) | |
| dTmp_3 = pd.DataFrame({'gexp': (min(stdGexp1) -0.25), | |
| 'day': pheno1, | |
| 'treatment': treatment}) | |
| dTmp = pd.concat([dTmp,dTmp_2], 0, sort=True) | |
| plt.subplot(grid[counts, count]) | |
| count += 1 | |
| plt.title(treatment,fontdict={'fontsize':8}) | |
| sns.lineplot(x='day',y='gexp', hue='TF', legend = False, palette = dict(zip(list(G.nodes()),node_colors)), err_style='bars', ci=95, data=dTmp) # pallete = 'colorblind' | |
| #sns.lineplot(x='day', y='gexp', hue='treatment', ci = 95, data = dTmp_3) | |
| counts += 1 | |
| # Build the main heatmap with states v GSE | |
| cName = [] | |
| for gsms in binExp[gse1]: | |
| cName.append(gsms) | |
| cName.append('attractorsAnd') | |
| cName.append('attractorsOr') | |
| statesHM = ['000', '100','010','001','110','011','101','111'] | |
| # statesHM = ['000','001','010','100','011','101','110','111'] | |
| hmdt = pd.DataFrame(0, index= [i for i in statesHM], columns=cName) | |
| # Filling the main heatmap matrix | |
| for gsm in binExp[gse1]: | |
| state1 = '' | |
| for i in geneSets[set1]: | |
| state1 += str(int(binExp[gse1][gsm].loc[int(symbol2entrez[i])])) | |
| #state1 = ''.join([str(int(binExp[gse1][gsm].loc[int(symbol2entrez[i])])) for i in geneSets[set1]]) | |
| hmdt.loc[state1,gsm] = 1 | |
| orAttr = '' | |
| for attrO in list(att_or_graph.nodes_with_selfloops()): | |
| for ind in hmdt.index: | |
| if ind == attrO: | |
| hmdt.loc[ind,'attractorsOr'] = 1 | |
| for attrA in list(att_and_graph.nodes_with_selfloops()): | |
| for ind in hmdt.index: | |
| if ind == attrA: | |
| hmdt.loc[ind, 'attractorsAnd'] = 1 | |
| count = 0 | |
| # Plot the heatmap with GSE title | |
| plt.subplot(grid[counts,count]) | |
| plt.title('General HeatMap') | |
| g = sns.heatmap(hmdt, cmap = 'BuPu', cbar = False, annot = False, xticklabels = False, yticklabels=True) | |
| g.set_yticklabels(g.get_yticklabels(), rotation=0,fontsize=8) | |
| count += 1 | |
| # Building the heatmap for each attractor possibility | |
| treatment = series[gse1].keys()[0] | |
| orgDay = pheno['Day'].loc[series[gse1][treatment]] | |
| htmp_and = pd.DataFrame(0, index = [i for i in series[gse1] if not i == 'all'], columns = orgDay.sort_values()) | |
| htmp_or = pd.DataFrame(0, index = [i for i in series[gse1] if not i == 'all'], columns = orgDay.sort_values()) | |
| # Creates the dictionary used to determine the color of the states and/or | |
| atAnd = list(att_and_graph.nodes_with_selfloops()) | |
| atOr = list(att_or_graph.nodes_with_selfloops()) | |
| atAndDict = {} | |
| atOrDict = {} | |
| counting = 1 | |
| for ant in atAnd: | |
| atAndDict[ant] = counting | |
| counting += 1 | |
| counting = 1 | |
| for ort in atOr: | |
| atOrDict[ort] = counting | |
| counting += 1 | |
| # Fills the dataframe for coloring and/or heatmaps | |
| cMA = {} # Used for the colormap labels | |
| cMO = {} # Used for the colormap labels | |
| for treat1 in series[gse1]: | |
| if not treat1 == 'all': | |
| orgDay = pheno['Day'].loc[series[gse1][treat1]] | |
| for day1 in range(len(orgDay)): | |
| state2 = hmdt.index[hmdt[orgDay.index[day1]]==1][0] | |
| if hmdt['attractorsAnd'][state2]==1: | |
| for state3 in atAndDict: | |
| if state3 == state2: | |
| cMA[atAndDict[state3]] = state3 # cMap label | |
| #print state3 | |
| htmp_and.loc[treat1].iloc[day1] = atAndDict[state3] | |
| if hmdt['attractorsOr'][state2]==1: | |
| for state3 in atOrDict: | |
| if state3 == state2: | |
| cMO[atOrDict[state3]] = state3 # cMap label | |
| #print state3 | |
| htmp_or.loc[treat1].iloc[day1] = atOrDict[state3] | |
| # Percent of the attractors in each GSE set for GSMs | |
| # Or | |
| percentOrs = {} | |
| for treats in htmp_or.index: | |
| orCount = 0 | |
| orTotal = 0 | |
| orTemp = 0 | |
| for value in htmp_or.loc[treats]: | |
| orTotal +=1 | |
| if not value == 0: | |
| orCount += 1 | |
| orTemp = round((float(orCount) / float(orTotal)),3) * 100 | |
| percentOrs[treats] = orTemp | |
| # And | |
| percentAnds = {} | |
| for treats in htmp_and.index: | |
| andCount = 0 | |
| andTotal = 0 | |
| andTemp = 0 | |
| for value in htmp_and.loc[treats]: | |
| andTotal +=1 | |
| if not value == 0: | |
| andCount += 1 | |
| andTemp = round((float(andCount) / float(andTotal)),3) * 100 | |
| percentAnds[treats] = andTemp | |
| # Build the color bar label list | |
| cMA[0] = 'None' | |
| cMO[0] = 'None' | |
| colorList = ['white', 'magenta', 'red','cyan', 'blue','green' ] | |
| cMapA = colorList[0:len(cMA)] | |
| cMapO = colorList[0:len(cMO)] | |
| # Function to label the plot bars on tippy top | |
| def toplabel(treatments): | |
| '''creates a label to sit above the bars for the percents of attractors present in each treatment''' | |
| for perc in treatments: | |
| height = perc.get_height() | |
| plt.annotate('{}'.format(height) + '%', | |
| xy = (perc.get_x() + perc.get_width() / 2, height), | |
| xytext = (0,3), | |
| textcoords = 'offset points', | |
| ha = 'center', | |
| va = 'bottom') | |
| # Plot the AND attractors heatmap | |
| plt.subplot(grid[counts, count]) | |
| plt.yticks() | |
| plt.xticks([]) | |
| plt.title(gse1 + ' and attractors', ha = 'center') | |
| andPlot = sns.heatmap(htmp_and, cmap = cMapA, cbar = True, annot = False, xticklabels = False) | |
| count += 1 | |
| # Colorbar | |
| cbar = andPlot.collections[0].colorbar | |
| cbar.set_ticks([i for i in cMA]) #1,2,3,4,5 | |
| cbar.set_ticklabels([i for i in cMA.values()]) | |
| # Create And percent bar plot | |
| plt.subplot(grid[counts, count]) | |
| plt.title('And Attr / Treatment', ha = 'center') | |
| plt.xticks(range(len(percentAnds)), list(percentAnds.keys())) | |
| plt.ylim([0, 100]) | |
| andBar = plt.bar(range(len(percentAnds)), list(percentAnds.values()), align = 'center') | |
| toplabel(andBar) | |
| count += 1 | |
| # Plot the OR attractors heatmap | |
| plt.subplot(grid[counts,count]) | |
| plt.yticks() | |
| plt.xticks([]) | |
| plt.title(' or attractors', ha = 'center') | |
| orPlot = sns.heatmap(htmp_or, cmap = cMapO, cbar = True, annot = False, xticklabels = False) | |
| count += 1 | |
| # Colorbar | |
| cbar = orPlot.collections[0].colorbar | |
| cbar.set_ticks([i for i in cMO]) #1,2,3,4,5 | |
| cbar.set_ticklabels([i for i in cMO.values()]) | |
| # Create Or percent bar plot | |
| plt.subplot(grid[counts, count]) | |
| plt.title('Or Attr / Treatment', ha = 'center') | |
| plt.xticks(range(len(percentOrs)), list(percentOrs.keys())) | |
| plt.ylim([0, 100]) | |
| orBar = plt.bar(range(len(percentOrs)), list(percentOrs.values()), align = 'center') | |
| toplabel(orBar) | |
| count +=1 | |
| counts +=1 # Start again | |
| count = 0 | |
| # Plot for the legends [0,3] | |
| legendDict = dict(zip(list(G.nodes()),node_colors)) | |
| plt.subplot(grid[0,3]) | |
| legendary = [1,1,1] | |
| legendBar = plt.bar(range(len(nodes)), legendary, color = ['k','r','g']) | |
| # sidelabel(legendBar) | |
| plt.xticks(range(len(nodes)), nodes) | |
| plt.yticks(range(len(nodes)), []) | |
| # for s in legendBar.patches: | |
| # legendBar.text(s.get_width() + 0.1, s.get_y()+ 0.31, str(round((s.get_width()), 2)), fontsize = 15, color = 'dimgrey') | |
| #break | |
| #break | |
| #kiiiiiiioy6 | |
| #if set1==6: | |
| # plt.show() | |
| #break | |
| pp.savefig(fig) | |
| #plt.show() | |
| #break | |
| pp.close() | |
| outFile.close() | |
| ''' | |
| 1. Treatment lines plt.arrow() | |
| 4. fix labels across tabels | |
| ''' | |
| # # Match the Or state, followed by And state | |
| # for treat3 in htmp_and.index: | |
| # for states in hmdt.index: | |
| # if hmdt.loc[states, 'attractorsOr'] == 1: | |
| # orT += sum(hmdt.loc[states,hmdt.columns[0:hmdt.shape[1]-2]]) | |
| # elif hmdt.loc[states, 'attractorsAnd'] == 1: | |
| # andT += sum(hmdt.loc[states,hmdt.columns[0:hmdt.shape[1]-2]]) | |
| # | |
| # # GSMs in the hmdt minus the last two | |
| # gsmCounter = 0 | |
| # for gsm in hmdt.columns[0:hmdt.shape[1]-2]: | |
| # gsmCounter += 1 | |
| # | |
| # # Percentage calculation | |
| # percentA = totalA / gsmCounter *100 | |
| # percentO = totalO / gsmCounter *100 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment