With survival analysis and network states.

plotNetworkMotifs.py

##########################################################
## 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
import matplotlib.pyplot as plt
from scipy.stats import pearsonr
from lifelines import KaplanMeierFitter
#import copy
#from multiprocessing import Pool, cpu_count, Manager

## Plot network information
##  ___________________
##  |        |        |
##  | Expr.  |  R^2   |
##  |        |        |
##  ___________________
##  |        |        |
##  | Net.   | Attra. |
##  |        |        |
##  ___________________
##  |        |        |
##  | Surv.  | State  |
##  |        | Dist.  |
##  ___________________
#def plotNetworks(geneSets, ids, consistencies, networks, rsq, exp, binExp, pheno, exp_lgg, binExp_lgg, pheno_lgg, symbol2entrez):
def plotNetworks(geneSets, ids, consistencies, networks, rsq, exp, binExp, pheno, symbol2entrez):
    nodeColors = ['k','r','g']
    for i in range(len(geneSets)):
        nodes = dict(zip(geneSets[i],nodeColors))
        fig = plt.figure(figsize=(11,8.5))
        grid = plt.GridSpec(3, 2, wspace=0.25, hspace=0.25, left=0.075, right=0.95, bottom=0.05, top=0.95)
        # Plot Gene Expresion [0,0]
        ax = plt.subplot(grid[0,0])
        gl1 = exp.loc[[int(symbol2entrez[str(geneSets[i][0])]),int(symbol2entrez[str(geneSets[i][1])]),int(symbol2entrez[str(geneSets[i][2])])]]
        #gl1 = exp.loc[[str(geneSets[i][0]),str(geneSets[i][1]),str(geneSets[i][2])]]
        cols1 = gl1.columns[gl1.mean().argsort()]
        gl2 = gl1[cols1]        
        plt.plot(gl2.loc[int(symbol2entrez[str(geneSets[i][0])])].tolist(),nodes[geneSets[i][0]],gl2.loc[int(symbol2entrez[str(geneSets[i][1])])].tolist(),nodes[geneSets[i][1]],gl2.loc[int(symbol2entrez[str(geneSets[i][2])])].tolist(),nodes[geneSets[i][2]])
        #plt.plot(gl2.loc[str(geneSets[i][0])].tolist(),nodes[geneSets[i][0]],gl2.loc[str(geneSets[i][1])].tolist(),nodes[geneSets[i][1]],gl2.loc[str(geneSets[i][2])].tolist(),nodes[geneSets[i][2]])
        ax.set_ylabel('Relative Expression')
        ax.set_xlabel('Patients')
        plt.tight_layout()
        
        # Plot R squared values [0,1]
        ax = plt.subplot(grid[0,1])
        index = np.arange(len(geneSets[i]))
        print index
        plt.xticks(index, geneSets[i])
        g1, g2, g3 = plt.bar(index, [consistencies[i]['all'][j] for j in geneSets[i]])
        print consistencies[i]['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[geneSets[i][0]])
        g2.set_facecolor(nodes[geneSets[i][1]])
        g3.set_facecolor(nodes[geneSets[i][2]])
        plt.tight_layout()
        
        # Network plot [1,0]
        plt.subplot(grid[1,0])
        G = networks[i]
        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)

        # Attractors [1,1]

        
        # 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']))
        tmp = binExp.loc[geneSets[1]].transpose()
        #print tmp
        tmp1 = dict(zip(tmp.index,[''.join([str(k) for k in tmp.loc[j].values]) for j in tmp.index]))
        tmp2 = []
        for j in pheno.index:
            if not j in tmp1:
                tmp2.append(np.nan)
            else:
                tmp2.append(tmp1[j])
        print len(tmp2)
        pheno2 = pheno.assign(binSet=tmp2)
    
        # Survival [2,0]
        ax = plt.subplot(grid[2,0])
        kmf = KaplanMeierFitter()
        pheno3 = pheno2[['SURVIVAL','DEAD','binSet']].dropna(axis='rows')
        #pheno3 = pheno2[['survival','vital_status','binSet']].dropna(axis='rows')
        T = pheno3['SURVIVAL']
        #T = pheno3['survival']
        E = pheno3['DEAD']=='DEAD'
        #E = pheno3['vital_status']=='dead'
        groups = pheno3['binSet']
        dist1 = groups.value_counts()
        for k in states:
            if k in dist1:
                ix = (groups==k)
                kmf.fit(T[ix], E[ix], label=k)
                kmf.plot(ax=ax, ci_show=False, color=states_colors[k])
        #print dir(ax)
        #print ax._get_lines().get_next_color()
        
    
        # Survival grouped by attractor states [1,1]
        ax = plt.subplot(grid[1,1])
        kmf = KaplanMeierFitter()
        pheno3 = pheno2[['SURVIVAL','DEAD','binSet']].dropna(axis='rows')
        #pheno3 = pheno2[['survival','vital_status','binSet']].dropna(axis='rows')
        T = pheno3['SURVIVAL']
        #T = pheno3['survival']
        E = pheno3['DEAD']=='DEAD'
        #E = pheno3['vital_status']=='dead'
        groups = pheno3['binSet']
        dist1 = groups.value_counts()
        onOff_colors = {'off':'k','on':'r'}
        onOff = {'off':['010','100','000','110'], 'on':['011','101','001','111']}
        for k in onOff:
            #if k in dist1:
            ix = []
            for i in groups:
                if not i in onOff[k]:
                    ix.append(False)
                else:
                    ix.append(True)
            print ix
            kmf.fit(T[ix], E[ix], label=k)
            kmf.plot(ax=ax, ci_show=False, color=onOff_colors[k])
        #print dir(ax)
        #print ax._get_lines().get_next_color()

        # Distribution of states [2,1]
        ax = plt.subplot(grid[2,1])
        groups = pheno2['binSet']
        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)
        index = np.arange(len(states))
        plt.xticks(index, states)
        g1, g2, g3, g4, g5, g6, g7, g8 = plt.bar(index, dist2)
        g1.set_facecolor(states_colors[states[0]])
        g2.set_facecolor(states_colors[states[1]])
        g3.set_facecolor(states_colors[states[2]])
        g4.set_facecolor(states_colors[states[3]])
        g5.set_facecolor(states_colors[states[4]])
        g6.set_facecolor(states_colors[states[5]])
        g7.set_facecolor(states_colors[states[6]])
        g8.set_facecolor(states_colors[states[7]])
        plt.ylim((0,1))
        plt.ylabel('% Patients')
        plt.xlabel('States')
        #ax.axhline(rsq,color='k',linestyle='--',alpha=0.5)
        #g1.set_facecolor(nodes[geneSets[i][0]])
        #g2.set_facecolor(nodes[geneSets[i][1]])
        #g3.set_facecolor(nodes[geneSets[i][2]])
        plt.tight_layout()

        plt.show()
        break

# 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.9
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 expression data
exp = pd.read_csv('tfExp.csv', header=0, index_col=0)#.transpose()
exp_lgg = pd.read_csv('gbmlgg.csv', header=0, index_col=0)#.transpose()

# Read in binarized data
binExp = pd.read_csv('tfBin_TCGA_GBM.csv', header=0, index_col=0)#.transpose()
binExp_lgg = pd.read_csv('binGeneExp_VkMeans_lgg.csv', header=0, index_col=0)#.transpose()
binExp_lgg.columns = [i.replace('.','-') for i in binExp_lgg.columns.values]

# Phenotypes
pheno = pd.read_csv('phenotypes.csv', header=0, index_col=0)
pheno_lgg = pd.read_csv('phenotypes_lgg.csv', header=0, index_col=0)

# Gather data
geneSets = []
ids = []
netMatrices = []
consistencies = []
networks = []
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])
                #netMatrices.append([i for i in 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()
                for gene1 in genes:
                    if gene1 in inEdges:
                        for gene2 in inEdges[gene1]:
                            if gene2 in genes:
                                if inEdges[gene1][gene2]=='positive':
                                    G.add_edge(gene2,gene1,color='g')
                                elif inEdges[gene1][gene2]=='negative':
                                    G.add_edge(gene2,gene1,color='r')
                networks.append(G)

# Plot them
#plotGenes(geneSets, probeSets, ids, netMatrices, consistencies, networks, rsq)
#plotGenes(geneSets, ids, netMatrices, consistencies, networks, rsq)
#plotNetworks(geneSets, ids, consistencies, networks, rsq, exp, binExp, pheno, exp_lgg, binExp_lgg, pheno_lgg, symbol2entrez)
plotNetworks(geneSets, ids, consistencies, networks, rsq, exp, binExp, pheno, symbol2entrez)