winnie_old.py

"""
WiNNie - Weighted motif eNrichmeNt analysis uses the weights received from TOBIAS in bigwig files to find enriched motifs
@author: Anastasiia Petrova
@contact: anastasiia.petrova(at)mpi-bn.mpg.de
"""

import argparse
import sys
import os
import re
import time
import multiprocessing
import logging
import subprocess
from Bio import SeqIO
import Bio.SeqIO.FastaIO as bio
import numpy as np
from collections import defaultdict
from scipy import stats
import pyBigWig
from statsmodels.sandbox.stats.multicomp import multipletests #for bonfferoni
import matplotlib.pyplot as plt
import random

import textwrap


import MOODS.scan
import MOODS.tools
import MOODS.parsers

logger = logging.getLogger('my_ame')
logger.setLevel(logging.INFO)

formatter = logging.Formatter("%(asctime)s : %(message)s", "%Y-%m-%d %H:%M")

fh = logging.FileHandler('my_ame_log.txt')
fh.setLevel(logging.INFO)
fh.setFormatter(formatter)
logger.addHandler(fh)

ch = logging.StreamHandler()
ch.setLevel(logging.INFO)
ch.setFormatter(formatter)
logger.addHandler(ch)

"""
#catch all the information about input and output files as well as information on the used tool (fimo or moods)
def parse_args():

	#formatter_class = argparse.RawDescriptionHelpFormatter(prog='winnie', description=textwrap.dedent('''\
	#	Please do not mess ip this text!
	#	--------------------------
	#	'''
	#	))
	
	parser = argparse.ArgumentParser(prog = 'winnie', description = textwrap.dedent('''
		                                  

		This script takes a list of motifs loaded from jaspar.genereg.net as a combined text file in MEME or .PFM format, a genome file in FASTA format and optionaly a .bed file (the one you want to be merged with the whole genome file) as input. If you want to merge a .bed file with the whole genome file, please enter --bed_file or -b bevor your .bed file. The tool will provide a file output_merge.fa, which you can also use for your research later on. If you already have a merged file, please give this one as genome file input. If there are several motifs in the input file, the tool will create a separate output file for each motif. Choose if you want to use fimo or moods with --use, this script uses by default fimo. Please note that the tool can not provide the calculation of q values with fimo due to the text mode that fimo needs to use. The tool sends merged genome file and motifs to fimo or moods, saves the sorted output for each of the given motifs as moods/fimo_output_[alternate name and id of the motif].txt in the output directory, then calculates the start and the end as real positions on the chromosom and writes this information in the ouput files. The columns in the output file are: chromosom, start, end, the name and score of TF. If a .bed file was given as input, the tool will also add the additional columns from it to the output. If the output file is empty, there were no machtes within given genome regions. Please note, if you want to have all intermediate output files, enter --clean nothing
		'''), epilog='That is what you need to make this script work for you. Enjoy it')
	
	#required_arguments = parser.add_argument_group('required arguments')
	#required_arguments.add_argument('-m', '--motifs', help='file in MEME format with mofits loaded from jaspar.genereg.net')
	#required_arguments.add_argument('-g', '--genome', help='a whole genome file or regions of interest in FASTA format to be scanned with motifs')
	#all other arguments are optional
	#parser.add_argument('-o', '--output_directory',  default='output', const='output', nargs='?', help='output directory, default ./output/')
	#parser.add_argument('-b', '--bed_file',  nargs='?', help='a .bed file to be merged with the whole genome file to find regions of interest')
	#parser.add_argument('--use', '--use_tool', default='fimo', const='fimo', nargs='?', choices=['fimo', 'moods'], help='choose the tool to work with, default tool is fimo')
	#parser.add_argument('--clean', nargs='*', choices=['nothing', 'all', 'cut_motifs', 'fimo_output', 'merge_output', 'moods_output'], dest='cleans', help='choose the files you want to delete from the output directory, the default is deleting all the temporary files from the directory')
	#parser.add_argument('--fimo', help='enter additional options for fimo using = inside "", for example fimo="--norc" to not score the reverse complement DNA strand. By default the --text mode is used and the calculation of the q values due to the --text mode is not possible')
	#parser.add_argument('--cores', type=int, help='number of cores allowed to use by this tool, by default the tool uses 2 cores', default=2)
	#parser.add_argument('-p', '--p_value', type=float, help='enter the p value, the default p value is 0.0001. Please note that if you enter the p value using --fimo="--thresh ..." as well, the one within --fimo call will be used', default=0.0001)
	#parser.add_argument('--resolve_overlaps', action='store_true', help='delete overlaps with greater p value, by default no overlaps are deleted')
	#parser.add_argument('--hide_info', action='store_true', help='while working with data write the information only into ./log.txt')
	#parser.add_argument('--moods_bg', nargs='+', type=float, help='set the bg for moods, by default moods uses the bg is 0.25 0.25 0.25 0.25')
	args = parser.parse_args()
	return args
"""

def check_directory(directory):
	if not os.path.exists(directory):
		os.makedirs(directory)
		logger.info('a new directory ' + directory + ' was created')
					
#merge the whole genome with the regions mentioned in .bed file
def merge(genome, bed_file, output_directory):
	print("entered merging")
	logger.info('the merging of files ' + genome + ' and ' + bed_file + ' will end soon, the result file is output_merged.fa')
	output_merge = os.path.join(output_directory, "output_merged.fa")
	#if os.path.isfile(output_merge):
	#	output_merge = os.path.join(output_directory, "custom_output_merge.fa")
	os.system("bedtools getfasta -fi " + genome + " -bed " + bed_file + " -fo " + output_merge)
	print("merging is finished")
	return output_merge

#split the motifs each in other file
def split_motifs(motifs, output_directory, usage):
	logger.info("the file with motifs " + motifs + " will be checked for motifs and if needed splitted in files each containing only one motif")
	
	first_line = subprocess.getoutput("head -1 " + motifs) #find the first line of the input file
	
	if usage == "moods":
		if first_line.startswith(">"):
			#the motif file probably has the .pfm format, try to read and split it
			splitted_motifs = read_pfm(motifs, output_directory)
		else: #maybe the file with motifs is in MEME format, so try to convert it
			logger.info("the input file has not the expected format, I will try to convert it to .pfm format")
			splitted_motifs = convert_meme_to_pfm(motifs, output_directory)

	elif usage == "fimo":
		if first_line.startswith("MEME version"):
			#the motifs file has probably the MEME format, try to read and split it
			splitted_motifs = read_meme(motifs, output_directory)
			
			#if the there was a convertion before, delete all the .pfm files as we don't need them
			for filename in os.listdir(output_directory):
				if filename.endswith(".pfm"):
					remove_file(os.path.join(output_directory, filename))

		else: #maybe the file with motifs is in .pfm format, so try to convert is
			logger.info("the input file has not the expected format, I will try to convert it to MEME format")
			splitted_motifs = convert_pfm_to_meme(motifs, output_directory)

	return splitted_motifs

