#!/usr/bin/python

import argparse, os
import MySQLdb

# Required
import rpy2.robjects as robjects
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import math

# Import sonicLength package (Rpy2 2.1+)
from rpy2.robjects.packages import importr
sonicLength = importr("sonicLength")

from scipy.stats import gaussian_kde
from collections import Counter

# Uncomment for Rpy2 2.0
# from rpy2.robjects import r
# r.library("sonicLength")


header = """

+--------------------------------------------------+
 Author:	Stefano Brasca, Giulio Spinozzi
 Date:		October 2013
 Contact:	brasca.stefano, spinozzi.giulio @hsr.it
 Version:	2.0
+--------------------------------------------------+

 Description
  - This is a wrapper with rpy2 for simulate the Berry's model in R.
  
 Note:
  - The script needs the dataset file (LTR60.LC90.shearsites.uniq.txt), a valid db_schema and db_table

 Steps
	1. Loads the input file
	2. Calls the Berry's Model in R
	3. Saves the output of the model in Python structures
	4. Gives the output files with Phi and Theta
""" 

description = "This is a wrapper with rpy2 for simulate the Berry's model in R."

usage_example = """
Examples of usage:
APP --dataset LTR60.LC90.shearsites.uniq.txt --db_schema sequence_qlam --db_table osr_p5n6_p20
"""

# print header

parser = argparse.ArgumentParser(usage = usage_example, epilog = "[ hSR-TIGET - Vector Integration Core - Bioinformatics ] \n", description = description)
parser.add_argument('--dataset', dest="dataset_file", help="Dataset file to process in R. No default option.", action="store", required=True)
parser.add_argument('--db_schema', dest="db_schema", help="Database schema. No default option.", action="store", required=True)
parser.add_argument('--db_table', dest="db_table", help="Database table. No default option.", action="store", required=True)
args = parser.parse_args()



#########################################################################################
### FUNCTIONS
#########################################################################################

def checkArgs(args):
	"""
	Check file path
	"""
	if not os.path.isfile(args.dataset_file):
		print "\n[AP]\tError while reading files: no valid paths.\n\tExit\n"
		sys.exit()


def generateDataset(data):
	# Load dataset
	dataset_file = open(data, "r")
	dataset = dataset_file.readlines()
	dataset_file.close()
	locations_list = []
	length_list = []
	for row in dataset:
		row = row.rstrip('\n')
		row_split = row.split("\t")
		locations_list.append(row_split[0])
		length_list.append(row_split[1])
	return locations_list, length_list


def fragmentsLengthPlot(length_phi,freq_phi,length_list,nameFile,dataset):
	length_phi_numbers = []
	for num in length_phi:
		length_phi_numbers.append(float(num))
	plt.figure()
	# Phi plot
	plt.plot(length_phi_numbers, freq_phi, 'r', hold=True, label="estimated distribution")
	length_list_numbers = []
	for num in length_list:
		length_list_numbers.append(float(num))
	binning = int(math.ceil(len(set(length_list_numbers))/2))
	# Plot (length-frequency) of input data
	plt.hist(length_list_numbers, bins=binning, normed=True, facecolor='green', alpha=0.3, hold=True, label="real distribution - histogram")
	density = gaussian_kde(length_list_numbers)
	xs = np.linspace(0,max(length_list_numbers)+25,len(set(length_list_numbers))*10+250)
	# Gaussian kde plot
	plt.plot(xs,density(xs), hold=True, label="real distribution - gaussian kde")
	fileName = dataset + "." + nameFile + ".fragmentsLengthPlot.pdf"
   	plt.legend()
	# Labels
	plt.xlabel('fragments length')
	plt.ylabel('probability')
	plt.title(dataset + "." + nameFile + ' Fragments length data')
   	plt.savefig(fileName, format="pdf")
   	return length_phi_numbers


def printThetaInfo(estimations_theta,locations_theta,nameFile):
	f_out = open(nameFile+".theta.tsv","w")
	dict_theta = {}
	for l in range(0,len(locations_theta)):
		dict_theta[locations_theta[l]] = estimations_theta[l]
		f_out.write(str(locations_theta[l]) + '\t' + str(estimations_theta[l]) + '\n')


def queryDataset(host,user,passwd,db,db_table,destfile,nameFile):
	# System Call
	query = "mysql -h %(host)s -u %(user)s --password=%(passwd)s %(db)s --skip-column-names -e \"SELECT DISTINCT sample, tissue, treatment, vector, enzyme FROM %(db_table)s WHERE tag='%(nameFile)s'\" > %(destfile)s" %{
     'host': host,
     'user': user,
     'passwd': passwd,
     'db': db,
     'db_table': db_table,
     'nameFile': nameFile,
     'destfile': destfile,
    }
	os.system(query)
	f_in = open(destfile, "r")
	value = None
	for line in f_in:
		string_splitted = line.split('\t')
		value = str(string_splitted[0]) + "." + str(string_splitted[1]) + "." + str(string_splitted[2]) + "." + str(string_splitted[3]) + "." + str(string_splitted[4])
	f_in.close()
	os.remove(destfile)
	return value

def box_plot (real_data_list, estimated_data_list, nameFile,dataset):
	"""
	nameFile - string: dataset ID
	real_data_list - list of sequence count for 'nameFile'
	estimated_data_list - list of estimated abundance (theta) for 'nameFile'
	"""
	#normalizing data (Z-score)
	normalized_real_data_list = stats.zscore(real_data_list)
	normalized_estimated_data_list = stats.zscore(estimated_data_list)
	#Set-up data
	data = [normalized_real_data_list, normalized_estimated_data_list]
	#Creating Boxplot
	plt.figure()
	ax = plt.axes()
	plt.hold(True)
	bp = plt.boxplot(data)
	#Set-up boxplot apparence
	plt.setp(bp['boxes'][0], color='red')
	plt.setp(bp['caps'][0], color='red')
	plt.setp(bp['caps'][1], color='red')
	plt.setp(bp['whiskers'][0], color='red')
	plt.setp(bp['whiskers'][1], color='red')
	plt.setp(bp['fliers'][0], color='red')
	plt.setp(bp['fliers'][1], color='red')
	plt.setp(bp['medians'][0], color='red')
	plt.setp(bp['boxes'][1], color='green')
	plt.setp(bp['caps'][2], color='green')
	plt.setp(bp['caps'][3], color='green')
	plt.setp(bp['whiskers'][2], color='green')
	plt.setp(bp['whiskers'][3], color='green')
	plt.setp(bp['fliers'][2], color='green')
	plt.setp(bp['fliers'][3], color='green')
	plt.setp(bp['medians'][1], color='green')
	ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', alpha=0.5)
	ax.set_axisbelow(True)
	ax.set_title('Boxplot comparison - {0}\n- n={1} redundant IS -'.format(dataset + "." + nameFile, len(real_data_list)))
	ax.set_ylabel('Abundance / SeqCount (Z-score normalized)')
	# draw temporary green and red lines and use them to create a legend
	hR, = plt.plot([1,1],'r-') #estimated
	hG, = plt.plot([1,1],'g-') #real
	plt.legend((hR, hG),('1) Estimated Data - Abundance', '2) Real Data - SequenceCount'))
	hR.set_visible(False)
	hG.set_visible(False)
	#save figure and #show
	plt.savefig(dataset + "." + nameFile + '.comparative_boxplot' + '.pdf', format='pdf')
	#Return
	return 0


def expected_lengths_given_theta (length_phi, freq_phi, theta):
	expected_lenght = 0
	i=0
	for leng in length_phi:
		expected_lenght = expected_lenght + (1 - np.exp(-1*theta*freq_phi[i]))
		i+=1
	return expected_lenght


def phi_VS_theta (length_phi, freq_phi,nameFile,dataset):
	# Retrieving points to plot
	x_expected_unique_lenghts = []
	y_number_of_parent_fragments = []
	N = 1000
	for theta in np.linspace(0, 10000, N):
		y_number_of_parent_fragments.append(theta)
		x_expected_unique_lenghts.append(expected_lengths_given_theta(length_phi, freq_phi, theta))
	#Plot
	plt.figure()
	plt.plot(x_expected_unique_lenghts, y_number_of_parent_fragments, 'b', hold=True)
	plt.xlabel('expected unique lengths (phi-i)')
	plt.ylabel('number of parent fragments (theta-i)')
	plt.title('phi VS theta')
	plt.savefig(dataset + "." + nameFile + '.phiVStheta' + '.pdf', format='pdf')
	return 0


def from_file_to_list (unique_file_path, file_type):
	unique_file_path_split = unique_file_path.split('.uniq.txt')
	source_file_path = unique_file_path_split[0] + file_type
	dataset_file = open(source_file_path, "r")
	file_as_list = dataset_file.readlines()
	dataset_file.close()
	return file_as_list


def redundant_reads_count (file_as_list):
	
	list_of_reads =[]
	for line in file_as_list:
		line_split = line.split('\t')
		line_split[-1] = line_split[-1].rstrip('\n')
		list_of_reads.append(" ".join([line_split[1], line_split[2], line_split[5]]))
	list_of_reads.sort()

	list_of_redundant_reads_count = []
	dic_of_redundant_reads_count = {}
	i=0
	count=1
	for line in list_of_reads[1:]:
		i+=1
		if (line == list_of_reads[i-1]):
			count = count + 1
		else:
			dic_of_redundant_reads_count.update({list_of_reads[i-1]:count})
			count = 1
	dic_of_redundant_reads_count.update({list_of_reads[i]:count})


	#Convert dic to list
	keys = dic_of_redundant_reads_count.keys()
	keys.sort()

	for key in keys:
		list_of_redundant_reads_count.append(str(dic_of_redundant_reads_count[key]))

	return dic_of_redundant_reads_count, list_of_redundant_reads_count


def corrected_reads_count (dic_of_redundant_reads_count, dic_of_theta):

	# Get Keys
	red_keys = dic_of_redundant_reads_count.keys()
	#theta_keys = dic_of_theta.keys()

	# Convert items in numbers
	for key,value in dic_of_redundant_reads_count.iteritems():
		dic_of_redundant_reads_count[key] = float(value)
	for key,value in dic_of_theta.iteritems():
		dic_of_theta[key] = float(value)	

	# Get normalization factors
	total_redundant = sum(dic_of_redundant_reads_count.values())
	total_theta = sum(dic_of_theta.values())

	# dic_of_relative_abundance
	dic_of_relative_abundance = {}
	for key,value in dic_of_theta.iteritems():
		dic_of_relative_abundance.update({key:(value/total_theta)})

	# dic_of_corrected_reads_count
	dic_of_corrected_reads_count = {}
	for key,value in dic_of_relative_abundance.iteritems():
		dic_of_corrected_reads_count.update({key:(value*total_redundant)})

	# dic_of_percentage_difference
	dic_of_percentage_difference = {}
	for key in red_keys:
		dic_of_percentage_difference.update({key:(round(((dic_of_corrected_reads_count[key]-dic_of_redundant_reads_count[key])/dic_of_redundant_reads_count[key]*100),1))})

	return dic_of_relative_abundance, dic_of_corrected_reads_count, dic_of_percentage_difference


