
## Author Qi Li (qli8@nd.edu)  ###
## If use this framework, please cite the following papers ##
# “Margolin, Adam A., etal. “ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context.” BMC bioinformatics. Vol. 7. No. 1. 
# Vân Anh Huynh-Thu, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. "Inferring regulatory networks from expression data using tree-based methods." PloS one 5.9 (2010).

#######################################################################################################

import numpy as np 
import pandas as pd 
import sys
import random
from scipy import stats as st 
import networkx as nx 
import os
import argparse
from function import *
from GENIE3 import *
from os import path
#######################################################################################################

parser = argparse.ArgumentParser()
parser.add_argument("-d", "--output directory", dest = 'd', type = str)
parser.add_argument("-k", "--number of samples", dest = 'k', type = int)
parser.add_argument("-b0", "--start bootstraps ID", dest = 'b0', type = int)
parser.add_argument("-b1", "--end bootstraps ID", dest = 'b1', type = int)
parser.add_argument("-f", "--input file", dest = "f", type = str )
args=parser.parse_args()

print (len(sys.argv))
if len(sys.argv) != 11:
	print ("ERROR MESSAGE: ")
	print ("comand should be: python RF_ARACNe.py -d <output directory> -k < number of columns you want to select each bootstrap[> -b0 <start bootstrap ID> -b1 <end bootstrap ID> -f <input gene expression file>")
	print ("example would be: python RF_ARACNe.py -d out-RF -k 183 -b0 1 -b1 100 -f ./examples/expression-data/Hu_data_renamed_nocontrol.txt")
	sys.exit()


### create output directory  ###
dir_path = os.getcwd()
out_dir = dir_path + '/' + args.d
if not os.path.exists(out_dir):
	os.makedirs(out_dir)


#######################################################################################################

## read in expression data
df = pd.read_csv(args.f, sep = "\t")

gene_ID = df.iloc[:, 0].tolist()
df = df.iloc[:, 1:]



#######################################################################################################
######## calculate the thresholds #########
###########################################

print ("Threshold process start")

## check if file exists:
if path.exists("./" + args.d +"/threshold_0.txt"):
	f = open("./" + args.d +"/threshold_0.txt", "r")
	thre = f.readline()
	thre = float(thre)
	print (thre)
else:
	## one random case 
	df_rand = Rand_EXP(df, axis=0)
	df_t = np.array(df_rand.T)
	edgemat_rand = GENIE3(df_t, nthreads = 4)
	MatrixToList(gene_ID, edgemat_rand,  args.d, 0, "thresholdraw")

	edgelist_rand = pd.read_csv("./" + args.d + "/thresholdraw_"+str(0)+".txt", header = None, sep = "\t")
	print (edgelist_rand.shape)
	print (edgelist_rand.columns)
	edgelist_rand = edgelist_rand.ix[(edgelist_rand[0] != edgelist_rand[1])]
	print (edgelist_rand.shape)

	threshold1(0, args.d, "threshold", edgelist_rand)

print ("Threshold process finished")

#######################################################################################################
########       100 bootstrap      #########
###########################################


for x in range(args.b0, args.b1):
	print ("bootstraping start")
	df_boot = Bootstrap(df, args.k, x)
	print ("bootstraping data ready: " + str(x))

	df_t = np.array(df_boot.T)
	edge_mat = GENIE3(df_t, nthreads = 4)

	start_time = time.time()

	MatrixToList(gene_ID, edge_mat, args.d, x, "bootstrapNetwork")

	edgelist = pd.read_csv("./" + args.d + "/bootstrapNetwork_"+str(x)+".txt", header = None, sep = "\t")
	edgelist.columns = ["V1", "V2", "E"]

	edgelist = edgelist.loc[edgelist['E'] > thre]
	print (edgelist.shape)

	edgelist_DPI = DPI(edgelist)
	edgelist_DPI.to_csv("./" + args.d + "/bootstrapNetwork_"+str(x)+".txt", index = False, sep = "\t")

print ("---%s seconds --" % (time.time() - start_time))