def read_pfm(motifs, output_directory):
	splitted_motifs = [] #to save the names of files after splitting
	motif = [] #to save the motif itself, which will be written to the file

	with open(motifs) as read_file:
		lines = 0
		for line in read_file:
			#as the motif has first line with the name and 4 lines with information, if the 5th line is something else than the name of the next motif, the exit will be forced
			if lines == 5 and not line.startswith(">"):
				logger.info('please make sure that the file with motifs has a right format and the number of lines is right in the motif file')
				sys.exit()
			else:
				if line.startswith(">"):
					if 'written_file' in locals():
						written_file.write(''.join(motif))
						motif = []
						lines = 0
						written_file.close()

					motif_alternate_name = check_name(re.split(' ', line)[1].rstrip())
					motif_id = re.split(' ', line[1:])[0] #[1:] meands do not use the first character
					motif_name = os.path.join(output_directory, motif_alternate_name + '_' + motif_id + '.pfm')
							
					splitted_motifs.append(motif_name)
					written_file = open(motif_name, 'w')
						
			if lines >= 1 and lines <= 4: #one motif has 5 lines, the first consists the name, the next 4 - the information we need to proceed the data within moods
				motif.append(line)
					
			lines = lines + 1
	written_file.write(''.join(motif))
	written_file.close()

	return splitted_motifs

def read_meme(motifs, output_directory):
	splitted_motifs = [] #to save the names of files after splitting
	motif = [] #to save the motif itself, which will be written to the file
	head = [] #define a list for header, as fimo needs a header in each motif file it proceedes

	with open(motifs) as read_file:
		lines = 0
		for line in read_file:
			#make the head part
			if lines <= 8:
				if lines == 0 and not line.startswith("MEME version"):
					logger.info('please make sure that the file with motifs has a right format and the number of lines is right in the motif file')
					sys.exit()

				head.append(line)
			else:
				#search for motifs and save each to another file
				if line.startswith("MOTIF"):
					if 'written_file' in locals():
						written_file.write(''.join(motif))
						motif = []
						written_file.close()

					#the alternate name will be checked for validity and the invalid chars will be replaced with '_'
					if len(re.split(' ', line.rstrip())) == 3: #in the input motif file the motif id and the alternate name are splitted using the tab, otherwise they are splitted using _, but we do not want to change it if so
						motif_alternate_name = check_name(re.split(' ', line)[2].rstrip())
						motif_id = re.split(' ', line)[1]
						motif_name = os.path.join(output_directory, motif_alternate_name + '_' + motif_id + '.meme')
					else: 
						motif_alternate_name = check_name(re.split(' ', line)[1].rstrip())
						motif_name = os.path.join(output_directory, motif_alternate_name + '.meme')

					#make a list with all the motif names to know which files to iterate when fimo is called
					splitted_motifs.append(motif_name)

					written_file = open(motif_name, 'w')
					written_file.write(''.join(head))

				motif.append(line)

			lines = lines + 1
					
		#write the last motif
		written_file.write(''.join(motif))
		written_file.close()

	read_file.close()

	return splitted_motifs
	
def convert_meme_to_pfm(motifs, output_directory):
	#i can only convert the file to pfm if the motifs file is in MEME format

	splitted_motifs = [] #to save the names of files after splitting
	rows = [[] for row in range(4)]

	with open(motifs) as read_file:
		lines = 0
		for line in read_file:
			if lines == 0 and not line.startswith("MEME version"):
				logger.info('please make sure that the file with motifs has a right format and the number of lines is right in the motif file')
				sys.exit()
			else:
				#search for motifs and save each to another file
				if line.startswith("MOTIF"):
					
					if 'written_file' in locals():
						for row in rows:
							written_file.write('\t'.join(row) + '\n')

						rows = [[] for row in range(4)]

						written_file.close()
					
					#the alternate name will be checked for validity and the invalid chars will be replaced with '_'
					if len(re.split(' ', line.rstrip())) == 3: #in the input motif file the motif id and the alternate name are splitted using the tab, otherwise they are splitted using _, but we do not want to change it if so
						motif_alternate_name = check_name(re.split(' ', line)[2].rstrip())
						motif_id = re.split(' ', line)[1]
						motif_name = os.path.join(output_directory, motif_alternate_name + '_' + motif_id + '.pfm')

					else: 
						motif_alternate_name = check_name(re.split(' ', line)[1].rstrip())
						motif_name = os.path.join(output_directory, motif_alternate_name + '.pfm')
					
					#make a list with all the motif names to know which files to iterate when fimo is called
					splitted_motifs.append(motif_name)

					written_file = open(motif_name, 'w')
					
				elif line.startswith("letter-probability matrix"):
					columns = int(re.split(' ', re.split('w= ', line)[1])[0]) #find the number of columns from the line out of the motifs file
					nsites = int(re.split(' ', re.split('nsites= ', line)[1])[0]) #find the nsites to count the frequency count for .pfm file

				elif line.startswith(' '): #each line with information about frequency starts in MEME format with ' '
					for i in range(len(rows)):
						rows[i].append(str(round(float(re.findall(r'\S+', line)[i])*nsites))) #split the line, do not mention how much whitespaces are in between, multiply it with nsites and save it to the corresponding row

			lines = lines + 1
					
		#write the last motif
		for row in rows:
			written_file.write('\t'.join(row) + '\n')		

		written_file.close()
	read_file.close()

	return splitted_motifs

def convert_pfm_to_meme(motifs, output_directory):
	#i can only convert the file to meme, if motifs file is in .pfm format

	#first we need to split the pfm motifs as the jaspar2meme does not work with the files containing several motifs, but with the directory consisting files each with only one motif in pfm format
	pfm_motifs = read_pfm(motifs, output_directory)

	converted_name = os.path.join(output_directory, "converted_motifs.meme")

	os.system("jaspar2meme -pfm " + output_directory + " > " + converted_name)

	#need to call split motifs for meme file
	splitted_motifs = split_motifs(converted_name, output_directory, "fimo")

	remove_file(converted_name)
	return splitted_motifs