def fragment_lengths_statistics (unique_file_path):

	# Function MAD calculates Median Absolute Deviation
	def MAD (median, list_of_values):
		list_of_absolute_deviations = []
		for value in list_of_values:
			list_of_absolute_deviations.append(abs(median-value))
		mad = round(np.median(list_of_absolute_deviations),1)
		return mad

	#Acquiring data
	dataset_file = open(unique_file_path, "r")
	line_list = dataset_file.read().splitlines()

	#Filling dic_of_unique_lengths, dic_of_unique_lengths_number
	dic_of_unique_lengths = {}
	dic_of_unique_lengths_number = {}

	#For cycle, step 1
	previous_line_split = line_list[0].split('\t')
	previous_line_split[1] = int(previous_line_split[1])
	tmp_list_of_lenghts = [previous_line_split[1]]

	# Whole for cycle
	for line in line_list[1:]:
		line_split = line.split('\t')
		line_split[1] = int(line_split[1])
		if (line_split[0] == previous_line_split[0]):
			tmp_list_of_lenghts.append(line_split[1])
		else:
			dic_of_unique_lengths.update({previous_line_split[0]:tmp_list_of_lenghts})
			dic_of_unique_lengths_number.update({previous_line_split[0]:len(tmp_list_of_lenghts)})
			previous_line_split = line_split
			tmp_list_of_lenghts = [previous_line_split[1]]

	# Store last entry in dic_of_unique_lengths
	dic_of_unique_lengths.update({previous_line_split[0]:tmp_list_of_lenghts})
	dic_of_unique_lengths_number.update({previous_line_split[0]:len(tmp_list_of_lenghts)})

	#Filling dic_of_median_of_unique_lengths, dic_of_MAD
	dic_of_median_of_unique_lengths = {}
	dic_of_MAD = {}
	for key,value in dic_of_unique_lengths.iteritems():
		dic_of_median_of_unique_lengths.update({key:(round(np.median(value),1))})
		dic_of_MAD.update({key:(MAD(dic_of_median_of_unique_lengths[key], value))})

	return dic_of_unique_lengths, dic_of_unique_lengths_number, dic_of_median_of_unique_lengths, dic_of_MAD


def is_CEM (genome_location_string):
	response = False
	cem_symbol = None
	cem_coordinates = None
	cem_locations = [('chr11',64537165,64537171), ('chr17',2032349,2032355), ('chr17',47732336,47732342), ('chr2',24546568,24546574), ('chr2',73762395,73762401), ('chr16',28497495,28497501), ('chr3', 52306691, 52306697), ('chr8', 8866483, 8866489)]
	#cem_locations.append(('chr1',1776600,1776700))	#Fake to test: 1776639
	cem_symbols_list = ['CEM6_37_1','CEM6_37_2','CEM6_37_3','CEM6_37_4','CEM6_37_5','CEM6_37_6', 'CEM1_25', 'CEM1_41_illumina']
	cem_coordinates_list = ['chr11:64537168','chr17:2032352','chr17:47732339','chr2:24546571','chr2:73762398','chr16:28497498', 'chr3:52306694', 'chr8:8866486']
	
	splitted_location = genome_location_string.split(' ') # 0 -> chr; 1 -> locus; 2 -> strand
	splitted_location[1] = int(splitted_location[1])

	i=0;
	for cem_location_tupla in cem_locations:
		if (splitted_location[0] in cem_location_tupla[0]):
			if (splitted_location[1] in range(cem_location_tupla[1],cem_location_tupla[2]+1)):
				response = True
				cem_symbol = cem_symbols_list[i]
				cem_coordinates = cem_coordinates_list[i]
				return response, cem_symbol, cem_coordinates
		i+=1

	return response, cem_symbol, cem_coordinates


def lengths_explicit_list (file_as_list): # from_file_to_list (unique_file_path, '.txt')

	dic_of_lengths = {}
	line_split = file_as_list[0].split('\t')
	dic_of_lengths.update({line_split[0]:[int(line_split[1])]})
	
	for line in file_as_list[1:]:
		line_split = line.split('\t')
		if (dic_of_lengths.has_key(line_split[0])):
			dic_of_lengths[line_split[0]].append(int(line_split[1]))
		else:
			dic_of_lengths.update({line_split[0]:[int(line_split[1])]})

	for key, value in dic_of_lengths.iteritems():
		count = Counter(value)
		tmp_list = []
		for key_2 in sorted(count.keys()):
			tmp_list.append(int(count[key_2]))
		dic_of_lengths.update({key:tmp_list})

	return dic_of_lengths



#########################################################################################
### MAIN
#########################################################################################

def main():
	"""
	Main part of the program.
	"""
	#first check args and file paths
	checkArgs(args)
	
	data = args.dataset_file
	f_name = data.split(".")
	print "\n[AP]\t"+"######## "+f_name[0] + '.' + f_name[1]+" ########"
	print "\n[AP]\tChecked inputs, now acquiring data"

	host = "localhost"
	user = "readonly"
	passwd = "readonlypswd"
	db = args.db_schema
	db_table = args.db_table

	nameFile = data[0:-20]
	dataset = queryDataset(host,user,passwd,db,db_table,"tmpFile.txt",nameFile)
	if dataset is not None:
		dataset = dataset.rstrip('\n')
		dataset = dataset.replace("/","-")

		locations_list, length_list = generateDataset(data)

		if len(locations_list) < 2:
			print "\n[SKIP]\t{dataset} has only one unique line! Can't estimate anything.\n\tSKIP THIS FILE!\n".format(dataset=str(dataset))
			return 0

		# Alias for estAbund calling
		estAbund = sonicLength.estAbund

		# Call estAbund and store returned object in results
		results = estAbund(robjects.StrVector(locations_list), robjects.FloatVector(length_list))

		# Put estimation for theta in estimations_theta and associated locations in locations_theta; then organize data in dic_of_theta
		theta = results.rx2("theta")
		estimations_theta = tuple(theta)
		locations_theta = tuple(theta.names)
		# dic_of_theta
		dic_of_theta = {}
		for i in range(len(locations_theta)):
			dic_of_theta.update({locations_theta[i]:estimations_theta[i]})

		# Put different fragment lengths in length_phi and associated frequencies in freq_phi
		phi = results.rx2("phi")
		freq_phi = tuple(phi)
		length_phi = tuple(phi.names)

		length_phi_numbers = fragmentsLengthPlot(length_phi,freq_phi,length_list,nameFile,dataset)

		printThetaInfo(estimations_theta,locations_theta,nameFile)

		# Retrieving redundant reads data
		dic_of_redundant_reads_count, sequence_count_list = redundant_reads_count(from_file_to_list(data,'.tsv'))

		# Box Plot
		sequence_count = []
		for v in sequence_count_list:
			sequence_count.append(int(v))
		box_plot(sequence_count, estimations_theta, nameFile,dataset)

		# Plot: unique lengths retrieved for a genomic location VS expected number of parent fragment for the same location
		phi_VS_theta(length_phi, freq_phi, nameFile, dataset)


		#######################################################################################################
		# Produce .tsv output about measured redundant reads count, abundance-corrected redundant reads count # 
		# and some descriptive of unique fragments lengths                                                    #
		#######################################################################################################

		# Retrieving data
		dic_of_relative_abundance, dic_of_corrected_reads_count, dic_of_percentage_difference = corrected_reads_count(dic_of_redundant_reads_count, dic_of_theta)
		dic_of_unique_lengths, dic_of_unique_lengths_number, dic_of_median_of_unique_lengths, dic_of_MAD = fragment_lengths_statistics(data)
		dic_of_lengths = lengths_explicit_list(from_file_to_list(data,'.txt'))

		# Writing File
		corrected_file = open(dataset + "." + nameFile+".outcomes"+".tsv", 'w')
		corrected_file.write("Chromosome\tIntegration_locus\tStrand\tSequence_Count\tEstimated_Relative_Abundance\tCorrected_Sequence_Count\tPercentage_Variation\tNumber_of_fragments_of_unique_lengths\tLength_Min\tLength_Max\tLenght_Median\tRounded_Lenght_Median\tMAD\tUnique_Lengths_List\tUnique_Lengths_Amount\tCEM_region_?") ## ! NB ! ## \tCEM_region_?" has to remain the last!!!
		genome_locations = dic_of_redundant_reads_count.keys()
		genome_locations.sort()
		for key in genome_locations:
			splitted_location = key.split(' ')
			corrected_file.write("\n" + splitted_location[0] + "\t" + splitted_location[1] + "\t" + splitted_location[2] + "\t" + str(dic_of_redundant_reads_count[key]) + "\t" + str(round(dic_of_relative_abundance[key],5)) + "\t" + str(round(dic_of_corrected_reads_count[key],0)) + "\t" + str(dic_of_percentage_difference[key]) + "\t" + str(dic_of_unique_lengths_number[key]) + "\t" + str(min(dic_of_unique_lengths[key])) + "\t" + str(max(dic_of_unique_lengths[key])) + "\t" + str(dic_of_median_of_unique_lengths[key]) + "\t" + str(math.ceil(dic_of_median_of_unique_lengths[key]))+ "\t" + str(dic_of_MAD[key]) + "\t" + str(dic_of_unique_lengths[key]) + "\t" + str(dic_of_lengths[key]))
			response, cem_symbol, cem_coordinates = is_CEM(key)
			if (response == True):
				corrected_file.write("\t" + cem_symbol)

		# Write database file - Like corrected_file with more field appended in the end
		db_file = open(dataset + "." + nameFile+".db_file"+".tsv", 'w')
		genome_locations = dic_of_redundant_reads_count.keys()
		genome_locations.sort()
		dataset_split = dataset.split('.')
		dataset_label = '_'.join(dataset_split)
		for key in genome_locations:
			splitted_location = key.split(' ')
			db_file.write(splitted_location[0] + "\t" + splitted_location[1] + "\t" + splitted_location[2] + "\t" + str(dic_of_redundant_reads_count[key]) + "\t" + str(round(dic_of_relative_abundance[key],5)) + "\t" + str(round(dic_of_corrected_reads_count[key],0)) + "\t" + str(dic_of_percentage_difference[key]) + "\t" + str(dic_of_unique_lengths_number[key]) + "\t" + str(min(dic_of_unique_lengths[key])) + "\t" + str(max(dic_of_unique_lengths[key])) + "\t" + str(dic_of_median_of_unique_lengths[key]) + "\t" + str(math.ceil(dic_of_median_of_unique_lengths[key]))+ "\t" + str(dic_of_MAD[key]) + "\t" + str(dic_of_unique_lengths[key])[1:-1] + "\t" + str(dic_of_lengths[key])[1:-1] + "\t")
			db_file.write("\t".join(dataset_split) + "\t" + dataset_label + "\t")
			response, cem_symbol, cem_coordinates = is_CEM(key)
			if (response == True):
				db_file.write(cem_symbol + "\t" + cem_coordinates)
			else:
				db_file.write("\t")
			db_file.write("\n")

		db_file.close()

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

		# Last print for user
		print "\n[AP]\tTask Finished, closing.\n"
	else:
		print "\n[AP]\tThe dataset is not in the reference DB. Skipped.\n"

	return 0


# sentinel
if __name__ == "__main__":
    main()
