# 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 from log import write_log_file 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_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.')