#if there are chars that are not allowed, they will be replaced with '_', to the possibly invalid names there will be added '_' at the beginning of the name
def check_name(name_to_test):
	badchars= re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$')
	badnames= re.compile(r'(aux|com[1-9]|con|lpt[1-9]|prn)(\.|$)')

	#replace all the chars that are not allowed with '_'
	name = badchars.sub('_', name_to_test)

	#check for the reserved by the os names
	if badnames.match(name):
		name = '_' + name
	return name

def call_moods(one_motif, genome, output_directory, p_value, moods_bg, fisher_dict_bg, fisher_dict, control, overexpression, control_dict, overexpression_dict, differences):

	#check if this is a bigwig file
	bw_control = pyBigWig.open(control)
	bw_overexpression = pyBigWig.open(overexpression)
	if not bw_control.isBigWig() or not bw_overexpression.isBigWig():
		logger.info("please provide the bigwig file!")
		sys.exit()
	else:
		# prepare everything for moods
		# setting standard parameters for moods
		pseudocount = 0.0001

		if moods_bg == None:
			bg = MOODS.tools.flat_bg(4)
		else:
			bg = tuple(moods_bg)

		logger.info("moods will work with the p_value " + str(p_value) + " and the bg " + str(bg))

		motif_name = os.path.basename(one_motif)

		matrix_names = [os.path.basename(one_motif)]

		matrices = []
		matrices_rc = []

		valid, matrix = pfm_to_log_odds(one_motif, bg, pseudocount)

		key_for_bed_dict = ''

		if valid:

			logger.info("please be patient, moods is working on the data")

			matrices.append(matrix)
			matrices_rc.append(MOODS.tools.reverse_complement(matrix,4))
			matrices_all = matrices + matrices_rc
			thresholds = [MOODS.tools.threshold_from_p(m, bg, p_value, 4) for m in matrices_all]

			scanner = MOODS.scan.Scanner(7)
			scanner.set_motifs(matrices_all, bg, thresholds)

			with open(genome) as handle:

				seq_iterator = bio.SimpleFastaParser(handle)

				for header, seq in seq_iterator:				

					header_splitted = re.split(r':', header)

					if len(header_splitted) == 1: #if there are no positions given
						header = header + ":0-" #set the first position as 0 and split it once more
						header_splitted = re.split(r':', header)
						logger.info("moods works with " + header)
					else: #the given genome file is a file with peaks, so use the header of the peak as a key to search in the bed dictionary for additional information later on
						key_for_bed_dict = header

						#------------------------delete?-------------------------------------
						if not key_for_bed_dict in fisher_dict_bg.keys():
							fisher_dict_bg[key_for_bed_dict] = 0
							fisher_dict[key_for_bed_dict] = 0
							#peak_length = int(re.split(r'-', header_splitted[1])[1]) - int(re.split(r'-', header_splitted[1])[0]) #find the length of the peak
							#fisher_dict[key_for_bed_dict] = np.zeros(peak_length) #initialize array of 0.0 with the length of the peak
						#----------------------------------------------------------------

					chromosom = header_splitted[0]
					positions = re.split(r'-', header_splitted[-1])

					#compute the background difference for this peak
					#bw_global_score_control = np.mean(np.nan_to_num(np.array(list(bw_control.values(chromosom, int(positions[0]), int(positions[1]))))))
					#bw_global_score_overexpression = np.mean(np.nan_to_num(np.array(list(bw_overexpression.values(chromosom, int(positions[0]), int(positions[1]))))))
					#
					#bw_global_difference = bw_global_score_control - bw_global_score_overexpression
					#print(bw_global_difference)

					results = scanner.scan(seq)

					fr = results[:len(matrix_names)] #forward strand
					rr = results[len(matrix_names):] #reverse strand

					results = [[(r.pos, r.score, '+', ()) for r in fr[i]] + 
						[(r.pos, r.score, '-', ()) for r in rr[i]] for i in range(len(matrix_names))] #use + and - to indicate strand

					for (matrix, matrix_name, result) in zip(matrices, matrix_names, results):

						motif_id = re.split(r'_', matrix_name)[-1].replace(".pfm", '') #find the id of the given morif
						motif_alternate_name = matrix_name.replace(motif_id, '')[:-1] #the alternate name of the motif is the name of the file without id and with cutted last character, that is _

						if len(matrix) == 4:
							l = len(matrix[0])
						if len(matrix) == 16:
							l = len(matrix[0] + 1)
						for r in sorted(result, key=lambda r: r[0]):
							strand = r[2]
							pos = r[0]
							hitseq = seq[pos:pos+l] #sequence
							score = format(r[1], '.15f') #round to 15 digits after floating point, already type str

							if key_for_bed_dict != '':
								start = pos + 1
								end = pos + len(hitseq)
								#chromosom = key_for_bed_dict #instead of only the name of chromosom write the key to search in the bed_file					
							else:
								start = int(positions[0]) + pos + 1
								end = start + len(hitseq) - 1

							"""# ---------------------------------------delete??-------------------------------
							if key_for_bed_dict not in scores_dict: #if this peak has no matches yet
								#fill the peak with 0
								peak_length = int(positions[1]) - int(positions[0]) + 1
								scores_dict[key_for_bed_dict] = np.full(peak_length, 0.0)
							"""#-------------------------------------------------------------------------------

							#find the real start and end positions on the chromosom		
							real_start = int(positions[0]) + int(start) #start of the peak + start of the motif within the peak, do not add 1, as bigwig is 0-based
							real_end = real_start + len(hitseq)

							#count this match to the fisher_dict_bg
							fisher_dict_bg[key_for_bed_dict] += 1

							#get the values from bw file 
							bw_scores_control = np.mean(np.nan_to_num(np.array(list(bw_control.values(chromosom, real_start, real_end)))))
							bw_scores_overexpression = np.mean(np.nan_to_num(np.array(list(bw_overexpression.values(chromosom, real_start, real_end)))))

							control_dict = save_bw_score(one_motif, key_for_bed_dict, control_dict, bw_scores_control, hitseq, score)
							overexpression_dict = save_bw_score(one_motif, key_for_bed_dict, overexpression_dict, bw_scores_overexpression, hitseq, score)

							bw_difference = abs(bw_scores_overexpression - bw_scores_control)
							
							if not np.isnan(bw_difference) and bw_difference != 0.0: #do not need to check for nan
								differences.append(bw_difference)

							"""#--------------------------------------------------delete????------------------------------
								if not key_for_bed_dict in matches_dict: 
								# and not start in matches_dict[key_for_bed_dict]: #and float(matches_dict[key_for_bed_dict][start]['difference']) < bw_difference: #if there is no information saved about this match yet and if so, the saved difference is smaller that the one we have just found
									matches_dict[key_for_bed_dict] = matches_dict.get(key_for_bed_dict, {})
									matches_dict[key_for_bed_dict][start] = {'moods_score': score, 'hitseq': hitseq, 'difference': bw_difference}
							"""

							"""
							for i in range(len(hitseq)):
								#check for the score, if there is already some score saved on this position, save the greatest (better) score
								if scores_dict[key_for_bed_dict][int(start) + i] != 0.0:
									#print(scores_dict[key_for_bed_dict][int(start) + i])
									scores_dict[key_for_bed_dict][int(start) + i] = max(bw_scores[i], scores_dict[key_for_bed_dict][int(start) + i])
									#print(scores_dict[key_for_bed_dict][int(start) + i])
								else:
									scores_dict[key_for_bed_dict][int(start) + i] = bw_scores[i]

							#if we already have background, use it to make the dictionary for fisher exact test
							if background_dict:
								#print(scores_dict[key_for_bed_dict][int(start):int(start) + len(hitseq)])
								overexpression_mean = np.mean(np.array(scores_dict[key_for_bed_dict][int(start):int(start) + len(hitseq)]))
								if not np.isnan(overexpression_mean) and overexpression_mean > background_dict[key_for_bed_dict]:
										motifs_dict[key_for_bed_dict] += 1
							else:
								motifs_dict[key_for_bed_dict] += 1
							"""#-----------------------------------------------------------------------------------------------------

			#one doesnt need to close file that was opened like so, as python does it on itself. file.closed says True			
			return fisher_dict_bg, fisher_dict, control_dict, overexpression_dict, differences
		else:
			logger.info("The input for moods was not validated by the MOODS.parsers.pfm. Please check if it has the right format (note that the MOODS accepts only the old version of .pfm files, that is one without the header containing the name and id of the motif)")
			sys.exit()

