path: root/Code/local_network.py

# 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.')