pathway_analysis_score_nodes.py

# import necessary modules
import argparse as argparse
import operator
import networkx as nx
import pickle
from ctypes import *
from simulation import paramClass, modelClass, NPsync
from utils import genInitValueList, setupEmptyKOKI, writeModel
from GA import GAsearchModel, localSearch

# calculate importance scores
def calcImportance(individual,params,model, sss,knockoutLists, knockinLists, boolC):
	importanceScores=[] # holder for impact scores
	# loop over all nodes
	for node in range(len(model.nodeList)):
		SSEs=[] # add SSE across samples
		nodeValues=[sss[j][model.nodeList[node]]for j in range(0,len(sss))]
		for j in range(0,len(sss)): # knock each node out and in to observe differences
			ss=sss[j]
			initValues=list(model.initValueList[j])
			knockerOuter=list(knockoutLists[j])
			knockerOuter.append(node) # knock out node
			boolValues1=NPsync(individual, model, params.cells, initValues, params, knockerOuter, knockinLists[j], boolC)
			knockerInner=list(knockinLists[j])
			knockerInner.append(node) # knock in node
			boolValues2=NPsync(individual, model, params.cells, initValues, params, knockoutLists[j], knockerInner, boolC)
			# find difference between knockout and knockin
			SSE=0
			for i in range(0, len(model.nodeList)):
				SSE+=(boolValues1[i]-boolValues2[i])**2
			SSEs.append(SSE)
		importanceScores.append(sum(SSEs))
	return importanceScores
if __name__ == '__main__':
	import time
	start_time = time.time()
	
	# read in arguments from shell scripts
	parser = argparse.ArgumentParser()
	parser.add_argument("graph")
	parser.add_argument("iterNum")
	results = parser.parse_args()
	graphName=results.graph
	iterNum=int(results.iterNum)
	name=graphName[:-8]+'_'+results.iterNum
	graph = nx.read_gpickle(graphName)
	
	# read in C function to run simulations
	updateBooler=cdll.LoadLibrary('./simulator.so')
	boolC=updateBooler.syncBool 

	# load data
	sampleList=pickle.Unpickler(open( graphName[:-8]+'_sss.pickle', "rb" )).load()
	
	# set up parameters of run, model
	params=paramClass()
	model=modelClass(graph,sampleList, False)
	model.updateCpointers()

	storeModel=[(model.size), list(model.nodeList), list(model.individualParse), list(model.andNodeList) , list(model.andNodeInvertList), list(model.andLenList),	list(model.nodeList), dict(model.nodeDict), list(model.initValueList)]
	
	# put lack of KOs, initial values into correct format
	knockoutLists, knockinLists= setupEmptyKOKI(len(sampleList))
	newInitValueList=genInitValueList(sampleList,model)
	model.initValueList=newInitValueList

	# find rules by doing GA then local search
	model1, dev, bruteOut =GAsearchModel(model, sampleList, params, knockoutLists, knockinLists, name, boolC) # run GA
	bruteOut1, equivalents, dev2 = localSearch(model1, bruteOut, sampleList, params, knockoutLists, knockinLists, boolC) # run local search
	
	# output results
	storeModel3=[(model.size), list(model.nodeList), list(model.individualParse), list(model.andNodeList) , list(model.andNodeInvertList), list(model.andLenList),	list(model.nodeList), dict(model.nodeDict), list(model.initValueList)]
	outputList=[bruteOut1,dev,storeModel, storeModel3, equivalents, dev2]
	pickle.dump( outputList, open( name+"_local1.pickle", "wb" ) ) # output rules

	# calculate importance scores and output
	scores1=calcImportance(bruteOut1,params,model1, sampleList,knockoutLists, knockinLists, boolC)
	pickle.dump( scores1, open( name+"_scores1.pickle", "wb" ) )

	# write rules
	with open(name+"_rules.txt", "w") as text_file:
		text_file.write(writeModel(bruteOut1, model1))
	print("--- %s seconds ---" % (time.time() - start_time))