#check if the bw score is already saved, if so check if it is bigger than the new one <--------------------------- do i need to save motif name in this matrix???
def save_bw_score(motif, key_for_bed_dict, matches_dict, bw_score, hitseq, score):

	if np.isnan(bw_score): bw_score = 0.0

	#bw_score = float(bw_score) * float(score) #apply the moods score

	if key_for_bed_dict in matches_dict:

		""" #---------------------------delete or leave as a feature??-----------
		if float(matches_dict[key_for_bed_dict]['bw_score']) < bw_score: #save the match with the bigger score
			#print(matches_dict[key_for_bed_dict]['bw_score'])
			matches_dict[key_for_bed_dict] = {'bw_score': bw_score, 'hitseq': hitseq, 'moods_score': score}
			#print(matches_dict[key_for_bed_dict]['bw_score'])
		"""
		#save the mean of the boths scores
		#matches_dict[key_for_bed_dict] = np.mean([matches_dict[key_for_bed_dict], bw_score])
		#----------------------------------------------------------------------

		#save the biggest of scores
		if matches_dict[key_for_bed_dict] < bw_score:
			matches_dict[key_for_bed_dict] = bw_score
	else:
		matches_dict[key_for_bed_dict] = bw_score

	return matches_dict

#help function for the moods call, convert pfm to log odds
def pfm_to_log_odds(filename, bg, pseudocount):
	if pfm(filename):
		mat = MOODS.parsers.pfm_to_log_odds(filename, bg, pseudocount)
		if len(mat) != 4: #if something went wrong, the empty list will be returned
			return False, mat
		else:
			return True, mat
	else:
		logger.info('please make sure the motif file has a .pfm format needed for moods')
		sys.exit()

#help function for the moods call, check if the file is in a pfm format using moods
def pfm(filename):
	mat = MOODS.parsers.pfm(filename)
	if len(mat) != 4:
		return False
	else:
		return True

def remove_file(file):
	if os.path.isfile(file):
		os.remove(file)

def clean_directory(cleans, output_directory, motif, tool_output_file): #<-----------------------check if it is working

	moods_output_unsorted = os.path.join(tool_output_file.replace("moods_output", "moods_output_unsorted"))
	
	for clean in cleans:
		if clean == 'all':
			remove_file(motif)
			remove_file(tool_output_file)
			remove_file(moods_output_unsorted)
		elif clean == 'cut_motifs':
			remove_file(motif)
		elif clean == 'moods_output':
			if os.path.basename(tool_output_file).startswith("moods"):
				remove_file(tool_output_file)