e69a633d867ff2c27fef56c9a08bb6482171aedb
#!/usr/bin/python

import argparse, os
import MySQLdb

# Required
import rpy2.robjects as robjects
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import math

# Import sonicLength package (Rpy2 2.1+)
from rpy2.robjects.packages import importr
sonicLength = importr("sonicLength")

from scipy.stats import gaussian_kde
from collections import Counter

# Uncomment for Rpy2 2.0
# from rpy2.robjects import r
# r.library("sonicLength")


header = """

+--------------------------------------------------+
 Author:	Stefano Brasca, Giulio Spinozzi
 Date:		October 2013
 Contact:	brasca.stefano, spinozzi.giulio @hsr.it
 Version:	2.0
+--------------------------------------------------+

 Description
  - This is a wrapper with rpy2 for simulate the Berry's model in R.
  
 Note:
  - The script needs the dataset file (LTR60.LC90.shearsites.uniq.txt), a valid db_schema and db_table

 Steps
	1. Loads the input file
	2. Calls the Berry's Model in R
	3. Saves the output of the model in Python structures
	4. Gives the output files with Phi and Theta
""" 

description = "This is a wrapper with rpy2 for simulate the Berry's model in R."

usage_example = """
Examples of usage:
APP --dataset LTR60.LC90.shearsites.uniq.txt --db_schema sequence_qlam --db_table osr_p5n6_p20
"""

# print header

parser = argparse.ArgumentParser(usage = usage_example, epilog = "[ hSR-TIGET - Vector Integration Core - Bioinformatics ] \n", description = description)
parser.add_argument('--dataset', dest="dataset_file", help="Dataset file to process in R. No default option.", action="store", required=True)
parser.add_argument('--db_schema', dest="db_schema", help="Database schema. No default option.", action="store", required=True)
parser.add_argument('--db_table', dest="db_table", help="Database table. No default option.", action="store", required=True)
args = parser.parse_args()



#########################################################################################
### FUNCTIONS
#########################################################################################

def checkArgs(args):
	"""
	Check file path
	"""
	if not os.path.isfile(args.dataset_file):
		print "\n[AP]\tError while reading files: no valid paths.\n\tExit\n"
		sys.exit()


def generateDataset(data):
	# Load dataset
	dataset_file = open(data, "r")
	dataset = dataset_file.readlines()
	dataset_file.close()
	locations_list = []
	length_list = []
	for row in dataset:
		row = row.rstrip('\n')
		row_split = row.split("\t")
		locations_list.append(row_split[0])
		length_list.append(row_split[1])
	return locations_list, length_list


def fragmentsLengthPlot(length_phi,freq_phi,length_list,nameFile,dataset):
	length_phi_numbers = []
	for num in length_phi:
		length_phi_numbers.append(float(num))
	plt.figure()
	# Phi plot
	plt.plot(length_phi_numbers, freq_phi, 'r', hold=True, label="estimated distribution")
	length_list_numbers = []
	for num in length_list:
		length_list_numbers.append(float(num))
	binning = int(math.ceil(len(set(length_list_numbers))/2))
	# Plot (length-frequency) of input data
	plt.hist(length_list_numbers, bins=binning, normed=True, facecolor='green', alpha=0.3, hold=True, label="real distribution - histogram")
	density = gaussian_kde(length_list_numbers)
	xs = np.linspace(0,max(length_list_numbers)+25,len(set(length_list_numbers))*10+250)
	# Gaussian kde plot
	plt.plot(xs,density(xs), hold=True, label="real distribution - gaussian kde")
	fileName = dataset + "." + nameFile + ".fragmentsLengthPlot.pdf"
   	plt.legend()
	# Labels
	plt.xlabel('fragments length')
	plt.ylabel('probability')
	plt.title(dataset + "." + nameFile + ' Fragments length data')
   	plt.savefig(fileName, format="pdf")
   	return length_phi_numbers


def printThetaInfo(estimations_theta,locations_theta,nameFile):
	f_out = open(nameFile+".theta.tsv","w")
	dict_theta = {}
	for l in range(0,len(locations_theta)):
		dict_theta[locations_theta[l]] = estimations_theta[l]
		f_out.write(str(locations_theta[l]) + '\t' + str(estimations_theta[l]) + '\n')


def queryDataset(host,user,passwd,db,db_table,destfile,nameFile):
	# System Call
	query = "mysql -h %(host)s -u %(user)s --password=%(passwd)s %(db)s --skip-column-names -e \"SELECT DISTINCT sample, tissue, treatment, vector, enzyme FROM %(db_table)s WHERE tag='%(nameFile)s'\" > %(destfile)s" %{
     'host': host,
     'user': user,
     'passwd': passwd,
     'db': db,
     'db_table': db_table,
     'nameFile': nameFile,
     'destfile': destfile,
    }
	os.system(query)
	f_in = open(destfile, "r")
	value = None
	for line in f_in:
		string_splitted = line.split('\t')
		value = str(string_splitted[0]) + "." + str(string_splitted[1]) + "." + str(string_splitted[2]) + "." + str(string_splitted[3]) + "." + str(string_splitted[4])
	f_in.close()
	os.remove(destfile)
	return value

