Skip to content

Instantly share code, notes, and snippets.

@cplaisier

cplaisier/consistilatorV2.py

Created June 5, 2018 18:56
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save cplaisier/098420499ce7643fe6df33c908d7dca3 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save cplaisier/098420499ce7643fe6df33c908d7dca3 to your computer and use it in GitHub Desktop.
Download ZIP
Code to do consistency checking and plotting of network motifs.
Raw
consistilatorV2.py
##########################################################
## Consistilator: consistilatorV2.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 *
from multiprocessing import Pool, cpu_count, Manager
from scipy.stats import pearsonr
import numpy as np
import pandas as pd
import statsmodels.api as sm
# Get a correlation p-value from R
def correlation(a1, a2):
"""
Calculate the correlation coefficient and p-value between two variables.
Input: Two arrays of float or integers.
Returns: Corrleation coefficient and p-value.
"""
"""
# Fire up R
rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
runMe = []
# Make the data into an R matrix
runMe.append('c1 = cor.test(c('+','.join([str(i) for i in a1])+'),c('+','.join([str(i) for i in a2])+'))')
runMe.append('c1$estimate')
runMe.append('c1$p.value')
runMe = '\n'.join(runMe)+'\n'
out = rProc.communicate(runMe)
# Process output
splitUp = out[0].strip().split('\n')
rho = float(splitUp[1])
pValue = float((splitUp[2].split(' '))[1])
"""
r1 = pearsonr(a1,a2)
return [r1[0], r1[1]]
# Send model for comparison versus
def regression_R(response, predictors):
"""
Fit a linear regression model of all terms.
Input: A response variable and a specified number of predictor variables.
- response = vector of floats
- predictors = dictionary of vectors of floats hashed on variable names
Returns: Overall model significance and term significance.
"""
# Fire up R
rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
runMe = []
# Make the data into an R matrix
tmp = []
names = {}
revNames = {}
for i in predictors:
names[i] = i.replace('-','_')
revNames[i.replace('-','_')] = i
tmp.append(names[i]+' = c('+','.join([str(i) for i in predictors[i]])+')')
runMe.append('d1 = data.frame(response = c('+','.join([str(i) for i in response])+')'+','+','.join(tmp)+')')
runMe.append('slm1 = summary(lm(response ~ .,data=d1))')
runMe.append('slm1$coefficients')
runMe.append('slm1$adj.r.squared')
runMe = '\n'.join(runMe)+'\n'
out = rProc.communicate(runMe)
splitUp = out[0].strip().split('\n')
adj_r_squared = splitUp.pop(-1).split(' ')[1]
splitUp.pop(0)
splitUp.pop(0)
splitUp = [[j for j in i.split(' ') if j] for i in splitUp]
tmp = dict(zip([revNames[i[0]] for i in splitUp],[{'estimate':i[1],'stdErr':i[2],'t':i[3],'p_value':i[4]} for i in splitUp]))
tmp['adj_r_squared'] = adj_r_squared
return tmp
# Send model for comparison versus
def regression(response, predictors):
"""
Fit a linear regression model of all terms.
Input: A response variable and a specified number of predictor variables.
- response = a Pandas dataframe of the response
- predictors = a Pandas dataframe of predictors
Returns: Overall model significance and term significance.
"""
# Run linear regression in Python
tmp = {}
model = sm.OLS(response, predictors).fit()
for p1 in list(predictors.columns.values):
tmp[p1] = {'estimate':model.params[p1],'stdErr':model.bse[p1],'t':model.tvalues[p1],'p_value':model.pvalues[p1]}
tmp['adj_r_squared'] = model.rsquared_adj
return tmp
### Load sample subsets ###
subsets = { 'all':range(0,12) }
### Load sub-networks ###
# Load up FANMOD subgraph enumeration results
inFile = open('mdraw_CHIR_curve_mdraw.txt.OUT','r')
# Get rid of headers
while 1:
line = inFile.readline()
if line.strip()=='Result overview:':
inFile.readline()
inFile.readline()
inFile.readline()
inFile.readline()
break
subnetworks = []
filter_pv = 0.05
filter_z = 2
while 1:
line = inFile.readline()
if not line:
break
line2 = inFile.readline().strip() # Get rid of second adjacency matrix line
line3 = inFile.readline().strip() # Get rid of third adjacency matrix line
#line4 = inFile.readline().strip() # Get rid of fourth adjacency matrix line
#print line4
inFile.readline() # Get rid of blank line
splitUp = [i for i in line.strip().split(' ') if i]
id = splitUp[1]+line2+line3 #+line4
if float(splitUp[6]) <= filter_pv and float(splitUp[5]) >= filter_z:
subnetworks.append(id)
inFile.close()
# Take out extrodinarily large subnetworks
subnetworks = [i for i in subnetworks if not i in [74, 2184]]
# Create dictionary to convert
id2gene = {}
inFile = open('mdraw_CHIR_curve_mdraw_LIST.txt','r')
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(' ')
id2gene[splitUp[1]] = splitUp[0]
inFile.close()
# Load up network motifs
networkMotifs = {}
inFile = open('mdraw_CHIR_curve_mdraw.txt.OUT.dump','r')
inFile.readline() # Get rid of header
inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(',')
#numMotif = int(splitUp.pop(0).replace('2','1'), 2)
motif = splitUp.pop(0)
if motif in subnetworks:
if not motif in networkMotifs:
networkMotifs[motif] = []
networkMotifs[motif].append([id2gene[i] for i in splitUp])
inFile.close()
# Read in Biotapesty file
outEdges = {}
inEdges = {}
inFile = open('biotapestry_CHIR_curve.csv','r')
counts = 0
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('"')
if not node1 in outEdges:
outEdges[node1] = []
if not node2 in outEdges[node1]:
outEdges[node1].append(node2)
if not node2 in inEdges:
inEdges[node2] = []
if not node1 in inEdges[node2]:
inEdges[node2].append(node1)
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
"""expression = {}
inFile = open('allExpression.csv','r')
inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(',')
gene = splitUp.pop(0)
splitUp = splitUp[2:]
expression[gene] = splitUp
inFile.close()
"""
expression = pd.read_csv('../../chir_gradient_EntrezID.csv', header=0, index_col=0).transpose()
# To test the survival function
#print regression(expression['Cluster-10'], {'SP8':expression['SP8'], 'SIM2':expression['SIM2'],'NKX2-1':expression['NKX2-1']})
def consistentInstance(instance, inEdges_sm, expression_sm, writeMe, cur, symbol2entrez):
print ' '+' '.join(instance)+' - '+cur['subset']+' - '+str(cur['subnetwork'])
entry = {}
for node in instance:
# A. Determine inputs for each node
inputs = []
if node in inEdges_sm:
inputs = [i for i in list(set(inEdges_sm[node]).intersection(instance)) if not i==node]
# B. Test model consistency for inputs to each node
if len(inputs)>0:
#res1 = regression(expression_sm[node],dict(zip(inputs,[expression_sm[i] for i in inputs])))
res1 = regression(expression_sm[int(symbol2entrez[node])],expression_sm[[int(symbol2entrez[i]) for i in inputs]])
consistent = 'Inconsistent'
if len([i for i in inputs if float(res1[int(symbol2entrez[i])]['p_value'])<=0.05])==len(inputs):
consistent = res1['adj_r_squared']
entry[node] = {'inputs':inputs,'consistency':consistent}
else:
entry[node] = {'inputs':inputs,'consistency':'NA'}
writeMe.append(cur['subset']+','+str(int(cur['subnetwork'].replace('2','1'), 2))+','+str(cur['subnetwork'])+','+';'.join(instance)+','+','.join([i+','+';'.join(entry[i]['inputs'])+','+str(entry[i]['consistency']) for i in instance]))
# Make shared memory objects
#cpus = cpu_count()
#mgr = Manager()
#inEdges_sm = mgr.dict(inEdges)
inEdges_sm = inEdges
#expression_sm = mgr.dict(expression)
expression_sm = expression
### For each subset of samples including all ###
# Goal is to write out a file with this header:
# Data Subset,Netowrk Motif,Instance,Node1,Node1.Inputs,Node1.Consistency,Node2,Node2.Inputs,Node2.Consistency,Node3,Node3.Inputs,Node3.Consistency,Node4,Node4.Inputs,Node4.Consistency
for subset in subsets:
#for subset in ['all']:
writeMe = [] #mgr.list()
writeMe.append('Data Subset,Netowrk Motif,Adjacency Matrix,Instance,Node1,Node1.Inputs,Node1.Consistency,Node2,Node2.Inputs,Node2.Consistency,Node3,Node3.Inputs,Node3.Consistency,Node4,Node4.Inputs,Node4.Consistency')
print 'Working on '+subset+' data subset.'
### For each network motif enriched in network ###
for subnetwork in subnetworks:
print ' Working on '+str(subnetwork)+' subnetwork.'
### For each subnetwork instance ###
for instance in networkMotifs[subnetwork]:
cur = {'subset':subset, 'subnetwork':subnetwork}
consistentInstance(instance, inEdges_sm, expression_sm, writeMe, cur, symbol2entrez)
"""cur = mgr.dict()
cur['subset'] = subset
cur['subnetwork'] = subnetwork
pool = Pool(processes=cpus)
res1 = pool.map(consistentInstance,networkMotifs[subnetwork])
pool.close()
pool.join()
"""
"""
for instance in networkMotifs[subnetwork]:
print ' '+' '.join(instance)+' - '+subset+' - '+str(subnetwork)
entry = {}
for node in instance:
# A. Determine inputs for each node
inputs = []
if node in inEdges:
inputs = [i for i in list(set(inEdges[node]).intersection(instance)) if not i==node]
# B. Test model consistency for inputs to each node
if len(inputs)>0:
res1 = regression(expression[node],dict(zip(inputs,[expression[i] for i in inputs])))
consistent = 'Inconsistent'
if len([i for i in inputs if float(res1[i]['p_value'])<=0.05])==len(inputs):
consistent = res1['adj_r_squared']
entry[node] = {'inputs':inputs,'consistency':consistent}
else:
entry[node] = {'inputs':inputs,'consistency':'NA'}
writeMe.append(subset+','+str(subnetwork)+','+';'.join(instance)+','+','.join([i+','+';'.join(entry[i]['inputs'])+','+str(entry[i]['consistency']) for i in instance]))
"""
outFile = open('results_'+subset+'.csv','w')
outFile.write('\n'.join(writeMe))
outFile.close()
Raw
plotNetworkMotifs.py
#################################################################
# @Program: plotNetworkMotifs.py #
# @Version: 1 #
# @Author: Chris Plaisier #
# @Sponsored by: #
# Nitin Baliga, ISB #
# Institute for Systems Biology #
# 1441 North 34th Street #
# Seattle, Washington 98103-8904 #
# (216) 732-2139 #
# #
# If this program is used in your analysis please mention who #
# built it. Thanks. :-) #
# #
# Copyrighted by Chris Plaisier 12/8/2014 #
#################################################################
from subprocess import *
#from multiprocessing import Pool, cpu_count, Manager
# Make descriptive plot for genes
#def plotGenes(geneSets, probeSets, ids, netMatrices, consistencies, networks, rsq):
def plotGenes(geneSets, ids, netMatrices, consistencies, networks, rsq):
runMe = []
runMe += ['pdf(\'networkMotifPlots_3_R_'+str(rsq)+'.pdf\',width=10.5,height=8)']
runMe += ['library(network)']
runMe += ['library(Hmisc)']
runMe += ['library(org.Hs.eg.db)']
runMe += ['x = org.Hs.egSYMBOL2EG',
'mapped_genes = mappedkeys(x)',
'xx = as.list(x[mapped_genes])']
# Load up data for plots
runMe += [
# Load up miRNA expression data
'dTF = read.csv(\'../../CHIR_curve_geneExp_all.csv\',header=T, row.names=1)',
#'dGSE45223 = read.csv(\'C:/Users/plais/Dropbox (ASU)/StemCellDifferentiation/GSE45223/GSE45223_genesExpMatrix_all.csv\',header=1,row.names=1)',
#'dGSE32658 = read.csv(\'C:/Users/plais/Dropbox (ASU)/StemCellDifferentiation/GSE32658/GSE32658_genesExpMatrix_all.csv\',header=1,row.names=1)',
#'dGSE26867 = read.csv(\'C:/Users/plais/Dropbox (ASU)/StemCellDifferentiation/GSE26867/GSE26867_genesExpMatrix_all.csv\',header=1,row.names=1)',
#'dGSE45223 = read.csv(\'C:/Users/cplaisie/Dropbox (ASU)/StemCellDifferentiation/GSE45223/GSE45223_genesExpMatrix_all.csv\',header=1,row.names=1)',
#'dGSE32658 = read.csv(\'C:/Users/cplaisie/Dropbox (ASU)/StemCellDifferentiation/GSE32658/GSE32658_genesExpMatrix_all.csv\',header=1,row.names=1)',
#'dGSE26867 = read.csv(\'C:/Users/cplaisie/Dropbox (ASU)/StemCellDifferentiation/GSE26867/GSE26867_genesExpMatrix_all.csv\',header=1,row.names=1)',
'dGSE45223 = read.csv(\'~/Dropbox (ASU)/StemCellDifferentiation/GSE45223/GSE45223_genesExpMatrix_all.csv\',header=1,row.names=1)',
'dGSE32658 = read.csv(\'~/Dropbox (ASU)/StemCellDifferentiation/GSE32658/GSE32658_genesExpMatrix_all.csv\',header=1,row.names=1)',
'dGSE26867 = read.csv(\'~/Dropbox (ASU)/StemCellDifferentiation/GSE26867/GSE26867_genesExpMatrix_all.csv\',header=1,row.names=1)',
'']
# Iterate through all geneSets
for i in range(0,len(geneSets)):
genes = geneSets[i]
print genes
#probes = probeSets[i]
id = ids[i]
netMatrix = netMatrices[i]
consistency = consistencies[i]
network = networks[i]
print genes,netMatrix,id
# Select probes for genes
runMe += [
'tmp = as.character(sapply(c('+','.join(['\''+gene+'\'' for gene in genes])+'), function(tf1) { xx[[tf1]] }))',
'tfsExp = as.matrix(dTF[tmp,])']
# Transform data for plotting
runMe += [
#'tfsExp_sw = sweep(tfsExp, 1, tfsExp[,1])',
'tfsExp_sw = tfsExp',
'ylim1 = c(floor(min(tfsExp_sw)),ceiling(max(tfsExp_sw)))'
]
runMe += [
# Plot CHIR curve (Brafman et al.)
'layout(rbind(c(1,1,5),c(1,1,6),c(2,3,4)))',
#'par(mar=c(2,2,2,1))',
'plot(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[1,],type=\'l\',col=1,lwd=2,lty=1,main=\'CHIR Curve\',ylim=ylim1,xaxt=\'n\',ylab=\'TF Expression\',xlab=\'[CHIR]\')',
'points(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[1,],col=1,pch=19,cex=1.5)',
'lines(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),rep(tfsExp_sw[1,1],12),col=1,lwd=1,lty=2)',
'lines(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[2,],col=2,lwd=2,lty=1)',
'lines(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),rep(tfsExp_sw[2,1],12),col=2,lwd=1,lty=2)',
'points(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[2,],col=2,pch=19,cex=1.5)',
'lines(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[3,],col=3,lwd=2,lty=1)',
'lines(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),rep(tfsExp_sw[3,1],12),col=3,lwd=1,lty=2)',
'points(c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),tfsExp_sw[3,],col=3,pch=19,cex=1.5)',
'abline(h=0,lwd=1,col=\'black\')',
'axis(1,at=c(0,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1))',
# GSE45223 nc
'sets = list(c(1:3),c(4:6),c(7:9),c(10:12),c(13:15),c(16:18))',
'tfsExp_nc = as.matrix(dGSE45223[tmp,c(1,13,25,2,14,26,4,16,28,6,18,30,8,20,32,11,23,35)])',
'ylim1 = c(min(tfsExp_nc,na.rm=T),max(tfsExp_nc,na.rm=T))',
'tfsExp_nc_med = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'tfsExp_nc_sd = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'for(i in 1:length(sets)) {',
' for(j in 1:nrow(tfsExp_nc)) {',
#' cat(sets[[i]])',
' tfsExp_nc_med[j,i] = median(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' tfsExp_nc_sd[j,i] = sd(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' }',
'}',
'time1 = c(0, 1, 3, 6, 8, 11)',
'plot(time1,tfsExp_nc_med[1,],type=\'l\',col=1,lwd=2,lty=1,main=\'Neural Crest (NC): GSE45223\',ylim=ylim1,xaxt=\'n\',ylab=\'TF Expression\',xlab=\'Days\')',
'lines(time1,rep(tfsExp_nc_med[1,1],length(time1)),type=\'l\',col=1,lwd=1,lty=2,ylim=ylim1)',
'axis(1,at=time1)',
'lines(c(2,11),rep(mean(ylim1),2),col=rgb(1,0,0,0.10),lwd=10000,lend=1)',
'errbar(time1,tfsExp_nc_med[1,], tfsExp_nc_med[1,]+tfsExp_nc_sd[1,], tfsExp_nc_med[1,]-tfsExp_nc_sd[1,],col=1,errbar.col=1,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[2,],type=\'l\',col=2,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[2,1],length(time1)),type=\'l\',col=2,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[2,], tfsExp_nc_med[2,]+tfsExp_nc_sd[2,], tfsExp_nc_med[2,]-tfsExp_nc_sd[2,],col=2,errbar.col=2,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[3,],type=\'l\',col=3,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[3,1],length(time1)),type=\'l\',col=3,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[3,], tfsExp_nc_med[3,]+tfsExp_nc_sd[3,], tfsExp_nc_med[3,]-tfsExp_nc_sd[3,],col=3,errbar.col=3,ylim=ylim1,add=T)',
'abline(h=0,lwd=1,col=\'black\')',
# GSE32658 lsfc
#'par(mar=c(2,2,2,1))',
'sets = list(c(1:3),c(4:6),c(7:9),c(10:12),c(13:15),c(16:18),c(19:25))',
'tfsExp_nc = as.matrix(dGSE32658[tmp,c(1,2,3,4,5,6,10,11,12,19,20,21,28,29,30,37,38,39,46,47,48,57,58,59,60)])',
'ylim1 = c(min(tfsExp_nc,na.rm=T),max(tfsExp_nc,na.rm=T))',
'tfsExp_nc_med = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'tfsExp_nc_sd = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'for(i in 1:length(sets)) {',
' for(j in 1:nrow(tfsExp_nc)) {',
#' cat(sets[[i]])',
' tfsExp_nc_med[j,i] = median(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' tfsExp_nc_sd[j,i] = sd(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' }',
'}',
'time1 = c(0,1,3,5,7,11,13)',
'plot(time1,tfsExp_nc_med[1,],type=\'l\',col=1,lwd=2,lty=1,main=\'DA Midbrain (LSFC): GSE32658\',ylim=ylim1,xaxt=\'n\',ylab=\'Expression\',xlab=\'Days\')',
'lines(time1,rep(tfsExp_nc_med[1,1],length(time1)),type=\'l\',col=1,lwd=1,lty=2,ylim=ylim1)',
'axis(1,at=time1)',
'lines(c(3,11),rep(mean(ylim1),2),col=rgb(1,0,0,0.1),lwd=10000,lend=1)',
'errbar(time1,tfsExp_nc_med[1,], tfsExp_nc_med[1,]+tfsExp_nc_sd[1,], tfsExp_nc_med[1,]-tfsExp_nc_sd[1,],col=1,errbar.col=1,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[2,],type=\'l\',col=2,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[2,1],length(time1)),type=\'l\',col=2,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[2,], tfsExp_nc_med[2,]+tfsExp_nc_sd[2,], tfsExp_nc_med[2,]-tfsExp_nc_sd[2,],col=2,errbar.col=2,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[3,],type=\'l\',col=3,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[3,1],length(time1)),type=\'l\',col=3,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[3,], tfsExp_nc_med[3,]+tfsExp_nc_sd[3,], tfsExp_nc_med[3,]-tfsExp_nc_sd[3,],col=3,errbar.col=3,ylim=ylim1,add=T)',
'abline(h=0,lwd=1,col=\'black\')',
# GSE26867 lsb3i
#'par(mar=c(2,2,2,1))',
#'sets = c(c(4:6),c(7:9),c(10:12),c(13:15),c(16:18),c(19:21),c(22:24),c(25:27),c(28:30),c(31:33))',
'sets = list(c(1:3),c(4:6),c(7:9),c(9:12),c(13:15),c(15:18))',
'tfsExp_nc = as.matrix(dGSE26867[tmp,c(1:3,19:33)])',
'ylim1 = c(min(tfsExp_nc,na.rm=T),max(tfsExp_nc,na.rm=T))',
'tfsExp_nc_med = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'tfsExp_nc_sd = matrix(nrow=nrow(tfsExp_nc),ncol=length(sets))',
'for(i in 1:length(sets)) {',
' for(j in 1:nrow(tfsExp_nc)) {',
#' cat(sets[[i]])',
' tfsExp_nc_med[j,i] = median(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' tfsExp_nc_sd[j,i] = sd(tfsExp_nc[j,sets[[i]]],na.rm=T)',
' }',
'}',
'time1 = c(2,3,5,7,9,15)',
'plot(time1,tfsExp_nc_med[1,],type=\'l\',col=1,lwd=2,lty=1,main=\'Noci. PNS (LSB3i): GSE26867\',ylim=ylim1,xaxt=\'n\',ylab=\'Expression\',xlab=\'Days\')',
'lines(time1,rep(tfsExp_nc_med[1,1],length(time1)),type=\'l\',col=1,lwd=1,lty=2,ylim=ylim1)',
'axis(1,at=time1)',
'lines(c(2,10),rep(mean(ylim1),2),col=rgb(1,0,0,0.1),lwd=10000,lend=1)',
'errbar(time1,tfsExp_nc_med[1,], tfsExp_nc_med[1,]+tfsExp_nc_sd[1,], tfsExp_nc_med[1,]-tfsExp_nc_sd[1,],col=1,errbar.col=1,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[2,],type=\'l\',col=2,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[2,1],length(time1)),type=\'l\',col=2,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[2,], tfsExp_nc_med[2,]+tfsExp_nc_sd[2,], tfsExp_nc_med[2,]-tfsExp_nc_sd[2,],col=2,errbar.col=2,ylim=ylim1,add=T)',
'lines(time1,tfsExp_nc_med[3,],type=\'l\',col=3,lwd=2,lty=1,ylim=ylim1)',
'lines(time1,rep(tfsExp_nc_med[3,1],length(time1)),type=\'l\',col=3,lwd=1,lty=2,ylim=ylim1)',
'errbar(time1,tfsExp_nc_med[3,], tfsExp_nc_med[3,]+tfsExp_nc_sd[3,], tfsExp_nc_med[3,]-tfsExp_nc_sd[3,],col=3,errbar.col=3,ylim=ylim1,add=T)',
'abline(h=0,lwd=1,col=\'black\')',
# Plot network and adjusted R^2
'par(mar=c(0.1,0.1,3,0.1))',
'm = t(matrix(c('+network+'),3))',
'plot(network(m),edge.col=c(\'green\',\'red\')[m],label=c('+','.join(['\''+gene+'\'' for gene in genes])+'),vertex.col=c(1,2,3),vertex.cex=4,vertex.border=0,arrowhead.cex=3,main=\''+id+'\')', # Expects a string
'par(mar=c(3,5,2,1))',
'ylab1 = expression(paste(\'Adjusted \',R^2,sep=\'\'))',
'barplot(t(rbind(\'All\'=c('+','.join(consistency['all'])+'))),col=c(1,2,3),beside=T,ylab=ylab1,ylim=c(-0.25,1))',
'abline(h=0,lwd=1,col=\'black\')',
'abline(h='+str(rsq)+',lwd=1,col=rgb(1,0,0,0.8),lty=2)'] # Expects a list of dictionaries indexed on 'all','l','lsf','lsfc'
# Close plotting device and write out PDF
runMe += ['dev.off()']
outFile = open('test.R','w')
outFile.write('\n'.join(runMe)+'\n')
outFile.close()
# Run in R
rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
runMe = '\n'.join(runMe)+'\n'
out = rProc.communicate(runMe)
print out
# Read in Biotapesty file
inEdges = {}
inFile = open('biotapestry_CHIR_curve.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()
# Gather data
geneSets = []
#probeSets = []
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)
"""gene2probe = {}
tmp = data[netMotif][instance]['all'][4:]
while 1:
if not tmp[1]=='NA':
if not tmp[0] in gene2probe:
gene2probe[tmp[0]] = []
if not tmp[1] in gene2probe[tmp[0]]:
gene2probe[tmp[0]].append(tmp[1])
for i in tmp[3].split(' '):
if not i=='':
tmp1 = i.split(':')
if not tmp1[0] in gene2probe:
gene2probe[tmp1[0]] = []
if not tmp1[1] in gene2probe[tmp1[0]]:
gene2probe[tmp1[0]].append(tmp1[1])
if len(tmp)==5:
break
tmp = tmp[5:]
probeSets.append([gene2probe[gene][0] for gene in 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':[data[netMotif][instance]['all'][i] for i in [6,9,12]]}
for i in tmp:
for j in range(0,len(tmp[i])):
if tmp[i][j]=='Inconsistent':
tmp[i][j] = '-0.25'
consistencies.append(tmp)
tmp = dict(zip(genes,[dict(zip(genes,[0 for gene in genes])) for gene in genes]))
for gene1 in genes:
if gene1 in inEdges:
for gene2 in inEdges[gene1]:
if gene2 in genes:
if inEdges[gene1][gene2]=='positive':
tmp[gene2][gene1] = 1
elif inEdges[gene1][gene2]=='negative':
tmp[gene2][gene1] = 2
networks.append(','.join([','.join([str(tmp[gene2][gene1]) for gene1 in genes]) for gene2 in genes]))
# Plot them
#plotGenes(geneSets, probeSets, ids, netMatrices, consistencies, networks, rsq)
plotGenes(geneSets, ids, netMatrices, consistencies, networks, rsq)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment