-
Notifications
You must be signed in to change notification settings - Fork 8
/
runTE.py
executable file
·79 lines (70 loc) · 3 KB
/
runTE.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
from jpype import *
import random
import math
import numpy
import os
import csv
import sys
import datetime
# Change location of jar to match yours:
abspath = os.path.abspath(sys.argv[0])
dname = os.path.dirname(abspath)
os.chdir(dname)
jarLocation = "infodynamics.jar"
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation,"-Xmx16G")
##ifile = open(os.environ.get('trajectory_file'))
ifile = open(sys.argv[3])
trajectory1=[]
for line in ifile:
trajectory1.append(float(line.split()[0]))
ifile = open(sys.argv[4])
##ifile = open(os.environ.get('cell_select_file'))
branch=[]
for line in ifile:
branch.append(int(line.split()[0]))
branch=numpy.array(branch)
historyLength = int(sys.argv[5])
trajectory1=numpy.array(trajectory1)
trajectory1=trajectory1[branch==1]
trajectory1SortIndex=numpy.argsort(trajectory1)
cell_gene_all = numpy.genfromtxt('cell_gene_trsps.csv',delimiter=',')
#cell_gene_all = cell_gene_all.astype(int)
expression_data = []
input_fname = sys.argv[1]
#input_fname='pair_jobs/pair_list_aaa'
list_pairs = numpy.genfromtxt(input_fname,delimiter=',')
list_pairs = list_pairs.astype(int)
#cg_file = open('list_genefiles')
#lines = cg_file.readlines()
TEresult=[None] * len(list_pairs)
for num_pair in range(len(list_pairs)):
expression_data = cell_gene_all[list_pairs[num_pair,0]-1][numpy.newaxis]
expression_data = numpy.append(expression_data, cell_gene_all[list_pairs[num_pair,1]-1][numpy.newaxis],axis=0)
expression_data1=[]
for i in range(len(expression_data)):
data_temp1 = numpy.array(expression_data[i])
data_temp1 = data_temp1[branch==1]
data_temp1 = data_temp1[trajectory1SortIndex]
data_temp1 = list(data_temp1)
expression_data1.append(data_temp1)
expression_data = expression_data1
del expression_data1
# Create a TE calculator and run it:
teCalcClass = JPackage("infodynamics.measures.continuous.kernel").TransferEntropyCalculatorKernel
teCalc = teCalcClass()
teCalc.setProperty("NORMALISE", "true") # Normalise the individual variables
teCalc.initialise(historyLength, 0.5) # Use history length 1 (Schreiber k=1), kernel width of 0.5 normalised units
resultTemp=[]
teCalc.setObservations(JArray(JDouble, 1)(expression_data[0]), JArray(JDouble, 1)(expression_data[1]))
resultTemp.append(teCalc.computeAverageLocalOfObservations())
teCalc.initialise(historyLength, 0.5) # Use history length 1 (Schreiber k=1), kernel width of 0.5 normalised units
teCalc.setObservations(JArray(JDouble, 1)(expression_data[1]), JArray(JDouble, 1)(expression_data[0]))
resultTemp.append(teCalc.computeAverageLocalOfObservations())
TEresult[num_pair] = numpy.ndarray.tolist(list_pairs[num_pair,:]) + resultTemp
if (num_pair % int(len(list_pairs)/2)) == 0:
print(datetime.datetime.now())
output_fname = sys.argv[2]
with open(output_fname, 'w') as f:
writer = csv.writer(f)
writer.writerows(TEresult)