def box_plot (real_data_list, estimated_data_list, nameFile,dataset):
	"""
	nameFile - string: dataset ID
	real_data_list - list of sequence count for 'nameFile'
	estimated_data_list - list of estimated abundance (theta) for 'nameFile'
	"""
	#normalizing data (Z-score)
	normalized_real_data_list = stats.zscore(real_data_list)
	normalized_estimated_data_list = stats.zscore(estimated_data_list)
	#Set-up data
	data = [normalized_real_data_list, normalized_estimated_data_list]
	#Creating Boxplot
	plt.figure()
	ax = plt.axes()
	plt.hold(True)
	bp = plt.boxplot(data)
	#Set-up boxplot apparence
	plt.setp(bp['boxes'][0], color='red')
	plt.setp(bp['caps'][0], color='red')
	plt.setp(bp['caps'][1], color='red')
	plt.setp(bp['whiskers'][0], color='red')
	plt.setp(bp['whiskers'][1], color='red')
	plt.setp(bp['fliers'][0], color='red')
	plt.setp(bp['fliers'][1], color='red')
	plt.setp(bp['medians'][0], color='red')
	plt.setp(bp['boxes'][1], color='green')
	plt.setp(bp['caps'][2], color='green')
	plt.setp(bp['caps'][3], color='green')
	plt.setp(bp['whiskers'][2], color='green')
	plt.setp(bp['whiskers'][3], color='green')
	plt.setp(bp['fliers'][2], color='green')
	plt.setp(bp['fliers'][3], color='green')
	plt.setp(bp['medians'][1], color='green')
	ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', alpha=0.5)
	ax.set_axisbelow(True)
	ax.set_title('Boxplot comparison - {0}\n- n={1} redundant IS -'.format(dataset + "." + nameFile, len(real_data_list)))
	ax.set_ylabel('Abundance / SeqCount (Z-score normalized)')
	# draw temporary green and red lines and use them to create a legend
	hR, = plt.plot([1,1],'r-') #estimated
	hG, = plt.plot([1,1],'g-') #real
	plt.legend((hR, hG),('1) Estimated Data - Abundance', '2) Real Data - SequenceCount'))
	hR.set_visible(False)
	hG.set_visible(False)
	#save figure and #show
	plt.savefig(dataset + "." + nameFile + '.comparative_boxplot' + '.pdf', format='pdf')
	#Return
	return 0


def expected_lengths_given_theta (length_phi, freq_phi, theta):
	expected_lenght = 0
	i=0
	for leng in length_phi:
		expected_lenght = expected_lenght + (1 - np.exp(-1*theta*freq_phi[i]))
		i+=1
	return expected_lenght


def phi_VS_theta (length_phi, freq_phi,nameFile,dataset):
	# Retrieving points to plot
	x_expected_unique_lenghts = []
	y_number_of_parent_fragments = []
	N = 1000
	for theta in np.linspace(0, 10000, N):
		y_number_of_parent_fragments.append(theta)
		x_expected_unique_lenghts.append(expected_lengths_given_theta(length_phi, freq_phi, theta))
	#Plot
	plt.figure()
	plt.plot(x_expected_unique_lenghts, y_number_of_parent_fragments, 'b', hold=True)
	plt.xlabel('expected unique lengths (phi-i)')
	plt.ylabel('number of parent fragments (theta-i)')
	plt.title('phi VS theta')
	plt.savefig(dataset + "." + nameFile + '.phiVStheta' + '.pdf', format='pdf')
	return 0


def from_file_to_list (unique_file_path, file_type):
	unique_file_path_split = unique_file_path.split('.uniq.txt')
	source_file_path = unique_file_path_split[0] + file_type
	dataset_file = open(source_file_path, "r")
	file_as_list = dataset_file.readlines()
	dataset_file.close()
	return file_as_list


def redundant_reads_count (file_as_list):
	
	list_of_reads =[]
	for line in file_as_list:
		line_split = line.split('\t')
		line_split[-1] = line_split[-1].rstrip('\n')
		list_of_reads.append(" ".join([line_split[1], line_split[2], line_split[5]]))
	list_of_reads.sort()

	list_of_redundant_reads_count = []
	dic_of_redundant_reads_count = {}
	i=0
	count=1
	for line in list_of_reads[1:]:
		i+=1
		if (line == list_of_reads[i-1]):
			count = count + 1
		else:
			dic_of_redundant_reads_count.update({list_of_reads[i-1]:count})
			count = 1
	dic_of_redundant_reads_count.update({list_of_reads[i]:count})


	#Convert dic to list
	keys = dic_of_redundant_reads_count.keys()
	keys.sort()

	for key in keys:
		list_of_redundant_reads_count.append(str(dic_of_redundant_reads_count[key]))

	return dic_of_redundant_reads_count, list_of_redundant_reads_count


def corrected_reads_count (dic_of_redundant_reads_count, dic_of_theta):

	# Get Keys
	red_keys = dic_of_redundant_reads_count.keys()
	#theta_keys = dic_of_theta.keys()

	# Convert items in numbers
	for key,value in dic_of_redundant_reads_count.iteritems():
		dic_of_redundant_reads_count[key] = float(value)
	for key,value in dic_of_theta.iteritems():
		dic_of_theta[key] = float(value)	

	# Get normalization factors
	total_redundant = sum(dic_of_redundant_reads_count.values())
	total_theta = sum(dic_of_theta.values())

	# dic_of_relative_abundance
	dic_of_relative_abundance = {}
	for key,value in dic_of_theta.iteritems():
		dic_of_relative_abundance.update({key:(value/total_theta)})

	# dic_of_corrected_reads_count
	dic_of_corrected_reads_count = {}
	for key,value in dic_of_relative_abundance.iteritems():
		dic_of_corrected_reads_count.update({key:(value*total_redundant)})

	# dic_of_percentage_difference
	dic_of_percentage_difference = {}
	for key in red_keys:
		dic_of_percentage_difference.update({key:(round(((dic_of_corrected_reads_count[key]-dic_of_redundant_reads_count[key])/dic_of_redundant_reads_count[key]*100),1))})

	return dic_of_relative_abundance, dic_of_corrected_reads_count, dic_of_percentage_difference