def tool_make_output(usage, motif, genome, output_directory, cleans, p_value, bed_dictionary, fimo_data, resolve_overlaps, moods_bg, bw_overexpression, bw_control, mu_control, std_control, mu_overexpression, std_overexpression, global_mean, global_std):

	standard_moods_bg = 0.25, 0.25, 0.25, 0.25

	fisher_dict = defaultdict(int)	
	fisher_dict_bg = defaultdict(int)	

	control_dict = {}
	overexpression_dict = {}

	differences = []
	#--------------------do we need the  differences here?
	fisher_dict_bg, fisher_dict, control_dict, overexpression_dict, differences = call_moods(motif, genome, output_directory, p_value, standard_moods_bg, fisher_dict_bg, fisher_dict, bw_control, bw_overexpression, control_dict, overexpression_dict, differences) 
	
	#make arrays of the dictionaries
	control_array = []
	overexpression_array = []

	for key in control_dict:
		control_array.append(control_dict[key])
		overexpression_array.append(overexpression_dict[key])

	"""
	#find the mean of both arrays to know if the score is higher in control or in overexpression
	direction = "control"

	control_mean = np.mean(control_array)
	overexpression_mean = np.mean(overexpression_array)

	if control_mean - overexpression_mean <= 0:
		#the score of overexpression is higher
		direction = "overexpression"

	#print(direction)
	"""

	"""
	#normalize the arrays
	control_array = normalize(control_array)
	overexpression_array = normalize(overexpression_array)
	"""

	"""
	#first concatenate two arrays to normalize them at once
	one_array_len = len(control_array)
	all_scores = control_array + overexpression_array
	#scale all the scores together
	all_scores = scale(all_scores)
	#cut them in two again
	control_array = all_scores[:one_array_len]
	overexpression_array = all_scores[one_array_len:]
	"""

	motif_name = motif.replace(output_directory, '')
	motif_name = motif_name.replace('/', '')

	#make the wilcoxon signed-rank test
	my_statistic, my_wilcoxon_pvalue, direction, differences, motif_std = my_wilcoxon(control_array, overexpression_array, mu_control, std_control, mu_overexpression, std_overexpression, global_mean, global_std, motif_name, correction = False)

	positive_count_bg, negative_count_bg = count_motifs_number(fisher_dict_bg) #<--------------delete this one?

	# ------------------do we need this as extra features?
	#find the pearson correlation
	#r_row, pearson_pvalue = stats.pearsonr(control_array, overexpression_array)
	#print(r_row, pearson_pvalue)

	#compute the ranksum statistic
	#statistic2, ranksum_pvalue = stats.ranksums(control_array, overexpression_array)
	#print(statistic2, ranksum_pvalue)
	
	""" #--------------------do we need fisher test???
	positive_count, negative_count = count_motifs_number(fisher_dict)
	positive_count_bg, negative_count_bg = count_motifs_number(fisher_dict_bg)

	#find the p value of the motif
	oddsratio, fisher_p_value = stats.fisher_exact([[positive_count, negative_count],
										[positive_count_bg, negative_count_bg]]) #, alternative='greater')
	"""
	#oddsratio, fisher_p_value = stats.fisher_exact([[positive_count, positive_count_bg],
	#										[negative_count, negative_count_bg]])
	#correlation, fisher_p_value = stats.spearmanr([positive_count], [negative_count]) # liefert nan als pvalue :(

	
	return my_wilcoxon_pvalue, fisher_dict, direction, positive_count_bg, differences, motif_std

def make_normal_distribution_plot(input_array, motif_name):
	figure_name = motif_name.replace("pfm", "png")

	output_directory = "./plots3"

	mu = np.mean(input_array)
	std = np.std(input_array, ddof = 1)

	#axes = plt.gca() #gca = get the current axes
	#axes.set_ylim([0, 80])

	fig, ax = plt.subplots()

	plt.hist(input_array, bins = len(input_array), color = 'g', normed=True)
	xmin, xmax = plt.xlim()
	x = np.linspace(xmin, xmax, 100)
	p = stats.norm.pdf(x, mu, std)
	plt.plot(x, p, 'r', linewidth = 3)
	motif_name = motif_name.replace(".pfm", "")
	title = motif_name + " fit results: mu = %.4f, std = %.4f" % (mu, std)

	plt.title(title)
	plt.grid()

	plt.ylim(0, 50)

	fig.savefig(os.path.join(output_directory, figure_name))
	#plt.show()

#modified from scipy.wilcoxon
def my_wilcoxon(x, y, mu_control, std_control, mu_overexpression, std_overexpression, global_mean, global_std, motif_name, correction = False):

	x, y = map(np.asarray, (x, y)) #apply np.asarray to both input arrays

	if len(x) != len(y):
		raise ValueError("The length of both arrays in Wilcoxon test should be the same. Aborting")

	#normalize both arrays considering the global mu and std for overexpression and control
	
	#make_normal_distribution_plot(x, "control_no_normalization")
	#make_normal_distribution_plot(y, "overexpression_no_normalization")

	#x = (x - mu_control) / std_control
	#y = (y - mu_overexpression) / std_overexpression

	#make_normal_distribution_plot(x, "control_with_normalization")
	#make_normal_distribution_plot(y, "overexpression_with_normalization")

	d = x - y #find the difference

	#make_normal_distribution_plot(d, "differences_with_normalization")

	#keep all non-zero differences
	d = np.compress(np.not_equal(d, 0), d) #in scipy axis = -1, in my case it does not matter, as i have flattened array

	#make_normal_distribution_plot(d, motif_name)

	#correct the differences according to the normal distribution
	#mu = np.mean(d)
	#std = np.std(d, ddof = 1) #ddof = 1 calculates corrected sample sd which is sqrt(N/(N-1)) times the population sd where N is the number of points, interpretes the data as samples, estimates true variance
	#https://stackoverflow.com/questions/15389768/standard-deviation-of-a-list

	#correct the differences according to the global mean and std
	d = (d - global_mean) / global_std

	#motif_name = motif_name.replace(".pfm", "_normalized.pfm")
	#make_normal_distribution_plot(d, motif_name)

	count = len(d)

	if count < 10:
		logger.info("The sampe size is too small for normal approximation")

	r = stats.rankdata(abs(d)) #assign ranks to data, dealing with ties appropriately

	r_plus = np.sum((d > 0) * r, axis = 0)
	r_minus = np.sum((d < 0) * r, axis = 0)

	"""
	#direction towards first array, plus
	first_array = max(d)
	#direction towards second array, minus
	second_array = abs(min(d))

	if first_array > second_array:
		#x was bigger
		direction = "first_array"
	else:
		#y was bigger
		direction = "second_array"

	#print(r_plus, r_minus)
	"""

	T = min(r_plus, r_minus)
	mn = count * (count + 1.) * 0.25
	se = count * (count + 1.) * (2. * count + 1.)

	replist, repnum = stats.find_repeats(r)

	if repnum.size != 0:
		#correction for repeated elements
		se -= 0.5 * (repnum * (repnum * repnum -1)).sum()

	se = np.sqrt(se / 24)
	correction = 0.5 * int(bool(correction)) * np.sign(T - mn)
	z = (T - mn - correction) / se
	#scale down
	prob = 2. * stats.norm.sf(abs(z), scale = 2.5)

	motif_std = np.std(d, ddof = 1)
	motif_mu = np.mean(d)

	direction = "control" #first_array
	if motif_mu < 0:
		direction = "overexpression" #second_array

	return T, prob, direction, d, motif_std

def scale(input_array):
	output_array = np.interp(input_array, (min(input_array), max(input_array)), (0, 1))
	return output_array

def normalize(input_array):
	norm = np.linalg.norm(input_array)
	if norm == 0:
		return input_array
	return input_array / norm

def count_motifs_number(fisher_dict):
	positive_count = 0
	negative_count = 0

	for k, v in fisher_dict.items():
		if v > 0: #the number of found motifs in one peak
			positive_count += 1
		else: #the number of found motifs in this peak is 0
			negative_count += 1
	return positive_count, negative_count

