From 97fdefab064f63642fa3ece05b807d29b459df31 Mon Sep 17 00:00:00 2001 From: Hui Lan Date: Wed, 4 Dec 2019 19:03:19 +0800 Subject: brain: add python and R code to local repository. --- Code/local_network.py | 1114 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1114 insertions(+) create mode 100644 Code/local_network.py (limited to 'Code/local_network.py') diff --git a/Code/local_network.py b/Code/local_network.py new file mode 100644 index 0000000..216832f --- /dev/null +++ b/Code/local_network.py @@ -0,0 +1,1114 @@ +# 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.') -- cgit v1.2.1