def fragment_lengths_statistics (unique_file_path):

	# Function MAD calculates Median Absolute Deviation
	def MAD (median, list_of_values):
		list_of_absolute_deviations = []
		for value in list_of_values:
			list_of_absolute_deviations.append(abs(median-value))
		mad = round(np.median(list_of_absolute_deviations),1)
		return mad

	#Acquiring data
	dataset_file = open(unique_file_path, "r")
	line_list = dataset_file.read().splitlines()

	#Filling dic_of_unique_lengths, dic_of_unique_lengths_number
	dic_of_unique_lengths = {}
	dic_of_unique_lengths_number = {}

	#For cycle, step 1
	previous_line_split = line_list[0].split('\t')
	previous_line_split[1] = int(previous_line_split[1])
	tmp_list_of_lenghts = [previous_line_split[1]]

	# Whole for cycle
	for line in line_list[1:]:
		line_split = line.split('\t')
		line_split[1] = int(line_split[1])
		if (line_split[0] == previous_line_split[0]):
			tmp_list_of_lenghts.append(line_split[1])
		else:
			dic_of_unique_lengths.update({previous_line_split[0]:tmp_list_of_lenghts})
			dic_of_unique_lengths_number.update({previous_line_split[0]:len(tmp_list_of_lenghts)})
			previous_line_split = line_split
			tmp_list_of_lenghts = [previous_line_split[1]]

	# Store last entry in dic_of_unique_lengths
	dic_of_unique_lengths.update({previous_line_split[0]:tmp_list_of_lenghts})
	dic_of_unique_lengths_number.update({previous_line_split[0]:len(tmp_list_of_lenghts)})

	#Filling dic_of_median_of_unique_lengths, dic_of_MAD
	dic_of_median_of_unique_lengths = {}
	dic_of_MAD = {}
	for key,value in dic_of_unique_lengths.iteritems():
		dic_of_median_of_unique_lengths.update({key:(round(np.median(value),1))})
		dic_of_MAD.update({key:(MAD(dic_of_median_of_unique_lengths[key], value))})

	return dic_of_unique_lengths, dic_of_unique_lengths_number, dic_of_median_of_unique_lengths, dic_of_MAD


def is_CEM (genome_location_string):
	response = False
	cem_symbol = None
	cem_coordinates = None
	cem_locations = [('chr11',64537165,64537171), ('chr17',2032349,2032355), ('chr17',47732336,47732342), ('chr2',24546568,24546574), ('chr2',73762395,73762401), ('chr16',28497495,28497501), ('chr3', 52306691, 52306697), ('chr8', 8866483, 8866489)]
	#cem_locations.append(('chr1',1776600,1776700))	#Fake to test: 1776639
	cem_symbols_list = ['CEM6_37_1','CEM6_37_2','CEM6_37_3','CEM6_37_4','CEM6_37_5','CEM6_37_6', 'CEM1_25', 'CEM1_41_illumina']
	cem_coordinates_list = ['chr11:64537168','chr17:2032352','chr17:47732339','chr2:24546571','chr2:73762398','chr16:28497498', 'chr3:52306694', 'chr8:8866486']
	
	splitted_location = genome_location_string.split(' ') # 0 -> chr; 1 -> locus; 2 -> strand
	splitted_location[1] = int(splitted_location[1])

	i=0;
	for cem_location_tupla in cem_locations:
		if (splitted_location[0] in cem_location_tupla[0]):
			if (splitted_location[1] in range(cem_location_tupla[1],cem_location_tupla[2]+1)):
				response = True
				cem_symbol = cem_symbols_list[i]
				cem_coordinates = cem_coordinates_list[i]
				return response, cem_symbol, cem_coordinates
		i+=1

	return response, cem_symbol, cem_coordinates


def lengths_explicit_list (file_as_list): # from_file_to_list (unique_file_path, '.txt')

	dic_of_lengths = {}
	line_split = file_as_list[0].split('\t')
	dic_of_lengths.update({line_split[0]:[int(line_split[1])]})
	
	for line in file_as_list[1:]:
		line_split = line.split('\t')
		if (dic_of_lengths.has_key(line_split[0])):
			dic_of_lengths[line_split[0]].append(int(line_split[1]))
		else:
			dic_of_lengths.update({line_split[0]:[int(line_split[1])]})

	for key, value in dic_of_lengths.iteritems():
		count = Counter(value)
		tmp_list = []
		for key_2 in sorted(count.keys()):
			tmp_list.append(int(count[key_2]))
		dic_of_lengths.update({key:tmp_list})

	return dic_of_lengths



#########################################################################################
### MAIN
#########################################################################################