def multiprocess(motifs, genome, output_directory, cleans, fimo_data, p_value, bed_dictionary, cpu_count, resolve_overlaps, usage, moods_bg, bw_overexpression, bw_control):
	print("entered multiprocess")

	#check_directory("./plots3")
	cpu_count = 30

	if cleans == None:
		cleans = ['all']

	pool = multiprocessing.Pool(cpu_count) #by default is cpu_count 2

	length = len(motifs) #the number of the motifs to find the percentage of the job that was done
	step = 100/length #the percentage that should be added after the job with each single motif is done

	#tasks = [] #here the jobs done by processes are saved

	matrix = {}
	fisher_p_values = []
	list_of_motifs = []
	fisher_dict = dict()

	#pick randomly 10000 scores of bigwig files and find the global score
	#global_mean, global_std = find_global_score(bw_overexpression, bw_control)
	global_mean, global_std, mu_control, std_control, mu_overexpression, std_overexpression = compute_differences(bed_dictionary, bw_overexpression, bw_control)
	#print(global_mean, global_std, mu_control, std_control, mu_overexpression, std_overexpression)
	#print(nastia)

	#find global difference: for each peak and the one for all peaks
	#global_differences, global_difference = compute_differences(bed_dictionary, bw_overexpression, bw_control)
	#print("global difference", global_difference)

	all_differences = []
	differences_file = open("test_difference12.txt", 'w')
	motifs_p_values = open("test_motifs_p_values12.txt", 'w')
	std_file = open("test_stds12.txt", 'w')

	all_stds = {}

	#first find motifs using standard background
	for motif in motifs:
		fisher_p_value, fisher_dict, direction, positive_count_bg, differences, motif_std = pool.apply_async(tool_make_output, args = (usage, motif, genome, output_directory, cleans, p_value, bed_dictionary, fimo_data, resolve_overlaps, moods_bg, bw_overexpression, bw_control, mu_control, std_control, mu_overexpression, std_overexpression, global_mean, global_std, )).get()

		motif_name = motif.replace(output_directory, "")
		motif_name = motif_name.replace("/", "")
		motif_name = motif_name.replace(".pfm", "")

		all_stds[motif_name] = motif_std
		all_differences.extend(differences)

		print(motif_name, fisher_p_value, direction, positive_count_bg)
		motifs_p_values.write('\t'.join([str(motif_name), str(fisher_p_value), direction, str(positive_count_bg)]) + '\n')
		differences_file.write(','.join(map(str, differences)) + ',')
		std_file.write('\t'.join([str(motif_name), str(motif_std)]) + '\n')

	#differences_file.write(','.join(map(str, (np.ndarray.tolist(global_differences)))))
	#differences_file.write(','.join(map(str, global_differences)))
	#differences_file.write("\n" + str(global_difference))

	#sort in descending order - greatest first
	#all_stds_sorted = sorted(all_stds, key = all_stds.get, reverse=True)

	#for key in all_stds_sorted:
	#	print(key)

	differences_file.close()
	motifs_p_values.close()
	std_file.close()

	#print(np.count_nonzero(all_differences), len(all_differences))

	#mu, std = stats.norm.fit(all_differences)
	mu = np.mean(all_differences)
	std = np.std(all_differences, ddof = 1) #ddof = 1 calculates corrected sample sd which is sqrt(N/(N-1)) times the population sd where N is the number of points, interpretes the data as samples, estimates true variance
	#https://stackoverflow.com/questions/15389768/standard-deviation-of-a-list
	plt.hist(all_differences, bins = 1500, color = 'g', normed=True)
	xmin, xmax = plt.xlim()
	x = np.linspace(xmin, xmax, 25000)
	#p = stats.norm.pdf(x, mu, std)
	p = stats.norm.pdf(x)
	plt.plot(x, p, 'r', linewidth = 3)
	title = "fit results: mu = %.4f, std = %.4f" % (mu, std)
	plt.title(title)
	plt.grid()
	plt.show()

	if mu != 0.0 and std != 1.0:
		all_differences_corrected = (all_differences - mu ) / std

		mu_corrected = np.mean(all_differences_corrected)
		std_corrected = np.std(all_differences_corrected, ddof = 1)
		plt.hist(all_differences_corrected, bins = 1500, color = 'b', normed=True)
		xmin, xmax = plt.xlim()
		x = np.linspace(xmin, xmax, 25000)
		#p = stats.norm.pdf(x, mu, std)
		p = stats.norm.pdf(x)
		plt.plot(x, p, 'r', linewidth = 3)
		title = "corrected results: mu = %.4f, std = %.4f" % (mu_corrected, std_corrected)
		plt.title(title)
		plt.grid()
		plt.show()




	"""
	tasks.append(pool.apply_async(tool_make_output, args = (usage, motif, genome, output_directory, cleans, p_value, bed_dictionary, fimo_data, resolve_overlaps, moods_bg, )))
	tasks_done = sum([task.ready() for task in tasks]) #the number of the processes that ended their job
	#check the number of the processes that are ready till the number of them reaches the number of processes started in the pool
	while tasks_done < len(tasks):
		#if the number of ready processes has changed, save the new number and print the percentage of the job done
		if sum([task.ready() for task in tasks]) != tasks_done:
			tasks_done = sum([task.ready() for task in tasks])
			sys.stdout.write("%-100s| %d%% \r" % (''*tasks_done, step*tasks_done))
			sys.stdout.flush()
			sys.stdout.write("\n")
		#update the number of ready processes each 0.05 seconds
		time.sleep(0.05)
	"""

	pool.close()
	pool.join() #make sure all the processes are done and exit
		
	#the processes should not delete the merged genome file
	#so make sure if this file is needed, otherwise delete it
	for clean in cleans:
		if clean == 'all' or clean == 'merge_output':
			for filename in os.listdir(output_directory):
				if filename == "output_merge.fa":
					remove_file(genome)

		if clean != 'nothing':
			logger.info('the directory ' + output_directory + ' was cleaned, only the required files were left')

