# Usage: python local_network.py
#
# Put this file under directory Code
# Prepare a paramter_for_buildCmatrix.txt and put it under dictory Data/parameter
# Execute the above command regularly.
#
# Edit Webapp/start_webapp.py, uncomment app.run(debug=True) and comment out the previous app.run(). To display the network, cd Webpap && python start_webapp.py.
# Enter http://127.0.0.1:5000 in net browser.
#
# Note:
# Tested on a 32-core server at slcu running Ubuntun. The program may slow down a personal computer.
# The program will check that the following required packages are installed.
# Required python packages: numpy, networkx, flask.
# Required R packages: rjson, mixtools, Rtsne.
# To install a python package, use the command pip install numpy.
# To install an R package, issue this command in R: install.packages('Rtsne', dependencies=TRUE, repos='http://cran.rstudio.com/')
# In a Mac, use export PYTHONPATH="/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages" to make installed python packages usable.
#
# Purpose: make a gene regulatory network locally. Periodically (e.g., per day) run this script to see if the network needs update. If yes, update it.
# You need to prepare and edit parameter_for_buildCmatrix.txt manually.
#
# Created 1 July 2017, hui.lan@slcu.cam.ac.uk, slcu
# Last modified 4 July 2017, hui lan, slcu
import os, sys
import numpy as np
import glob
import time
import subprocess
from datetime import datetime
from param4net import make_global_param_dict, get_key_value
FORCE_MAKE_EDGES = 'NO'
CODE_DIR = os.getcwd() # Get current working directory. It is important that you execute this script under it directory.
# DON'T CHANGE THE FOLLOWING PATHS AND NAMES
HISTORY_DIR = '../Data/history/edges/many_targets' # each edge file contains edges for many targets
HISTORY_DIR2 = '../Data/history/edges/one_target' # edges.txt.* files are here, all edge files have the name edges.txt.*, the leading string 'edges.txt' must be present.
FILE_TIMESTAMP = '../Data/log/file_timestamp.txt'
SAMPLE_SIZE_FILE = '../Data/log/total.samples.txt' # each line contains a date and the number of samples on and after that date
TEMP_DIR = '../Data/temp'
PARAMETER_FOR_BUILDCMATRIX = '../Data/parameter/parameter_for_buildCmatrix.txt'
PARAMETER_FOR_BUILDRMATRIX = '../Data/parameter/parameter_for_buildRmatrix.txt'
PARAMETER_FOR_NET = '../Data/parameter/parameter_for_net.txt'
EDGE_FILE = '../Data/history/edges/edges.txt'
BINDING_FILE = '../Data/history/bind/binding.txt'
TPM_FILE = '../Data/history/expr/TPM.txt'
LOG_FILE = '../Data/log/update.network.log.txt'
NEW_OR_UPDATED_CHIP_FILE = '../Data/log/new.or.updated.chip.file.txt'
RNA_SEQ_INFO_DATABASE = '../Data/information/rnaseq_info_database.txt'
RNA_SEQ_INFO_DATABASE_JSON = '../Data/information/rnaseq_info_database.json'
FILE_LIST_TO_CHECK = [PARAMETER_FOR_BUILDCMATRIX,
PARAMETER_FOR_BUILDRMATRIX,
PARAMETER_FOR_NET,
EDGE_FILE,
BINDING_FILE,
TPM_FILE]
## help functions
def ok_webapp_dir(para_for_net):
''' we are now under Code '''
glb_param_dict = make_global_param_dict(para_for_net)
# genes.json is not here, create one
if not os.path.exists('../Webapp/static/json/genes.json'):
print('genes.json not here, make one')
cmd = 'python text2json.py %s > ../Webapp/static/json/genes.json' % (glb_param_dict['GENE_ID_AND_GENE_NAME'])
os.system(cmd)
def make_paths(s):
if not os.path.isdir(s):
os.makedirs(s)
def make_important_dirs():
make_paths('../Data/history/edges/many_targets')
make_paths('../Data/history/edges/one_target')
make_paths('../Data/log')
make_paths('../Data/information')
make_paths('../Data/temp')
make_paths('../Data/parameter')
make_paths('../Data/R/Mapped')
make_paths('../Data/R/Mapped/public')
make_paths('../Data/R/Mapped/inhouse')
make_paths('../Data/R/Mapped/other')
make_paths('../Data/R/Raw')
make_paths('../Data/C/Mapped')
make_paths('../Data/C/Mapped/Columns')
make_paths('../Data/C/Raw')
make_paths('../Data/history/edges')
make_paths('../Data/history/bind')
make_paths('../Data/history/expr')
make_paths('../Webapp/static/json')
make_paths('../Webapp/static/edges')
make_paths('../Webapp/templates')
def num_line(fname):
f = open(fname)
lines = f.readlines()
f.close()
return len(lines)
def record_new_or_updated_chip_ids(lst, fname):
f = open(fname, 'a')
curr_time = datetime.now().strftime('%Y%m%d')
for x in lst:
f.write('%s\t%s\n' % (x, curr_time))
f.close()
def write_log_file(s, fname):
f = open(fname, 'a')
print(s)
curr_time = datetime.now().strftime('%Y-%m-%d %H:%M')
s = '[' + curr_time + ']: ' + s
if not '\n' in s:
s += '\n'
f.write(s)
f.close()
def write_sample_size_file(sample_size_file, curr_date, tpm_sample_size):
if not os.path.exists(sample_size_file):
f = open(sample_size_file, 'w')
else:
f = open(sample_size_file, 'a')
f.write('%s\t%s\n' % (curr_date, tpm_sample_size))
f.close()
def age_of_file(fname):
st = os.stat(fname)
days = (time.time() - st.st_mtime)/(3600*24.0)
return days
def hold_on(fname):
f = open(fname)
lines = f.readlines()
f.close()
for line in lines[:100]: # check the first 100 lines for HOLDON
line = line.strip()
if line.startswith('%%HOLDON=YES'):
return True
return False
def all_files_present(lst):
missing_file_lst = []
for path in lst: # lst is a list of file names to check
if not os.path.exists(path):
if 'edges.txt' in path:
write_log_file('WARNING: must have %s to update network.' % (path), LOG_FILE)
missing_file_lst.append(path)
return missing_file_lst
def record_file_time(lst):
f = open(FILE_TIMESTAMP, 'w')
s = ''
for x in lst:
if os.path.exists(x):
s += '%s\t%d\n' % (os.path.basename(x), int(os.stat(x).st_mtime))
else:
s += '%s\t%d\n' % (os.path.basename(x), 0)
f.write(s)
f.close()
def read_file_timestamp(ftimestamp):
d = {}
f = open(ftimestamp)
for line in f:
line = line.strip()
lst = line.split()
fname = lst[0]
t = lst[1]
d[fname] = int(t)
f.close()
return d
def file_updated(fname, d):
ft = int(os.stat(fname).st_mtime)
k = os.path.basename(fname)
return ft > d[k]
def get_updated_files(lst, d):
result = []
for x in lst:
if file_updated(x, d):
result.append(os.path.basename(x))
return result
def not_bad_line(s):
if s.strip() == '':
return False
if 'WARNING' in s:
return False
if 'number' in s:
return False
if 'Need' in s:
return False
if 'Error' in s:
return False
if 'Too' in s:
return False
if not s.startswith('AT'): # comment out this test if the organism is not Arabidopsis CHANGE
return False
return True
def get_rcond_string(s):
s = s.strip()
if s.startswith('R0000DRR') or s.startswith('R0000ERR') or s.startswith('R0000SRR'):
s = s.replace('R0000DRR', 'R0DRR').replace('R0000SRR', 'R0SRR').replace('R0000ERR', 'R0ERR') # remove extra 0's introduced in earlier edges.txt
return s
def compute_metric(d):
'''
d has the form {'freq':0, 'total_RNAseq_ID':0, 'sum_abs_R':0,
'most_recent_date':'20161201'}.
Metric is a combination of sevaral quantities, average absolute
correlation coefficient, average number of RNA-seq data used,
frequency of this edge, and recentness. It is used to rank edges.
So larger correlation cofficients based on more RNA-seq data that
frequenly appear and recent will be ranked on top.
Formula: avg.abs(r)*log(avg.RN)*log(F+1)*(CurrentDate-MostRecentDate)^-0.2
'''
avg_abs_r = 1.0 * d['sum_abs_R'] / d['freq']
log_avg_RN = np.log10(1.0 * d['total_RNAseq_ID'] / d['freq'])
log_freq = np.log2(d['freq'] + 1)
S = 200.0 # strength of memory, larger value means better memory
recentness = np.exp(1.0*(int(d['most_recent_date']) - int(datetime.now().strftime('%Y%m%d')))/S) # a monotonic decreasing function exp(-t/S), exponential nature of forgetting
return avg_abs_r * log_avg_RN * log_freq * recentness
def make_sample_size_dict(fname):
if not os.path.exists(fname):
return {}
f = open(fname)
lines = f.readlines()
f.close()
d = {}
for line in lines:
line = line.strip()
if line != '' and not line.startswith('#'):
lst = line.split('\t')
d[lst[0]] = int(lst[1])
return d, sorted(d.keys())
def get_sample_size(d, sorted_keys, day):
if len(d) == 0:
return 1200 # a default number of sample size, CHANGE
for x in sorted_keys:
if x >= day:
return d[x]
k = sorted_keys[-1] # last key, latest date
return d[k]
def merge_edges(file_lst, edge_fname, sample_size_file):
'''
Write to edge_fname.
file_lst -- a list of edges.txt.*, computed using different sets
of data. Each edges.txt file has the same format.
This function merges all edges together. Correlation based on
larger number of RNA-seq experiments is favoured. Each edge has
the following format: target_gene tf_gene value type
RNA_experiments ChIP_experiments loglik date, and metric. metric is newly computed.
'''
sample_size_dict, sample_size_keys = make_sample_size_dict(sample_size_file)
ll_dict = { # loglikehood to tissue or method
'-999.0':'seedling',
'-998.0':'meristem',
'-997.0':'root',
'-996.0':'leaf',
'-995.0':'flower',
'-994.0':'shoot',
'-993.0':'seed',
'-992.0':'stem',
'-990.0':'aerial',
'-991.0':'hclust',
'-1000.0':'wedge (post.translation.3)',
'-1001.0':'wedge (post.translation.4)',
'.':'all'
}
d = {} # hold edges, best will be kept for mix pos, mix neg and all.
d_mix_pos = {} # for computing rank metric 'freq':0, 'total_RNAseq_ID':0, 'sum_abs_R':0, 'most_recent_date':'20161201'
d_mix_neg = {}
d_all = {}
d_tissue = {} # hold subset information (tissue or method)
max_rcond_size = 10
for fname in file_lst: # check each edge file
if not os.path.exists(fname):
write_log_file('%s missing.' % (fname), LOG_FILE)
continue
f = open(fname)
lines = f.readlines()
f.close()
for line in lines:
line = line.strip()
lst = line.split('\t') # fields are separated by TABs
if not_bad_line(line) and len(lst) >= 7: # an edge line must have at least 7 fields
target_id = lst[0].split()[0] # AGI id only, no gene name
tf_id = lst[1].split()[0] # same as above
k = target_id + '_' + tf_id # key
date = '20161201' # edge creation date. if there's no date field, use 20161201
rcond_s = get_rcond_string(lst[4]) # remove extra 0's from old (obsolete) RNA-seq experiment ids
loglik = lst[6]
t = lst[3] # current line edge type, mix or all
if '=' in loglik:
loglik = loglik.split('=')[1]
if t == 'mix' and loglik == '.': # post translation case, i.e., post.translation.3, see create_edegs4.py
loglik = '-1000.0'
tissue_or_method_name = ll_dict[loglik] if loglik in ll_dict else 'MixReg'
if len(lst) == 8: # has a date field
date = lst[7]
# add for each edge an tissue or method information. For example, the edge is based on a certain tissue, or is derived by a certain method. ll_dict contains the mapping.
if not k in d_tissue:
d_tissue[k] = [tissue_or_method_name]
else:
d_tissue[k].append(tissue_or_method_name)
if not k in d: # first edge to be added between two nodes: target_id and tf_id
d[k] = [{'target':lst[0], 'tf':lst[1], 'value':lst[2], 'type':lst[3], 'rcond':rcond_s, 'ccond':lst[5], 'loglik':loglik, 'date':date}] # a list of dicitionaries, {target, tf, value, type, rcond, ccond, loglik, date}
else: # two nodes can have multiple edges. an edge already exists, determine whether or not add this new edge.
v = lst[2] # current line value
rcond = rcond_s
ccond = lst[5]
len_r = len(rcond.split()) # number of RNA-seq experiment ids in current line. if it is a dot, then length is 1
ignore = False # assume this edge is to be included. will be set False if otherwise.
for i in range(len(d[k])): # search each of already added edges. If current one is better, replace the old one
xd = d[k][i] # xd is a dictionary, representing an edge {target, tf, type, value, rcond, ccond, loglik, date}
fv = float(v) # value (confidence) of the candiadte edge
fx = float(xd['value'])
len_rx = len(xd['rcond'].split())
if xd['type'] == t and t == 'all': # edge type is all, use most recent result as the edge is based on all RNA-seq experiments.
ignore = True # either ignore it or replace the old with this one
if date > xd['date']:
d[k][i]['value'] = v
d[k][i]['date'] = date
break
if xd['type'] == t and t == 'mix': # current line represents a better edge, i.e., larger r value and based on more RNA-seq experiments
if fv*fx > 0 and abs(fv*np.log10(len_r)) > abs(fx*np.log10(len_rx)): # fv*fx > 0 means they have the same sign.
d[k][i]['value'] = v
d[k][i]['rcond'] = rcond
d[k][i]['ccond'] = ccond
d[k][i]['loglik'] = loglik
d[k][i]['date'] = date
ignore = True
break
elif fv*fx > 0: # curr line has same sign, but based on less RNA-seq experiment, ignore it.
ignore = True
break
if v == xd['value'] and len_r == len_rx: # ChIPs are updated, but based on same number of rna-seq experiments
if xd['ccond'] != ccond:
merged_cond = xd['ccond'] + ' ' + ccond
merged_cond_str = ' '.join(sorted(list(set(merged_cond.split()))))
d[k][i]['ccond'] = merged_cond_str
if xd['date'] < date:
d[k][i]['date'] = date
d[k][i]['loglik'] = loglik
ignore = True
break
if ignore == False:
d[k].append({'target':lst[0], 'tf':lst[1], 'value':lst[2], 'type':lst[3], 'rcond':rcond_s, 'ccond':lst[5], 'loglik':loglik, 'date':date})
# fill d_mix_pos, d_mix_neg and d_all
curr_rcond_size = len(rcond_s.split())
if t == 'mix':
if float(lst[2]) >= 0: # lst[2] is value
if not k in d_mix_pos:
d_mix_pos[k] = {'freq':1, 'total_RNAseq_ID':curr_rcond_size, 'sum_abs_R':abs(float(lst[2])), 'most_recent_date':date}
else:
d_mix_pos[k]['freq'] += 1
d_mix_pos[k]['total_RNAseq_ID'] += curr_rcond_size
d_mix_pos[k]['sum_abs_R'] += abs(float(lst[2]))
if date > d_mix_pos[k]['most_recent_date']:
d_mix_pos[k]['most_recent_date'] = date
else:
if not k in d_mix_neg:
d_mix_neg[k] = {'freq':1, 'total_RNAseq_ID':curr_rcond_size, 'sum_abs_R':abs(float(lst[2])), 'most_recent_date':date}
else:
d_mix_neg[k]['freq'] += 1
d_mix_neg[k]['total_RNAseq_ID'] += curr_rcond_size
d_mix_neg[k]['sum_abs_R'] += abs(float(lst[2]))
if date > d_mix_neg[k]['most_recent_date']:
d_mix_neg[k]['most_recent_date'] = date
if curr_rcond_size > max_rcond_size:
max_rcond_size = curr_rcond_size
if t == 'all':
all_rcond_size = get_sample_size(sample_size_dict, sample_size_keys, date)
if not k in d_all:
d_all[k] = {'freq':1, 'total_RNAseq_ID':all_rcond_size, 'sum_abs_R':abs(float(lst[2])), 'most_recent_date':date}
else:
d_all[k]['freq'] += 1
d_all[k]['total_RNAseq_ID'] += all_rcond_size
d_all[k]['sum_abs_R'] += abs(float(lst[2]))
if date > d_all[k]['most_recent_date']:
d_all[k]['most_recent_date'] = date
# rewrite all total_RNAseq_ID value, the total RNAseq ID is actually an overestimate.
if len(sample_size_dict) == 0: # sample size file is not available
for k in d_all:
d_all[k]['total_RNAseq_ID'] = d_all[k]['freq'] * max_rcond_size
f = open(edge_fname, 'w')
for k in sorted(d.keys()):
lst = d[k] # a list of dictionaries
tissue_or_method = ', '.join(list(set(d_tissue[k])))
for xd in lst:
if xd['type'] == 'all':
metric = '%4.2f' % ( compute_metric(d_all[k]) )
s = '\t'.join([xd['target'], xd['tf'], xd['value'], xd['type'], xd['rcond'], xd['ccond'], xd['loglik'], xd['date'], metric, tissue_or_method]) + '\n'
if xd['type'] == 'mix':
if float(xd['value']) >= 0:
metric = '%4.2f' % ( compute_metric(d_mix_pos[k]) )
else:
metric = '%4.2f' % ( compute_metric(d_mix_neg[k]) )
s = '\t'.join([xd['target'], xd['tf'], xd['value'], xd['type'], xd['rcond'], xd['ccond'], xd['loglik'], xd['date'], metric, tissue_or_method]) + '\n'
f.write(s)
f.close()
def get_value(s, delimit):
lst = s.split(delimit, 1) # only split at the first delimit
return lst[1].strip()
def make_data_dict(fname):
d = {} # keep a list of chip id's, such as C0001100007100
f = open(fname)
lines = f.readlines()
f.close()
for line in lines:
line = line.strip()
if line == '' or line.startswith('#'):
continue
if line.startswith('@'):
s = line[line.rfind('@')+1:]
s = s.strip()
if s in d:
write_log_file('In make_data_dict: ID %s duplicated' % (s), LOG_FILE)
sys.exit()
d[s] = {'PROTEIN_ID':'', 'PROTEN_NAME':'', 'DATA_NAME':'', 'DATA_FORMAT':'', 'DESCRIPTION':'', 'LOCATION':'', 'NOTE':''}
if line.startswith('DESCRIPTION:'):
d[s]['DESCRIPTION'] = get_value(line, ':')
elif line.startswith('PROTEN_NAME:'):
d[s]['PROTEN_NAME'] = get_value(line, ':')
elif line.startswith('PROTEIN_ID:'):
d[s]['PROTEIN_ID'] = get_value(line, ':')
elif line.startswith('DATA_NAME:'):
d[s]['DATA_NAME'] = get_value(line, ':')
elif line.startswith('DATA_FORMAT:'):
d[s]['DATA_FORMAT'] = get_value(line, ':')
elif line.startswith('LOCATION:'):
d[s]['LOCATION'] = get_value(line, ':')
elif line.startswith('NOTE:'):
d[s]['NOTE'] = get_value(line, ':')
return d
def get_bad_chip_ids(d):
''' a id chip is bad if its Note field has obsolete '''
lst = []
for k in d:
note = d[k]['NOTE'].lower()
if 'obsolete' in note:
lst.append(k)
return lst
def get_update_date_chip_ids(d):
''' Return a list of ChIP ids with update in the NOTE: field. '''
ud = {}
for k in d:
note = d[k]['NOTE'].lower()
if 'update' in note:
if 'update:' in note: # has a specific date, e.g., 20170101
idx = note.find('update:')
udate = note[idx+7:idx+15] # get date string yyyymmdd
else: # if only update, but no specific date, assume it is 20161201
udate = '20161201'
ud[k] = udate
return ud
def get_chip_ids_from_edge_file(fname):
lst = []
f = open(fname)
for line in f:
line = line.strip()
l = line.split('\t')
c = l[5] # ids for ChIP experiments
l2 = c.split()
lst.extend(l2)
f.close()
return list(set(lst))
def get_nonexistent_chip_ids(d, fname):
''' fname -- edges.txt. Get ids that are in fname, but not in d. '''
lst = []
ids = get_chip_ids_from_edge_file(fname)
for k in ids:
if not k in d:
lst.append(k)
return lst
def created_after_update(x, created_date, update_date_chip_ids):
''' the edge is created after the ChIP is updated.'''
if not x in update_date_chip_ids: # update_date_chip_ids dose not contain id x, which means x is not updated since ChIP id is created
return True
else:
return created_date >= update_date_chip_ids[x] # check whether edge is created after x is updated.
def rm_chip_ids_from_edge_file(fname, bad_id_lst, update_date_chip_ids):
''' fname -- edges.txt '''
f = open(fname)
lines = f.readlines()
f.close()
f = open(fname, 'w')
for line in lines:
line = line.strip()
lst = line.split('\t')
c = lst[5] # chip experiment ids
l = c.split() # list of chip experiment ids
created_date = lst[7]
l2 = [] # for keeping good chip experiment ids
for x in l:
if not x in bad_id_lst and created_after_update(x, created_date, update_date_chip_ids): # if an edge is created before this ChIP experiment is updated, then the binding information and thus the edge may be no longer valid. If this ChIP experiment is the sole evidence of binding, then this edge is ignored.
l2.append(x)
if l2 != []: # still have some chips
lst[5] = ' '.join(l2)
f.write('\t'.join(lst) + '\n')
f.close()
def make_file(fname, s):
f = open(fname, 'w')
f.write(s)
f.close()
def get_chip_ids(fname):
if not os.path.exists(fname):
print('ERROR: %s not exists.' % (fname))
sys.exit()
f = open(fname)
lines = f.readlines()
f.close()
head_line = lines[0].strip()
lst = head_line.split('\t')
return lst[1:]
def is_recent(note, ndays): # less than 2-days old
idx = note.find('update:')
if idx == -1: # not found
return False
udate = note[idx+7:idx+15] # get date string yyyymmdd
curr_date = datetime.now().strftime('%Y%m%d')
d1 = datetime.strptime(udate, "%Y%m%d")
d2 = datetime.strptime(curr_date, "%Y%m%d")
return (d2 - d1).days <= ndays
def get_note_date(s):
''' s has this format: NOTE: update:20170101 '''
if not 'update:' in s: # if people forgot to put an update date in Note:, then it is 20161201
return '20161201'
index = s.rfind('update:')
if len(s[index+7:]) < 8:
return '20161201'
result = s[index+7:index+15].strip()
if result <= '20170101' or len(result) < 8:
return '20161201'
return result
def is_recently_updated(note, chip_id, fname):
''' fname keeps track of chip ids and their most recent update date '''
if not is_recent(note, 15): # if update happened 15 days ago, just ignore it. no further action will be taken for this chip id.
return False
note_date = get_note_date(note)
log_date = '20161201'
if os.path.exists(fname): # search update date of chip_id, if it has been incorperated, then log date should be no less than note's update date
f = open(fname)
for line in f:
line = line.strip()
if line != '':
lst = line.split()
if len(lst) >= 2 and lst[0] == chip_id:
log_date = lst[1]
f.close()
return note_date > log_date
def get_new_chip_ids(old_lst, para_c_dict, new_or_updated_chip_fname):
result = []
for k in para_c_dict:
note = para_c_dict[k]['NOTE'].lower()
if 'obsolete' in note:
continue
if (not 'obsolete' in note and not k in old_lst) or (k in old_lst and 'update:' in note and is_recently_updated(note, k, new_or_updated_chip_fname)): # First case: the ChIP-seq ID is not in old list and not marked obsolte. Definitely add it. A second case is more subtle. The narrowPeak of an old ChIP-seq ID has recently been updated for whatever purposes. This update should happen reasonably recently. Why is_recently_updated is important? Because we don't want to treat ChIP-seq that is updated weeks ago as new. new_or_updated_chip_fname keeps track of newly updated ChIP-seq and their update dates.
result.append(k)
return sorted(result)
def edit_file(fname, line_starter, new_str):
f = open(fname)
lines = f.readlines()
f.close()
f = open(fname, 'w')
for line in lines:
if line.startswith(line_starter):
s = line_starter
s += new_str
f.write(s + '\n')
else:
f.write(line)
f.close()
def number_rnaseq_id(tpm_file):
f = open(tpm_file)
first_line = f.readlines()[0]
f.close()
first_line = first_line.strip()
return len(first_line.split()) - 1
def number_rnaseq_diff(para_file, tpm_file):
''' count the number @ in para_file, and count the number of columns in tpm_file, return their difference '''
a = 0
f = open(para_file)
for line in f:
line = line.strip()
if line.startswith('@'):
a += 1
f.close()
b = number_rnaseq_id(tpm_file)
return a - b
def get_key_value(s):
lst = s.split('=')
k, v = lst[0], lst[1]
return (k.strip(), v.strip())
def get_value(s, delimit):
lst = s.split(delimit, 1)
return lst[1].strip()
def validate_gene_file(fname):
f = open(fname)
lines = f.readlines()
f.close()
for line in lines: # check all lines
line = line.strip()
lst = line.split('\t')
if len(lst) < 6:
print('Not enought fields: %s. Only %d are given. Each line must have gene_id, gene_name, chr, start, end, strand, description (optional). See prepare_gene_file.py in the documentation on how to prepare this file.' % (line, len(lst)))
sys.exit()
def validate_parameter_for_buildcmatrix(fname):
# first the file must exist
if not os.path.exists(fname):
print('CANNOT FIND %s.' % (fname))
sys.exit()
f = open(fname)
lines = f.readlines()
f.close()
d = {}
location_count = 0
for line in lines:
line = line.strip()
if line.startswith('%%'):
k, v = get_key_value(line[2:])
d[k] = v
if k == 'GENE_FILE' or k == 'CHR_INFO':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'GENE_FILE':
validate_gene_file(v)
if k == 'DESTINATION':
if not os.path.isdir(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'TARGET_RANGE':
if int(v) <= 0:
print('Target range (%d) must be greater than 0.' % (v))
sys.exit()
if line.startswith('LOCATION:'):
v = get_value(line, ':')
location_count += 1
if not os.path.exists(v):
print('Location %s does not exists.' % (v))
sys.exit()
if not 'GENE_FILE' in d:
print('Must specify GENE_FILE.')
sys.exit()
if not 'DESTINATION' in d:
print('Must specify DESTINATION.')
sys.exit()
if not 'CHR_INFO' in d:
print('Must specify CHR_INFO.')
sys.exit()
if location_count == 0:
print('Must contain at least one ChIP-seq.')
sys.exit()
def validate_parameter_for_buildrmatrix(fname):
# first the file must exist
if not os.path.exists(fname):
print('CANNOT FIND %s.' % (fname))
sys.exit()
f = open(fname)
lines = f.readlines()
f.close()
d = {}
location_count = 0
for line in lines:
line = line.strip()
if line.startswith('%%'):
k, v = get_key_value(line[2:])
d[k] = v
if k == 'GENE_LIST':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if line.startswith('LOCATION:'):
v = get_value(line, ':')
location_count += 1
if not os.path.exists(v):
print('Location %s does not exists.' % (v))
sys.exit()
if not 'GENE_LIST' in d:
print('Must specify GENE_LIST.')
sys.exit()
if location_count == 0:
print('Must contain at least one RNA-seq.')
sys.exit()
def validate_parameter_for_net(fname):
# first the file must exist
if not os.path.exists(fname):
print('CANNOT FIND %s.' % (fname))
sys.exit()
f = open(fname)
lines = f.readlines()
f.close()
d = {}
location_count = 0
for line in lines:
line = line.strip()
if line.startswith('%%'):
k, v = get_key_value(line[2:])
d[k] = v
if k == 'GENE_LIST':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'GENE_ID_AND_GENE_NAME':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'BINDING_INFO':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'EXPRESSION_INFO':
if not os.path.exists(v):
print('%s not exists.' % (v))
sys.exit()
if k == 'BINDING_MATRIX':
if not os.path.exists(v):
print('%s not exists.' % (v))
#print('Use python buildCmatrix.py parameter_for_buildCmatrix.txt > binding.txt to create binding.txt.')
if k == 'EXPRESSION_MATRIX':
if not os.path.exists(v):
print('%s not exists.' % (v))
print('Use python buildRmatrix.py parameter_for_buildRmatrix.txt to create TPM.txt.')
if not 'GENE_LIST' in d:
print('Must specify GENE_FILE.')
sys.exit()
if not 'GENE_ID_AND_GENE_NAME' in d:
print('Must specify GENE_ID_AND_GENE_NAME.')
sys.exit()
if not 'BINDING_INFO' in d:
print('Must specify BINDING_INFO.')
sys.exit()
if not 'EXPRESSION_INFO' in d:
print('Must specify EXPRESSION_INFO.')
sys.exit()
if not 'BINDING_MATRIX' in d:
print('%s not exists.' % (v))
print('Use python buildCmatrix.py paramter_for_buildCmatrix.txt > binding.txt to create binding.txt.')
if not 'EXPRESSION_MATRIX' in d:
print('%s not exists.' % (v))
print('Use python buildRmatrix.py paramter_for_buildRmatrix.txt to create TPM.txt.')
def file_contains(fname, s):
if not os.path.exists(fname):
return False
f = open(fname)
for line in f:
if s in line:
return True
f.close()
return False
def check_required_packages():
# mixtools, networkx, rjson, flask
# python libraries
try:
import numpy
except ImportError, e:
print('numpy not available. Install it via: pip install numpy')
sys.exit()
try:
import networkx
except ImportError, e:
print('networkx not available. Install it via: pip install networkx')
sys.exit()
try:
import flask
except ImportError, e:
print('flask not available. Install it via: pip install flask')
sys.exit()
# R libraries
cmd = 'echo \"is.installed <- function(mypkg) is.element(mypkg, installed.packages()[,1]); is.installed(\'mixtools\')\" > ../Data/temp/check.r.library.R && Rscript ../Data/temp/check.r.library.R > ../Data/temp/check.r.library.result'
os.system(cmd)
if not file_contains('../Data/temp/check.r.library.result', 'TRUE'):
print('R package mixtools not available. Install it first.')
sys.exit()
os.remove('../Data/temp/check.r.library.result')
cmd = 'echo \"is.installed <- function(mypkg) is.element(mypkg, installed.packages()[,1]); is.installed(\'rjson\')\" > ../Data/temp/check.r.library.R && Rscript ../Data/temp/check.r.library.R > ../Data/temp/check.r.library.result'
os.system(cmd)
if not file_contains('../Data/temp/check.r.library.result', 'TRUE'):
print('R package rjson not available. Install it first.')
sys.exit()
os.remove('../Data/temp/check.r.library.result')
cmd = 'echo \"is.installed <- function(mypkg) is.element(mypkg, installed.packages()[,1]); is.installed(\'Rtsne\')\" > ../Data/temp/check.r.library.R && Rscript ../Data/temp/check.r.library.R > ../Data/temp/check.r.library.result'
os.system(cmd)
if not file_contains('../Data/temp/check.r.library.result', 'TRUE'):
print('R package Rtsne not available. Install it first.')
sys.exit()
os.remove('../Data/temp/check.r.library.result')
## main
## Shipped with this distribution: TPM.txt.gz, experiment.and.tissue.txt, rnaseq_info_database.json
## BINDING_FILE='../Data/history/bind/binding.txt' is not given, as users should provide a BED file to create it.
if not os.path.isdir(CODE_DIR): # make sure that CODE_DIR exists
print('ERROR: %s does not exists.' % (CODE_DIR))
sys.exit()
os.chdir(CODE_DIR) # run this file at the Code directory
if not os.path.isdir('../Data/C/Mapped'):
make_important_dirs() # make important directories (if non-existentt) for holding all kinds of files, must be after os.chdir(CODE_DIR)
write_log_file('Copy BED files to Data/C/Mapped. Update the LOCATION field in Data/parameter/parameter_for_buildCmatrix.txt. Run local_network.py again.' , LOG_FILE)
sys.exit()
ok_webapp_dir(PARAMETER_FOR_NET) # make sure Webapp contains necessary files
# BEFORE WE START, WE SHOULD CHECK REQUIRED SOFTWARES ARE INSTALLED.
# For example, numpy, mixtools, networkx, rjson, flask (TBD)
# make sure required packages are available
check_required_packages()
# check rnaseq_info_database.txt and rnaseq_info_database.json
if not os.path.exists(RNA_SEQ_INFO_DATABASE):
write_log_file('NEED CREATE %s.' % (RNA_SEQ_INFO_DATABASE), LOG_FILE)
sys.exit()
if not os.path.exists(RNA_SEQ_INFO_DATABASE_JSON):
write_log_file('NEED CREATE %s.' % (RNA_SEQ_INFO_DATABASE_JSON), LOG_FILE)
sys.exit()
# make sure parameter files are present and valid (very rudimentary check but important)
validate_parameter_for_buildcmatrix(PARAMETER_FOR_BUILDCMATRIX)
#validate_parameter_for_buildrmatrix(PARAMETER_FOR_BUILDRMATRIX)
validate_parameter_for_net(PARAMETER_FOR_NET)
# remove binding file, if any
if os.path.exists(BINDING_FILE):
os.remove(BINDING_FILE)
# update edges.txt, a merged file from several sources, HISTORY_DIR and HISTORY_DIR2.
edge_file_lst = [] # collect edge files.
for fname in glob.glob(os.path.join(HISTORY_DIR, 'edges.txt.*')): # edges.txt.* are to be merged
edge_file_lst.append(fname)
for fname in glob.glob(os.path.join(HISTORY_DIR2, 'edges.txt.*')): # edges.txt.* are to be merged
edge_file_lst.append(fname)
# make sure all needed files are present, if not, make them if possible
miss_lst = all_files_present(FILE_LIST_TO_CHECK) # check if any of them are missing
if (miss_lst != [] and edge_file_lst == []) or FORCE_MAKE_EDGES == 'YES':
write_log_file('Cannot find these required files: %s. The program will prepare them.' % (', '.join(miss_lst)), LOG_FILE)
# initially, we only have three parameter files, but not binding.txt
important_miss_number = 0
if PARAMETER_FOR_BUILDCMATRIX in miss_lst:
print('Must prepare %s first' % (PARAMETER_FOR_BUILDCMATRIX))
important_miss_number += 1
if PARAMETER_FOR_BUILDRMATRIX in miss_lst:
print('Must prepare %s first' % (PARAMETER_FOR_BUILDRMATRIX))
important_miss_number += 1
if PARAMETER_FOR_NET in miss_lst:
print('Must prepare %s first' % (PARAMETER_FOR_NET))
important_miss_number += 1
if important_miss_number > 0:
sys.exit() # need to provide all the above three files; otherwise cannot proceed
target_tf_fname = '../Data/information/target_tf.txt'
if os.path.exists(target_tf_fname):
os.remove(target_tf_fname)
if BINDING_FILE in miss_lst:
write_log_file('Make initial binding.txt ...', LOG_FILE)
cmd = 'python get_binding.py %s' % (PARAMETER_FOR_BUILDCMATRIX)
os.system(cmd)
cmd = 'python buildCmatrix.py %s > %s' % (PARAMETER_FOR_BUILDCMATRIX, BINDING_FILE)
os.system(cmd)
#print('IMPORATNT:Now check BINDING_MATRIX in %s is set %s and rerun local_network.py.' % (PARAMETER_FOR_NET, BINDING_FILE))
#sys.exit()
if TPM_FILE in miss_lst:
if not os.path.exists(TPM_FILE + '.gz'):
write_log_file('Cannot find %s. Try %s.' % (TPM_FILE + '.gz', TPM_FILE), LOG_FILE)
if not os.path.exists(TPM_FILE):
sys.exit()
else:
write_log_file('Unzip initial TPM.txt', LOG_FILE)
cmd = 'gunzip %s' % (TPM_FILE + '.gz')
os.system(cmd)
#print('IMPORTANT:Now check EXPRESSION_MATRIX in %s is set %s and rerun local_network.py.' % (PARAMETER_FOR_NET, TPM_FILE))
#sys.exit()
miss_lst2 = all_files_present(FILE_LIST_TO_CHECK) # check files again
if (len(miss_lst2) == 1 and miss_lst2[0] == EDGE_FILE) or (len(miss_lst2) == 0 and os.path.getmtime(EDGE_FILE) < os.path.getmtime(BINDING_FILE)): # all other files are ready except edges.txt, make one.
# should assert the files needed by the following scripts are present
# big correlation matrix may not eat all computer's memory, so need to modify the code. or test the user's memory first.
print('Make some edges, wait ... Change MAX_PROCESS in Code/create_edges4.py to change the number of processes. Default=10.')
cmd = 'nohup python create_edges4.py %s &' % (PARAMETER_FOR_NET) # this will create target_tf.txt needed by the following scripts
#os.system(cmd)
time.sleep(10)
# wait or make target_tf.txt
wait_sec = 0
WAIT_SECONDS = 120
while not os.path.exists(target_tf_fname):
time.sleep(WAIT_SECONDS)
wait_sec += WAIT_SECONDS
if wait_sec > WAIT_SECONDS * 5:
write_log_file('Make Data/information/target_tf.txt', LOG_FILE)
cmd = 'python make_target_tf.py %s > %s' % (PARAMETER_FOR_NET , target_tf_fname) # make target_tf.txt CHANGE better to make a temperory copy for this program
os.system(cmd)
break
time.sleep(5)
write_log_file('Create group-specific (fixed) edges.txt using new TPM.txt (size=%d).' % (number_rnaseq_id(TPM_FILE)), LOG_FILE)
cmd = 'Rscript correlation_per_group_fixed_number.R &'
os.system(cmd)
time.sleep(5)
write_log_file('Create SIMPLE edges.txt using new TPM.txt (size=%d). SIMPLE means using all RNA-seq samples.' % (number_rnaseq_id(TPM_FILE)), LOG_FILE)
cmd = 'python create_edges0.py %s &' % (PARAMETER_FOR_NET) # use all samples in TPM.txt results will be written to history/edges.txt.simple.correlation.all.conditions
os.system(cmd)
time.sleep(5)
write_log_file('Create tissue-specific edges.txt using new TPM.txt (size=%d).' % (number_rnaseq_id(TPM_FILE)), LOG_FILE)
cmd = 'python create_edges0B.py %s &' % (PARAMETER_FOR_NET) # call correlation_per_tissue.R
os.system(cmd)
time.sleep(5)
write_log_file('Create group-specific edges.txt using new TPM.txt (size=%d).' % (number_rnaseq_id(TPM_FILE)), LOG_FILE)
cmd = 'Rscript correlation_per_group.R &'
os.system(cmd)
time.sleep(5)
write_log_file('Create wedge-shape edges.txt using new TPM.txt (size=%d).' % (number_rnaseq_id(TPM_FILE)), LOG_FILE)
cmd = 'Rscript wedge.R &'
os.system(cmd)
# make json (sliced TPM.txt) and json2 (sliced binding.txt) if they don't exist
if os.path.exists(TPM_FILE):
if not os.path.isdir('../Data/history/expr/json'):
write_log_file('Make directory ../Data/history/expr/json', LOG_FILE)
cmd = 'python slice_TPM_to_JSON.py %s' % (PARAMETER_FOR_NET)
os.system(cmd)
if os.path.exists(BINDING_FILE):
if not os.path.isdir('../Data/history/bind/json2'):
write_log_file('Make directory ../Data/history/bind/json2', LOG_FILE)
cmd = 'python slice_binding_to_JSON.py %s' % (PARAMETER_FOR_NET)
os.system(cmd)
elif os.path.getmtime('../Data/history/bind/json2') < os.path.getmtime(BINDING_FILE):
write_log_file('Make directory ../Data/history/bind/json2', LOG_FILE)
cmd = 'python slice_binding_to_JSON.py %s' % (PARAMETER_FOR_NET)
os.system(cmd)
# if the file timestamp does not exist, create one
if not os.path.exists(FILE_TIMESTAMP):
record_file_time(FILE_LIST_TO_CHECK)
# get update time of must-have files
timestamp_dict = read_file_timestamp(FILE_TIMESTAMP)
# update edges.txt, a merged file from several sources, HISTORY_DIR and HISTORY_DIR2.
edge_file_lst = [] # collect edge files.
for fname in glob.glob(os.path.join(HISTORY_DIR, 'edges.txt.*')): # edges.txt.* are to be merged
edge_file_lst.append(fname)
for fname in glob.glob(os.path.join(HISTORY_DIR2, 'edges.txt.*')): # edges.txt.* are to be merged
edge_file_lst.append(fname)
if edge_file_lst == []:
write_log_file('No files to merge. Run this script again a few hours later.', LOG_FILE)
sys.exit()
# merge edge files
write_log_file('Merge edge files ...', LOG_FILE)
merge_edges(edge_file_lst, EDGE_FILE, SAMPLE_SIZE_FILE) # merge individual files to EDGE_FILE, so EDGE_FILE is always updated. A new field, metric, is appended.
# delete edges if their supporting ChIP expriments become obsolete
#write_log_file('Delete edges ...', LOG_FILE)
para_c_dict = make_data_dict(PARAMETER_FOR_BUILDCMATRIX)
bad_chip_ids = get_bad_chip_ids(para_c_dict) # e.g., those marked with obsolete in parameter_for_buildCmatrix.txt
nonexistent_chip_ids = get_nonexistent_chip_ids(para_c_dict, EDGE_FILE) # in edges.txt but not in parameter_for_buildCmatrix.txt, either because it is removed or commented out
bad_chip_ids.extend(nonexistent_chip_ids)
bad_chip_ids = list(set(bad_chip_ids)) # unique elements
update_date_chip_ids = get_update_date_chip_ids(para_c_dict) # a dictionary, get the update date of each ChIP experiment with update in Note field.
rm_chip_ids_from_edge_file(EDGE_FILE, bad_chip_ids, update_date_chip_ids) # edges.txt is updated, with bad chip ids removed.
nlines = num_line(EDGE_FILE)
write_log_file('Number of total edges %d.' % (nlines), LOG_FILE)
if nlines == 0:
write_log_file('Empty edges.txt.', LOG_FILE)
sys.exit()
# get a list of updated files
# updated_file_list = get_updated_files(FILE_LIST_TO_CHECK, timestamp_dict)
updated_file_list = ['edges.txt']
# check edges.txt, if updated, re-make static html summary
if 'edges.txt' in updated_file_list: # if edges.txt is updated
write_log_file('Rebuild html pages ...', LOG_FILE)
cmd = 'python html_network.py -f %s -r %s -c %s -n %s' % (EDGE_FILE, PARAMETER_FOR_BUILDRMATRIX, PARAMETER_FOR_BUILDCMATRIX, PARAMETER_FOR_NET) # produce edges and summary.html
os.system(cmd)
# update webapp folder
cmd = 'cp %s ../Webapp/static/edges/' % (EDGE_FILE)
os.system(cmd)
# kill and restart start_webapp.py
# write_log_file('Terminate start_webapp ...', LOG_FILE)
# cmd = 'kill $(ps aux | grep \'[p]ython start_webapp\' | awk \'{print $2}\')'
# os.system(cmd)
# write_log_file('Restart start_webapp ...', LOG_FILE)
# os.chdir('../Webapp')
# os.system('python start_webapp.py &')
# os.chdir(CODE_DIR) # return to CODE_DIR
# update time stamp file
record_file_time(FILE_LIST_TO_CHECK)
# remove .R files in Data/temp, files older than 3 days will be removed
cmd = 'find %s -mtime +1 -name \"*.R\" -delete' % (TEMP_DIR)
os.system(cmd)
write_log_file('Network update done at %s.\n\n' % (datetime.now().strftime('%Y-%m-%d %H:%M:%S')), LOG_FILE)
print('Creating edges... Wait one hour to have edges ready.\nUncomment app.run(debug=True) in Webapp/start_webapp.py and comment out the previous app.run().\nTo display the network, cd Webpap && python start_webapp.py. Enter http://127.0.0.1:5000 in the address bar in your browser.')