def main():
	"""
	Main part of the program.
	"""
	#first check args and file paths
	checkArgs(args)
	
	data = args.dataset_file
	f_name = data.split(".")
	print "\n[AP]\t"+"######## "+f_name[0] + '.' + f_name[1]+" ########"
	print "\n[AP]\tChecked inputs, now acquiring data"

	host = "localhost"
	user = "readonly"
	passwd = "readonlypswd"
	db = args.db_schema
	db_table = args.db_table

	nameFile = data[0:-20]
	dataset = queryDataset(host,user,passwd,db,db_table,"tmpFile.txt",nameFile)
	if dataset is not None:
		dataset = dataset.rstrip('\n')
		dataset = dataset.replace("/","-")

		locations_list, length_list = generateDataset(data)

		if len(locations_list) < 2:
			print "\n[SKIP]\t{dataset} has only one unique line! Can't estimate anything.\n\tSKIP THIS FILE!\n".format(dataset=str(dataset))
			return 0

		# Alias for estAbund calling
		estAbund = sonicLength.estAbund

		# Call estAbund and store returned object in results
		results = estAbund(robjects.StrVector(locations_list), robjects.FloatVector(length_list))

		# Put estimation for theta in estimations_theta and associated locations in locations_theta; then organize data in dic_of_theta
		theta = results.rx2("theta")
		estimations_theta = tuple(theta)
		locations_theta = tuple(theta.names)
		# dic_of_theta
		dic_of_theta = {}
		for i in range(len(locations_theta)):
			dic_of_theta.update({locations_theta[i]:estimations_theta[i]})

		# Put different fragment lengths in length_phi and associated frequencies in freq_phi
		phi = results.rx2("phi")
		freq_phi = tuple(phi)
		length_phi = tuple(phi.names)

		length_phi_numbers = fragmentsLengthPlot(length_phi,freq_phi,length_list,nameFile,dataset)

		printThetaInfo(estimations_theta,locations_theta,nameFile)

		# Retrieving redundant reads data
		dic_of_redundant_reads_count, sequence_count_list = redundant_reads_count(from_file_to_list(data,'.tsv'))

		# Box Plot
		sequence_count = []
		for v in sequence_count_list:
			sequence_count.append(int(v))
		box_plot(sequence_count, estimations_theta, nameFile,dataset)

		# Plot: unique lengths retrieved for a genomic location VS expected number of parent fragment for the same location
		phi_VS_theta(length_phi, freq_phi, nameFile, dataset)


		#######################################################################################################
		# Produce .tsv output about measured redundant reads count, abundance-corrected redundant reads count # 
		# and some descriptive of unique fragments lengths                                                    #
		#######################################################################################################

		# Retrieving data
		dic_of_relative_abundance, dic_of_corrected_reads_count, dic_of_percentage_difference = corrected_reads_count(dic_of_redundant_reads_count, dic_of_theta)
		dic_of_unique_lengths, dic_of_unique_lengths_number, dic_of_median_of_unique_lengths, dic_of_MAD = fragment_lengths_statistics(data)
		dic_of_lengths = lengths_explicit_list(from_file_to_list(data,'.txt'))

		# Writing File
		corrected_file = open(dataset + "." + nameFile+".outcomes"+".tsv", 'w')
		corrected_file.write("Chromosome\tIntegration_locus\tStrand\tSequence_Count\tEstimated_Relative_Abundance\tCorrected_Sequence_Count\tPercentage_Variation\tNumber_of_fragments_of_unique_lengths\tLength_Min\tLength_Max\tLenght_Median\tRounded_Lenght_Median\tMAD\tUnique_Lengths_List\tUnique_Lengths_Amount\tCEM_region_?") ## ! NB ! ## \tCEM_region_?" has to remain the last!!!
		genome_locations = dic_of_redundant_reads_count.keys()
		genome_locations.sort()
		for key in genome_locations:
			splitted_location = key.split(' ')
			corrected_file.write("\n" + splitted_location[0] + "\t" + splitted_location[1] + "\t" + splitted_location[2] + "\t" + str(dic_of_redundant_reads_count[key]) + "\t" + str(round(dic_of_relative_abundance[key],5)) + "\t" + str(round(dic_of_corrected_reads_count[key],0)) + "\t" + str(dic_of_percentage_difference[key]) + "\t" + str(dic_of_unique_lengths_number[key]) + "\t" + str(min(dic_of_unique_lengths[key])) + "\t" + str(max(dic_of_unique_lengths[key])) + "\t" + str(dic_of_median_of_unique_lengths[key]) + "\t" + str(math.ceil(dic_of_median_of_unique_lengths[key]))+ "\t" + str(dic_of_MAD[key]) + "\t" + str(dic_of_unique_lengths[key]) + "\t" + str(dic_of_lengths[key]))
			response, cem_symbol, cem_coordinates = is_CEM(key)
			if (response == True):
				corrected_file.write("\t" + cem_symbol)

		# Write database file - Like corrected_file with more field appended in the end
		db_file = open(dataset + "." + nameFile+".db_file"+".tsv", 'w')
		genome_locations = dic_of_redundant_reads_count.keys()
		genome_locations.sort()
		dataset_split = dataset.split('.')
		dataset_label = '_'.join(dataset_split)
		for key in genome_locations:
			splitted_location = key.split(' ')
			db_file.write(splitted_location[0] + "\t" + splitted_location[1] + "\t" + splitted_location[2] + "\t" + str(dic_of_redundant_reads_count[key]) + "\t" + str(round(dic_of_relative_abundance[key],5)) + "\t" + str(round(dic_of_corrected_reads_count[key],0)) + "\t" + str(dic_of_percentage_difference[key]) + "\t" + str(dic_of_unique_lengths_number[key]) + "\t" + str(min(dic_of_unique_lengths[key])) + "\t" + str(max(dic_of_unique_lengths[key])) + "\t" + str(dic_of_median_of_unique_lengths[key]) + "\t" + str(math.ceil(dic_of_median_of_unique_lengths[key]))+ "\t" + str(dic_of_MAD[key]) + "\t" + str(dic_of_unique_lengths[key])[1:-1] + "\t" + str(dic_of_lengths[key])[1:-1] + "\t")
			db_file.write("\t".join(dataset_split) + "\t" + dataset_label + "\t")
			response, cem_symbol, cem_coordinates = is_CEM(key)
			if (response == True):
				db_file.write(cem_symbol + "\t" + cem_coordinates)
			else:
				db_file.write("\t")
			db_file.write("\n")

		db_file.close()

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

		# Last print for user
		print "\n[AP]\tTask Finished, closing.\n"
	else:
		print "\n[AP]\tThe dataset is not in the reference DB. Skipped.\n"

	return 0


# sentinel
if __name__ == "__main__":
    main()