def compute_differences(bed_dictionary, overexpression, control):
	print("entered compute_differences")

	bw_control = pyBigWig.open(control)
	bw_overexpression = pyBigWig.open(overexpression)

	global_differences = {} #dict
	differences_array = [] #to compute the mean at the end
	control_array = []
	overexpression_array = []

	for header in bed_dictionary:
		header_splitted = re.split(r':', header)
		chromosom = header_splitted[0]
		positions = re.split(r'-', header_splitted[-1])

		#compute the background difference for this peak
		bw_global_score_control = np.mean(np.nan_to_num(np.array(list(bw_control.values(chromosom, int(positions[0]), int(positions[1]))))))
		bw_global_score_overexpression = np.mean(np.nan_to_num(np.array(list(bw_overexpression.values(chromosom, int(positions[0]), int(positions[1]))))))
		bw_global_difference = bw_global_score_control - bw_global_score_overexpression

		global_differences[header] = bw_global_difference

		control_array.append(bw_global_score_control)
		overexpression_array.append(bw_global_score_overexpression)

		differences_array.append(bw_global_difference)

	mu = np.mean(differences_array)
	std = np.std(differences_array, ddof = 1)

	mu_control = np.mean(control_array)
	std_control = np.std(control_array, ddof = 1)

	mu_overexpression = np.mean(overexpression_array)
	std_overexpression = np.std(overexpression_array, ddof = 1)

	#make_normal_distribution_plot(control_array, "global_control_peaks")
	#make_normal_distribution_plot(overexpression_array, "global_overexpression_peaks")

	if (mu_control != mu_overexpression):
		print("not the same! The mean of overexpression_array is " + str(mu_overexpression) + ", the mean of control is " + str(mu_control))
		print("need to normalize!!!")
		if abs(mu_control) > mu_overexpression:
			print("control is bigger")
		else:
			print("overexpression is bigger")

	return mu, std, mu_control, std_control, mu_overexpression, std_overexpression

def find_global_score(overexpression, control):
	print("entered find_global_score")

	differences = []

	bw_control = pyBigWig.open(control)
	bw_overexpression = pyBigWig.open(overexpression)

	if not bw_control.isBigWig() or not bw_overexpression.isBigWig():
		logger.info("please provide the bigwig file!")
		sys.exit()

	chr_list = bw_control.chroms()

	i = 0
	while i < 10000:
		#find a random chromosome from the list
		chromosom = random.choice(list(chr_list))
		#find a random start
		start = random.randint(0, bw_control.chroms(chromosom) - 30)
		length = random.randint(7, 30) #we are considering that the smallest TF has length of 7,vf and the biggest - of 30 bp
		end = start + length

		bw_scores_control = np.mean(np.array(list(bw_control.values(chromosom, start, end))))
		bw_scores_overexpression = np.mean(np.array(list(bw_overexpression.values(chromosom, start, end))))

		bw_difference = bw_scores_control - bw_scores_overexpression

		if not np.isnan(bw_difference) and bw_difference != 0.0:
			differences.append(bw_difference)
			i += 1

	print(len(differences))
	difference = np.mean(differences)
	std = np.std(differences, ddof = 1)

	return difference, std

def is_fasta(check_fasta):
	if not os.path.isfile(check_fasta):
		logger.info('there is no file with genome, the exit is forced')
		sys.exit()
	else:
		# modified code from https://stackoverflow.com/questions/44293407/how-can-i-check-whether-a-given-file-is-fasta
	    with open(check_fasta, "r") as handle:
	        fasta = SeqIO.parse(handle, "fasta")
	        return any(fasta)  # False when `fasta` is empty, i.e. wasn't a FASTA file

def check_existing_input_files(args): #---------------------add check for bigwig files
	if not args.motifs or not args.genome:
		logger.info('there is no satisfied input, please enter --help or -h to learn how to use this tool')
		sys.exit()
	elif not is_fasta(args.genome):
		logger.info('please make sure the input genome file has a fasta format')
		sys.exit()
	#check if the file with motifs exists
	elif not os.path.isfile(args.motifs):
		logger.info('there is no file with motifs, the exit is forced')
		sys.exit()
	#if the bedfile was given as input, check if this file exists
	elif args.bed_file:
		if not os.path.isfile(args.bed_file):
			logger.info('there is no such bed file ' + args.bed_file + ', the exit is forced')
			sys.exit()

def find_window_length(bed_file):

	window_length = 0
	lengths = list()

	with open(bed_file) as read_bed_file:
		for bed_line in read_bed_file:
			bed_line_array = re.split(r'\t', bed_line.rstrip('\n'))
			if bed_line_array[1].isdigit() and bed_line_array[2].isdigit() and int(bed_line_array[1]) <= int(bed_line_array[2]): #in the real bedfile the second column is a start position, and the third column is an end position, so we are checking if these are integers and if the start position is smaller than the end one
				lengths.append(int(bed_line_array[2]) - int(bed_line_array[1])) #find the length of the peak
			else: #this is not a bed file, force exit
				logger.info('please make sure the input bed file has a right format, the problem occured on the line ' + bed_line)
				sys.exit()

	read_bed_file.close()

	window_length = int(np.mean(lengths))

	if window_length % 2 == 1: #if found window_length is an uneven number
		window_length += 1 #add 1 to make it even for the later division and work
	
	return window_length#-------------we don't need this function

def make_own_peaks(output_directory, bed_file, internal_window_length):
	custom_peaks_name = os.path.join(output_directory, "custom_peaks.bed")
	custom_peaks_file = open(custom_peaks_name, 'w')
	radius = internal_window_length / 2

	with open(bed_file) as read_bed_file:
		for bed_line in read_bed_file:
			bed_line_array = re.split(r'\t', bed_line.rstrip('\n'))
			if bed_line_array[1].isdigit() and bed_line_array[2].isdigit() and int(bed_line_array[1]) <= int(bed_line_array[2]): #in the real bedfile the second column is a start position, and the third column is an end position, so we are checking if these are integers and if the start position is smaller than the end one
				center = int((int(bed_line_array[2]) - int(bed_line_array[1])) / 2) + int(bed_line_array[1])
				custom_peaks_file.write('\t'.join([bed_line_array[0], str(int(center - radius)), str(int(center + radius))]) + '\n')

	custom_peaks_file.close()
	return custom_peaks_name#-------------we don't need this function

def make_bed_array(bed_file): #--------------------------do we need this function??

	bed_array = list()

	with open(bed_file) as read_bed_file:
		for bed_line in read_bed_file:
			bed_line_array = re.split(r'\t', bed_line.rstrip('\n'))
			if bed_line_array[1].isdigit() and bed_line_array[2].isdigit() and int(bed_line_array[1]) <= int(bed_line_array[2]): #in the real bedfile the second column is a start position, and the third column is an end position, so we are checking if these are integers and if the start position is smaller than the end one
				peak = bed_line_array[0] + ":" + bed_line_array[1] + "-" + bed_line_array[2]
				bed_array.append(peak)
			else: #this is not a bed file, force exit
				logger.info('please make sure the input bed file has a right format, the problem occured on the line ' + bed_line)
				sys.exit()

	read_bed_file.close()

	return bed_array

def make_bed_dictionary(bed_file):

	bed_dictionary = {}

	with open(bed_file) as read_bed_file:
		for bed_line in read_bed_file:
			bed_line_array = re.split(r'\t', bed_line.rstrip('\n'))
			if bed_line_array[1].isdigit() and bed_line_array[2].isdigit() and int(bed_line_array[1]) <= int(bed_line_array[2]): #in the real bedfile the second column is a start position, and the third column is an end position, so we are checking if these are integers and if the start position is smaller than the end one
				key = bed_line_array[0] + ":" + bed_line_array[1] + "-" + bed_line_array[2]
				value = []
				for i in range(3, len(bed_line_array)):
					value.append(bed_line_array[i]) 

				bed_dictionary[key] = value
			else: #this is not a bed file, force exit
				logger.info('please make sure the input bed file has a right format, the problem occured on the line ' + bed_line)
				sys.exit()

	read_bed_file.close()

	return bed_dictionary

def print_small_logo():
	winnie_logo = """\
	 __      ___ _  _ _  _ _     
	 \ \    / (_) \| | \| (_)___ 
	  \ \/\/ /| | .` | .` | / -_)
	   \_/\_/ |_|_|\_|_|\_|_\___|
	                             
	"""
	print(winnie_logo)

def print_big_logo():
	winnie_big_logo = """\


     
           ▄▒▀▀▀▒▄             ▄▄▀▀▒▄
        ▄▄▒▀       ▀▒▄▄     ▄▄▒▀       ▀▒▄
      ▐▒               ▐▒ ▒▀               ▒
      ▐▌               ▐▒ ▒                ▐▌
      ▐▌               ▐▒ ▒                ▐▌
      ▐▌               ▐▒ ▒          ▒█▌   ▐▌                      ▒█▌
      ▐▌        ▐█▄    ▐▒ ▒      █▀   ▀▀   ▐█▄   ▐█     █▄    █     ▀▀
        ▀▒▄▄      █▄▄▒▀ ▐█  ▀▒▄▄█      ▄▄▒▀▐██▄  ▐█    ▐██▌   █            ▄█▀▌
            ▀▒▄▄▄▒▀█▌▄▒▐█▀█▒▄▄██▀▄▄▄▄▐█    ▐█ ▀▌ ▐█    ▐█ ▀█  █     █     ▓█▄▄▐█▄
                 ▄▄▒▀█ █▌ █▌ █▀░▄▄   ▐█    ▐█  ▀█▄█    ▐█  ▀█▄█     █     ▀█▀▀▀▀
                ▐▌   ▀██   ██    ▐▌  ▐█    ▐█    ██    ▐█    ██     █       █▄
                ▐                ▐▌
                ▐                ▐▌
                ▐                ▐▌
                ▐▒▄             ▄▒▌
                   ▀▒▄▄     ▄▄▒▀
                       ▀▒▄▒▀
     
    

	"""
	print(winnie_big_logo)

def main():

	print_big_logo()

	#args = parse_args()
	logger.removeHandler(ch)

	#motifs = "small_database.meme"
	#genome = "small_hg19.fasta"
	#bed_file = "small_peaks.bed"

	#motifs = "../../PaperInPrep/TOBIAS/buenrostro_analysis/data/buenrostro_motifs.meme"
	#genome = "../analysis_my_tool/chipseq/hg19.fasta"
	#bed_file = "../analysis_my_tool/chipseq/hg19_peaks.bed"

	#bed_file = "../TOBIAS/tobias/TFBScan/peaks.bed"
	#output_directory = "./"

	

	#bw_overexpression = "overexpression_footprints.bw"
	#bw_control = "control_footprints.bw"

	#bw_overexpression = "../../mette.bentsen/to_anastasiia/duxbl_footprints.bw"
	#bw_control = "../../mette.bentsen/to_anastasiia/gfp_footprints.bw"

	bw_control = "../../mette.bentsen/to_anastasiia/duxbl_footprints.bw"
	bw_overexpression = "../../mette.bentsen/to_anastasiia/gfp_footprints.bw"

	#to check with ame
	motifs = "../my_ame/JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.meme"

	genome = "../my_ame/mm10.fa"
	#genome = "../my_ame/overexpression_peaks.fa"
	#genome = "/new_test/merged.fa"

	bed_file = "../../mette.bentsen/to_anastasiia/all_merged.bed"
	#bed_file = "../my_ame/overexpression_peaks.bed"

	output_directory_moods = "new_test"
	
	cleans = "nothing"
	fimo = ""
	p_value = 0.0001
	cores = 31
	resolve_overlaps = False
	moods_bg = 0.2955, 0.2045, 0.2045, 0.2955

	#check if there is an existing directory that user gave as input, otherwise create this directory from the path provided from the user
	check_directory(output_directory_moods)
	check_directory("./plots3")

	#merge the genome file with the normal peaks file
	genome = merge(genome, bed_file, output_directory_moods)

	#internal_window_length = find_window_length(bed_file)
	#print("internal window length: ", internal_window_length)
	#window_length = internal_window_length * 5

	#make_bed_array(bed_file)

	#make custom peaks file with all peaks of the same length
	#custom_peaks = make_own_peaks(output_directory, bed_file, internal_window_length)

	#merge thetool_make_output genome file with the custom peaks file
	#custom_genome = merge(genome, custom_peaks, output_directory)	
	
	splitted_motifs_moods = split_motifs(motifs, output_directory_moods, "moods")
	
	bed_dictionary = make_bed_dictionary(bed_file)

	moods_start = time.time()
	multiprocess(splitted_motifs_moods, genome, output_directory_moods, cleans, fimo, p_value, bed_dictionary, cores, resolve_overlaps, "moods", moods_bg, bw_overexpression, bw_control) # <---------------------- put the last as custom_genome

	print("----- MOODS needed %s seconds -----" % (time.time() - moods_start))
	
	for handler in logger.handlers:
		handler.close()
		logger.removeFilter(handler)
	
if __name__ == "__main__":
	main()