Repository: broadinstitute/genepy (language: python)

Issue:
async
some functions take a lot of time and one might want to not await theem in some cases (or run them in parallel) this can now be done using async/await


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 3f5355a..89a9e0e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -233,4 +233,7 @@ bh_unicode_properties.cache
 # https://packagecontrol.io/packages/sublime-github
 GitHub.sublime-settings
 
-# End of https://www.gitignore.io/api/r,macos,python,sublimetext
\ No newline at end of file
+# tmp files
+tmp.py
+
+# End of https://www.gitignore.io/api/r,macos,python,sublimetext
diff --git a/.gitmodules b/.gitmodules
index e69de29..ffaf141 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "genepy/epigenetics/rose"]
+	path = genepy/epigenetics/rose
+	url = https://github.com/jkobject/rose.git
diff --git a/README.md b/README.md
index 287444a..e8dac5b 100644
--- a/README.md
+++ b/README.md
@@ -1,14 +1,14 @@
-# GenePy
+# genepy
 
 _what is [genepy](https://en.wikipedia.org/wiki/G%C3%A9n%C3%A9pi)?_
 
-A set of awesome functions & tools for Computational Genomists
+A set of awesome functions & tools for Computational Geneticists
 
 ![long genome](documentation/genome.jpg)
 
 ## Content
 
-- **utils**: where a bunch of helper functions and usefull general scripts are stoed
+- **utils**: where a bunch of helper functions and usefull general scripts are stored
   - **plots**: a set of plotting tools based on [matplotlib]() and [bokeh]() to make volcano plots / CNV maps etc..
   - **helper**: and additional helper functions to save data, do merging of dataframes...
 - **terra**: contains a set of functions that uses [dalmatian]() to interact with the [GCP]() powered genomics HPC platform: [Terra](). 
@@ -18,67 +18,53 @@ A set of awesome functions & tools for Computational Genomists
 - **mutations**: a set of functions to work with maf files, vcf files etc..
 - **google**: functions and packages linked to google's apis
   - **google_sheet**: function to upload a df as a google sheet
-  - **gcp**: sets of functions to interact with google storage (relies on gsutil)
+  - **gcp**: sets of functions to interact with google storage (relies on `gsutil`)
 - **epigenetics**: where we have things related to epigenomics
-  - **chipseq**: has functions to read, merge, denoise, ChIP seq data, it contains a lot of functions required for the AML paper.
+  - **chipseq**: has functions to read, merge, denoise, ChIP seq data.
+  - **plot**: has functions to plot ChIP seq data.
 
 ### Helper tools
 
-_tools that you do not need to use directly as they have binding functions in GenePy._ 
+_tools that you do not need to use directly as they have binding functions in genepy._ 
 
 - **epigenetics/rose:**: where an updated version of the rose algorithm is stored (as a git submodule) 
-- **cell_line_mapping**: a set of functions to map cell line ids to other cell line ids based on an up to date google spreadsheet. 
+- **cell_line_mapping-master/python/cell_line_mapper**: a set of functions to map cell line ids to other cell line ids based on an up to date google spreadsheet. 
 
 
 ## Install
 
-### with pip (WIP)
+### with pip
 
 `pip install broad-genepy`
-### dev mode (better for now)
 
-```bash
-git clone git://github.com/BroadInstitute/GenePy.git
-cd GenePy
-git submodule update --init
-```
-
-then you can import files in python with e.g:
-```python
-from GenePy import TerraFunction as terra
-```
-
-if GenePy is not in your path, first do:
-
-```python
-import sys
-sys.path.append(RELATIVE_PATH_TO_GenePy)
-```
+and then use with `from genepy.utils/epigenetics/... import ...` 
 
-now you can install the necessary python packages:
+### dev mode
 
 ```bash
-pip install requirements.txt
-pip install rpy2-bioconductor-extensions gseapy macs2 deeptools
+git clone git://github.com/BroadInstitute/genepy.git
+pip install -e genepy
 ```
 
-or if not using the requirements.txt (computation results might change):
+then you can import files in python with e.g:
+```python
+from genepy import terra
+from genepy.utils import helper as h
+from genepy.google import gcp
+from genepy.utils import plot
+from genepy.epigenetics import chipseq
 
-```bash
-pip install numpy pandas
 ```
 
-```bash
-pip install bokeh dalmatian firecloud_dalmatian google_api_python_client gsheets gspread ipdb ipython matplotlib Pillow pybedtools pyBigWig pysam pytest requests rpy2 scikit_learn scipy seaborn setuptools taigapy taigapy typing venn rpy2-bioconductor-extensions gseapy macs2 deeptools
-```
+## installation: to get access to all bindings
 
-then install the following tools:
+Install the following tools:
 - [htslib/samtools](http://www.htslib.org/)
 - [bwa](https://github.com/lh3/bwa)
 just used once:
 - [bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)
 
-finaly you can install R packages (GSEABase, erccdashboard, GSVA, DESeq2):
+Finaly you can install R packages (GSEABase, erccdashboard, GSVA, DESeq2):
 
 ```bash
 R -e 'if(!requireNamespace("BiocManager", quietly = TRUE)){install.packages("BiocManager")};BiocManager::install(c("GSEABase", "erccdashboard", "GSVA", "DESeq2"));'
@@ -89,6 +75,6 @@ please do contribute, we do not have time to fix all issues or work on feature r
 
 Jeremie Kalfon jkalfon@broadinstitute.org jkobject@gmail.com https://jkobject.com
 
-
+Javad Noorbakhsh jnoorbak@broadinstitute.org
 
 Apache license 2.0.
diff --git a/epigenetics/README.md b/epigenetics/README.md
deleted file mode 100644
index 769ea1b..0000000
--- a/epigenetics/README.md
+++ /dev/null
@@ -1,47 +0,0 @@
-# epigenomics
-
-Especially targeted to functions related to the analysis of epigenomics data. It has functions to read, merge, denoise, ChIP seq data.
-
-## Available functions:
-
-_in ./chipseq.py_
-
-- bigWigFrom
-- ReadRoseSuperEnhancers
-- loadPeaks
-- pysam_getPeaksAt
-- bedtools_getPeaksAt
-- makeProfiles
-- getPeaksAt
-- computeMeanCov
-- substractPeaks
-- simpleMergePeaks
-- findpeakpath
-- findBestPeak
-- mergeReplicatePeaks
-- col_nan_scatter
-- col_nan_kde_histo
-- findAdditionalPeaks
-- negLogLikelihood
-- poissonFit
-- putInBed
-- pairwiseOverlap
-- enrichment
-- findAdditionalCobindingSignal
-- fullDiffPeak
-- diffPeak
-- MakeSuperEnhancers
-- runChromHMM
-- loadMEMEmotifs
-- simpleMergeMotifs
-- substractPeaksTo
-
-## highly recommended packages
-
-*This package won't contain anything that overlap with those and might use those packages for what it is doing.*
-- Bedtools
-- deepTools
-- MACS2
-- ROSE
-- MEME
-- ChromHMM
diff --git a/epigenetics/chipseq.py b/epigenetics/chipseq.py
deleted file mode 100644
index 096f710..0000000
--- a/epigenetics/chipseq.py
+++ /dev/null
@@ -1,1431 +0,0 @@
-import os
-import pandas as pd
-import pysam
-import sys
-import numpy as np
-from GenePy.utils import helper as h
-from GenePy.utils import plot
-import re
-from pybedtools import BedTool
-import seaborn as sns
-import pyBigWig
-import matplotlib.pyplot as plt
-import pdb
-import signal
-from scipy.optimize import curve_fit,minimize
-from sklearn.preprocessing import normalize
-from scipy.stats import poisson, zscore, boxcox, binom, fisher_exact
-from scipy.special import factorial
-import subprocess
-from pandas.io.parsers import ParserError, EmptyDataError
-import warnings
-import itertools
-from statsmodels.stats.multitest import multipletests
-
-size = {"GRCh37": 2864785220,
-		"GRCh38": 2913022398}
-
-cmaps = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
-		 'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
-		 'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']
-
-chroms = {'chr1','chr10','chr11','chr12','chr13','chr14','chr15','chr16','chr17','chr18','chr19',
-'chr2','chr20','chr21','chr22','chr3','chr4','chr5','chr6','chr7','chr8','chr9','chrX',
- 'chrY','1','10','11','12','13','14','15','16','17','18','19','2','20','21','22','3','4','5','6',
- '7','8','9','X','Y'}
-
-
-def bigWigFrom(bams, folder="", numthreads=8, genome='GRCh37', scaling=None, verbose=1):
-	"""
-	run the bigwig command line for a set of bam files in a folder
-	"""
-	if "bigwig" not in os.listdir(folder if folder else '.'):
-		os.mkdir(folder + "bigwig")
-	for i, bam in enumerate(bams):
-		in1 = folder + bam
-		out1 = folder + "bigwig/" + bam.split('/')[-1].split('.')[0] + '.bw'
-		cmd = "bamCoverage --effectiveGenomeSize " + str(size[genome]) + " -p " + str(numthreads) +\
-			" -b " + in1 + " -o " + out1
-		if scaling is not None:
-			cmd += ' --scaleFactor ' + str(scaling[i])
-		if verbose == 0:
-			cmd += ' 2> ' + bam + '.error.log'
-		res = subprocess.run(cmd, capture_output=True, shell=True)
-		if res.returncode != 0:
-			raise ValueError('issue with the command: ' + str(res.stderr))
-
-
-def ReadRoseSuperEnhancers(roseFolder, containsINPUT=True, name="MV411"):
-	beds = os.listdir(roseFolder)
-	superenhan = pd.DataFrame()
-	for i in beds:
-		if i.endswith('_SuperEnhancers.table.txt'):
-			superenhan = superenhan.append(pd.read_csv(roseFolder+i, sep='\t', skiprows=5)).drop(columns=['enhancerRank','isSuper'])
-	data = [superenhan.columns[6]] + superenhan.columns[8:].tolist()
-	if containsINPUT:
-		inputd = superenhan.columns[7]
-	superenhan['name'] = [i.split(name)[1].split('_')[2] for i in superenhan['REGION_ID']]
-	superenhan = superenhan.rename(columns={'CONSTITUENT_SIZE':'size','REGION_ID':"id",'CHROM':'chrom','START':'start','STOP':'end',"NUM_LOCI":'num'}).replace(np.nan,0)
-	superenhan['foldchange'] = superenhan[data].sum(1)/superenhan[inputd]
-	superenhan = superenhan.drop(columns=data+[inputd])
-	return superenhan.sort_values(by=['chrom','start','end'])
-
-
-def loadPeaks(peakFile=None,peakfolder=None, isMacs=True, CTFlist=[], skiprows=0):
-	# for each data peak type file listed in a MACS2 way in a given MACS2 output folder, will merge them
-	# all into one dataframe and output the dataframe
-	if peakfolder:
-		bindings = pd.DataFrame()
-		for folder in os.listdir(peakfolder):
-			if isMacs:
-				if any(tf in folder for tf in CTFlist) or not CTFlist:
-					binding = pd.read_csv(peakfolder + folder + "/NA_peaks.narrowPeak", sep='\t', header=None) if\
-					os.exists(peakfolder + folder + "/NA_peaks.narrowPeak") else peakfolder + folder + "/NA_peaks.broadPeak"
-					binding['name'] = folder.replace('.narrowPeak', '').replace('.broadPeak','')
-					bindings = bindings.append(binding)
-			else:
-				file = folder
-				if file[-10:] in ["narrowPeak",".broadPeak"] and (any(tf in file for tf in CTFlist) or not CTFlist):
-					print('reading: '+file)
-					binding = pd.read_csv(peakfolder + file, sep='\t', header=None, skiprows=skiprows)
-					binding['name'] = file.replace('.narrowPeak', '').replace('.broadPeak','')
-					bindings = bindings.append(binding)
-	elif peakFile:
-		bindings = pd.read_csv(peakFile, sep='\t', header=None, skiprows=skiprows)
-		bindings['name'] = peakFile.split('/')[-1].split('.')[0]
-	else:
-		raise ValueError("need to provide one of peakFile or peakfolder")
-	bindings = bindings.drop(5, 1).drop(4, 1)
-	bindings = bindings.rename(columns={
-		0: "chrom",
-		1: 'start',
-		2: 'end',
-		3: 'peak_number',
-		6: "foldchange",
-		7: "-log10pvalue",
-		8: "-log10qvalue",
-		9: 'relative_summit_pos'})
-	bindings = bindings.sort_values(by=["chrom", "start", "end"], axis=0)
-	bindings.start = bindings.start.astype('int')
-	bindings.end = bindings.end.astype('int')
-	bindings['relative_summit_pos'] = bindings.relative_summit_pos.astype(
-		float) if 'relative_summit_pos' in bindings.columns else bindings.end - bindings.start
-	bindings.foldchange = bindings.foldchange.astype('float')
-	bindings["-log10pvalue"] = bindings["-log10pvalue"].astype('float')
-	bindings['-log10qvalue'] = bindings['-log10qvalue'].astype('float')
-	bindings = bindings.reset_index(drop=True)
-	loc = bindings['relative_summit_pos'].isna()
-	bindings.loc[bindings[loc].index, 'relative_summit_pos'] = bindings[loc].end - bindings[loc].start
-	bindings.relative_summit_pos = bindings.relative_summit_pos.astype(int)
-	return bindings
-
-
-def pysam_getPeaksAt(peaks, bams, folder='data/seqs/', window=1000, numpeaks=1000, numthreads=8):
-
-	# get pysam data
-	# ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
-	# for each counts, do a rolling average (or a convolving of the data) with numpy
-	# append to an array
-	# return array, normalized
-	loaded = {}
-	res = {i: np.zeros((len(peaks), window * 2)) for i in bams}
-	peaks = peaks.sort_values(by="foldchange", ascending=False).iloc[:numpeaks]
-	peaks.chrom = peaks.chrom.astype(str)
-	for val in bams:
-		loaded.update({val: pysam.AlignmentFile(folder + val, 'rb', threads=numthreads)})
-	for k, bam in loaded.items():
-		for num, (i, val) in enumerate(peaks.iterrows()):
-			print(int(num / len(peaks)), end='\r')
-			center = int((val['start'] + val['end']) / 2)
-			for pileupcolumn in bam.pileup(val['chrom'], start=center - window,
-										   stop=center + window, truncate=True):
-				res[k][num][pileupcolumn.pos - (center - window)] = pileupcolumn.n
-	fig, ax = plt.subplots(1, len(res))
-	for i, (k, val) in enumerate(res.items()):
-		sns.heatmap(val, ax=ax[i])
-		ax[i].set_title(k.split('.')[0])
-	fig.show()
-	return res, fig
-
-
-def bedtools_getPeaksAt(peaks, bams, folder='data/seqs/', window=1000, numpeaks=1000, numthreads=8):
-	"""
-	get pysam data
-	ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
-	for each counts, do a rolling average (or a convolving of the data) with numpy
-	append to an array
-	return array, normalized
-	"""
-	loaded = {}
-	center = [int((val['start'] + val['end']) / 2) for k, val in peaks.iterrows()]
-	peaks['start'] = [c - window for c in center]
-	peaks['end'] = [c + window - 1 for c in center]
-	peaks[peaks.columns[:3]].sort_values(by=['chrom', 'start']).to_csv('temp/peaks.bed', sep='\t', index=False, header=False)
-	bedpeaks = BedTool('temp/peaks.bed')
-
-	fig, ax = plt.subplots(1, len(bams))
-	peakset = peaks["foldchange"].values.argsort()[::-1][:numpeaks]
-	for i, val in enumerate(bams):
-		coverage = BedTool(folder + val).intersect(bedpeaks).genome_coverage(bga=True, split=True)\
-			.intersect(bedpeaks).to_dataframe(names=['chrom', 'start', 'end', 'coverage'])
-		cov = np.zeros((len(peaks), window * 2), dtype=int)
-		j = 0
-		pdb.set_trace()
-		for i, (k, val) in enumerate(peaks.iterrows()):
-			print(i / len(peaks), end='\r')
-			while coverage.iloc[j].start > val.start:
-				j -= 1
-			while coverage.iloc[j].start < val.end:
-				cov[i][coverage.iloc[j].start - val.start:coverage.iloc[j].end - val.start] =\
-					coverage.iloc[j].coverage
-				j += 1
-		sns.heatmap(coverage, ax=ax[i])
-		ax[i].set_title(val.split('.')[0])
-	fig.show()
-	return None, fig
-
-
-def makeProfiles(matx=[], folder='', matnames=[], title='',
-				   name='temp/peaksat.pdf', refpoint="TSS", scale=None,
-				   sort=False, withDeeptools=True, cluster=1, vmax=None, vmin=None, overlap=False,
-				   legendLoc=None):
-	if withDeeptools:
-		if not (len(matnames) == 2 and len(matx) == 2):
-			raise ValueError('you need two mat.gz files and two names')
-		h.createFoldersFor(name)
-		cmd = 'computeMatrixOperations relabel -m '
-		cmd += matx[0] + ' -o '+matx[0]+' --groupLabels '+matnames[0]
-		cmd += ' && computeMatrixOperations relabel -m '
-		cmd += matx[1] + ' -o '+matx[1]+' --groupLabels '+matnames[1]
-		cmd += ' && computeMatrixOperations rbind -m '
-		cmd += matx[0] + ' ' + matx[1] + " -o " + '.'.join(name.split('.')[:-1]) + ".gz"
-		cmd += ' && plotProfile'
-		cmd += " --matrixFile " + '.'.join(name.split('.')[:-1]) + ".gz"
-		cmd += " --outFileName " + name
-		cmd += " --refPointLabel " + refpoint
-		if vmax is not None:
-			cmd += " -max "+str(vmax)
-		if vmin is not None:
-			cmd += " -min "+str(vmin)
-		if cluster >1:
-			cmd += " --perGroup --kmeans "+str(cluster)
-		if legendLoc:
-			cmd += " --legendLocation "+legendLoc
-		if title:
-			cmd += " --plotTitle " + title
-		data = subprocess.run(cmd, shell=True, capture_output=True)
-		print(data)
-
-def getPeaksAt(peaks, bigwigs, folder='', bigwignames=[], peaknames=[], window=1000, title='', numpeaks=4000, numthreads=8,
-				   width=5, length=10,torecompute=False, name='temp/peaksat.pdf', refpoint="TSS", scale=None,
-				   sort=False, withDeeptools=True, onlyProfile=False, cluster=1, vmax=None, vmin=None, overlap=False,
-				   legendLoc=None):
-	"""
-	get pysam data
-	ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
-	for each counts, do a rolling average (or a convolving of the data) with numpy
-	append to an array
-	return array, normalized
-	"""
-	if withDeeptools:
-		if isinstance(peaks, pd.DataFrame):
-			peaks = 'peaks.bed '
-			peaks.to_csv('peaks.bed', sep='\t', index=False, header=False)
-		elif type(peaks) == list:
-			pe = ''
-			i=0
-			for n, p in enumerate(peaks):
-				if 20 < int(os.popen('wc -l ' + p).read().split(' ')[0]):
-					pe += p + ' '
-				elif len(peaknames) > 0:
-					peaknames.pop(n-i)
-					i+=1
-			peaks = pe
-		elif type(peaks) == str:
-			peaks += ' '
-		else:
-			raise ValueError(' we dont know this filetype')
-		if type(bigwigs) is list:
-			pe = ''
-			for val in bigwigs:
-				pe += folder + val + ' '
-			bigwigs = pe
-		else:
-			bigwigs = folder + bigwigs + ' '
-		h.createFoldersFor(name)
-		cmd= ''
-		if not os.path.exists('.'.join(name.split('.')[:-1]) + ".gz") or torecompute:
-			cmd += "computeMatrix reference-point -S "
-			cmd += bigwigs
-			cmd += " --referencePoint "+refpoint
-			cmd += " --regionsFileName " + peaks
-			cmd += " --missingDataAsZero"
-			cmd += " --outFileName " + '.'.join(name.split('.')[:-1]) + ".gz"
-			cmd += " --upstream " + str(window) + " --downstream " + str(window)
-			cmd += " --numberOfProcessors " + str(numthreads) + ' && '
-		if type(name) is list:
-			cmd+= " --matrixFile " + '.gz '.join(name) + ".gz"
-		cmd += "plotHeatmap" if not onlyProfile else 'plotProfile'
-		cmd += " --matrixFile " + '.'.join(name.split('.')[:-1]) + ".gz"
-		cmd += " --outFileName " + name
-		cmd += " --refPointLabel "+ refpoint
-		if vmax is not None:
-			cmd += " -max "+str(vmax)
-		if vmin is not None:
-			cmd += " -min "+str(vmin)
-		if cluster>1:
-			cmd += " --perGroup --kmeans "+str(cluster)
-		if overlap:
-			if onlyProfile:
-				cmd += " --plotType overlapped_lines"
-			else:
-				raise ValueError("overlap only works when onlyProfile is set")
-		if legendLoc:
-			cmd+=" --legendLocation "+legendLoc
-
-		if len(peaknames) > 0:
-			pe = ''
-			for i in peaknames:
-				pe += ' ' + i
-			peaknames = pe
-			cmd += " --regionsLabel" + peaknames
-		if type(bigwigs) is list:
-			if len(bigwignames) > 0:
-				pe = ''
-				for i in bigwignames:
-					pe += ' ' + i
-				bigwignames = pe
-				cmd += " --samplesLabel" + bigwignames
-		if title:
-			cmd += " --plotTitle " + title
-		data = subprocess.run(cmd, shell=True, capture_output=True)
-		print(data)
-	else:
-		if 'relative_summit_pos' in peaks.columns:
-			center = [int((val['start'] + val['relative_summit_pos'])) for k, val in peaks.iterrows()]
-		else:
-			center = [int((val['start'] + val['end']) / 2) for k, val in peaks.iterrows()]
-		pd.set_option('mode.chained_assignment', None)
-		peaks['start'] = [c - window for c in center]
-		peaks['end'] = [c + window for c in center]
-		fig, ax = plt.subplots(1, len(bigwigs), figsize=[width, length], title=title if title else 'Chip Heatmap')
-		if sort:
-			peaks = peaks.sort_values(by=["foldchange"], ascending=False)
-		if numpeaks > len(peaks):
-			numpeaks = len(peaks) - 1
-		cov = {}
-		maxs = []
-		for num, bigwig in enumerate(bigwigs):
-			bw = pyBigWig.open(folder + bigwig)
-			co = np.zeros((numpeaks, window * 2), dtype=int)
-			scale = scale[bigwig] if scale is dict else 1
-			for i, (k, val) in enumerate(peaks.iloc[:numpeaks].iterrows()):
-				try:
-					co[i] = np.nan_to_num(bw.values(str(val.chrom), val.start, val.end), 0)
-				except RuntimeError as e:
-					print(str(val.chrom), val.start, val.end)
-					pass
-			cov[bigwig] = co
-			maxs.append(co.max())
-		for num, bigwig in enumerate(bigwigs):
-			sns.heatmap(cov[bigwig] * scale, ax=ax[num], vmax=max(maxs), yticklabels=[], cmap=cmaps[num],
-						cbar=True)
-			ax[num].set_title(bigwig.split('.')[0])
-		fig.subplots_adjust(wspace=0.1)
-		fig.show()
-		fig.savefig(name)
-		return cov, fig
-
-
-"""
-def computeMeanCov(bigwigFolder, meanOnly=True, ref="GRCh38", averageFragSize=150, outputfolder='/data/coverage/'):
-	meancov = {}
-	for val in os.listdir(bigwigFolder):
-		bw = pyBigWig.open(folder + bigwig)
-		if bw:
-			meancov[val.split('.')[0]] = bw.
-	return meancov
-
-
-def substractPeaks(peaks1, to):
-	peaks1 = peaks1.sort_values(by=['chrom', 'start'])
-	peaks2 = to.sort_values(by=['chrom', 'start'])
-	j = 0
-	i = 0
-	newpeaks = pd.DataFrame(columns=peaks2.columns)
-	while i < len(peaks2):
-		if peaks1.iloc[j].chrom == peaks2.iloc[j].chrom:
-			else
-		intersection = h.intersection([peaks2.iloc[i]['start'], peaks2.iloc[i]['end']],
-										   [peaks1.iloc[j]['start'], peaks1.iloc[j]['start']])
-		if intersection:
-
-			j += 1
-		while peaks2.iloc[i]['end'] > peaks1.iloc[j]['start']:
-			if peaks2.iloc[i]['end'] > peaks1.iloc[j]['end']:
-				newpeaks.append(peaks2.iloc[i], ignore_index=True)
-			else:
-				newpeaks.append({'start': peaks2.iloc[i]['start'],
-								 'end': peaks1.iloc[j]['start'],
-								 'chromomsome': peaks1.iloc[j]['chrom']}, ignore_index=True)
-"""
-
-
-def simpleMergePeaks(peaks, window=0, totpeaknumber=0, maxp=True, mergedFold="mean"):
-	"""
-	simply merges bedfiles from peak callers. providing a concaneted dataframe of bed-like tables
-
-	will recompute pvalues and foldchange from that.
-
-	"""
-	peaks = peaks.sort_values(by=['chrom', 'start','end'])
-	tfs = list(set(peaks['name']))
-	mergedpeaksdict = {}
-	remove = []
-	peaknumber = 0
-	merged_bed = {
-		"chrom": [peaks.iloc[0]['chrom']],
-		"start": [peaks.iloc[0]['start']],
-		"end": [],
-		"peak_number": [peaknumber + totpeaknumber],
-		"foldchange": [],
-		"-log10pvalue": [],
-		"-log10qvalue": [],
-		"relative_summit_pos": []
-	}
-	foldchange = [peaks.iloc[0].get('foldchange', 1)]
-	log10pvalue = [peaks.iloc[0].get('-log10pvalue', 0)]
-	log10qvalue = [peaks.iloc[0].get('-log10qvalue', 0)]
-	relative_summit_pos = peaks.iloc[1].get('relative_summit_pos', peaks.iloc[0]['start'])
-	# computes overlap by extending a bit the window (100bp?) should be ~readsize
-	prev_end = peaks.iloc[0]['end']
-	prev_chrom = peaks.iloc[0]['chrom']
-	tfmerged = {a: [0] for a in tfs}
-	tfmerged[peaks.iloc[0]['name']][-1] = peaks.iloc[0].get('foldchange', 1)
-	for i, (pos, peak) in enumerate(peaks.iloc[1:].iterrows()):
-		print(str(i / len(peaks)), end="\r")
-		if prev_end + window > peak['start'] and prev_chrom == peak['chrom']:
-			# can be merged
-			if peak.get('foldchange', 1) > max(foldchange):
-				relative_summit_pos = peak.get('relative_summit_pos', peaks['start'])
-			foldchange.append(peak.get('foldchange', 1))
-			log10pvalue.append(peak.get('-log10pvalue', 0))
-			log10qvalue.append(peak.get('-log10qvalue', 0))
-
-		else:
-			# newpeak
-			for k, val in tfmerged.items():
-				val.append(0)
-			peaknumber += 1
-			merged_bed['chrom'].append(peak['chrom'])
-			merged_bed['start'].append(peak['start'])
-			merged_bed['end'].append(prev_end)
-			merged_bed['peak_number'].append(peaknumber + totpeaknumber)
-			if mergedFold=="mean":
-				merged_bed['foldchange'].append(np.mean(foldchange))
-			elif mergedFold=="max":
-				merged_bed['foldchange'].append(max(foldchange))
-			elif mergedFold=="sum":
-				merged_bed['foldchange'].append(sum(foldchange))
-			else:
-				raise ValueError("mergedFold needs to be one of:")
-			merged_bed['-log10pvalue'].append(max(log10pvalue) if maxp else np.prod(log10pvalue))
-			merged_bed['-log10qvalue'].append(max(log10qvalue) if maxp else np.prod(log10qvalue))
-			merged_bed['relative_summit_pos'].append(relative_summit_pos)
-			foldchange = [peak.get('foldchange', 1)]
-			log10pvalue = [peak.get('-log10pvalue', 0)]
-			log10qvalue = [peak.get('-log10qvalue', 0)]
-			relative_summit_pos = peak.get('relative_summit_pos', peak['start'])
-		prev_end = peak['end']
-		prev_chrom = peak['chrom']
-		tfmerged[peak['name']][-1] = peak.get('foldchange', 1)
-	merged_bed['end'].append(prev_end)
-	if mergedFold=="mean":
-		merged_bed['foldchange'].append(np.mean(foldchange))
-	elif mergedFold=="max":
-		merged_bed['foldchange'].append(max(foldchange))
-	elif mergedFold=="sum":
-		merged_bed['foldchange'].append(sum(foldchange))
-	else:
-		raise ValueError("mergedFold needs to be one of:")
-	merged_bed['-log10pvalue'].append(max(log10pvalue) if maxp else np.prod(log10pvalue))
-	merged_bed['-log10qvalue'].append(max(log10qvalue) if maxp else np.prod(log10qvalue))
-	merged_bed['relative_summit_pos'].append(relative_summit_pos)
-
-	merged_bed = pd.DataFrame(merged_bed)
-	tfmerged = pd.DataFrame(tfmerged)
-	return pd.concat([merged_bed, tfmerged], axis=1, sort=False)
-
-
-def findpeakpath(folder, proteiname):
-	"""
-	given a folder of bigwigs and a protein name, finds the right bigwig
-	"""
-	res = None
-	for val in os.listdir(folder):
-		if str(proteiname) in val:
-			if res:
-				raise ValueError('more than 1 bigwig file found')
-			res= val
-	if res:
-		return res
-	raise ValueError('no bigwig file found')
-
-
-def findBestPeak(presence):
-	"""
-	given a list of -sets of peak locations for each replicate- will return the best replicate given a simple metric
-	"""
-	tot = []
-	for ind, el in enumerate(presence):
-		val = len(el)
-		pres = [x for j,x in enumerate(presence) if j!=ind]
-		for jnd in range(1, len(pres)+1):
-			for comb in itertools.combinations(pres, jnd):
-				ov = el
-				for knd in range(jnd):
-					ov = ov & comb[knd]
-				val += len(ov)*(jnd+1)
-		tot.append(val)
-	return np.argsort(tot)[::-1]
-
-
-def mergeReplicatePeaks(peaks, bigwigfolder, markedasbad=None, window=100,
-						sampling=3000, mincov=4, doPlot=True, cov={}, minKL=8, use='max',
-						MINOVERLAP=0.3, lookeverywhere=True, only='', saveloc=''):
-	"""
-	
-	
-	/!/ should only be passed peaks with at least one good replicate
-	for each TFpeaksets,
-	1. find the replicate that have the most peaks
-	2. correlate peaks and get in highest correlation order with the replicate found in 1
-	3. find overlap of both and get size of second replicate
-	4. if small(er)-> use only to increase statistics
-	  1. if a lot of uncalled peaks in replicate 2 at replicate 1 peaks (flag for mergebam)
-	5. if similar size -> get only intersect
-	  2. add to intersect, find uncalled peaks in both replicates which are called in the other
-	6. repeat for all replicates
-	-------------------------
-	if full overlap of one of the peak replicate, only use the overlapped one to increase confidence on peak
-	if >80% average non overlap,
-	  print warning and percentage of overlap
-
-	if <20% average non overlap,
-	  take the overlap and increase confidence and avg logfold
-
-	if one is <20%:
-	  if other <40% average non overlap,
-		take the overlap and increase confidence and avg logfold
-	  else
-		take
-
-	gets the max cov at the genomic window and if above some threshold, accepts the peak.
-
-	extend peak by X bp if no TSS
-	remove TSS from peaks
-
-	create a new data frame containing merged peak size, reassembled peak data (p value etc..) and
-	a the value for presence of each TF listed in previous df
-	------------------------------------
-
-	args:
-	----
-	peaks: df[bed-like] all the peaks into the sameBam with a column containing the 'name'
-	being the id of the sample, the 'replicate' number of this sample, the 'tf' chiped here
-	bamfolder: str, foldername
-	avgCov: dict(filename:int) a dict where for each bam filename is given an averageCoverage
-	if use=='max':
-		window:
-		mincov:
-
-	if use=='max':
-
-
-	returns:
-	-------
-	mergedpeaks: dict{df-peakslike}
-	bamtomerge: [[bam1,bam2]]
-
-	"""
-	def col_nan_scatter(x, y, **kwargs):
-		df = pd.DataFrame({'x': x[:], 'y': y[:]})
-		df = df[df.sum(0) != 0]
-		x = df['x']
-		y = df['y']
-		plt.gca()
-		plt.scatter(x, y)
-	def col_nan_kde_histo(x, **kwargs):
-		df = pd.DataFrame({'x':x[:]})
-		df = df[df['x']!=0]
-		x = df['x']
-		plt.gca()
-		sns.kdeplot(x)
-	print("/!/ should only be passed peaks with at least one good replicate")
-	# for a df containing a set of peaks in bed format and an additional column of different TF
-	tfs = list(set(peaks['tf']))
-	totpeaknumber = 0
-	mergedpeaksdict = {}
-	remove = []
-	tomergebam = []
-	ratiosofunique = {}
-	h.createFoldersFor(saveloc)
-	f = open(saveloc+'results.txt', 'w')
-	warnings.simplefilter("ignore")
-	for tf in tfs:
-		if only and tf!=only:
-			continue
-		cpeaks = peaks[peaks.tf==tf]
-		print('_____________________________________________________')
-		f.write('_____________________________________________________' + '\n')
-		if len(set(cpeaks['replicate'])) == 1:
-			if cpeaks.name.tolist()[0] in markedasbad:
-				print('the only replicate is considered bad!')
-				f.write('the only replicate is considered bad!'+"\n")
-				print('wrong TF: '+tf)
-				f.write('wrong TF: '+tf+"\n")
-				mergedpeaksdict.update({tf: cpeaks})
-				remove.append(tf)
-				continue
-			print("we only have one replicate for " + tf + " .. pass")
-			f.write("we only have one replicate for " + tf + " .. pass"+"\n")
-			mergedpeaksdict.update({tf: cpeaks})
-			continue
-		print("merging " + tf + " peaks")
-		f.write("merging " + tf + " peaks"+"\n")
-		merged = simpleMergePeaks(cpeaks, window=window, maxp=False)
-		merged_bed = merged[merged.columns[8:]]
-		finalpeaks = merged[merged.columns[:8]]
-		print('--> finish first overlaps lookup')
-		f.write('--> finish first overlaps lookup'+"\n")
-		# flag when  biggest is <1000 peaks
-		if len(finalpeaks) < 1000:
-			print('!TF has less than 1000 PEAKS!')
-			f.write('!TF has less than 1000 PEAKS!'+"\n")
-		# for each TF (replicates), compute number of peaks
-		peakmatrix = merged_bed.values.astype(bool)
-
-		presence = []
-		for peakpres in peakmatrix.T:  # https://github.com/tctianchi/pyvenn
-			presence.append(set([i for i, val in enumerate(peakpres) if val == 1]))
-		# compute overlap matrix (venn?)
-		if peakmatrix.shape[1] < 7 and doPlot:
-			plot.venn(presence, [i+'_BAD' if i.split('-')[0] in markedasbad else i for i in merged_bed.columns], title=tf+"_before_venn", folder=saveloc)
-			plt.show()
-		else:
-			print('too many replicates for Venn: '+str(peakmatrix.shape[1]))
-			f.write('too many replicates for Venn: '+str(peakmatrix.shape[1])+"\n")
-		if doPlot:
-			fig = sns.pairplot(merged_bed,corner=True, diag_kind="kde", kind="reg", plot_kws={"scatter_kws":{"alpha":.05}})
-			#fig = fig.map_upper(col_nan_scatter)
-			#fig = fig.map_upper(col_nan_kde_histo)
-			plt.suptitle("correlation of peaks in each replicate", y=1.08)
-			if saveloc:
-				fig.savefig(saveloc+tf+"_before_pairplot.pdf")
-			plt.show()
-			for i, val in enumerate(merged_bed):
-				unique_inval = np.logical_and(np.delete(peakmatrix,i,axis=1).sum(1).astype(bool)==0, peakmatrix[:,i])
-				sns.kdeplot(merged_bed[val][unique_inval], legend=True).set(xlim=(0,None))
-			plt.title("distribution of unique peaks in each replicate")
-			if saveloc:
-				plt.savefig(saveloc+tf+"_before_unique_kdeplot.pdf")
-			plt.show()
-
-		bigwigs = os.listdir(bigwigfolder)
-
-		foundgood=False
-		sort = findBestPeak(presence)
-		for ib,sb in enumerate(sort):
-			if merged_bed.columns[sb].split('-')[0] not in markedasbad:
-				foundgood=True
-				break
-		if not foundgood:
-			print('no peaks were good enough quality')
-			f.write('no peaks were good enough quality'+"\n")
-			print('bad TF: '+tf)
-			f.write('bad TF: '+tf+"\n")
-			remove.append(tf)
-			ib = 0
-		# distplot
-		# correlation plot
-
-
-		biggest_ind = sort[ib]
-		peakmatrix = peakmatrix.T
-		biggest = merged_bed.columns[biggest_ind]
-		print('-> main rep is: '+str(biggest))
-		f.write('-> main rep is: '+str(biggest)+'\n')
-		tot = peakmatrix[biggest_ind].copy().astype(int)
-		# starts with highest similarity and go descending
-		j = 0
-		recovered = 0
-		additionalpeaksinbig = np.array([])
-		for i, val in enumerate(sort):
-			if i==ib:
-				continue
-			j+=1
-			# if avg non overlap > 60%, and first, and none small flag TF as unreliable.
-			overlap = len(presence[biggest_ind] & presence[val]) / len(presence[biggest_ind])
-			peakname = merged_bed.columns[val]
-			print('- '+peakname)
-			f.write('- '+peakname+'\n')
-			print('  overlap: ' + str(overlap*100)+"%")
-			f.write('  overlap: ' + str(overlap*100)+"%"+'\n')
-			if overlap < MINOVERLAP:
-				smallsupport = len(presence[biggest_ind] & presence[val]) / len(presence[val])
-				print(' --> not enough overlap')
-				f.write(' --> not enough overlap'+'\n')
-				if smallsupport < MINOVERLAP:
-					# if the secondary does not have itself the required support
-					if j == 1 and merged_bed.columns[val].split('-')[0] not in markedasbad:
-						print("  Wrong TF: "+tf)
-						f.write("  Wrong TF: "+tf+'\n')
-						remove.append(tf)
-						break
-					# if not first, throw the other replicate and continue
-					print("  not using this replicate from the peakmatrix")
-					f.write("  not using this replicate from the peakmatrix"+'\n')
-					continue
-			if lookeverywhere:
-				tolookfor = peakmatrix[val] == 0
-			else:
-				tolookfor = np.logical_and(peakmatrix[biggest_ind], peakmatrix[val] == 0)
-			# ones that we have in the Primary but not in the secondary
-			additionalpeaksinsec = findAdditionalPeaks(finalpeaks, tolookfor, bigwigfolder + findpeakpath(bigwigfolder, peakname), sampling=sampling, mincov=mincov, window=window, minKL=minKL, use=use)
-			if len(additionalpeaksinsec[additionalpeaksinsec>0])>0:
-				sns.kdeplot(additionalpeaksinsec[additionalpeaksinsec>0],label=peakname, legend=True).set(xlim=(0,None))
-				print('  min,max from newly found peaks: '+str((additionalpeaksinsec[additionalpeaksinsec>0].min(),additionalpeaksinsec[additionalpeaksinsec>0].max())))
-				f.write('  min,max from newly found peaks: '+str((additionalpeaksinsec[additionalpeaksinsec>0].min(),additionalpeaksinsec[additionalpeaksinsec>0].max()))+'\n')
-			# for testing purposes mainly
-			finalpeaks[additionalpeaksinsec.astype(bool)].to_csv('additionalpeaksinsec_mp'+merged_bed.columns[val]+'.bed',sep='\t',index=None,header=False)
-			peakmatrix[val] = np.logical_or(peakmatrix[val], additionalpeaksinsec.astype(bool))
-			overlap = np.sum(np.logical_and(peakmatrix[val],peakmatrix[biggest_ind]))/np.sum(peakmatrix[biggest_ind])
-			if overlap < MINOVERLAP:
-				newsmalloverlap = np.sum(np.logical_and(peakmatrix[val],peakmatrix[biggest_ind]))/np.sum(peakmatrix[val])
-				print("  we did not had enough initial overlap.")
-				f.write("  we did not had enough initial overlap."+'\n')
-				if newsmalloverlap < MINOVERLAP:
-					if merged_bed.columns[val].split('-')[0] in markedasbad:
-						print('  replicate ' + merged_bed.columns[val] + ' was too bad and had not enough overlap')
-						f.write('  replicate ' + merged_bed.columns[val] + ' was too bad and had not enough overlap'+'\n')
-						continue
-					elif h.askif("we have two good quality peaks that don't merge well at all: "+merged_bed.columns[val] +\
-							" and " +merged_bed.columns[biggest_ind]+ " can the first one be removed?:\n  \
-							overlap: "+str(overlap*100)+'%\n  smalloverlap: '+str(smalloverlap*100)+'%\n  new smalloverlap: '+str(newsmalloverlap*100)+"%"):
-						continue
-					else:
-						print("  enough from small overlaps")
-						f.write("  enough from small overlaps"+'\n')
-			print(' --> enough overlap')
-			f.write(' --> enough overlap'+'\n')
-			recovered += np.sum(additionalpeaksinsec.astype(bool))
-			if merged_bed.columns[val].split('-')[0] not in markedasbad:
-				tot += peakmatrix[val].astype(int)
-			# ones that we have in the Primary but not in the secondary
-			if not lookeverywhere or len(additionalpeaksinbig)==0:
-				tolookfor = peakmatrix[biggest_ind] == 0 if lookeverywhere else np.logical_and(peakmatrix[biggest_ind]==0, peakmatrix[val])
-				additionalpeaksinbig = findAdditionalPeaks(finalpeaks, tolookfor, bigwigfolder + findpeakpath(bigwigfolder, biggest), sampling=sampling, mincov=mincov, window=window, minKL=minKL, use=use)
-				if len(additionalpeaksinbig[additionalpeaksinbig>0])>0:
-					sns.kdeplot(additionalpeaksinbig[additionalpeaksinbig>0],label=biggest, legend=True).set(xlim=(0,None))
-					print('  min,max from newly found peaks: '+str((additionalpeaksinbig[additionalpeaksinbig>0].min(),additionalpeaksinbig[additionalpeaksinbig>0].max())))
-					f.write('  min,max from newly found peaks: '+str((additionalpeaksinbig[additionalpeaksinbig>0].min(),additionalpeaksinbig[additionalpeaksinbig>0].max()))+'\n')
-
-				peakmatrix[biggest_ind] = np.logical_or(peakmatrix[biggest_ind], additionalpeaksinbig)
-				tot += additionalpeaksinbig.astype(bool).astype(int)
-				recovered += np.sum(additionalpeaksinbig.astype(bool))
-			print('  we have recovered ' + str(recovered)+' peaks, equal to '+ str(100*recovered/np.sum(peakmatrix[biggest_ind]))+\
-				'% of the peaks in main replicate')
-			f.write('  we have recovered ' + str(recovered)+' peaks, equal to '+ str(100*recovered/np.sum(peakmatrix[biggest_ind]))+\
-				'% of the peaks in main replicate'+'\n')
-			if overlap < (MINOVERLAP+0.2)/1.2:
-				# we recompute to see if the overlap changed
-				newoverlap = np.sum(np.logical_and(peakmatrix[val],peakmatrix[biggest_ind]))/np.sum(peakmatrix[biggest_ind])
-				smalloverlap = np.sum(np.logical_and(peakmatrix[val],peakmatrix[biggest_ind]))/np.sum(peakmatrix[val])
-				if newoverlap < (MINOVERLAP+0.2)/1.2:
-					if smalloverlap < (2+MINOVERLAP)/3:
-						print("  not enough overlap to advice to merge the bams.\n  oldnew overlap: "+str(overlap*100)+'%\n  \
-							new overlap: '+str(newoverlap*100)+"%")
-						f.write("  not enough overlap to advice to merge the bams.\n  oldnew overlap: "+str(overlap*100)+'%\n  \
-							new overlap: '+str(newoverlap*100)+"%"+'\n')
-						continue
-					else:
-						print('  enough from small overlap to advice to merge the peaks')
-						f.write('  enough from small overlap to advice to merge the peaks'+'\n')
-			tomergebam.append([biggest, peakname])
-			#the quality is good enough in the end we can pop from the list if it exists
-			if tf in remove:
-				remove.remove(tf)
-		plt.title('distribution of new found peaks')
-		if saveloc:
-			plt.savefig(saveloc+tf+"_new_found_peaks_kdeplot.pdf")
-		plt.show()
-		# new distplot
-		# new correlation plot
-		ratiosofunique[tf] = len(np.argwhere(peakmatrix.sum(0)==1))/peakmatrix.shape[1]
-		if doPlot:
-			sns.pairplot(merged_bed,corner=True, diag_kind="kde", kind="reg", plot_kws={"scatter_kws":{"alpha":.05}})
-			#fig = fig.map_upper(col_nan_scatter)
-			#fig = fig.map_upper(col_nan_kde_histo)
-			plt.suptitle("correlation and distribution of peaks after recovery", y=1.08)
-			if saveloc:
-				fig.savefig(saveloc+tf+"_after_pairplot.pdf")
-			plt.show()
-			for i, val in enumerate(merged_bed):
-				unique_inval = np.logical_and(np.delete(peakmatrix,i,axis=0).sum(0).astype(bool)==0, peakmatrix[i])
-				sns.kdeplot(merged_bed[val][unique_inval], legend=True).set(xlim=(0,None))
-			plt.title("distribution of unique peaks in each replicate after recovery")
-			if saveloc:
-				plt.savefig(saveloc+tf+"_after_unique_kdeplot.pdf")
-			plt.show()
-		if len(peakmatrix.shape) > 1 and doPlot:
-			if peakmatrix.shape[0] < 7:
-				presence = []
-				for peakpres in peakmatrix:  # https://github.com/tctianchi/pyvenn
-					presence.append(set([i for i, val in enumerate(peakpres) if val == 1]))
-				title = tf + '_recovered (TOREMOVE)' if tf in remove else tf+'_recovered'
-				h.venn(presence, [i+'_BAD' if i.split('-')[0] in markedasbad else i for i in merged_bed.columns], title=title, folder=saveloc)
-				plt.show()
-			else:
-				print('too many replicates for Venn')
-				f.write('(too many replicates for Venn)'+'\n')
-			finalpeaks = finalpeaks[np.logical_or(tot>1,peakmatrix[biggest_ind])]
-		finalpeaks['name'] = biggest
-		finalpeaks['tf'] = tf
-		mergedpeaksdict.update({tf: finalpeaks})
-		print(str((tf,len(finalpeaks))))
-		f.write(str((tf,len(finalpeaks)))+'\n')
-	mergedpeak = pd.concat([peaks for _, peaks in mergedpeaksdict.items()]).reset_index(drop=True)
-	if doPlot:
-		df= pd.DataFrame(data=ratiosofunique,index=['percentage of unique'])
-		df['proteins'] = df.index
-		fig = sns.barplot(data=df)
-		plt.xticks(rotation=60,ha='right')
-		plt.title("ratios of unique in replicates across experiments")
-		if saveloc:
-			plt.savefig(saveloc+"All_ratios_unique.pdf")
-		plt.show()
-	f.close()
-	mergedpeak['name'] = mergedpeak.tf
-	return mergedpeak, tomergebam, remove, ratiosofunique
-
-
-def findAdditionalPeaks(peaks, tolookfor, filepath, sampling=1000, mincov=4,
-						window=100, cov={}, minKL=8, use='max'):
-	"""
-	findAdditionalPeaks: for all peaks in A and/or B find in coverage file if zone has relative cov
-	of more than thresh then add to peak
-	if B is small and > 20% of peaks in A are found back, increase confidence and
-	flag for mergeBams
-	if < 20% don't flag for merge bam
-	f B is big and now mean non overlap < 40%, take union and flag for mergeBam else, throw B.
-
-	Args:
-	-----
-		peaks
-		tolookfor
-		filepath
-		sampling
-		mincov
-		window
-		cov
-		minKL
-		use
-	returns:
-	-------
-		np.array(bool) for each peaks in peakset, returns a binary
-	"""
-	# def poisson(k, lamb, scale): return scale * (lamb**k / factorial(k)) * np.exp(-lamb)
-
-	def KLpoisson(lamb1, lamb2): return lamb1 * np.log(lamb1 / lamb2) + lamb2 - lamb1
-	def poisson(k, lamb): return (lamb**k/factorial(k)) * np.exp(-lamb)
-	def negLogLikelihood(params, data): return - np.sum(np.log(poisson(data, params[0])))
-	def poissonFit(data): return float(minimize(negLogLikelihood,x0=np.ones(1),args=(data,),method='Powell').x)
-	bw = pyBigWig.open(filepath)
-	res = np.zeros(len(peaks))
-	prevchrom = ''
-	lamb = {}
-	cov = {}
-	#ignore by message
-	warnings.filterwarnings("ignore", message="encountered in")
-	for i, has in enumerate(tolookfor):
-		if has:
-			val = peaks.iloc[i]
-			if val.chrom not in chroms:
-				continue
-			if val.chrom != prevchrom:
-				if val.chrom not in cov:
-					cov[val.chrom] = bw.stats(str(val.chrom))[0]
-					prevchrom = val.chrom
-					if use == 'poisson':
-						#TODO: compute on INPUT file instead
-						samples = np.zeros(window * sampling)
-						sam = np.random.rand(sampling)
-						sam = sam * (bw.chroms(str(val.chrom))-window)
-						for j, sample in enumerate(sam.astype(int)):
-							samples[j*window:(j + 1)*window] = np.nan_to_num(bw.values(str(val.chrom), sample, sample + window), 0)
-						scale = np.unique(samples)[1]
-						samples = (samples/scale).astype(int)
-						lamb[val.chrom] = (poissonFit(samples),scale)
-
-			start = max([val.start - window, 0])
-			end = min(val.end + window, bw.chroms(str(val.chrom)))
-			zone = np.nan_to_num(bw.values(str(val.chrom), start, end), 0)
-			if use == 'max':
-				if max(zone) / cov[val.chrom] > mincov*1.5 or sum(zone) / (cov[val.chrom] * (end - start)) > mincov:
-					res[i] = max(zone) / cov[val.chrom]
-			elif use == 'poisson':
-				#TODO: compute -log10pvalue
-				la = poissonFit((zone/lamb[val.chrom][1]).astype(int))
-				kl=KLpoisson(la, lamb[val.chrom][0])
-				if kl > minKL:
-					res[i] = max(zone) / cov[val.chrom] #foldchange from macs2
-
-	return res
-
-
-def putInBed(conscensus, value, window=10, mergetype='mean'):
-	""" given ordered dataframes
-
-	conscensus df[start,end,chrom]
-	value df[start, end, chrom,foldchange]
-
-	type: oneof: mean,first,last,
-	"""
-	conscensus = conscensus.sort_values(by=['chrom','start','end']).reset_index(drop=True)
-	value = value.sort_values(by=['chrom','start','end']).reset_index(drop=True)
-	locinvalue=0
-	loc=0
-	tot=0
-	num = []
-	res = np.zeros(len(conscensus))
-	not_end=True
-	def add(res,num,not_end,loc):
-		if len(num)>0:
-			if mergetype=='mean':
-				res[loc] = np.mean(num)
-			elif mergetype=='first':
-				res[loc]=num[0]
-			elif mergetype=='last':
-				res[loc]=num[-1]
-			else:
-				raise ValueError('must be one of')
-			num= []
-		loc+=1
-		if loc == len(conscensus):
-			not_end=False
-		return res, num, not_end,loc
-	while not_end:
-		print(loc/len(conscensus),end="\r")
-		a = conscensus.iloc[loc]
-		b = value.iloc[locinvalue]
-		if b.chrom < a.chrom:
-			locinvalue+=1
-			if locinvalue == len(value):
-				not_end=False
-		elif b.chrom > a.chrom:
-			loc+=1
-			if loc == len(conscensus):
-				not_end=False
-		elif b.start<a.start:
-			if b.end+window>a.start:
-				tot+=1
-				num.append(b.foldchange)
-				if b.end>a.end+window:
-					res,num,not_end,loc = add(res,num,not_end,loc)
-					continue
-			locinvalue+=1
-			if locinvalue == len(value):
-				not_end=False
-		elif b.start<a.end+window:
-			tot+=1
-			num.append(b.foldchange)
-			if b.end>a.end+window:
-				res,num,not_end,loc = add(res,num,not_end,loc)
-				continue
-			locinvalue+=1
-			if locinvalue == len(value):
-				not_end=False
-		else:
-			res,num,not_end,loc = add(res,num,not_end,loc)
-	print(str(tot)+' were merged into conscensus')
-	return res
-
-
-def pairwiseOverlap(bedfile, norm=True, bedcol=8, correct=True, docorrelation=True):
-	"""
-	considering a befile representing a conscensus set of peaks
-	with each columns after the 7th one representing the signal of a given ChIP experiment
-	over this conscensus
-
-	overlap of j in i
-	overlap of row values in col values
-	"""
-	if correct:
-		print("we will be correcting for fully similar lines/ columns by removing 1 on their last value")
-	dat = bedfile[bedfile.columns[bedcol:]].values
-	prob = dat.astype(bool).sum(0)/len(dat)
-	correlation = np.ones((dat.shape[1],dat.shape[1]))
-	overlap = np.ones((dat.shape[1],dat.shape[1]))
-	for i, col in enumerate(dat.T):
-		#pdb.set_trace()
-		overlapping = np.delete(dat,i,axis=1)[col!=0]
-		col = col[col!=0]
-		add=0
-		for j, val in enumerate(overlapping.T):
-			if j==i:
-				add=1
-			if docorrelation:
-				if norm and not np.isnan(zscore(col[val!=0]).sum()) and not np.isnan(zscore(val[val!=0]).sum()) or not correct:
-					correlation[i,j+add] = np.corrcoef(zscore(val[val!=0]),zscore(col)[val!=0])[0,1]
-				else:
-					tmp = np.corrcoef(val[val != 0], col[val != 0])[0, 1]
-					if np.isnan(tmp) and correct:
-						if len(col[val!=0]) == 0 or len(val[val!=0]) == 0:
-						# one has no overlap
-							correlation[i,j+add] =  0
-						else:
-							# one contains only the same value everywhere
-							col[-1]-=max(0.01,abs(np.mean(col)))
-							val[-1]-=max(0.01,abs(np.mean(val)))
-							correlation[i,j+add] = np.corrcoef(val,col)[0,1]
-			overlap[i,j+add]=len(val[val!=0])/len(col)
-	if docorrelation:
-		correlation = pd.DataFrame(data=correlation, index=bedfile.columns[bedcol:], columns=bedfile.columns[bedcol:])
-		correlation[correlation.isna()] = 0
-	overlap = pd.DataFrame(data=overlap, index=bedfile.columns[bedcol:], columns=bedfile.columns[bedcol:]).T
-	return overlap, correlation if docorrelation else None
-
-
-def enrichment(bedfile, bedcol=8, groups=None, okpval=10**-3):
-	"""
-	considering a befile representing a conscensus set of peaks
-	with each columns after the 7th one representing the signal of a given ChIP experiment
-	over this conscensus
-
-	enrichment of j in i
-	enrichment of row values in col values
-	"""
-	dat = bedfile[bedfile.columns[bedcol:]].values
-	# pdb.set_trace()
-	prob = dat.astype(bool).sum(0)/len(dat)
-	enrichment = np.zeros((dat.shape[1] if groups is None else len(set(groups)), dat.shape[1]))
-	pvals = np.zeros(
-		(dat.shape[1] if groups is None else len(set(groups)), dat.shape[1]))
-	if groups is not None:
-		for i in set(groups):
-			overlapping = dat[groups==i]
-			for j,val in enumerate(overlapping.T):
-				# enrichment of j in i
-				e, p = fisher_exact([
-					[len(val[val != 0]), len(val[val == 0])],
-					[prob[j]*len(dat), (1-prob[j])*len(dat)]])
-				enrichment[i, j] = np.log2(e)
-				pvals[i, j] = p
-	else:
-		for i, col in enumerate(dat.T):
-			overlapping = np.delete(dat,i,axis=1)[col!=0]
-			col = col[col!=0]
-			add=0
-			for j, val in enumerate(overlapping.T):
-				if j==i:
-					add=1
-					enrichment[i,i]=0
-				e, p = fisher_exact([[len(val[val != 0]), len(val[val == 0])], [
-					prob[j+add]*len(dat), (1-prob[j+add])*len(dat)]])
-				enrichment[i, j+add] = np.log2(e)
-				pvals[i, j+add] = p
-		enrichment[i,i]=0
-	enrichment = pd.DataFrame(data=enrichment, index=bedfile.columns[bedcol:] if groups is None else set(groups), columns=bedfile.columns[bedcol:]).T
-	enrichment[enrichment==-np.inf] = -1000
-	enrichment[enrichment.isna()] = 0
-	enrichment[enrichment == np.inf] = 1000
-	pvals = np.reshape(multipletests(pvals.ravel(),
-									 0.1, method="bonferroni")[1], pvals.shape)
-	pvals = pd.DataFrame(
-		data=pvals, index=bedfile.columns[bedcol:]  if groups is None else set(groups), columns=bedfile.columns[bedcol:]).T
-	enrichment[pvals>okpval] = 0
-	return enrichment, pvals
-
-
-def findAdditionalCobindingSignal(conscensus, known=None, bigwigs=[], window=100):
-	"""
-	somewhat similar concept to computePeaksAt
-
-	# get pysam data
-	# ask for counts only at specific locus based on peak windows from mergedpeakset
-	# append to an array
-	# return array, normalized
-	"""
-	if known:
-		print('getting '+ str(len(peaks.tf))+' peaks. Using the peaks values directly if \
-				available and using the bigwigs otherwise.')
-		res = known.values.astype(float)
-	elif len(bigwigs)>0:
-		print('getting '+str(len(bigwigs))+' bigwigs, no peaks passed. Will compute the cobinding values\
-			across the conscensus for each bigwigs.')
-		res = np.zeros((len(bigwigs), len(conscensus)), dtype=float)
-	else:
-		raise ValueError('you need to pass a list of path to bigwigs for each/some samples')
-	for i, bw in enumerate(bigwigs):
-		if known:
-			found = False
-			for j, val in enumerate(known.tf):
-				if val in bw:
-					if found:
-						raise ValueError('found two or more matching tf for bigwig: '+str(bw))
-					found = True
-					i=j
-					break
-				if not found:
-					print('no tf found in known for tf: '+bw)
-					raise ValueError('you need to have an amount of known columns equal to your bigwigs')
-		print('doing file ' + str(bw))
-		bw = pyBigWig.open(bw)
-		for k, val in conscensus.iterrows():
-			if known:
-				if res[i][k]!=0:
-					continue
-			start = max([val.start - window, 0])
-			end = min(val.end + window, bw.chroms(str(val.chrom)))
-			res[i][k] = bw.stats(str(val.chrom), start, end)[0]
-	res = np.nan_to_num(res, 0)
-	return conscensus.join(pd.Dataframe(data=(res.T / res.max(1)).T, columns=bigwigs if not known else known.columns))
-
-
-
-def annotatePeaks():
-	"""
-	get H3k27AC peaks and compute Rose
-	for each TF peak
-	assign super enhancer if within it and create a super enhancer TFgroup
-	for each peaks
-	apply similar to rose and merge TF together. (create new TFgroup)
-	for each TF groups say
-	  if within super enhancer
-	  if within high h3k27ac marks
-	  its nearest gene (distance to it, cis/trans)
-
-
-	FIND A WAY TO FILTER TO ONLY PLACES WITH H3K27AC marks??
-
-	TAD points where most contacts on one side happened on this side.
-	specific distance. zone of most concentration of contacts for a peak region.
-	"""
-	# if way == 'Closest':
-
-	# elif way == 'ClosestExpressed':
-
-	# elif way == 'ActivityByContact':
-
-	# else:
-	#     raise ValueError('needs to be oneof Closest ClosestExpressed ActivityByContact')
-
-
-# os.system()
-
-
-def getCoLocalization():
-	"""
-	for each annotations (super enhancer & TFgroups)
-	for each TF, find highest peak/meanCov. if above thresh add to localization
-	"""
-
-
-def refineGroupsWithHiC():
-	"""
-	given HiC data, for each loops (should be less than X Mb distance. create a Xkb zone around both sides
-	and find + merge each TF/TFgroup at this location)
-	"""
-
-
-def fullDiffPeak(bam1, bam2, control1, size=None, control2=None, scaling=None, directory='diffData/',
-				 res_directory="diffPeaks/", isTF=False, compute_size=True, pairedend=True):
-	"""
-	will use macs2 to call differential peak binding from two bam files and their control
-
-	one can also provide some spike in scaling information
-
-	Args:
-	-----
-	bam1
-	bam2()
-	control1
-	control2
-	scaling
-	"""
-	print("doing diff from " + bam1 + " and " + bam2)
-	if scaling is not None:
-		if max(scaling) > 1:
-			raise ValueError("scalings need to be between 0-1")
-	name1 = bam1.split('/')[-1].split('.')[0]
-	name2 = bam2.split('/')[-1].split('.')[0]
-	if size is None:
-		if isTF:
-			size = 147
-		else:
-			size = 200
-	if compute_size:
-		print('computing the fragment avg size')
-		cmd = "macs2 predictd -i " + bam1
-		ret = subprocess.run(cmd, capture_output=True, shell=True)
-		size = re.findall("# predicted fragment length is (\d+)", str(ret.stderr))[0]
-		print(size)
-	else:
-		print('using default|given size')
-	pairedend = "BAMPE" if pairedend else "BAM"
-	if control2 is None:
-		control2 = control1
-	cmd1 = "macs2 callpeak -B -t " + bam1 + " -c " + control1 + " --nomodel --extsize " + str(size) + " -n " + name1 + " --outdir " + directory + " -f " + pairedend
-	cmd2 = "macs2 callpeak -B -t " + bam2 + " -c " + control2 + " --nomodel --extsize " + str(size) + " -n " + name2 + " --outdir " + directory + " -f " + pairedend
-	print('computing the scaling values')
-	ret = subprocess.run(cmd1, capture_output=True, shell=True)
-	print(ret.stderr)
-	scaling1a = int(re.findall(" after filtering in treatment: (\d+)", str(ret.stderr))[0])
-	scaling1b = int(re.findall(" after filtering in control: (\d+)", str(ret.stderr))[0])
-	scaling1 = scaling1a if scaling1a <= scaling1b else scaling1b
-	ret = subprocess.run(cmd2, capture_output=True, shell=True)
-	print(ret.stderr)
-	scaling2a = int(re.findall(" after filtering in treatment: (\d+)", str(ret.stderr))[0])
-	scaling2b = int(re.findall(" after filtering in control: (\d+)", str(ret.stderr))[0])
-	scaling2 = scaling2a if scaling2a <= scaling2b else scaling2b
-	if scaling is not None:
-		scaling1 = int(scaling1/scaling[0])
-		scaling2 = int(scaling2/scaling[1])
-	print(scaling1, scaling2)
-	return diffPeak(directory+name1+"_treat_pileup.bdg", directory+name2+"_treat_pileup.bdg",
-		directory+name1+"_control_lambda.bdg", directory+name2+"_control_lambda.bdg",
-		res_directory, scaling1, scaling2, size)
-
-
-def diffPeak(name1, name2, control1, control2, res_directory, scaling1, scaling2, size):
-	"""
-	calls MACS2 bdgdiff given the parameters
-	"""
-	print("doing differential peak binding")
-	cmd = "macs2 bdgdiff --t1 " + name1 + " --c1 "
-	cmd += control1+" --t2 " + name2 +" --c2 " + control2
-	cmd += " --d1 " + str(scaling1) + " --d2 " + str(scaling2) + " -g 60 "
-	cmd += "-l " + str(size) + " --o-prefix " + name1.split('/')[-1].split('.')[0] + "_vs_"
-	cmd += name2.split('/')[-1].split('.')[0] + " --outdir " + res_directory
-	res = subprocess.run(cmd, capture_output=True, shell=True)
-	return res
-
-
-def AssignToClosestExpressed(bed,countFile,genelocFile):
-	print("the bed file and genelocFile should use the same assembly")
-	genelocFile = pd.read_csv(genelocFile,sep="\t", compression="", columns=['chrom','start','end',"name"])
-	#for val in bed.iterrows():
-
-
-def MakeSuperEnhancers(MACS2bed, bamFile, outdir, baiFile=None, rosePath=".",
-	stitching_distance=None, TSS_EXCLUSION_ZONE_SIZE="2500", assembly="hg38",controlBam=None,controlBai=None):
-	"""
-	Calls super enhancer from H3K27ac with the ROSE algorithm
-
-	Args:
-	----
-		MACS2GFF
-		bamFile
-		outdir
-		baiFile
-		rosePath
-		stitching_distance
-		TSS_EXCLUSION_ZONE_SIZE
-		assembly
-		controlBam
-		controlBai
-
-	Returns:
-	--------
-		a bed-like dataframe with the superenhancers
-
-	outdir has to be an absolute path or a ~/path
-	"""
-	print("we are going to move your input files to "+rosePath)
-	cmd = 'mv '+MACS2bed+' '+rosePath+MACS2bed.split('/')[-1]+'.bed'
-	cmd += ' && mv '+bamFile+' '+rosePath
-	baiFile = baiFile if baiFile else bamFile[:-1]+'i'
-	cmd += ' && mv '+baiFile +' '+rosePath
-	res = subprocess.run(cmd, capture_output=True, shell=True)
-	if res.returncode != 0:
-		raise SystemError('failed moving files: '+str(res))
-	cmd = "cd "+rosePath+" && python ROSE_main.py -g "+assembly+" -i "+MACS2bed.split('/')[-1]+'.bed' + " -r " + bamFile.split('/')[-1] + " -o " + outdir
-	if TSS_EXCLUSION_ZONE_SIZE:
-		cmd+=" -t "+TSS_EXCLUSION_ZONE_SIZE
-	if controlBam:
-		os.system('mv '+controlBam+' '+rosePath)
-		controlBai = controlBai if controlBai else controlBam[:-1]+'i'
-		os.system('mv '+controlBai+' '+rosePath)
-		cmd+=" -c "+controlBam.split('/')[-1]
-	if stitching_distance:
-		cmd+=" -s "+ stitching_distance
-
-	res = subprocess.run(cmd, capture_output=True, shell=True)
-	fail = False
-	if res.returncode != 0:
-		v = 'ROSE failed:' +str(res)
-		fail = True
-	print('finished.. moving them back to their original folder')
-	cmd = 'mv '+rosePath+MACS2bed.split('/')[-1]+'.bed ' + MACS2bed
-	cmd += ' && mv '+rosePath+bamFile.split('/')[-1]+' '+bamFile
-	cmd+= ' && mv '+rosePath+baiFile.split('/')[-1]+' '+baiFile
-	if controlBam:
-		cmd += ' && mv '+rosePath+controlBam.split('/')[-1]+' '+controlBam
-		cmd += ' && mv '+rosePath+controlBai.split('/')[-1]+' '+controlBai
-	res = subprocess.run(cmd, capture_output=True, shell=True)
-	if res.returncode != 0:
-		raise SystemError('failed moving files: '+str(res))
-	if fail:
-		raise SystemError(v)
-	print('worked')
-	print(res)
-	return ReadRoseSuperEnhancers(outdir ,bool(controlBam))
-
-
-
-def runChromHMM(outdir, data, numstates=15, datatype='bed', folderPath=".", chromHMMFolderpath="~/ChromHMM/", assembly="hg38",control_bam_dir=None):
-	"""
-	runs chromHMM algorithm
-
-	Args:
-	-----
-		outdir str: an existing dir where the results should be saved
-		data: a df[cellname,markname,markbed|bam|bigwig, ?controlbed|bam|bigwig]
-		folderpath
-		numstates
-		datatype
-		folderPath
-		chromHMMFolderpath
-		assembly
-		control_bam_dir
-
-	Returns:
-	-------
-		A dict of bed like dataframes containing the regions of the different states
-	"""
-	print("you need to have ChromHMM")
-	chromHMM = "java -mx8000M -jar "+chromHMMFolderpath+"ChromHMM.jar "
-	h.createFoldersFor(outdir+'binarized/')
-	data.to_csv(outdir+"input_data.tsv", sep='\t',index=None,header=None)
-	cmd = chromHMM
-	if datatype=="bed":
-		cmd+="BinarizeBed "
-	elif datatype=="bigwig":
-		cmd+="BinarizeSignal "
-	elif datatype=="bam":
-		cmd+="BinarizeBam "
-	else:
-		raise ValueError('you need to provide one of bam, bigwig, bed')
-	cmd+= chromHMMFolderpath+"CHROMSIZES/"+assembly+".txt "+ folderPath+" "+outdir+"input_data.tsv "+outdir+"binarized"
-	if control_bam_dir:
-		cmd+=" -c "+control_bam_dir
-	res1 = subprocess.run(cmd, capture_output=True, shell=True)
-	print(res1)
-	if res1.returncode!=0:
-		raise ValueError(str(res1.stderr))
-	cmd = chromHMM + "LearnModel -printposterior -noautoopen "
-	if len(data)<10:
-		cmd += '-init load -m '+chromHMMFolderpath+'model_15_coreMarks.txt '
-	cmd += outdir+"binarized "+outdir+" "+str(numstates)+" "+assembly
-	res2 = subprocess.run(cmd, capture_output=True, shell=True)
-	print(res2)
-	if res2.returncode!=0:
-		raise ValueError(res2.stderr)
-	ret = {}
-	for v in set(data[0]):
-		ret[v] = pd.read_csv(outdir+v+'_'+str(numstates)+'_dense.bed', sep='\t', header=None,
-			skiprows=1).drop(columns=[4,5,6,7]).rename(columns=
-			{0:'chrom',1:'start',2:'end',3:'state',8:"color"})
-	return ret
-
-def loadMEMEmotifs(file, tfsubset=[],motifspecies='HUMAN'):
-	if file.endswith('.gff'):
-		print('converting to bed, you need to have "gfftobed" installed')
-		cmd = 'gff2bed < '+file+' > '+file+'.bed'
-		file = file+'.bed'
-		res = subprocess.run(cmd,capture_output=True, shell=True)
-		if res.returncode != 0:
-			raise ValueError('issue with the command: ' + str(res.stderr))
-		else:
-			print(res.stdout.decode("utf-8"))
-	## What are the motifs of our CRC members in ATACseq but not in our matrix
-	merged_motif = pd.read_csv(file, sep='\t',skiprows=0,index_col=None, names=['pos',"fimo",
-		"nucleotide_motif","relStart","relEnd","pval","strand",".","data"])
-	merged_motif['tf']=[i[5:].split("_"+motifspecies)[0] for i in merged_motif.data]
-	if tfsubset:
-		merged_motif = merged_motif[merged_motif.tf.isin(tfsubset)]
-	merged_motif['chrom'] = [i.split(':')[0][3:] for i in merged_motif.index]
-	merged_motif['start'] = [i.split(':')[1].split('-')[0] for i in merged_motif.index]
-	merged_motif['end'] = [i.split(':')[1].split('-')[1] for i in merged_motif.index]
-	merged_motif = merged_motif.reset_index(drop=True)
-	merged_motif = merged_motif.rename(columns={'pos':'relStart','fimo':'relEnd','nucleotide_motif':'pos',
-		'relStart':'pval','relEnd':'strand','pval':'fimo','strand':'motif'})
-	merged_motif['motif'] = [i.split('sequence=')[1].split(';')[0] for i in merged_motif.data]
-	merged_motif['p_val'] = [i.split('pvalue=')[1].split(';')[0] for i in merged_motif.data]
-	merged_motif['q_val'] = [i.split('qvalue=')[1].split(';')[0] for i in merged_motif.data]
-	merged_motif = merged_motif.drop(columns=['pos','.','fimo','data'])
-	merged_motif = merged_motif[merged_motif.columns[[6,7,8,0,1,3,2,9,10,5,4]]]
-	merged_motif = merged_motif.sort_values(by=['chrom','start','end']).reset_index(drop=True)
-	return merged_motif
-
-
-def simpleMergeMotifs(motifs, window=0):
-	if type(motifs) is list:
-		motifs = pd.concat(motifs)
-	motifs = motifs.sort_values(by=['chrom', 'start'])
-	toremove = []
-	issues = []
-	prevmotif = motifs.iloc[0]
-	for i, (pos, motif) in enumerate(motifs.iloc[1:].iterrows()):
-		print(str(i / len(motifs)), end="\r")
-		if prevmotif['end'] + window > motif['start'] and prevmotif['chrom'] == motif['chrom']:
-			# can be merged
-			if motif['tf']!= prevmotif['tf'] or motif['motif'] != prevmotif['motif']:
-				print('found different motifs overlapping')
-				issues.extend([motif,prevmotif])
-			else:
-				toremove.append(pos)
-		prevmotif = motif
-	motifs = motifs.drop(index=toremove).reset_index(drop=True)
-	issues = pd.concat(issues)
-	return motifs, issues
-
-
-def substractPeaksTo(peaks,loci, bp=50):
-	"""
-	removes all peaks that are not within a bp distance to a set of loci
-
-	Args:
-	----
-		peaks: a bed file df with a chrom,start, end column at least
-		loci: a df witth a chrom & loci column
-		bp: the max allowed distance to the loci
-
-	Returns:
-	-------
-		all the peaks that are within this distance
-	"""
-	i=0
-	j=0
-	keep=[]
-	bp=50
-	while j<len(peaks) and i<len(loci):
-		h.showcount(j,len(peaks))
-		if peaks.loc[j].chrom > loci.loc[i].chrom:
-			i+=1
-			continue
-		if peaks.loc[j].chrom < loci.loc[i].chrom:
-			j+=1
-			continue
-		if peaks.loc[j].start - bp > loci.loc[i].loci:
-			i+=1
-			continue
-		if peaks.loc[j].end + bp< loci.loc[i].loci:
-			j+=1
-			continue
-		if peaks.loc[j].end + bp >= loci.loc[i].loci and peaks.loc[j].start - bp <= loci.loc[i].loci:
-			keep.append(j)
-			j+=1
-	return peaks.loc[set(keep)]
-
-	
diff --git a/epigenetics/rose b/epigenetics/rose
deleted file mode 160000
index 028ee65..0000000
--- a/epigenetics/rose
+++ /dev/null
@@ -1,1 +0,0 @@
-Subproject commit 028ee650a3006c757b32328b37f7644823571441
diff --git a/__init__.py b/genepy/__init__.py
similarity index 100%
rename from __init__.py
rename to genepy/__init__.py
diff --git a/cell_line_mapping-master/.gitignore b/genepy/cell_line_mapping-master/.gitignore
similarity index 100%
rename from cell_line_mapping-master/.gitignore
rename to genepy/cell_line_mapping-master/.gitignore
diff --git a/cell_line_mapping-master/.travis.yml b/genepy/cell_line_mapping-master/.travis.yml
similarity index 100%
rename from cell_line_mapping-master/.travis.yml
rename to genepy/cell_line_mapping-master/.travis.yml
diff --git a/cell_line_mapping-master/README.md b/genepy/cell_line_mapping-master/README.md
similarity index 100%
rename from cell_line_mapping-master/README.md
rename to genepy/cell_line_mapping-master/README.md
diff --git a/cell_line_mapping-master/celllinemapr/.Rbuildignore b/genepy/cell_line_mapping-master/celllinemapr/.Rbuildignore
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/.Rbuildignore
rename to genepy/cell_line_mapping-master/celllinemapr/.Rbuildignore
diff --git a/cell_line_mapping-master/celllinemapr/.celllinemapr.Rds b/genepy/cell_line_mapping-master/celllinemapr/.celllinemapr.Rds
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/.celllinemapr.Rds
rename to genepy/cell_line_mapping-master/celllinemapr/.celllinemapr.Rds
diff --git a/cell_line_mapping-master/celllinemapr/DESCRIPTION b/genepy/cell_line_mapping-master/celllinemapr/DESCRIPTION
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/DESCRIPTION
rename to genepy/cell_line_mapping-master/celllinemapr/DESCRIPTION
diff --git a/cell_line_mapping-master/celllinemapr/NAMESPACE b/genepy/cell_line_mapping-master/celllinemapr/NAMESPACE
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/NAMESPACE
rename to genepy/cell_line_mapping-master/celllinemapr/NAMESPACE
diff --git a/cell_line_mapping-master/celllinemapr/R/cell_line_mapping.R b/genepy/cell_line_mapping-master/celllinemapr/R/cell_line_mapping.R
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/R/cell_line_mapping.R
rename to genepy/cell_line_mapping-master/celllinemapr/R/cell_line_mapping.R
diff --git a/cell_line_mapping-master/celllinemapr/SOP.Rmd b/genepy/cell_line_mapping-master/celllinemapr/SOP.Rmd
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/SOP.Rmd
rename to genepy/cell_line_mapping-master/celllinemapr/SOP.Rmd
diff --git a/cell_line_mapping-master/celllinemapr/celllinemapr.Rproj b/genepy/cell_line_mapping-master/celllinemapr/celllinemapr.Rproj
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/celllinemapr.Rproj
rename to genepy/cell_line_mapping-master/celllinemapr/celllinemapr.Rproj
diff --git a/cell_line_mapping-master/celllinemapr/man/arxspan.to.ccle.Rd b/genepy/cell_line_mapping-master/celllinemapr/man/arxspan.to.ccle.Rd
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/man/arxspan.to.ccle.Rd
rename to genepy/cell_line_mapping-master/celllinemapr/man/arxspan.to.ccle.Rd
diff --git a/cell_line_mapping-master/celllinemapr/man/ccle.to.arxspan.Rd b/genepy/cell_line_mapping-master/celllinemapr/man/ccle.to.arxspan.Rd
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/man/ccle.to.arxspan.Rd
rename to genepy/cell_line_mapping-master/celllinemapr/man/ccle.to.arxspan.Rd
diff --git a/cell_line_mapping-master/celllinemapr/naming.csv b/genepy/cell_line_mapping-master/celllinemapr/naming.csv
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/naming.csv
rename to genepy/cell_line_mapping-master/celllinemapr/naming.csv
diff --git a/cell_line_mapping-master/python/cell_line_mapper/__init__.py b/genepy/cell_line_mapping-master/python/cell_line_mapper/__init__.py
similarity index 100%
rename from cell_line_mapping-master/python/cell_line_mapper/__init__.py
rename to genepy/cell_line_mapping-master/python/cell_line_mapper/__init__.py
diff --git a/cell_line_mapping-master/python/name_mapping.csv b/genepy/cell_line_mapping-master/python/name_mapping.csv
similarity index 100%
rename from cell_line_mapping-master/python/name_mapping.csv
rename to genepy/cell_line_mapping-master/python/name_mapping.csv
diff --git a/cell_line_mapping-master/python/requirements.txt b/genepy/cell_line_mapping-master/python/requirements.txt
similarity index 100%
rename from cell_line_mapping-master/python/requirements.txt
rename to genepy/cell_line_mapping-master/python/requirements.txt
diff --git a/cell_line_mapping-master/python/setup.py b/genepy/cell_line_mapping-master/python/setup.py
similarity index 100%
rename from cell_line_mapping-master/python/setup.py
rename to genepy/cell_line_mapping-master/python/setup.py
diff --git a/genepy/epigenetics/CREME.md b/genepy/epigenetics/CREME.md
new file mode 100644
index 0000000..6f4b81a
--- /dev/null
+++ b/genepy/epigenetics/CREME.md
@@ -0,0 +1,92 @@
+# genepy/CREME: ChIP REplicate MErger
+
+For Introduction we will link to the [article](https://ro-che.info/articles/2018-07-11-chip-seq-consensus) by Roman Cheplyaka on the subject.
+
+We built this tool noticing the lack of publicly available simple Chip Merger working for [MACS2](https://github.com/macs3-project/MACS)'s output, with replicates of broadly different quality. We wanted a 1 function tool that would work in python.
+
+CREME is part of the [genepy](https://github.com/broadinstitute/GenePy) package.
+
+We will although note tools such as:
+- [PePr](https://pubmed.ncbi.nlm.nih.gov/24894502/) [code](https://github.com/shawnzhangyx/PePr) which can substitute itself from MACS2 by ccalling on mutliple bam files at the same time. It will work by counting reads and looking at the peak shape.
+- [multiGPS](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003501), [code](https://github.com/seqcode/multigps) which is mostly  for differential binding chipseq but can work with replicates and work in java + R. 
+- [MSPC](https://academic.oup.com/bioinformatics/article/31/17/2761/183989), [code](https://github.com/Genometric/MSPC) in .NET, which is very well documented, simple and provide some QC by the user.
+
+- [genoGAM](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2238-7) in R, [code](https://github.com/gstricker/GenoGAM). Which calls peaks by itself as well and seem to handle replicates.
+
+- [sierra Platinum](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5025614/), [code](https://github.com/sierraplatinum/sierra), which does not seem to be maintained.
+
+## Our tool
+
+The goal of creme is to be a simple, 1 function tool. It works with 1 to many sets of replicates for each pulled protein/mark.
+
+CREME takes as input a pandas dataframe. This dataframe is the concatenation of each replicates' bed files and can be loaded from a set of MACS2 bedfiles using genepy's loadPeaks function.
+
+CREME will output, amongst other thing, a dataframe representing a concatenation f bedfiles of merged replicates.
+
+## Process
+
+### Selecion: Finding the best replicate
+
+A first goal of CREME was to find the best replicate, to do so, it can take manual annotation of bad replicates. These can be provided by visual inspection of bigwig tracks + bed files on IGV, from thresholding on QC results such as FRiP scores.
+
+![plot igv]()
+
+Given all available replicates, CREME will compute a conscensus, considering any peaks at most 150 bp from another peak, to be in overlap. We have noticed that changing this parameter from 0 to 150 decreased the total number of peaks found by only 8%.
+
+Non overlapping peaks are kept in the conscensus. When we have an overlap we take the mean of signals and the product of pvalues across overlapping replicates.
+
+![plot venn](MED1_before_venn_venn.png)
+
+Then, CREME will look at their overlap and select the one that has the best overlap score:
+
+$O_{score}(A) = \sum{i from 0 to m} \sum{K in comb(i, G)} i * \sum {j from 0 to n}  AND(A[j],...K[j])$
+
+Where:
+- $G$ is a binary matrix of size (row/col) $m*n$ of $m$ replicates with $n$ conscensus peaks and a value of 1 if replicate $m_i$ has peak on conscensus peak $n_i$.
+- $comb(i, G)$ is a list of all possible matrices made from taking $i$ elements (row) from matrix $G$ without replacement.
+- $AND$ is a binary operation returning 1 if all passed elements are 1 else 0.
+
+The non-bad quality replicate with the best score will be selected as the __main replicate__.
+
+In addition to the venn diagram, correlation between each replicate's peak signal is computed and displayed to the user.
+
+![pairplot of replicates](MED1_before_pairplot.png)
+
+### Validation: Finding new peaks
+
+For each additional replicate S, we will now look for new peaks. 
+First, if we find that the second best replicate and the first best replicate have both less than 30% of their peaks in common we __discard__ that protein/mark and only return the main replicate.
+
+Taking peaks that are found in the main replicate, we will call peaks, using S's bigwig and a lower threshold than what MACS2 is using by default. We then do the same for peaks in S that were not in the main replicate.
+
+If after calling new peaks we get less than 30% overlap in both replicates, we discard the replicate.
+
+Else, we finalize the merging of overlapping peaks and update the __main replicate__ with this overlap.
+
+### Calling Peaks
+
+The process of calling peaks is loosely based on MACS2's peak calling algorithm:
+
+We compute a distance: the [KL divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence), between two poisson distributions. One is representing the distribution of signal from a bigwig file under a region. The other is the representing the same signal under the entire chromosome where that region lies. The region here is the peak in the other sample that we want to look for in the current sample.
+
+If that distance is above a threshold: here 8, we validate the region as being a peak.
+
+### Output and QC
+
+The output of our tool is a dataframe of concatenated merged replicates. The pipeline also outputs a set of bad quality replicates and bad quality proteins/marks.
+
+Additionally, inforrmation on distribution of peak signal across replicates and number of peaks found is provided to the user.
+
+![kdeplot of new found peaks](MED1_new_found_peaks_kdeplot.png)
+
+## WIP and current issues
+
+1. For now, we are not using the exact same algorithm as MACS2, as we are comparing the peak's read distribution to overall reads in the chromosome using KL divergence. But MACS2 is comparing 4 terms: the distribution in the likely region of the sample BAM, the distribution in the likely region of the INPUT BAM, the distribution in the sample BAM's chromosome, the distribution in the INPUT BAM's chromosome. Moreover, MACS is comparing them using something like a fisher's exact test and corrects for FDR using the BH method.
+
+2. For now, we are not computing a perfect overall replicate quality ourselves. Our  scoring method did not work in 5% of cases. We might want to mitigate it by adding peaks' Qvalues and the replicate's Frip score and total read count in our analysis.
+
+3. For now, we do not compute pvalue when we compute new peaks.
+
+4. For now, we do not integrate the pvalue/ signal of newly found peaks in the conscensus merger.  
+
+5. More long term: We would hope to do something more akin to joint calling across replicates using graphical models to call peaks.
diff --git a/genepy/epigenetics/CREME.py b/genepy/epigenetics/CREME.py
new file mode 100644
index 0000000..55a0ad3
--- /dev/null
+++ b/genepy/epigenetics/CREME.py
@@ -0,0 +1,474 @@
+import os
+import pandas as pd
+import numpy as np
+from genepy.utils import helper as h
+from genepy.utils import plot
+from genepy.epigenetics.chipseq import *
+import seaborn as sns
+import pyBigWig
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize
+from scipy.special import factorial
+import warnings
+import itertools
+
+
+def findpeakpath(folder, proteiname):
+	"""
+	given a folder of bigwigs and a protein name, finds the right bigwig
+	"""
+	res = None
+	for val in os.listdir(folder):
+		if str(proteiname) in val:
+			if res:
+				raise ValueError('more than 1 bigwig file found')
+			res = val
+	if res:
+		return res
+	raise ValueError('no bigwig file found')
+
+
+def findBestPeak(presence):
+	"""
+	given a list of -sets of peak locations for each replicate- will return the best replicate given a simple metric
+	"""
+	tot = []
+	for ind, el in enumerate(presence):
+		val = len(el)
+		pres = [x for j, x in enumerate(presence) if j != ind]
+		for jnd in range(1, len(pres)+1):
+			for comb in itertools.combinations(pres, jnd):
+				ov = el
+				for knd in range(jnd):
+					ov = ov & comb[knd]
+				val += len(ov)*(jnd+1)
+		tot.append(val)
+	return np.argsort(tot)[::-1]
+
+
+def mergeReplicatePeaks(peaks, bigwigfolder, markedasbad=None, window=100,
+                        sampling=3000, mincov=4, doPlot=True, cov={}, minKL=8, use='max',
+                        MINOVERLAP=0.3, lookeverywhere=True, only='', saveloc=''):
+  """
+	
+	
+	/!/ should only be passed peaks with at least one good replicate
+	for each TFpeaksets,
+	1. find the replicate that have the most peaks
+	2. correlate peaks and get in highest correlation order with the replicate found in 1
+	3. find overlap of both and get size of second replicate
+	4. if small(er)-> use only to increase statistics
+	  1. if a lot of uncalled peaks in replicate 2 at replicate 1 peaks (flag for mergebam)
+	5. if similar size -> get only intersect
+	  2. add to intersect, find uncalled peaks in both replicates which are called in the other
+	6. repeat for all replicates
+	-------------------------
+	if full overlap of one of the peak replicate, only use the overlapped one to increase confidence on peak
+	if >80% average non overlap,
+	  print warning and percentage of overlap
+
+	if <20% average non overlap,
+	  take the overlap and increase confidence and avg logfold
+
+	if one is <20%:
+	  if other <40% average non overlap,
+		take the overlap and increase confidence and avg logfold
+	  else
+		take
+
+	gets the max cov at the genomic window and if above some threshold, accepts the peak.
+
+	extend peak by X bp if no TSS
+	remove TSS from peaks
+
+  create a new data frame containing merged peak size, reassembled peak data (p value etc..) and
+  a the value for presence of each TF listed in previous df
+  ------------------------------------
+
+  args:
+  ----
+  peaks: df[bed-like] all the peaks into the sameBam with a column containing the 'name'
+  being the id of the sample, the 'replicate' number of this sample, the 'tf' chiped here
+  bamfolder: str, foldername
+  avgCov: dict(filename:int) a dict where for each bam filename is given an averageCoverage
+  if use=='max':
+    window:
+    mincov:
+
+  if use=='max':
+
+
+  returns:
+  -------
+  mergedpeaks: dict{df-peakslike}
+  bamtomerge: [[bam1,bam2]]
+
+  """
+  def col_nan_scatter(x, y, **kwargs):
+    df = pd.DataFrame({'x': x[:], 'y': y[:]})
+    df = df[df.sum(0) != 0]
+    x = df['x']
+    y = df['y']
+    plt.gca()
+		plt.scatter(x, y)
+
+	def col_nan_kde_histo(x, **kwargs):
+		df = pd.DataFrame({'x': x[:]})
+		df = df[df['x'] != 0]
+		x = df['x']
+		plt.gca()
+		sns.kdeplot(x)
+	print("/!/ should only be passed peaks with at least one good replicate")
+	# for a df containing a set of peaks in bed format and an additional column of different TF
+	tfs = list(set(peaks['tf']))
+	totpeaknumber = 0
+	mergedpeaksdict = {}
+	remove = []
+	tomergebam = []
+	ratiosofunique = {}
+	h.createFoldersFor(saveloc)
+	f = open(saveloc+'results.txt', 'w')
+	warnings.simplefilter("ignore")
+	for tf in tfs:
+		if only and tf != only:
+			continue
+		cpeaks = peaks[peaks.tf == tf]
+		print('_____________________________________________________')
+		f.write('_____________________________________________________' + '\n')
+		if len(set(cpeaks['replicate'])) == 1:
+			if cpeaks.name.tolist()[0] in markedasbad:
+				print('the only replicate is considered bad!')
+				f.write('the only replicate is considered bad!'+"\n")
+				print('wrong TF: '+tf)
+				f.write('wrong TF: '+tf+"\n")
+				mergedpeaksdict.update({tf: cpeaks})
+				remove.append(tf)
+				continue
+			print("we only have one replicate for " + tf + " .. pass")
+			f.write("we only have one replicate for " + tf + " .. pass"+"\n")
+			mergedpeaksdict.update({tf: cpeaks})
+			continue
+		print("merging " + tf + " peaks")
+		f.write("merging " + tf + " peaks"+"\n")
+		merged = simpleMergePeaks(cpeaks, window=window, maxp=False)
+		merged_bed = merged[merged.columns[8:]]
+		finalpeaks = merged[merged.columns[:8]]
+		print('--> finish first overlaps lookup')
+		f.write('--> finish first overlaps lookup'+"\n")
+		# flag when  biggest is <1000 peaks
+		if len(finalpeaks) < 1000:
+			print('!TF has less than 1000 PEAKS!')
+			f.write('!TF has less than 1000 PEAKS!'+"\n")
+		# for each TF (replicates), compute number of peaks
+		peakmatrix = merged_bed.values.astype(bool)
+
+		presence = []
+		for peakpres in peakmatrix.T:  # https://github.com/tctianchi/pyvenn
+			presence.append(set([i for i, val in enumerate(peakpres) if val == 1]))
+		# compute overlap matrix (venn?)
+		if peakmatrix.shape[1] < 7 and doPlot:
+		plot.venn(presence, [i+'_BAD' if i.split('-')[0]
+                       in markedasbad else i for i in merged_bed.columns], title=tf+"_before_venn", folder=saveloc)
+		plt.show()
+		else:
+			print('too many replicates for Venn: '+str(peakmatrix.shape[1]))
+			f.write('too many replicates for Venn: '+str(peakmatrix.shape[1])+"\n")
+		if doPlot:
+			fig = sns.pairplot(merged_bed, corner=True, diag_kind="kde",
+			                   kind="reg", plot_kws={"scatter_kws": {"alpha": .05}})
+			#fig = fig.map_upper(col_nan_scatter)
+			#fig = fig.map_upper(col_nan_kde_histo)
+			plt.suptitle("correlation of peaks in each replicate", y=1.08)
+			if saveloc:
+				fig.savefig(saveloc+tf+"_before_pairplot.pdf")
+			plt.show()
+			for i, val in enumerate(merged_bed):
+				unique_inval = np.logical_and(
+					np.delete(peakmatrix, i, axis=1).sum(1).astype(bool) == 0, peakmatrix[:, i])
+				sns.kdeplot(merged_bed[val][unique_inval], legend=True).set(xlim=(0, None))
+			plt.title("distribution of unique peaks in each replicate")
+			if saveloc:
+				plt.savefig(saveloc+tf+"_before_unique_kdeplot.pdf")
+			plt.show()
+
+		bigwigs = os.listdir(bigwigfolder)
+
+		foundgood = False
+		sort = findBestPeak(presence)
+		for ib, sb in enumerate(sort):
+			if merged_bed.columns[sb].split('-')[0] not in markedasbad:
+				foundgood = True
+				break
+		if not foundgood:
+			print('no peaks were good enough quality')
+			f.write('no peaks were good enough quality'+"\n")
+			print('bad TF: '+tf)
+			f.write('bad TF: '+tf+"\n")
+			remove.append(tf)
+			ib = 0
+		# distplot
+		# correlation plot
+
+		biggest_ind = sort[ib]
+		peakmatrix = peakmatrix.T
+		biggest = merged_bed.columns[biggest_ind]
+		print('-> main rep is: '+str(biggest))
+		f.write('-> main rep is: '+str(biggest)+'\n')
+		tot = peakmatrix[biggest_ind].copy().astype(int)
+		# starts with highest similarity and go descending
+		j = 0
+		recovered = 0
+		additionalpeaksinbig = np.array([])
+		for i, val in enumerate(sort):
+			if i == ib:
+				continue
+			j += 1
+			# if avg non overlap > 60%, and first, and none small flag TF as unreliable.
+			overlap = len(presence[biggest_ind] & presence[val]
+			              ) / len(presence[biggest_ind])
+			peakname = merged_bed.columns[val]
+			print('- '+peakname)
+			f.write('- '+peakname+'\n')
+			print('  overlap: ' + str(overlap*100)+"%")
+			f.write('  overlap: ' + str(overlap*100)+"%"+'\n')
+			if overlap < MINOVERLAP:
+				smallsupport = len(presence[biggest_ind] &
+				                   presence[val]) / len(presence[val])
+				print(' --> not enough overlap')
+				f.write(' --> not enough overlap'+'\n')
+				if smallsupport < MINOVERLAP:
+					# if the secondary does not have itself the required support
+					if j == 1 and merged_bed.columns[val].split('-')[0] not in markedasbad:
+						print("  Wrong TF: "+tf)
+						f.write("  Wrong TF: "+tf+'\n')
+						remove.append(tf)
+						break
+					# if not first, throw the other replicate and continue
+					print("  not using this replicate from the peakmatrix")
+					f.write("  not using this replicate from the peakmatrix"+'\n')
+					continue
+			if lookeverywhere:
+				tolookfor = peakmatrix[val] == 0
+			else:
+				tolookfor = np.logical_and(peakmatrix[biggest_ind], peakmatrix[val] == 0)
+			# ones that we have in the Primary but not in the secondary
+			additionalpeaksinsec = findAdditionalPeaks(finalpeaks, tolookfor, bigwigfolder + findpeakpath(
+				bigwigfolder, peakname), sampling=sampling, mincov=mincov, window=window, minKL=minKL, use=use)
+			if len(additionalpeaksinsec[additionalpeaksinsec > 0]) > 0:
+				sns.kdeplot(additionalpeaksinsec[additionalpeaksinsec > 0],
+				            label=peakname, legend=True).set(xlim=(0, None))
+				print('  min,max from newly found peaks: ' +
+				      str((additionalpeaksinsec[additionalpeaksinsec > 0].min(), additionalpeaksinsec[additionalpeaksinsec > 0].max())))
+				f.write('  min,max from newly found peaks: '+str((additionalpeaksinsec[additionalpeaksinsec > 0].min(
+				), additionalpeaksinsec[additionalpeaksinsec > 0].max()))+'\n')
+			# for testing purposes mainly
+			finalpeaks[additionalpeaksinsec.astype(bool)].to_csv(
+				'additionalpeaksinsec_mp'+merged_bed.columns[val]+'.bed', sep='\t', index=None, header=False)
+			peakmatrix[val] = np.logical_or(
+				peakmatrix[val], additionalpeaksinsec.astype(bool))
+			overlap = np.sum(np.logical_and(
+				peakmatrix[val], peakmatrix[biggest_ind]))/np.sum(peakmatrix[biggest_ind])
+			if overlap < MINOVERLAP:
+				newsmalloverlap = np.sum(np.logical_and(
+					peakmatrix[val], peakmatrix[biggest_ind]))/np.sum(peakmatrix[val])
+				print("  we did not had enough initial overlap.")
+				f.write("  we did not had enough initial overlap."+'\n')
+				if newsmalloverlap < MINOVERLAP:
+					if merged_bed.columns[val].split('-')[0] in markedasbad:
+						print('  replicate ' +
+						      merged_bed.columns[val] + ' was too bad and had not enough overlap')
+						f.write('  replicate ' +
+						        merged_bed.columns[val] + ' was too bad and had not enough overlap'+'\n')
+						continue
+					elif h.askif("we have two good quality peaks that don't merge well at all: "+merged_bed.columns[val] +
+                                            " and " + merged_bed.columns[biggest_ind] + " can the first one be removed?:\n  \
+							overlap: "+str(overlap*100)+'%\n  smalloverlap: '+str(smalloverlap*100)+'%\n  new smalloverlap: '+str(newsmalloverlap*100)+"%"):
+						continue
+					else:
+						print("  enough from small overlaps")
+						f.write("  enough from small overlaps"+'\n')
+			print(' --> enough overlap')
+			f.write(' --> enough overlap'+'\n')
+			recovered += np.sum(additionalpeaksinsec.astype(bool))
+			if merged_bed.columns[val].split('-')[0] not in markedasbad:
+				tot += peakmatrix[val].astype(int)
+			# ones that we have in the Primary but not in the secondary
+			if not lookeverywhere or len(additionalpeaksinbig) == 0:
+				tolookfor = peakmatrix[biggest_ind] == 0 if lookeverywhere else np.logical_and(
+					peakmatrix[biggest_ind] == 0, peakmatrix[val])
+				additionalpeaksinbig = findAdditionalPeaks(finalpeaks, tolookfor, bigwigfolder + findpeakpath(
+					bigwigfolder, biggest), sampling=sampling, mincov=mincov, window=window, minKL=minKL, use=use)
+				if len(additionalpeaksinbig[additionalpeaksinbig > 0]) > 0:
+					sns.kdeplot(additionalpeaksinbig[additionalpeaksinbig > 0],
+					            label=biggest, legend=True).set(xlim=(0, None))
+					print('  min,max from newly found peaks: ' +
+					      str((additionalpeaksinbig[additionalpeaksinbig > 0].min(), additionalpeaksinbig[additionalpeaksinbig > 0].max())))
+					f.write('  min,max from newly found peaks: '+str((additionalpeaksinbig[additionalpeaksinbig > 0].min(
+					), additionalpeaksinbig[additionalpeaksinbig > 0].max()))+'\n')
+
+				peakmatrix[biggest_ind] = np.logical_or(
+					peakmatrix[biggest_ind], additionalpeaksinbig)
+				tot += additionalpeaksinbig.astype(bool).astype(int)
+				recovered += np.sum(additionalpeaksinbig.astype(bool))
+			print('  we have recovered ' + str(recovered)+' peaks, equal to ' + str(100*recovered/np.sum(peakmatrix[biggest_ind])) +
+                            '% of the peaks in main replicate')
+			f.write('  we have recovered ' + str(recovered)+' peaks, equal to ' + str(100*recovered/np.sum(peakmatrix[biggest_ind])) +
+                            '% of the peaks in main replicate'+'\n')
+			if overlap < (MINOVERLAP+0.2)/1.2:
+				# we recompute to see if the overlap changed
+				newoverlap = np.sum(np.logical_and(
+					peakmatrix[val], peakmatrix[biggest_ind]))/np.sum(peakmatrix[biggest_ind])
+				smalloverlap = np.sum(np.logical_and(
+					peakmatrix[val], peakmatrix[biggest_ind]))/np.sum(peakmatrix[val])
+				if newoverlap < (MINOVERLAP+0.2)/1.2:
+					if smalloverlap < (2+MINOVERLAP)/3:
+						print("  not enough overlap to advice to merge the bams.\n  oldnew overlap: "+str(overlap*100)+'%\n  \
+							new overlap: '+str(newoverlap*100)+"%")
+						f.write("  not enough overlap to advice to merge the bams.\n  oldnew overlap: "+str(overlap*100)+'%\n  \
+							new overlap: '+str(newoverlap*100)+"%"+'\n')
+						continue
+					else:
+						print('  enough from small overlap to advice to merge the peaks')
+						f.write('  enough from small overlap to advice to merge the peaks'+'\n')
+			tomergebam.append([biggest, peakname])
+			#the quality is good enough in the end we can pop from the list if it exists
+			if tf in remove:
+				remove.remove(tf)
+		plt.title('distribution of new found peaks')
+		if saveloc:
+			plt.savefig(saveloc+tf+"_new_found_peaks_kdeplot.pdf")
+		plt.show()
+		# new distplot
+		# new correlation plot
+		ratiosofunique[tf] = len(np.argwhere(
+			peakmatrix.sum(0) == 1))/peakmatrix.shape[1]
+		if doPlot:
+			sns.pairplot(merged_bed, corner=True, diag_kind="kde",
+			             kind="reg", plot_kws={"scatter_kws": {"alpha": .05}})
+			#fig = fig.map_upper(col_nan_scatter)
+			#fig = fig.map_upper(col_nan_kde_histo)
+			plt.suptitle("correlation and distribution of peaks after recovery", y=1.08)
+			if saveloc:
+				fig.savefig(saveloc+tf+"_after_pairplot.pdf")
+			plt.show()
+			for i, val in enumerate(merged_bed):
+				unique_inval = np.logical_and(
+					np.delete(peakmatrix, i, axis=0).sum(0).astype(bool) == 0, peakmatrix[i])
+				sns.kdeplot(merged_bed[val][unique_inval], legend=True).set(xlim=(0, None))
+			plt.title("distribution of unique peaks in each replicate after recovery")
+			if saveloc:
+				plt.savefig(saveloc+tf+"_after_unique_kdeplot.pdf")
+			plt.show()
+		if len(peakmatrix.shape) > 1 and doPlot:
+			if peakmatrix.shape[0] < 7:
+				presence = []
+				for peakpres in peakmatrix:  # https://github.com/tctianchi/pyvenn
+					presence.append(set([i for i, val in enumerate(peakpres) if val == 1]))
+				title = tf + '_recovered (TOREMOVE)' if tf in remove else tf+'_recovered'
+				h.venn(presence, [i+'_BAD' if i.split('-')[0]
+                                    in markedasbad else i for i in merged_bed.columns], title=title, folder=saveloc)
+				plt.show()
+			else:
+				print('too many replicates for Venn')
+				f.write('(too many replicates for Venn)'+'\n')
+			finalpeaks = finalpeaks[np.logical_or(tot > 1, peakmatrix[biggest_ind])]
+		finalpeaks['name'] = biggest
+		finalpeaks['tf'] = tf
+		mergedpeaksdict.update({tf: finalpeaks})
+		print(str((tf, len(finalpeaks))))
+		f.write(str((tf, len(finalpeaks)))+'\n')
+	mergedpeak = pd.concat(
+		[peaks for _, peaks in mergedpeaksdict.items()]).reset_index(drop=True)
+	if doPlot:
+		df = pd.DataFrame(data=ratiosofunique, index=['percentage of unique'])
+		df['proteins'] = df.index
+		fig = sns.barplot(data=df)
+		plt.xticks(rotation=60, ha='right')
+		plt.title("ratios of unique in replicates across experiments")
+		if saveloc:
+			plt.savefig(saveloc+"All_ratios_unique.pdf")
+		plt.show()
+	f.close()
+	mergedpeak['name'] = mergedpeak.tf
+	return mergedpeak, tomergebam, remove, ratiosofunique
+
+
+def findAdditionalPeaks(peaks, tolookfor, filepath, sampling=1000, mincov=4,
+                        window=100, cov={}, minKL=8, use='max'):
+	"""
+	findAdditionalPeaks: for all peaks in A and/or B find in coverage file if zone has relative cov
+	of more than thresh then add to peak
+	if B is small and > 20% of peaks in A are found back, increase confidence and
+	flag for mergeBams
+	if < 20% don't flag for merge bam
+	f B is big and now mean non overlap < 40%, take union and flag for mergeBam else, throw B.
+
+	Args:
+	-----
+		peaks
+		tolookfor
+		filepath
+		sampling
+		mincov
+		window
+		cov
+		minKL
+		use
+	returns:
+	-------
+		np.array(bool) for each peaks in peakset, returns a binary
+	"""
+	# def poisson(k, lamb, scale): return scale * (lamb**k / factorial(k)) * np.exp(-lamb)
+
+	def KLpoisson(lamb1, lamb2): return lamb1 * \
+            np.log(lamb1 / lamb2) + lamb2 - lamb1
+
+	def poisson(k, lamb): return (lamb**k/factorial(k)) * np.exp(-lamb)
+
+	def negLogLikelihood(params, data): return - \
+            np.sum(np.log(poisson(data, params[0])))
+
+	def poissonFit(data): return float(
+		minimize(negLogLikelihood, x0=np.ones(1), args=(data,), method='Powell').x)
+	bw = pyBigWig.open(filepath)
+	res = np.zeros(len(peaks))
+	prevchrom = ''
+	lamb = {}
+	cov = {}
+	#ignore by message
+	warnings.filterwarnings("ignore", message="encountered in")
+	for i, has in enumerate(tolookfor):
+		if has:
+			val = peaks.iloc[i]
+			if val.chrom not in chroms:
+				continue
+			if val.chrom != prevchrom:
+				if val.chrom not in cov:
+					cov[val.chrom] = bw.stats(str(val.chrom))[0]
+					prevchrom = val.chrom
+					if use == 'poisson':
+						#TODO: compute on INPUT file instead
+						samples = np.zeros(window * sampling)
+						sam = np.random.rand(sampling)
+						sam = sam * (bw.chroms(str(val.chrom))-window)
+						for j, sample in enumerate(sam.astype(int)):
+							samples[j*window:(j + 1)*window] = np.nan_to_num(
+								bw.values(str(val.chrom), sample, sample + window), 0)
+						scale = np.unique(samples)[1]
+						samples = (samples/scale).astype(int)
+						lamb[val.chrom] = (poissonFit(samples), scale)
+
+			start = max([val.start - window, 0])
+			end = min(val.end + window, bw.chroms(str(val.chrom)))
+			zone = np.nan_to_num(bw.values(str(val.chrom), start, end), 0)
+			if use == 'max':
+				if max(zone) / cov[val.chrom] > mincov*1.5 or sum(zone) / (cov[val.chrom] * (end - start)) > mincov:
+					res[i] = max(zone) / cov[val.chrom]
+			elif use == 'poisson':
+				#TODO: compute -log10pvalue
+				la = poissonFit((zone/lamb[val.chrom][1]).astype(int))
+				kl = KLpoisson(la, lamb[val.chrom][0])
+				if kl > minKL:
+					res[i] = max(zone) / cov[val.chrom]  # foldchange from macs2
+
+	return res
diff --git a/genepy/epigenetics/README.md b/genepy/epigenetics/README.md
new file mode 100644
index 0000000..3ba5e91
--- /dev/null
+++ b/genepy/epigenetics/README.md
@@ -0,0 +1,42 @@
+# epigenomics
+
+Especially targeted to functions related to the analysis of epigenomics data. It has functions to read, merge, denoise, ChIP seq data.
+
+## Available functions:
+
+### chipseq.py
+
+- bigWigFrom: run the bigwig command line for a set of bam files in a folder
+- ReadRoseSuperEnhancers: reads ROSE2's output and returns its superenhancer bedfile as a pd dataframe. 
+- loadPeaks: loads 1 to many peak bedfile into one pandas dataframe.
+- simpleMergePeaks: simply merges bedfiles from peak callers. providing a concaneted dataframe of bed-like tables
+- putInBed: given a conscensus bed-like dataframe and another one, will merge the second one into the first
+- pairwiseOverlap: compute pairwise overlap and correlation on this overlap for a set of peaks mappe to a conscensus 
+- enrichment: compute pairwise enrichment and correlation for a set of peaks mappe to a conscensus 
+- fullDiffPeak: will use macs2 to call differential peak binding from two bam files and their control
+- diffPeak: calls MACS2 bdgdiff given some parameters
+- MakeSuperEnhancers: Calls super enhancer from H3K27ac with the ROSE algorithm
+- runChromHMM: runs the chromHMM algorithm
+- loadMEMEmotifs: loads motif from the output file of MEME after running fimo.
+- simpleMergeMotifs: aggregates the motifs if they overlap, into one motif file
+- substractPeaksTo: removes all peaks that are not within a bp distance to a set of loci
+
+### CREME.py
+
+The goal of creme is to be a simple, 1 function tool. It works with 1 to many sets of replicates for each pulled protein/mark.
+
+CREME takes as input a pandas dataframe. This dataframe is the concatenation of each replicates' bed files and can be loaded from a set of MACS2 bedfiles using genepy's loadPeaks function.
+
+CREME will output, amongst other thing, a dataframe representing a concatenation f bedfiles of merged replicates.
+
+find out more at __CREME.md__
+
+## highly recommended packages
+
+*This package won't contain anything that overlap with those and might use those packages for what it is doing.*
+- Bedtools
+- deepTools
+- MACS2
+- ROSE
+- MEME
+- ChromHMM
diff --git a/epigenetics/__init__.py b/genepy/epigenetics/__init__.py
similarity index 100%
rename from epigenetics/__init__.py
rename to genepy/epigenetics/__init__.py
diff --git a/genepy/epigenetics/chipseq.py b/genepy/epigenetics/chipseq.py
new file mode 100644
index 0000000..e95ab93
--- /dev/null
+++ b/genepy/epigenetics/chipseq.py
@@ -0,0 +1,888 @@
+import os
+import pandas as pd
+import numpy as np
+from genepy.utils import helper as h
+import re
+import pyBigWig
+from scipy.stats import zscore, fisher_exact
+import subprocess
+from statsmodels.stats.multitest import multipletests
+
+size = {"GRCh37": 2864785220, "GRCh38": 2913022398}
+
+cmaps = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
+        'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
+        'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']
+
+chroms = {'chr1', 'chr10', 'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19',
+'chr2', 'chr20', 'chr21', 'chr22', 'chr3', 'chr4', 'chr5', 'chr6', 'chr7', 'chr8', 'chr9', 'chrX',
+ 'chrY', '1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6',
+ '7', '8', '9', 'X', 'Y'}
+
+
+def bigWigFrom(bams, folder="", numthreads=8, genome='GRCh37', scaling=None, verbose=1):
+  """
+  run the bigwig command line for a set of bam files in a folder
+
+  Can apply some scaling and process files in parallel
+
+  Args:
+  ----
+    bams: list[str] of filepath to bam files
+    folder: str folder where the bam files would be stored (will save the bigwigs there as well, else to current path ./)
+    numthreadds: int number of threads to process in parallel
+    genome: str genome version to use (for genome size) only GRCh37 GRCh38 (see size global variable in this file)
+    scaling: list[float] a list of scaling values to apply to the bigwigs. if provided, won't scale to read counts.
+    verbose: int verbose level
+
+  """
+  if "bigwig" not in os.listdir(folder if folder else '.'):
+    os.mkdir(folder + "bigwig")
+  for i, bam in enumerate(bams):
+    in1 = folder + bam
+    out1 = folder + "bigwig/" + bam.split('/')[-1].split('.')[0] + '.bw'
+    cmd = "bamCoverage --effectiveGenomeSize " + str(size[genome]) + " -p " + str(numthreads) +\
+      " -b " + in1 + " -o " + out1
+    if scaling is not None:
+      cmd += ' --scaleFactor ' + str(scaling[i])
+    if verbose == 0:
+      cmd += ' 2> ' + bam + '.error.log'
+    res = subprocess.run(cmd, capture_output=True, shell=True)
+    if res.returncode != 0:
+      raise ValueError('issue with the command: ' + str(res.stderr))
+
+
+def ReadRoseSuperEnhancers(roseFolder, containsINPUT=True, name="MV411"):
+  """
+  reads ROSE2's output and returns its superenhancer bedfile as a pd dataframe. 
+
+  Can be multiple superenhancers from a set of HRK27ac bam files
+
+  Args:
+  -----
+    roseFolder: str folderpath to ROSE's output
+    containsINPUT: bool whether the bedfile contains INPUT signal as well
+    name: str sample name from which we computed superenhancers
+
+  Returns:
+  --------
+    a dataframe of the bed representation of the superenhancers
+  """
+  beds = os.listdir(roseFolder)
+  superenhan = pd.DataFrame()
+  for i in beds:
+    if i.endswith('_SuperEnhancers.table.txt'):
+      superenhan = superenhan.append(pd.read_csv(roseFolder+i, sep='\t', skiprows=5)).drop(columns=['enhancerRank','isSuper'])
+  data = [superenhan.columns[6]] + superenhan.columns[8:].tolist()
+  if containsINPUT:
+    inputd = superenhan.columns[7]
+  superenhan['name'] = [i.split(name)[1].split('_')[2] for i in superenhan['REGION_ID']]
+  superenhan = superenhan.rename(columns={'CONSTITUENT_SIZE':'size','REGION_ID':"id",'CHROM':'chrom','START':'start','STOP':'end',"NUM_LOCI":'num'}).replace(np.nan,0)
+  superenhan['foldchange'] = superenhan[data].sum(1)/superenhan[inputd]
+  superenhan = superenhan.drop(columns=data+[inputd])
+  return superenhan.sort_values(by=['chrom','start','end'])
+
+
+def loadPeaks(peakFile=None, peakfolder=None, isNF=True, CTFlist=[], skiprows=0):
+  """
+  loads 1 to many peak bedfile into one pandas dataframe.
+
+  all og the peaks will be concatenated into one dataframe. this function can 
+  work with jkobject|nfcore/chipseq nextflow pipelines output (isNF)
+
+  Args:
+  -----
+    peakFile: str filepath to a peak bedfile
+    peakfolder: str folderpath to a folder congaining beddfile
+    isNF: bool if the peaks come from an nfcore/chipseq nextflow pipeline (with VA_names)
+    CTFlist: list[str] only load peaks corresponding from this list of names
+    skiprows: int if there is rows to skip in the file (can be 1, often)
+
+  Returns:
+  -------
+    a bed like dataframe containing the concatenated peaks of the bed files
+  """
+  if peakfolder:
+    bindings = pd.DataFrame()
+    for folder in os.listdir(peakfolder):
+      if isNF:
+        if any(tf in folder for tf in CTFlist) or not CTFlist:
+          binding = pd.read_csv(peakfolder + folder + "/NA_peaks.narrowPeak", sep='\t', header=None) if\
+          os.exists(peakfolder + folder + "/NA_peaks.narrowPeak") else peakfolder + folder + "/NA_peaks.broadPeak"
+          binding['name'] = folder.replace('.narrowPeak', '').replace('.broadPeak','')
+          bindings = bindings.append(binding)
+      else:
+        file = folder
+        if file[-10:] in ["narrowPeak",".broadPeak"] and (any(tf in file for tf in CTFlist) or not CTFlist):
+          print('reading: '+file)
+          binding = pd.read_csv(peakfolder + file, sep='\t', header=None, skiprows=skiprows)
+          binding['name'] = file.replace('.narrowPeak', '').replace('.broadPeak','')
+          bindings = bindings.append(binding)
+  elif peakFile:
+    bindings = pd.read_csv(peakFile, sep='\t', header=None, skiprows=skiprows)
+    bindings['name'] = peakFile.split('/')[-1].split('.')[0]
+  else:
+    raise ValueError("need to provide one of peakFile or peakfolder")
+  if not len(bindings.columns)==6:
+    bindings = bindings.drop(5, 1).drop(4, 1).rename(columns={6:4, 7:5, 8:6, 9:7,})
+  bindings = bindings.rename(columns={
+    0: "chrom",
+    1: 'start',
+    2: 'end',
+    3: 'peak_number',
+    4: "foldchange",
+    5: "-log10pvalue",
+    6: "-log10qvalue",
+    7: 'relative_summit_pos'})
+  bindings = bindings.sort_values(by=["chrom", "start", "end"], axis=0)
+  bindings.start = bindings.start.astype('int')
+  bindings.end = bindings.end.astype('int')
+  if not len(bindings.columns)==6:
+    bindings['relative_summit_pos'] = bindings.relative_summit_pos.astype(
+      float) if 'relative_summit_pos' in bindings.columns else bindings.end - bindings.start
+    bindings["-log10pvalue"] = bindings["-log10pvalue"].astype('float')
+    bindings['-log10qvalue'] = bindings['-log10qvalue'].astype('float')
+    loc = bindings['relative_summit_pos'].isna()
+    bindings.loc[bindings[loc].index, 'relative_summit_pos'] = bindings[loc].end - bindings[loc].start
+    bindings.relative_summit_pos = bindings.relative_summit_pos.astype(int)
+  bindings.foldchange = bindings.foldchange.astype('float')
+  bindings = bindings.reset_index(drop=True)
+  return bindings
+
+
+def simpleMergePeaks(peaks, window=0, totpeaknumber=0, maxp=True, mergedFold="mean"):
+  """
+  simply merges bedfiles from peak callers. providing a concaneted dataframe of bed-like tables
+
+  will recompute pvalues and foldchange from that.
+
+  Args:
+  -----
+    peaks: pd.df of concatenated bed-like dataframe with a name column differentiating each peak group
+      to merge
+    window: int the max size between each peaks under which to still merge the peaks together 
+    totpeaknumber: int, set to more than 0 for peaks to be counted from this value 
+      (e.g. 10, with 100 peaks will give peak_number of 10-110)
+    maxp: bool, set to true to merge pvalue uusing the max pvalue, else will take the product of all
+    mergedFold: flag "mean"|"max"|"sum", on how to merge foldchanges
+
+  Returns:
+  -------
+    df bed-like of the merged peaks
+  """
+  peaks = peaks.sort_values(by=['chrom', 'start','end'])
+  tfs = list(set(peaks['name']))
+  mergedpeaksdict = {}
+  remove = []
+  peaknumber = 0
+  merged_bed = {
+    "chrom": [peaks.iloc[0]['chrom']],
+    "start": [peaks.iloc[0]['start']],
+    "end": [],
+    "peak_number": [peaknumber + totpeaknumber],
+    "foldchange": [],
+    "-log10pvalue": [],
+    "-log10qvalue": [],
+    "relative_summit_pos": []
+  }
+  foldchange = [peaks.iloc[0].get('foldchange', 1)]
+  log10pvalue = [peaks.iloc[0].get('-log10pvalue', 0)]
+  log10qvalue = [peaks.iloc[0].get('-log10qvalue', 0)]
+  relative_summit_pos = peaks.iloc[1].get('relative_summit_pos', peaks.iloc[0]['start'])
+  # computes overlap by extending a bit the window (100bp?) should be ~readsize
+  prev_end = peaks.iloc[0]['end']
+  prev_chrom = peaks.iloc[0]['chrom']
+  tfmerged = {a: [0] for a in tfs}
+  tfmerged[peaks.iloc[0]['name']][-1] = peaks.iloc[0].get('foldchange', 1)
+  for i, (pos, peak) in enumerate(peaks.iloc[1:].iterrows()):
+    print(str(i / len(peaks)), end="\r")
+    if prev_end + window > peak['start'] and prev_chrom == peak['chrom']:
+      # can be merged
+      if peak.get('foldchange', 1) > max(foldchange):
+        relative_summit_pos = peak.get('relative_summit_pos', peaks['start'])
+      foldchange.append(peak.get('foldchange', 1))
+      log10pvalue.append(peak.get('-log10pvalue', 0))
+      log10qvalue.append(peak.get('-log10qvalue', 0))
+
+    else:
+      # newpeak
+      for k, val in tfmerged.items():
+        val.append(0)
+      peaknumber += 1
+      merged_bed['chrom'].append(peak['chrom'])
+      merged_bed['start'].append(peak['start'])
+      merged_bed['end'].append(prev_end)
+      merged_bed['peak_number'].append(peaknumber + totpeaknumber)
+      if mergedFold=="mean":
+        merged_bed['foldchange'].append(np.mean(foldchange))
+      elif mergedFold=="max":
+        merged_bed['foldchange'].append(max(foldchange))
+      elif mergedFold=="sum":
+        merged_bed['foldchange'].append(sum(foldchange))
+      else:
+        raise ValueError("mergedFold needs to be one of:")
+      merged_bed['-log10pvalue'].append(max(log10pvalue) if maxp else np.prod(log10pvalue))
+      merged_bed['-log10qvalue'].append(max(log10qvalue) if maxp else np.prod(log10qvalue))
+      merged_bed['relative_summit_pos'].append(relative_summit_pos)
+      foldchange = [peak.get('foldchange', 1)]
+      log10pvalue = [peak.get('-log10pvalue', 0)]
+      log10qvalue = [peak.get('-log10qvalue', 0)]
+      relative_summit_pos = peak.get('relative_summit_pos', peak['start'])
+    prev_end = peak['end']
+    prev_chrom = peak['chrom']
+    tfmerged[peak['name']][-1] = peak.get('foldchange', 1)
+  merged_bed['end'].append(prev_end)
+  if mergedFold=="mean":
+    merged_bed['foldchange'].append(np.mean(foldchange))
+  elif mergedFold=="max":
+    merged_bed['foldchange'].append(max(foldchange))
+  elif mergedFold=="sum":
+    merged_bed['foldchange'].append(sum(foldchange))
+  else:
+    raise ValueError("mergedFold needs to be one of:")
+  merged_bed['-log10pvalue'].append(max(log10pvalue) if maxp else np.prod(log10pvalue))
+  merged_bed['-log10qvalue'].append(max(log10qvalue) if maxp else np.prod(log10qvalue))
+  merged_bed['relative_summit_pos'].append(relative_summit_pos)
+
+  merged_bed = pd.DataFrame(merged_bed)
+  tfmerged = pd.DataFrame(tfmerged)
+  return pd.concat([merged_bed, tfmerged], axis=1, sort=False)
+
+
+def putInBed(conscensus, value, window=10, mergetype='mean'):
+  """ 
+  given a conscensus bed-like dataframe and another one, will merge the second one into the first
+
+  Args:
+  -----
+    conscensus df[start,end,chrom] the conscensus (first one)
+    value df[start, end, chrom,foldchange] the value one (second one)
+    mergetype: flag: mean,first,last, on how to merge ttwo peaks that would fall on the same one
+    window: int on max distance of second df peaks' to the first ones, to still merge them
+      re.g. 0 is hard overap, 10000000 is nearest peak)
+
+  Returns:
+  -------
+    np.array of values of the second dataframe to add to the first one 
+    (e.g. can do df1[newcol] = returned_array)
+  """
+  conscensus = conscensus.sort_values(by=['chrom','start','end']).reset_index(drop=True)
+  value = value.sort_values(by=['chrom','start','end']).reset_index(drop=True)
+  locinvalue=0
+  loc=0
+  tot=0
+  num = []
+  res = np.zeros(len(conscensus))
+  not_end=True
+  def add(res,num,not_end,loc):
+    if len(num)>0:
+      if mergetype=='mean':
+        res[loc] = np.mean(num)
+      elif mergetype=='first':
+        res[loc]=np.sum(num)
+      elif mergetype=='first':
+        res[loc]=num[0]
+      elif mergetype=='last':
+        res[loc]=num[-1]
+      else:
+        raise ValueError('must be one of')
+      num= []
+    loc+=1
+    if loc == len(conscensus):
+      not_end=False
+    return res, num, not_end,loc
+  while not_end:
+    print(loc/len(conscensus),end="\r")
+    a = conscensus.iloc[loc]
+    b = value.iloc[locinvalue]
+    if b.chrom < a.chrom:
+      locinvalue+=1
+      if locinvalue == len(value):
+        not_end=False
+    elif b.chrom > a.chrom:
+      loc+=1
+      if loc == len(conscensus):
+        not_end=False
+    elif b.start<a.start:
+      if b.end+window>a.start:
+        tot+=1
+        num.append(b.foldchange)
+        if b.end>a.end+window:
+          res,num,not_end,loc = add(res,num,not_end,loc)
+          continue
+      locinvalue+=1
+      if locinvalue == len(value):
+        not_end=False
+    elif b.start<a.end+window:
+      tot+=1
+      num.append(b.foldchange)
+      if b.end>a.end+window:
+        res,num,not_end,loc = add(res,num,not_end,loc)
+        continue
+      locinvalue+=1
+      if locinvalue == len(value):
+        not_end=False
+    else:
+      res,num,not_end,loc = add(res,num,not_end,loc)
+  print(str(tot)+' were merged into conscensus')
+  return res
+
+
+async def pairwiseOverlap(bedfile, norm=True, bedcol=8, correct=True, docorrelation=True):
+  """
+  compute pairwise overlap and correlation on this overlap for a set of peaks mappe to a conscensus 
+
+  with each columns after the 7th one representing the signal of a given ChIP experiment
+  over this conscensus. will present overlap of row values in col values
+
+  Args:
+  ----
+    bedfile: df bed-like representing a conscensus set of peaks, and a set of values/foldchanges
+      over it
+    norm: bool wether to zscore normalise the signal or not
+    bedcol: int col where the bed information ends and signal information columns start
+    correct: bool whether to correct for fully similar lines/ columns by removing 1 on their last value
+    docorrelation: bool whether or not to compute correlation as well as enrichment
+
+  Returns:
+  -------
+    a dataframe[values_name x values_name] of % overlap of row values in col values in the consensus
+    a dataframe[values_name x values_name] of correlation of values signal over the overlaps
+  """
+  if correct:
+    print("we will be correcting for fully similar lines/ columns by removing 1 on their last value")
+  dat = bedfile[bedfile.columns[bedcol:]].values
+  correlation = np.ones((dat.shape[1],dat.shape[1]))
+  overlap = np.ones((dat.shape[1],dat.shape[1]))
+  for i, col in enumerate(dat.T):
+    #pdb.set_trace()
+    overlapping = np.delete(dat,i,axis=1)[col!=0]
+    col = col[col!=0]
+    add=0
+    for j, val in enumerate(overlapping.T):
+      if j==i:
+        add=1
+      if docorrelation:
+        if norm and not np.isnan(zscore(col[val!=0]).sum()) and not np.isnan(zscore(val[val!=0]).sum()) or not correct:
+          correlation[i,j+add] = np.corrcoef(zscore(val[val!=0]),zscore(col)[val!=0])[0,1]
+        else:
+          tmp = np.corrcoef(val[val != 0], col[val != 0])[0, 1]
+          if np.isnan(tmp) and correct:
+            if len(col[val!=0]) == 0 or len(val[val!=0]) == 0:
+            # one has no overlap
+              correlation[i,j+add] =  0
+            else:
+              # one contains only the same value everywhere
+              col[-1]-=max(0.01,abs(np.mean(col)))
+              val[-1]-=max(0.01,abs(np.mean(val)))
+              correlation[i,j+add] = np.corrcoef(val,col)[0,1]
+      overlap[i,j+add]=len(val[val!=0])/len(col)
+  if docorrelation:
+    correlation = pd.DataFrame(data=correlation, index=bedfile.columns[bedcol:], columns=bedfile.columns[bedcol:])
+    correlation[correlation.isna()] = 0
+  overlap = pd.DataFrame(data=overlap, index=bedfile.columns[bedcol:], columns=bedfile.columns[bedcol:]).T
+  return overlap, correlation if docorrelation else None
+
+
+async def enrichment(bedfile, bedcol=8, groups=None, correct=True, okpval=10**-3):
+  """
+  compute pairwise enrichment and correlation for a set of peaks mappe to a conscensus 
+
+  with each columns after the 7th one representing the signal of a given ChIP experiment
+  over this conscensus. will present enrichment of row values in col values
+
+  Args:
+  ----
+    bedfile: df bed-like representing a conscensus set of peaks, and a set of values/foldchanges
+      over it
+    bedcol: int col where the bed information ends and signal information columns start
+    correct: bool whether to correct for multiple hypothesis testing or not
+    docorrelation: bool whether or not to compute correlation as well as enrichment
+    okpval: float max pvalue over which to set the enrichment to 0
+
+  Returns:
+  -------
+    a dataframe[values_name x values_name] of enrichment of row values in col values
+    a dataframe[values_name x values_name] of correlation of values signal over the overlaps
+  """
+  dat = bedfile[bedfile.columns[bedcol:]].values
+  prob = dat.astype(bool).sum(0)/len(dat)
+  enrichment = np.zeros((dat.shape[1] if groups is None else len(set(groups)), dat.shape[1]))
+  pvals = np.zeros(
+    (dat.shape[1] if groups is None else len(set(groups)), dat.shape[1]))
+  if groups is not None:
+    for i in set(groups):
+      overlapping = dat[groups==i]
+      for j,val in enumerate(overlapping.T):
+        # enrichment of j in i
+        e, p = fisher_exact([
+          [len(val[val != 0]), len(val[val == 0])],
+          [prob[j]*len(dat), (1-prob[j])*len(dat)]])
+        enrichment[i, j] = np.log2(e)
+        pvals[i, j] = p
+  else:
+    for i, col in enumerate(dat.T):
+      overlapping = np.delete(dat,i,axis=1)[col!=0]
+      col = col[col!=0]
+      add=0
+      for j, val in enumerate(overlapping.T):
+        if j==i:
+          add=1
+          enrichment[i,i]=0
+        e, p = fisher_exact([[len(val[val != 0]), len(val[val == 0])], [
+          prob[j+add]*len(dat), (1-prob[j+add])*len(dat)]])
+        enrichment[i, j+add] = np.log2(e)
+        pvals[i, j+add] = p
+    enrichment[i,i]=0
+  enrichment = pd.DataFrame(data=enrichment, index=bedfile.columns[bedcol:] if groups is None else set(groups), columns=bedfile.columns[bedcol:]).T
+  enrichment[enrichment==-np.inf] = -1000
+  enrichment[enrichment.isna()] = 0
+  enrichment[enrichment == np.inf] = 1000
+  if correct:
+    pvals = np.reshape(multipletests(pvals.ravel(),
+                    0.1, method="bonferroni")[1], pvals.shape)
+  pvals = pd.DataFrame(
+    data=pvals, index=bedfile.columns[bedcol:]  if groups is None else set(groups), columns=bedfile.columns[bedcol:]).T
+  enrichment[pvals>okpval] = 0
+  return enrichment, pvals
+
+
+def findAdditionalCobindingSignal(conscensus, known=None, bigwigs=[], window=100):
+  """
+  somewhat similar concept to computePeaksAt
+
+  # get pysam data
+  # ask for counts only at specific locus based on peak windows from mergedpeakset
+  # append to an array
+  # return array, normalized
+  """
+  raise NotImplementedError()
+  if known:
+    print('getting '+ str(len(peaks.tf))+' peaks. Using the peaks values directly if \
+        available and using the bigwigs otherwise.')
+    res = known.values.astype(float)
+  elif len(bigwigs)>0:
+    print('getting '+str(len(bigwigs))+' bigwigs, no peaks passed. Will compute the cobinding values\
+      across the conscensus for each bigwigs.')
+    res = np.zeros((len(bigwigs), len(conscensus)), dtype=float)
+  else:
+    raise ValueError('you need to pass a list of path to bigwigs for each/some samples')
+  for i, bw in enumerate(bigwigs):
+    if known:
+      found = False
+      for j, val in enumerate(known.tf):
+        if val in bw:
+          if found:
+            raise ValueError('found two or more matching tf for bigwig: '+str(bw))
+          found = True
+          i=j
+          break
+        if not found:
+          print('no tf found in known for tf: '+bw)
+          raise ValueError('you need to have an amount of known columns equal to your bigwigs')
+    print('doing file ' + str(bw))
+    bw = pyBigWig.open(bw)
+    for k, val in conscensus.iterrows():
+      if known:
+        if res[i][k]!=0:
+          continue
+      start = max([val.start - window, 0])
+      end = min(val.end + window, bw.chroms(str(val.chrom)))
+      res[i][k] = bw.stats(str(val.chrom), start, end)[0]
+  res = np.nan_to_num(res, 0)
+  return conscensus.join(pd.Dataframe(data=(res.T / res.max(1)).T, columns=bigwigs if not known else known.columns))
+
+
+
+def annotatePeaks():
+  """
+  get H3k27AC peaks and compute Rose
+  for each TF peak
+  assign super enhancer if within it and create a super enhancer TFgroup
+  for each peaks
+  apply similar to rose and merge TF together. (create new TFgroup)
+  for each TF groups say
+    if within super enhancer
+    if within high h3k27ac marks
+    its nearest gene (distance to it, cis/trans)
+
+
+  FIND A WAY TO FILTER TO ONLY PLACES WITH H3K27AC marks??
+
+  TAD points where most contacts on one side happened on this side.
+  specific distance. zone of most concentration of contacts for a peak region.
+  """
+  # if way == 'Closest':
+
+  # elif way == 'ClosestExpressed':
+
+  # elif way == 'ActivityByContact':
+
+  # else:
+  #     raise ValueError('needs to be oneof Closest ClosestExpressed ActivityByContact')
+  raise NotImplementedError()
+
+
+# os.system()
+
+
+def getCoLocalization():
+  """
+  for each annotations (super enhancer & TFgroups)
+  for each TF, find highest peak/meanCov. if above thresh add to localization
+  """
+  raise NotImplementedError()
+
+
+def refineGroupsWithHiC():
+  """
+  given HiC data, for each loops (should be less than X Mb distance. create a Xkb zone around both sides
+  and find + merge each TF/TFgroup at this location)
+  """
+  raise NotImplementedError()
+
+
+async def fullDiffPeak(bam1, bam2, control1, size=None, control2=None, scaling=None, directory='diffData/',
+         res_directory="diffPeaks/", isTF=False, compute_size=True, pairedend=True):
+  """
+  will use macs2 to call differential peak binding from two bam files and their control
+
+  one can also provide some spike in scaling information
+
+  Args:
+  -----
+  bam1: str bamfilepath
+  bam2: str bamfilepath
+  control1: str control bam filepath (INPUT/IGG)
+  control2: str control bam filepath (INPUT/IGG)
+  scaling: tuple(float1,float2) of scaling of each bams 
+  size: int the read length (extsize parameter)
+  directory: str the directory where to save the produced data
+  res_directory: str the directory where to save the results
+  isTF: bool true if TF else false (will change the size)
+  compute_size: bool whether to compute the extsize value with `macs2 predictd`
+  pairedend: bool if paired end
+
+  Returns:
+  -------
+    str the log of the macs2 bdgdiff command
+  """
+  print("doing diff from " + bam1 + " and " + bam2)
+  if scaling is not None:
+    if max(scaling) > 1:
+      raise ValueError("scalings need to be between 0-1")
+  name1 = bam1.split('/')[-1].split('.')[0]
+  name2 = bam2.split('/')[-1].split('.')[0]
+  if size is None:
+    if isTF:
+      size = 147
+    else:
+      size = 200
+  if compute_size:
+    print('computing the fragment avg size')
+    cmd = "macs2 predictd -i " + bam1
+    ret = subprocess.run(cmd, capture_output=True, shell=True)
+    size = re.findall("# predicted fragment length is (\d+)", str(ret.stderr))[0]
+    print(size)
+  else:
+    print('using default|given size')
+  pairedend = "BAMPE" if pairedend else "BAM"
+  if control2 is None:
+    control2 = control1
+  cmd1 = "macs2 callpeak -B -t " + bam1 + " -c " + control1 + " --nomodel --extsize " + str(size) + " -n " + name1 + " --outdir " + directory + " -f " + pairedend
+  cmd2 = "macs2 callpeak -B -t " + bam2 + " -c " + control2 + " --nomodel --extsize " + str(size) + " -n " + name2 + " --outdir " + directory + " -f " + pairedend
+  print('computing the scaling values')
+  ret = subprocess.run(cmd1, capture_output=True, shell=True)
+  print(ret.stderr)
+  scaling1a = int(re.findall(" after filtering in treatment: (\d+)", str(ret.stderr))[0])
+  scaling1b = int(re.findall(" after filtering in control: (\d+)", str(ret.stderr))[0])
+  scaling1 = scaling1a if scaling1a <= scaling1b else scaling1b
+  ret = subprocess.run(cmd2, capture_output=True, shell=True)
+  print(ret.stderr)
+  scaling2a = int(re.findall(" after filtering in treatment: (\d+)", str(ret.stderr))[0])
+  scaling2b = int(re.findall(" after filtering in control: (\d+)", str(ret.stderr))[0])
+  scaling2 = scaling2a if scaling2a <= scaling2b else scaling2b
+  if scaling is not None:
+    scaling1 = int(scaling1/scaling[0])
+    scaling2 = int(scaling2/scaling[1])
+  print(scaling1, scaling2)
+  return diffPeak(directory+name1+"_treat_pileup.bdg", directory+name2+"_treat_pileup.bdg",
+    directory+name1+"_control_lambda.bdg", directory+name2+"_control_lambda.bdg",
+    res_directory, scaling1, scaling2, size)
+
+
+async def diffPeak(name1, name2, control1, control2, res_directory, scaling1, scaling2, size):
+  """
+  calls MACS2 bdgdiff given some parameters
+
+  Args:
+  -----
+    name1: str bamfilepath
+    name2: str bamfilepath
+    control1: str control bam filepath (INPUT/IGG)
+    control2: str control bam filepath (INPUT/IGG)
+    res_directory: str the directory where to save the results
+    scaling1: float the scaling factor of the first bam
+    scaling2: float the scaling factor of the second bam
+    size: int the read length (extsize parameter)
+
+  Returns:
+  -------
+    str the log of the macs2 bdgdiff command
+  """
+  print("doing differential peak binding")
+  cmd = "macs2 bdgdiff --t1 " + name1 + " --c1 "
+  cmd += control1+" --t2 " + name2 +" --c2 " + control2
+  cmd += " --d1 " + str(scaling1) + " --d2 " + str(scaling2) + " -g 60 "
+  cmd += "-l " + str(size) + " --o-prefix " + name1.split('/')[-1].split('.')[0] + "_vs_"
+  cmd += name2.split('/')[-1].split('.')[0] + " --outdir " + res_directory
+  res = subprocess.run(cmd, capture_output=True, shell=True)
+  return res
+
+
+def AssignToClosestExpressed(bed,countFile,genelocFile):
+  print("the bed file and genelocFile should use the same assembly")
+  genelocFile = pd.read_csv(genelocFile,sep="\t", compression="", columns=['chrom','start','end',"name"])
+  #for val in bed.iterrows():
+  raise NotImplementedError()
+
+
+
+async def makeSuperEnhancers(MACS2bed, bamFile, outdir, baiFile=None, rosePath=".",
+  stitching_distance=None, TSS_EXCLUSION_ZONE_SIZE="2500", assembly="hg38",controlBam=None,controlBai=None):
+  """
+  Calls super enhancer from H3K27ac with the ROSE algorithm
+
+  Args:
+  ----
+    MACS2bed: str he bed filepath to a MACS2 called set of peaks for an H3K27ac bam file
+    bamFile: str bam filepath of the H3K27ac bam file
+    outdir: str folderpath of the output. has to be an absolute path or a ~/path
+    baiFile: str bam filepath of the H3K27ac bam file
+    rosePath: str path to the ROSE folder containing the rose algorithm
+    stitching_distance: int max stitching distance
+    TSS_EXCLUSION_ZONE_SIZE: int size around TSS to excllude from SEs
+    assembly: flag one of hg38, hg37...
+    controlBam: str filepath to input bam file
+    controlBai: str filepath to input bai file
+
+  Returns:
+  --------
+    a bed-like dataframe of the superenhancers
+
+  """
+  print("we are going to move your input files to "+rosePath)
+  cmd = 'cp '+MACS2bed+' '+rosePath+MACS2bed.split('/')[-1]+'.bed'
+  cmd += ' && mv '+bamFile+' '+rosePath
+  baiFile = baiFile if baiFile else bamFile[:-1]+'i'
+  cmd += ' && cp '+baiFile +' '+rosePath
+  res = subprocess.run(cmd, capture_output=True, shell=True)
+  if res.returncode != 0:
+    raise SystemError('failed moving files: '+str(res))
+  cmd = "cd "+rosePath+" && python ROSE_main.py -g "+assembly+" -i "+MACS2bed.split('/')[-1]+'.bed' + " -r " + bamFile.split('/')[-1] + " -o " + outdir
+  if TSS_EXCLUSION_ZONE_SIZE:
+    cmd+=" -t "+TSS_EXCLUSION_ZONE_SIZE
+  if controlBam:
+    os.system('mv '+controlBam+' '+rosePath)
+    controlBai = controlBai if controlBai else controlBam[:-1]+'i'
+    os.system('cp '+controlBai+' '+rosePath)
+    cmd+=" -c "+controlBam.split('/')[-1]
+  if stitching_distance:
+    cmd+=" -s "+ stitching_distance
+
+  res = subprocess.run(cmd, capture_output=True, shell=True)
+  fail = False
+  if res.returncode != 0:
+    v = 'ROSE failed:' +str(res)
+    fail = True
+  print('finished.. moving them back to their original folder')
+  cmd = 'rm '+rosePath+MACS2bed.split('/')[-1]+'.bed'
+  cmd += ' && mv '+rosePath+bamFile.split('/')[-1]+' '+bamFile
+  cmd+= ' && rm '+rosePath+baiFile.split('/')[-1]
+  if controlBam:
+    cmd += ' && mv '+rosePath+controlBam.split('/')[-1]+' '+controlBam
+    cmd += ' && rm '+rosePath+controlBai.split('/')[-1]
+  res = subprocess.run(cmd, capture_output=True, shell=True)
+  if res.returncode != 0:
+    raise SystemError('failed moving files: '+str(res))
+  if fail:
+    raise SystemError(v)
+  print('worked')
+  print(res)
+  return ReadRoseSuperEnhancers(outdir ,bool(controlBam))
+
+
+
+def runChromHMM(outdir, data, numstates=15, datatype='bed', folderPath=".", 
+chromHMMFolderpath="~/ChromHMM/", assembly="hg38",control_bam_dir=None):
+  """
+  runs the chromHMM algorithm
+
+  Args:
+  -----
+    outdir str: an existing dir where the results should be saved
+    data: df[cellname,markname,markbed|bam|bigwig, ?controlbed|bam|bigwig]
+    numstates: number of states to use
+    datatype: flag one of bed
+    folderPath: str folder where to save chromHMM's work
+    chromHMMFolderpath: str folderpath to chromHMM algorithm
+    assembly: flag one of hg38, hg37 ...
+    control_bam_dir: str directory where the control would be stored (if not given in the ddf)
+
+  Returns:
+  -------
+    A dict of bed like dataframes containing the regions of the different states
+  """
+  print("you need to have ChromHMM")
+  chromHMM = "java -mx8000M -jar "+chromHMMFolderpath+"ChromHMM.jar "
+  h.createFoldersFor(outdir+'binarized/')
+  data.to_csv(outdir+"input_data.tsv", sep='\t',index=None,header=None)
+  cmd = chromHMM
+  if datatype=="bed":
+    cmd+="BinarizeBed "
+  elif datatype=="bigwig":
+    cmd+="BinarizeSignal "
+  elif datatype=="bam":
+    cmd+="BinarizeBam "
+  else:
+    raise ValueError('you need to provide one of bam, bigwig, bed')
+  cmd+= chromHMMFolderpath+"CHROMSIZES/"+assembly+".txt "+ folderPath+" "+outdir+"input_data.tsv "+outdir+"binarized"
+  if control_bam_dir:
+    cmd+=" -c "+control_bam_dir
+  res1 = subprocess.run(cmd, capture_output=True, shell=True)
+  print(res1)
+  if res1.returncode!=0:
+    raise ValueError(str(res1.stderr))
+  cmd = chromHMM + "LearnModel -printposterior -noautoopen "
+  if len(data)<10:
+    cmd += '-init load -m '+chromHMMFolderpath+'model_15_coreMarks.txt '
+  cmd += outdir+"binarized "+outdir+" "+str(numstates)+" "+assembly
+  res2 = subprocess.run(cmd, capture_output=True, shell=True)
+  print(res2)
+  if res2.returncode!=0:
+    raise ValueError(res2.stderr)
+  ret = {}
+  for v in set(data[0]):
+    ret[v] = pd.read_csv(outdir+v+'_'+str(numstates)+'_dense.bed', sep='\t', header=None,
+      skiprows=1).drop(columns=[4,5,6,7]).rename(columns=
+      {0:'chrom',1:'start',2:'end',3:'state',8:"color"})
+  return ret
+
+
+def loadMEMEmotifs(file, tfsubset=[],motifspecies='HUMAN'):
+  """
+  loads motif from the output file of MEME after running fimo.
+
+  Args:
+  ----
+    file: str location of the fimo gff|bed file
+    tfsubset: list[str] subset of tf to load
+    motifspecies: str name of the specied concerned (should be the same one as provided to MEME)
+
+  Returns:
+  --------
+    df bed-like of each motifs across the genome
+  """
+  if file.endswith('.gff'):
+    print('converting to bed, you need to have "gfftobed" installed')
+    cmd = 'gff2bed < '+file+' > '+file+'.bed'
+    file = file+'.bed'
+    res = subprocess.run(cmd,capture_output=True, shell=True)
+    if res.returncode != 0:
+      raise ValueError('issue with the command: ' + str(res.stderr))
+    else:
+      print(res.stdout.decode("utf-8"))
+  ## What are the motifs of our CRC members in ATACseq but not in our matrix
+  merged_motif = pd.read_csv(file, sep='\t',skiprows=0,index_col=None, names=['pos',"fimo",
+    "nucleotide_motif","relStart","relEnd","pval","strand",".","data"])
+  merged_motif['tf']=[i[5:].split("_"+motifspecies)[0] for i in merged_motif.data]
+  if tfsubset:
+    merged_motif = merged_motif[merged_motif.tf.isin(tfsubset)]
+  merged_motif['chrom'] = [i.split(':')[0][3:] for i in merged_motif.index]
+  merged_motif['start'] = [i.split(':')[1].split('-')[0] for i in merged_motif.index]
+  merged_motif['end'] = [i.split(':')[1].split('-')[1] for i in merged_motif.index]
+  merged_motif = merged_motif.reset_index(drop=True)
+  merged_motif = merged_motif.rename(columns={'pos':'relStart','fimo':'relEnd','nucleotide_motif':'pos',
+    'relStart':'pval','relEnd':'strand','pval':'fimo','strand':'motif'})
+  merged_motif['motif'] = [i.split('sequence=')[1].split(';')[0] for i in merged_motif.data]
+  merged_motif['p_val'] = [i.split('pvalue=')[1].split(';')[0] for i in merged_motif.data]
+  merged_motif['q_val'] = [i.split('qvalue=')[1].split(';')[0] for i in merged_motif.data]
+  merged_motif = merged_motif.drop(columns=['pos','.','fimo','data'])
+  merged_motif = merged_motif[merged_motif.columns[[6,7,8,0,1,3,2,9,10,5,4]]]
+  merged_motif = merged_motif.sort_values(by=['chrom','start','end']).reset_index(drop=True)
+  return merged_motif
+
+
+def simpleMergeMotifs(motifs, window=0):
+  """
+  aggregates the motifs if they overlap, into one motif file
+
+  Args:
+  ----
+    motifs: df bed-like of motifs locations
+    window: int maxsize around motif for which to still merge
+
+  Returns:
+  --------
+    df bedlike of merged motif:
+    df bedlike of motifs that were not merged as they were different but still overlapping
+  """
+  if type(motifs) is list:
+    motifs = pd.concat(motifs)
+  motifs = motifs.sort_values(by=['chrom', 'start'])
+  toremove = []
+  issues = []
+  prevmotif = motifs.iloc[0]
+  for i, (pos, motif) in enumerate(motifs.iloc[1:].iterrows()):
+    print(str(i / len(motifs)), end="\r")
+    if prevmotif['end'] + window > motif['start'] and prevmotif['chrom'] == motif['chrom']:
+      # can be merged
+      if motif['tf']!= prevmotif['tf'] or motif['motif'] != prevmotif['motif']:
+        print('found different motifs overlapping')
+        issues.extend([motif,prevmotif])
+      else:
+        toremove.append(pos)
+    prevmotif = motif
+  motifs = motifs.drop(index=toremove).reset_index(drop=True)
+  issues = pd.concat(issues)
+  return motifs, issues
+
+
+def substractPeaksTo(peaks,loci, bp=50):
+  """
+  removes all peaks that are not within a bp distance to a set of loci
+
+  Args:
+  ----
+    peaks: a bed file df with a chrom,start, end column at least
+    loci: a df witth a chrom & loci column
+    bp: the max allowed distance to the loci
+
+  Returns:
+  -------
+    all the peaks that are within this distance
+  """
+  i=0
+  j=0
+  keep=[]
+  bp=50
+  while j<len(peaks) and i<len(loci):
+    h.showcount(j,len(peaks))
+    if peaks.loc[j].chrom > loci.loc[i].chrom:
+      i+=1
+      continue
+    if peaks.loc[j].chrom < loci.loc[i].chrom:
+      j+=1
+      continue
+    if peaks.loc[j].start - bp > loci.loc[i].loci:
+      i+=1
+      continue
+    if peaks.loc[j].end + bp< loci.loc[i].loci:
+      j+=1
+      continue
+    if peaks.loc[j].end + bp >= loci.loc[i].loci and peaks.loc[j].start - bp <= loci.loc[i].loci:
+      keep.append(j)
+      j+=1
+  return peaks.loc[set(keep)]
+
+  
diff --git a/genepy/epigenetics/plot.py b/genepy/epigenetics/plot.py
new file mode 100644
index 0000000..c3537c0
--- /dev/null
+++ b/genepy/epigenetics/plot.py
@@ -0,0 +1,285 @@
+import pandas as pd
+import numpy as np
+import seaborn as sns
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from genepy.epigenetics import chipseq as chip
+
+def plotAverageOfSamples(samples, folder="", showAll=False, maxv=None, minv=None):
+  res = [] 
+  plt.figure()
+  plt.ylim(minv,maxv)
+  for sample in samples:
+    data = pd.read_csv(sample, sep='\t', skiprows=1, header=None, names=['chr', 'start', 'end', 'name', "foldchange","."]+list(range(600)))
+    r = data[list(range(600))].mean().tolist()
+    res.append(r)
+    if showAll:
+      sns.lineplot(data=np.array(r), color="#BFBFFF")
+  sns.lineplot(data=np.array(res).mean(0))
+  if folder:
+    plt.savefig(folder+"_averageofsamples.pdf", color="#1F1FFF")
+  return res
+
+
+def pysam_getPeaksAt(peaks, bams, folder='data/seqs/', window=1000, numpeaks=1000, numthreads=8):
+
+	# get pysam data
+	# ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
+	# for each counts, do a rolling average (or a convolving of the data) with numpy
+	# append to an array
+	# return array, normalized
+	loaded = {}
+	res = {i: np.zeros((len(peaks), window * 2)) for i in bams}
+	peaks = peaks.sort_values(by="foldchange", ascending=False).iloc[:numpeaks]
+	peaks.chrom = peaks.chrom.astype(str)
+	for val in bams:
+		loaded.update({val: pysam.AlignmentFile(
+            folder + val, 'rb', threads=numthreads)})
+	for k, bam in loaded.items():
+		for num, (i, val) in enumerate(peaks.iterrows()):
+			print(int(num / len(peaks)), end='\r')
+			center = int((val['start'] + val['end']) / 2)
+			for pileupcolumn in bam.pileup(val['chrom'], start=center - window,
+                                  stop=center + window, truncate=True):
+				res[k][num][pileupcolumn.pos - (center - window)] = pileupcolumn.n
+	fig, ax = plt.subplots(1, len(res))
+	for i, (k, val) in enumerate(res.items()):
+		sns.heatmap(val, ax=ax[i])
+		ax[i].set_title(k.split('.')[0])
+	fig.show()
+	return res, fig
+
+
+def bedtools_getPeaksAt(peaks, bams, folder='data/seqs/', window=1000, numpeaks=1000, numthreads=8):
+	"""
+	get pysam data
+	ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
+	for each counts, do a rolling average (or a convolving of the data) with numpy
+	append to an array
+	return array, normalized
+	"""
+	loaded = {}
+	center = [int((val['start'] + val['end']) / 2) for k, val in peaks.iterrows()]
+	peaks['start'] = [c - window for c in center]
+	peaks['end'] = [c + window - 1 for c in center]
+	peaks[peaks.columns[:3]].sort_values(by=['chrom', 'start']).to_csv(
+            'temp/peaks.bed', sep='\t', index=False, header=False)
+	bedpeaks = BedTool('temp/peaks.bed')
+
+	fig, ax = plt.subplots(1, len(bams))
+	peakset = peaks["foldchange"].values.argsort()[::-1][:numpeaks]
+	for i, val in enumerate(bams):
+		coverage = BedTool(folder + val).intersect(bedpeaks).genome_coverage(bga=True, split=True)\
+			.intersect(bedpeaks).to_dataframe(names=['chrom', 'start', 'end', 'coverage'])
+		cov = np.zeros((len(peaks), window * 2), dtype=int)
+		j = 0
+		pdb.set_trace()
+		for i, (k, val) in enumerate(peaks.iterrows()):
+			print(i / len(peaks), end='\r')
+			while coverage.iloc[j].start > val.start:
+				j -= 1
+			while coverage.iloc[j].start < val.end:
+				cov[i][coverage.iloc[j].start - val.start:coverage.iloc[j].end - val.start] =\
+					coverage.iloc[j].coverage
+				j += 1
+		sns.heatmap(coverage, ax=ax[i])
+		ax[i].set_title(val.split('.')[0])
+	fig.show()
+	return None, fig
+
+
+def makeProfiles(matx=[], folder='', matnames=[], title='',
+                 name='temp/peaksat.pdf', refpoint="TSS", scale=None,
+                 sort=False, withDeeptools=True, cluster=1, vmax=None, vmin=None, overlap=False,
+                 legendLoc=None):
+	if withDeeptools:
+		if not (len(matnames) == 2 and len(matx) == 2):
+			raise ValueError('you need two mat.gz files and two names')
+		h.createFoldersFor(name)
+		cmd = 'computeMatrixOperations relabel -m '
+		cmd += matx[0] + ' -o '+matx[0]+' --groupLabels '+matnames[0]
+		cmd += ' && computeMatrixOperations relabel -m '
+		cmd += matx[1] + ' -o '+matx[1]+' --groupLabels '+matnames[1]
+		cmd += ' && computeMatrixOperations rbind -m '
+		cmd += matx[0] + ' ' + matx[1] + " -o " + \
+                    '.'.join(name.split('.')[:-1]) + ".gz"
+		cmd += ' && plotProfile'
+		cmd += " --matrixFile " + '.'.join(name.split('.')[:-1]) + ".gz"
+		cmd += " --outFileName " + name
+		cmd += " --refPointLabel " + refpoint
+		if vmax is not None:
+			cmd += " -max "+str(vmax)
+		if vmin is not None:
+			cmd += " -min "+str(vmin)
+		if cluster > 1:
+			cmd += " --perGroup --kmeans "+str(cluster)
+		if legendLoc:
+			cmd += " --legendLocation "+legendLoc
+		if title:
+			cmd += " --plotTitle " + title
+		data = subprocess.run(cmd, shell=True, capture_output=True)
+		print(data)
+
+
+def getPeaksAt(peaks, bigwigs, folder='', bigwignames=[], peaknames=[], window=1000, title='', numpeaks=4000, numthreads=8,
+               width=5, length=10, torecompute=False, name='temp/peaksat.pdf', refpoint="TSS", scale=None,
+               sort=False, withDeeptools=True, onlyProfile=False, cluster=1, vmax=None, vmin=None, overlap=False,
+               legendLoc=None):
+	"""
+	get pysam data
+	ask for counts only at specific locus based on windows from center+-size from sorted MYC peaks
+	for each counts, do a rolling average (or a convolving of the data) with numpy
+	append to an array
+	return array, normalized
+	"""
+	if withDeeptools:
+		if isinstance(peaks, pd.DataFrame):
+			peaks = 'peaks.bed '
+			peaks.to_csv('peaks.bed', sep='\t', index=False, header=False)
+		elif type(peaks) == list:
+			pe = ''
+			i = 0
+			for n, p in enumerate(peaks):
+				if 20 < int(os.popen('wc -l ' + p).read().split(' ')[0]):
+					pe += p + ' '
+				elif len(peaknames) > 0:
+					peaknames.pop(n-i)
+					i += 1
+			peaks = pe
+		elif type(peaks) == str:
+			peaks += ' '
+		else:
+			raise ValueError(' we dont know this filetype')
+		if type(bigwigs) is list:
+			pe = ''
+			for val in bigwigs:
+				pe += folder + val + ' '
+			bigwigs = pe
+		else:
+			bigwigs = folder + bigwigs + ' '
+		h.createFoldersFor(name)
+		cmd = ''
+		if not os.path.exists('.'.join(name.split('.')[:-1]) + ".gz") or torecompute:
+			cmd += "computeMatrix reference-point -S "
+			cmd += bigwigs
+			cmd += " --referencePoint "+refpoint
+			cmd += " --regionsFileName " + peaks
+			cmd += " --missingDataAsZero"
+			cmd += " --outFileName " + '.'.join(name.split('.')[:-1]) + ".gz"
+			cmd += " --upstream " + str(window) + " --downstream " + str(window)
+			cmd += " --numberOfProcessors " + str(numthreads) + ' && '
+		cmd += "plotHeatmap" if not onlyProfile else 'plotProfile'
+		if type(name) is list:
+			if not onlyProfile:
+				raise ValueError('needs to be set to True, can\'t average heatmaps')
+			cmd += " --matrixFile " + '.gz '.join(name) + ".gz"
+			if average:
+				cmd += "--averageType mean"
+		else:
+			cmd += " --matrixFile " + '.'.join(name.split('.')[:-1]) + ".gz"
+		cmd += " --outFileName " + name
+		cmd += " --refPointLabel " + refpoint
+		if vmax is not None:
+			cmd += " -max "+str(vmax)
+		if vmin is not None:
+			cmd += " -min "+str(vmin)
+		if cluster > 1:
+			cmd += " --perGroup --kmeans "+str(cluster)
+		if overlap:
+			if onlyProfile:
+				cmd += " --plotType overlapped_lines"
+			else:
+				raise ValueError("overlap only works when onlyProfile is set")
+		if legendLoc:
+			cmd += " --legendLocation "+legendLoc
+
+		if len(peaknames) > 0:
+			pe = ''
+			for i in peaknames:
+				pe += ' ' + i
+			cmd += " --regionsLabel" + pe
+		if type(bigwigs) is list:
+			if len(bigwignames) > 0:
+				pe = ''
+				for i in bigwignames:
+					pe += ' "' + i + '"'
+				cmd += " --samplesLabel" + pe
+		if title:
+			cmd += " --plotTitle '"+title+"'"
+		data = subprocess.run(cmd, shell=True, capture_output=True)
+		print(data)
+	else:
+		if 'relative_summit_pos' in peaks.columns:
+			center = [int((val['start'] + val['relative_summit_pos']))
+                    for k, val in peaks.iterrows()]
+		else:
+			center = [int((val['start'] + val['end']) / 2)
+                    for k, val in peaks.iterrows()]
+		pd.set_option('mode.chained_assignment', None)
+		peaks['start'] = [c - window for c in center]
+		peaks['end'] = [c + window for c in center]
+		fig, ax = plt.subplots(1, len(bigwigs), figsize=[
+		                       width, length], title=title if title else 'Chip Heatmap')
+		if sort:
+			peaks = peaks.sort_values(by=["foldchange"], ascending=False)
+		if numpeaks > len(peaks):
+			numpeaks = len(peaks) - 1
+		cov = {}
+		maxs = []
+		for num, bigwig in enumerate(bigwigs):
+			bw = pyBigWig.open(folder + bigwig)
+			co = np.zeros((numpeaks, window * 2), dtype=int)
+			scale = scale[bigwig] if scale is dict else 1
+			for i, (k, val) in enumerate(peaks.iloc[:numpeaks].iterrows()):
+				try:
+					co[i] = np.nan_to_num(bw.values(str(val.chrom), val.start, val.end), 0)
+				except RuntimeError as e:
+					print(str(val.chrom), val.start, val.end)
+					pass
+			cov[bigwig] = co
+			maxs.append(co.max())
+		for num, bigwig in enumerate(bigwigs):
+			sns.heatmap(cov[bigwig] * scale, ax=ax[num], vmax=max(maxs), yticklabels=[], cmap=cmaps[num],
+                    cbar=True)
+			ax[num].set_title(bigwig.split('.')[0])
+		fig.subplots_adjust(wspace=0.1)
+		fig.show()
+		fig.savefig(name)
+		return cov, fig
+
+
+def andrew(groups, merged, annot, enr=None, pvals=None, cols=8, precise=True, title = "sorted clustermap of cobindings clustered", folder="", rangeval=4, okpval=10**-3, size=(20,15),vmax=3, vmin=0):
+	if enr is None or pvals is None:
+		enr, pvals = chip.enrichment(merged, groups=groups)
+	rand = np.random.choice(merged.index,5000)
+	subgroups = groups[rand]
+	sorting = np.argsort(subgroups)
+	redblue = cm.get_cmap('RdBu_r',256)
+	subenr = enr.iloc[annot-cols:]
+	subenr[subenr>rangeval]=rangeval
+	subenr[subenr<-rangeval]=-rangeval
+	subenr = subenr/rangeval
+	data = []
+	#colors = []
+	impv = pvals.values
+	for i in subgroups[sorting]:
+		#colors.append(viridis(i))
+		a = redblue((128+(subenr[i]*128)).astype(int)).tolist()
+		for j in range(len(a)):
+			a[j] = [1.,1.,1.,1.] if impv[j,i] > okpval else a[j]
+		data.append(a)
+	data = pd.DataFrame(data=data,columns=list(subenr.index),index= rand[sorting])
+	#data["clusters"]  = colors
+	
+	a = np.log2(1.01+merged[merged.columns[cols:annot]].iloc[rand].iloc[sorting].T)
+	if not precise:
+		for i in set(groups):
+			e = a[a.columns[subgroups[sorting]==i]].mean(1)
+			e = pd.DataFrame([e for i in range((subgroups[sorting]==i).sum())]).T
+			a[a.columns[subgroups[sorting]==i]] = e
+	
+	fig = sns.clustermap(a, vmin=vmin, vmax=vmax, figsize=size, z_score=0, colors_ratio=0.01, col_cluster=False,col_colors=data, xticklabels=False)
+	fig.ax_col_dendrogram.set_visible(False)
+	fig.fig.suptitle(title)
+	fig.savefig(folder + str(len(set(groups))) + '_clustermap_cobinding_enrichment_andrewplot.pdf')
+	plt.show()
\ No newline at end of file
diff --git a/genepy/epigenetics/rose b/genepy/epigenetics/rose
new file mode 160000
index 0000000..3acc207
--- /dev/null
+++ b/genepy/epigenetics/rose
@@ -0,0 +1,1 @@
+Subproject commit 3acc2079498742598818687fcd50307a0d90783b
diff --git a/google/README.md b/genepy/google/README.md
similarity index 100%
rename from google/README.md
rename to genepy/google/README.md
diff --git a/google/__init__.py b/genepy/google/__init__.py
similarity index 100%
rename from google/__init__.py
rename to genepy/google/__init__.py
diff --git a/google/gcp.py b/genepy/google/gcp.py
similarity index 96%
rename from google/gcp.py
rename to genepy/google/gcp.py
index f470068..30b119d 100644
--- a/google/gcp.py
+++ b/genepy/google/gcp.py
@@ -7,12 +7,11 @@
 import dalmatian as dm
 import numpy as np
 import os
-import ipdb
 import pdb
 import subprocess
 import signal
 import re
-from GenePy.utils import helper as h
+from genepy.utils import helper as h
 
 
 def list_blobs_with_prefix(bucket_name, prefix, delimiter=None):
@@ -179,7 +178,7 @@ def rmFiles(files, group=50, add='', dryrun=True):
                 break
 
 
-def recoverFiles(files):
+def recoverFiles(files, cores=1):
     """
     recover a set of files in parallel that were erased 
 
@@ -194,7 +193,18 @@ def recoverFiles(files):
     h.parrun(cmd, cores=cores)
 
 
-def patternRN(rename_dict, location, wildcards, types=[], dryrun=True, check_dependencies=True, cores=1):
+def folderRN(gspath, newpath, cores=1):
+    """
+    """
+    lis = lsFiles([gspath])
+    if lis!=0:
+        h.parrun(['gsutil -m mv ' + val + " " + newpath for val in lis], cores=cores)
+    else:
+        raise ValueError('no such folder')
+        
+
+def patternRN(rename_dict, location, wildcards, types=[], dryrun=True, 
+check_dependencies=True, cores=1):
     """
     rename/move a bunch of GCP objects found in some specific places
 
@@ -211,7 +221,6 @@ def patternRN(rename_dict, location, wildcards, types=[], dryrun=True, check_dep
         test: if test, just shows the command but does not run it
         cores:  cores tells on how many processor to parallelize the tas#k
     """
-    r = 0
     val = []
     for k, v in rename_dict.items():
         val.append(v)
diff --git a/google/good-retention.json b/genepy/google/good-retention.json
similarity index 100%
rename from google/good-retention.json
rename to genepy/google/good-retention.json
diff --git a/google/google_sheet.py b/genepy/google/google_sheet.py
similarity index 100%
rename from google/google_sheet.py
rename to genepy/google/google_sheet.py
diff --git a/google/gsheet_upload.py b/genepy/google/gsheet_upload.py
similarity index 100%
rename from google/gsheet_upload.py
rename to genepy/google/gsheet_upload.py
diff --git a/mutations/README.md b/genepy/mutations/README.md
similarity index 100%
rename from mutations/README.md
rename to genepy/mutations/README.md
diff --git a/mutations/__init__.py b/genepy/mutations/__init__.py
similarity index 100%
rename from mutations/__init__.py
rename to genepy/mutations/__init__.py
diff --git a/rna/README.md b/genepy/rna/README.md
similarity index 100%
rename from rna/README.md
rename to genepy/rna/README.md
diff --git a/rna/__init__.py b/genepy/rna/__init__.py
similarity index 95%
rename from rna/__init__.py
rename to genepy/rna/__init__.py
index a1e509e..2b38300 100644
--- a/rna/__init__.py
+++ b/genepy/rna/__init__.py
@@ -7,8 +7,8 @@
 from matplotlib import pyplot as plt
 from bokeh.palettes import *
 from bokeh.plotting import *
-from GenePy.rna import pyDESeq2
-from GenePy.utils import helper as h
+from genepy.rna import pyDESeq2
+from genepy.utils import helper as h
 import pdb
 import os
 import seaborn as sns
@@ -97,7 +97,7 @@ def convertGenes(listofgenes, from_idtype="ensembl_gene_id", to_idtype="symbol")
     return(renamed, not_parsed)
 
 
-def getSpikeInControlScales(refgenome, fastq=None, fastQfolder='', mapper='bwa', pairedEnd=False, cores=1,
+async def getSpikeInControlScales(refgenome, fastq=None, fastQfolder='', mapper='bwa', pairedEnd=False, cores=1,
                             pathtosam='samtools', pathtotrim_galore='trim_galore', pathtobwa='bwa',
                             totrim=True, tomap=True, tofilter=True, results='res/', toremove=False):
     """
@@ -210,7 +210,7 @@ def getSpikeInControlScales(refgenome, fastq=None, fastQfolder='', mapper='bwa',
     return norm, mapped,  # unique_mapped
 
 
-def GSEAonExperiments(data, experiments, res={}, savename='', scaling=[], geneset='GO_Biological_Process_2015',
+async def GSEAonExperiments(data, experiments, res={}, savename='', scaling=[], geneset='GO_Biological_Process_2015',
                       cores=8, cleanfunc=lambda i: i.split('(GO')[0]):
     """
 
@@ -337,13 +337,15 @@ def runERCC(ERCC, experiments, featurename="Feature", issingle=False, dilution=1
         d.update({
             featurename: 'Feature'
         })
+        cont = {}
         for i in cols:
             if val + '-' in i:
                 e += 1
                 d.update({i: val.split('_')[-1] + '_' + str(e)})
             if control + "-" in i:
                 c += 1
-                d.update({i: control + "_" + str(c)})
+                cont.update({i: control + "_" + str(c)})
+        d.update(cont)
         a = ERCC[list(d.keys())].rename(columns=d)
         a.to_csv('/tmp/ERCC_estimation.csv', index=None)
         val = val.split('_')[-1]
@@ -377,9 +379,9 @@ def runERCC(ERCC, experiments, featurename="Feature", issingle=False, dilution=1
 
         ipython.magic("R -o rm rm <- exDat$Results$r_m.res$r_m.mn")
         ipython.magic("R -o se se <- exDat$Results$r_m.res$r_m.mnse")
-        ipython.magic("R write.csv(c(rm,se), file = '/tmp/GenePy_ercc.csv')")
+        ipython.magic("R write.csv(c(rm,se), file = '/tmp/genepy_ercc.csv')")
         ipython.magic("R print(se,rm)")
-        l = h.fileToList("/tmp/GenePy_ercc.csv")
+        l = h.fileToList("/tmp/genepy_ercc.csv")
         res[val] = (float(l[1][4:]), float(l[2][4:]))
     for i, v in res.items():
         if abs(v[0]) > v[1]:
@@ -515,18 +517,18 @@ def DESeqSamples(data, experiments, scaling=None, keep=True, rescaling=None, res
                               data=np.array([contr, cond]).T)
         design.index = design.index.astype(str).str.replace('-', '.')
         deseq = pyDESeq2.pyDESeq2(count_matrix=data, design_matrix=design,
-                                  design_formula='~DMSO + Target', gene_column="gene_id")
+                                  design_formula='~Target', gene_column="gene_id")
         if type(scaling) is bool:
-            print("scaling using ERCC")
+            print("  scaling using ERCC")
             if scaling:
                 deseq.run_estimate_size_factors(
                     controlGenes=data.gene_id.str.contains(spikecontrolscontain))
         elif type(scaling) is list or type(scaling) is set:
-            print("scaling using a gene set")
+            print("  scaling using a gene set")
             deseq.run_estimate_size_factors(controlGenes=data.gene_id.isin(scaling))
         elif type(scaling) is dict:
             if val in scaling:
-                print("auto scaling from ERCCdashboard mean/std values")
+                print("  auto scaling from ERCCdashboard mean/std values")
                 if abs(scaling[val][0]) > threshforscaling*scaling[val][1]:
                     print("  estimating sizeFactors for this one")
                     deseq.run_estimate_size_factors(
@@ -551,16 +553,16 @@ def DESeqSamples(data, experiments, scaling=None, keep=True, rescaling=None, res
     return results
 
 
-def gsva(data, geneset_file, pathtoGenePy, method='ssgsea'):
+async def gsva(data, geneset_file, pathtogenepy, method='ssgsea'):
   print('you need to have R installed with GSVA and GSEABase library installed')
-  data.to_csv('/tmp/data_GenePyhelper_gsva.csv')
-  cmd = "Rscript "+pathtoGenePy + \
-      "/rna/ssGSEA.R /tmp/data_GenePyhelper_gsva.csv " + geneset_file + " " + method
+  data.to_csv('/tmp/data_genepyhelper_gsva.csv')
+  cmd = "Rscript "+pathtogenepy + \
+      "/rna/ssGSEA.R /tmp/data_genepyhelper_gsva.csv " + geneset_file + " " + method
   res = subprocess.run(cmd, shell=True, capture_output=True)
   if res.returncode != 0:
     raise ValueError('issue with the command: ' + str(res))
   print(res)
-  res = pd.read_csv("/tmp/res_GenePy_ssGSEA.tsv", sep='\t')
+  res = pd.read_csv("/tmp/res_genepy_ssGSEA.tsv", sep='\t')
   return res
 
 
@@ -658,19 +660,22 @@ def filterRNAfromQC(rnaqc, folder='tempRNAQCplot/', plot=True, qant1=0.07, qant3
         plt.show()
         plot.get_figure().savefig(folder+'failed_qc.pdf')
 
-        for val in rnaqc.index:
+        num_cols = 10
+        num_rows = math.ceil(len(rnaqc)/num_cols)
+        _, axes = plt.subplots(num_rows, num_cols, figsize=(20, num_rows*2))
+        for val_idx, val in enumerate(rnaqc.index):
+            ax = axes.flatten()[val_idx]
             qc = rnaqc.loc[val]
-            boxplot = sns.violinplot(y=qc)
+            sns.violinplot(y=qc, ax=ax)
             q1 = qc.quantile(qant1)
             q3 = qc.quantile(qant3)
             outlier_top_lim = q3 + 1.5 * (q3 - q1)
             outlier_bottom_lim = q1 - 1.5 * (q3 - q1)
             for k, v in qc[(qc < outlier_bottom_lim) | (qc > outlier_top_lim)].iteritems():
-                plt.text(0.05, v, k, ha='left', va='center',
+                ax.text(0.05, v, k, ha='left', va='center',
                          color='red' if k in a else 'black')
-            plt.show()
-            plt.savefig(folder +
-                        val.replace(' ', '_').replace('/', '_')+'.pdf')
+        plt.tight_layout()
+        plt.savefig('{}/qc_metrics.pdf'.format(folder), bbox_inches='tight')
     return res
 
 
diff --git a/rna/pyDESeq2.py b/genepy/rna/pyDESeq2.py
similarity index 100%
rename from rna/pyDESeq2.py
rename to genepy/rna/pyDESeq2.py
diff --git a/rna/ssGSEA.R b/genepy/rna/ssGSEA.R
similarity index 87%
rename from rna/ssGSEA.R
rename to genepy/rna/ssGSEA.R
index 12e2dfc..bb5d45b 100644
--- a/rna/ssGSEA.R
+++ b/genepy/rna/ssGSEA.R
@@ -13,4 +13,4 @@ gsc_obj <- GSEABase::getGmt(gmtfile,
 collectionType = GSEABase::BroadCollection(),
 geneIdType = GSEABase::EntrezIdentifier())
 gsea <- GSVA::gsva(mat, gsc_obj, method = method)
-write.table(gsea, file = "/tmp/res_GenePy_ssGSEA.tsv", sep = '\t', quote = F)
+write.table(gsea, file = "/tmp/res_genepy_ssGSEA.tsv", sep = '\t', quote = F)
diff --git a/sequencing/README.md b/genepy/sequencing/README.md
similarity index 100%
rename from sequencing/README.md
rename to genepy/sequencing/README.md
diff --git a/sequencing/__init__.py b/genepy/sequencing/__init__.py
similarity index 97%
rename from sequencing/__init__.py
rename to genepy/sequencing/__init__.py
index 4d0c7dc..559f7e6 100644
--- a/sequencing/__init__.py
+++ b/genepy/sequencing/__init__.py
@@ -9,7 +9,8 @@
 
 import pandas as pd
 import numpy as np
-
+from genepy.utils import helper as h
+from genepy.google import gcp
 from taigapy import TaigaClient
 tc = TaigaClient()
 
@@ -48,7 +49,7 @@ def fromGTF2BED(gtfname, bedname, gtftype='geneAnnot'):
         prevname = ''
         newbed = {'chr': [], 'start': [], 'end': [], 'gene': []}
         for i, val in gtf.iterrows():
-            showcount(i, len(gtf))
+            h.showcount(i, len(gtf))
             if val['name'] == prevname:
                 newbed['end'][-1] = val['stop']
             else:
@@ -107,7 +108,7 @@ def getBamDate(bams, split='-', order="des", unknown='U'):
     return DTs
 
 
-def indexBams(bucketpath, cores=4):
+async def indexBams(bucketpath, cores=4):
     """
     given a bucket path, will index all .bam files without an associated index and return their paths
     """
@@ -191,7 +192,6 @@ def findReplicatesBams(folder, sep='-', namings='-r([0-9])', namepos=2):
 		if val[-4:] == '.bam':
 			match = re.search(namings, val)
 			if match:
-				number = match.groups()[0]
 				name = val.split(sep)[namepos]
 				if name in rep:
 					rep[name].append(val)
@@ -226,7 +226,7 @@ def pairedEnd(pairedend, folder="", numthreads=8, peaksFolder="peaks/",
 	# run the paired end pipeline
 	"""
 	print("you need to have bowtie2 installed: http://bowtie-bio.sourceforge.net/bowtie2/index.shtml")
-	for key, val in pairedend.items():
+	for _, val in pairedend.items():
 		out1 = folder + mappedFolder + val[0].split('.')[0] + ".mapped.sam"
 		in1 = folder + val[0]
 		in2 = folder + val[1]
@@ -241,10 +241,11 @@ def pairedEnd(pairedend, folder="", numthreads=8, peaksFolder="peaks/",
 		# it can take many TB so better delete
 
 
-def mergeBams(rep):
+async def mergeBams(rep):
 	"""
 	uses samtools to merge a set of replicates considered into one file
 	"""
+	in1 =''
 	for i, val in rep.items():
 		out1 = i + '.merged.bam'
 		for bam in val:
diff --git a/terra/README.md b/genepy/terra/README.md
similarity index 100%
rename from terra/README.md
rename to genepy/terra/README.md
diff --git a/terra/__init__.py b/genepy/terra/__init__.py
similarity index 98%
rename from terra/__init__.py
rename to genepy/terra/__init__.py
index 10eff85..90399af 100644
--- a/terra/__init__.py
+++ b/genepy/terra/__init__.py
@@ -10,8 +10,8 @@
 import os
 import re
 import signal
-from GenePy.utils import helper as h
-from GenePy.google import gcp
+from genepy.utils import helper as h
+from genepy.google import gcp
 import pdb
 import subprocess
 from gsheets import Sheets
@@ -41,7 +41,7 @@ def createManySubmissions(workspace, workflow, references, entity=None, expressi
   return submission_ids
 
 
-def waitForSubmission(workspace, submissions, raise_errors=True):
+async def waitForSubmission(workspace, submissions, raise_errors=True):
   """
   wrapper to create many submissions for a workflow
 
@@ -556,7 +556,7 @@ def renametsvs(workspace, wmto=None, index_func=None):
         wmto.update_sample_set(i, val.samples)
 
 
-def findBackErasedDuplicaBamteFromTerraBucket(workspace, gsfolder, bamcol="WES_bam", baicol="WES_bai"):
+async def findBackErasedDuplicaBamteFromTerraBucket(workspace, gsfolder, bamcol="WES_bam", baicol="WES_bai"):
   """
   If you have erased bam files in gcp with bai files still present and the bam files are stored elsewhere
   and their location is in a terra workspace.
@@ -614,7 +614,7 @@ def findBackErasedDuplicaBamteFromTerraBucket(workspace, gsfolder, bamcol="WES_b
       print("no data for " + str(k))
 
 
-def shareTerraBams(users, workspace, samples, bamcols=["internal_bam_filepath", "internal_bai_filepath"]):
+async def shareTerraBams(users, workspace, samples, bamcols=["internal_bam_filepath", "internal_bai_filepath"]):
   """
   will share some files from gcp with a set of users using terra as metadata repo.
 
@@ -652,7 +652,7 @@ def shareTerraBams(users, workspace, samples, bamcols=["internal_bam_filepath",
   return togiveaccess
 
 
-def shareCCLEbams(users, samples, groups=[], raise_error=True, arg_max_length=100000, bamcols=["internal_bam_filepath", "internal_bai_filepath"],
+async def shareCCLEbams(users, samples, groups=[], raise_error=True, arg_max_length=100000, bamcols=["internal_bam_filepath", "internal_bai_filepath"],
                   refsheet_url="https://docs.google.com/spreadsheets/d/1XkZypRuOEXzNLxVk9EOHeWRE98Z8_DBvL4PovyM01FE",
                   privacy_sheeturl="https://docs.google.com/spreadsheets/d/115TUgA1t_mD32SnWAGpW9OKmJ2W5WYAOs3SuSdedpX4"):
   """
@@ -682,7 +682,7 @@ def shareCCLEbams(users, samples, groups=[], raise_error=True, arg_max_length=10
   blacklisted = set(blacklist) & set(samples)
   print("we have " + str(len(blacklist)) + ' blacklisted files')
   if len(blacklisted):
-    print("these lines are blacklisted " + blacklisted)
+    print("these lines are blacklisted " + str(blacklisted))
     if raise_error:
       raise ValueError("blacklistedlines")
   if type(users) is str:
@@ -741,7 +741,7 @@ def saveConfigs(workspace, filepath):
   h.dictToFile(params, filepath + '.json')
 
 
-def cleanWorkspace(workspaceid, only=[], toleave=[], defaulttoleave=['workspace', 'scripts',
+async def cleanWorkspace(workspaceid, only=[], toleave=[], defaulttoleave=['workspace', 'scripts',
                                                                      'notebooks', 'files', 'data', 'hound', 'references', 'name', 'folder']):
   """
   removes all processing folder in a terra workspace easily
@@ -880,7 +880,7 @@ def removeFromFailedWorkflows(workspaceid, maxtime = '2020-06-10', everythingFor
                     delete_job(workspaceid,val.submission_id,a,dryrun=dryrun)
 
 
-def deleteHeavyFiles(workspaceid, unusedOnly=True):
+async def deleteHeavyFiles(workspaceid, unusedOnly=True):
     """
     deletes all files above a certain size in a workspace (that are used or unused)
 
diff --git a/utils/Datanalytics.py b/genepy/utils/Datanalytics.py
similarity index 100%
rename from utils/Datanalytics.py
rename to genepy/utils/Datanalytics.py
diff --git a/utils/README.md b/genepy/utils/README.md
similarity index 100%
rename from utils/README.md
rename to genepy/utils/README.md
diff --git a/utils/RScript.R b/genepy/utils/RScript.R
similarity index 100%
rename from utils/RScript.R
rename to genepy/utils/RScript.R
diff --git a/utils/__init__.py b/genepy/utils/__init__.py
similarity index 100%
rename from utils/__init__.py
rename to genepy/utils/__init__.py
diff --git a/utils/helper.py b/genepy/utils/helper.py
similarity index 86%
rename from utils/helper.py
rename to genepy/utils/helper.py
index b743e13..1d2977d 100644
--- a/utils/helper.py
+++ b/genepy/utils/helper.py
@@ -3,22 +3,17 @@
 # in 2019
 
 from __future__ import print_function
-from matplotlib import pyplot as plt
 import json
 import os
-import sys
 import string
 import subprocess
-
+import pandas as pd
 import pdb
-import ipdb
 import pandas as pd
 import numpy as np
 import itertools
 import random
-
-from taigapy import TaigaClient
-tc = TaigaClient()
+import xmltodict
 
 
 rename_mut = {'contig': 'chr', 'position': 'pos', 'Reference_Allele': 'ref', 'ref_allele': 'ref', 'alt_allele': 'alt',
@@ -357,38 +352,39 @@ def datetoint(dt, split='-', unknown='U', order="des"):
 
 
 def showcount(i, size):
-  """
-  pretty print of i/size%, to put in a for loop
-  """
-  print(str(1 + int(100 * (i / size))) + '%', end='\r')
+    """
+    pretty print of i/size%, to put in a for loop
+    """
+    print(str(1 + int(100 * (i / size))) + '%', end='\r')
 
 
 def combin(n, k):
-  """
-  Nombre de combinaisons de n objets pris k a k
-  outputs the number of comabination of n object taken k at a time
-  """
-  if k > n // 2:
-    k = n - k
-  x = 1
-  y = 1
-  i = n - k + 1
-  while i <= n:
-    x = (x * i) // y
-    y += 1
-    i += 1
-  return x
+    """
+    Nombre de combinaisons de n objets pris k a k
+    outputs the number of comabination of n object taken k at a time
+    """
+    if k > n // 2:
+        k = n - k
+    x = 1
+    y = 1
+    i = n - k + 1
+    while i <= n:
+        x = (x * i) // y
+        y += 1
+        i += 1
+    return x
 
 
 def dups(lst):
-  """
-	shows the duplicates in a list
-  """
-  seen = set()
-  # adds all elements it doesn't know yet to seen and all other to seen_twice
-  seen_twice = set(x for x in lst if x in seen or seen.add(x))
-  # turn the set into a list (as requested)
-  return list(seen_twice)
+    """
+        shows the duplicates in a list
+    """
+    seen = set()
+    # adds all elements it doesn't know yet to seen and all other to seen_twice
+    seen_twice = set(x for x in lst if x in seen or seen.add(x))
+    # turn the set into a list (as requested)
+    return list(seen_twice)
+
 
 def makeCombinations(size, proba):
     """
@@ -419,11 +415,12 @@ def makeCombinations(size, proba):
     sums[0] = 1-sum(list(sums.values()))
     return sums
 
+
 def closest(lst, K):
-  """
-  returns the index of the value closest to K in a lst
-  """
-  return lst[min(range(len(lst)), key = lambda i: abs(lst[i]-K))]
+    """
+    returns the index of the value closest to K in a lst
+    """
+    return lst[min(range(len(lst)), key = lambda i: abs(lst[i]-K))]
 
 
 def compareDfs(df1, df2):
@@ -446,14 +443,42 @@ def compareDfs(df1, df2):
     print('FOUND new 0s in df1: ' + str(new0s))
     return nmissmatchCols, omissmatchCols, nmissmatchInds, omissmatchInds, newNAs, new0s
 
+
 def stringifydict(res):
     """
     
     """
     a = {}
-    for k,v in res.items():
+    for k, v in res.items():
         if type(v) is dict:
             a[k] = stringifydict(v)
         else:
             a[str(k)] = v
-    return a
\ No newline at end of file
+    return a
+
+
+def readXMLs(folder=None, file=None, rename=None):
+    if file is not None:
+        if type(file) is str: 
+            print('reading 1 file')
+            files = [file]
+        else:
+            print('reading files')
+            files = file
+    if folder is not None:
+        print("reading from folder")
+        files = [i for i in os.listdir(folder) if i.endswith(".xml")]
+    df = pd.DataFrame()
+    for file in files:
+        res = []
+        a = open(file, "r").read()
+        xmldict = xmltodict.parse(a)
+        data = xmldict['Workbook']['Worksheet']['Table']["Row"]
+        for val in data:
+            res.append([v['Data']['#text'] if "#text" in v['Data']
+                        else None for v in val["Cell"]])
+        res = pd.DataFrame(data=res[2:], columns=res[0])
+        df = df.append(res)
+    if rename is not None:
+        df = df.rename(columns=rename)
+    return df
diff --git a/utils/plot.py b/genepy/utils/plot.py
similarity index 82%
rename from utils/plot.py
rename to genepy/utils/plot.py
index 39bd0ae..99ee774 100644
--- a/utils/plot.py
+++ b/genepy/utils/plot.py
@@ -20,7 +20,6 @@
 import colorcet as cc
 from PIL import Image, ImageDraw, ImageFont
 import seaborn as sns
-from GenePy.epigenetics import chipseq as chip 
 
 import pandas as pd
 from math import pi
@@ -502,6 +501,136 @@ def correlationMatrix(data, names, colors=None, pvals=None, maxokpval=10**-9, ot
         plt.show()
 
 
+def heatmap(data, colors=None, title="correlation Matrix", size=40, 
+            folder='', interactive=False, pvals=None, maxokpval=10**-9, maxval=None, minval=None):
+    """
+    Make an interactive heatmap from a dataframe using bokeh
+
+    Args:
+    -----
+      data: dataframe of int / float/ bool of size(names1*names2)
+      colors: list[int] of size(names) a color for each names (good to display clusters)
+      pvals: arraylike of int / float/ bool of size(names*val) or (names*names) with the corresponding pvalues
+      maxokpval: float threshold when pvalue is considered good. otherwise lowers the size of the square
+        until 10**-3 when it disappears
+      title: str the plot title
+      size: int the plot size
+      folder: str of folder location where to save the plot, won't save if empty
+      interactive: bool whether or not to make the plot interactive (else will use matplotlib)
+      maxval: float clamping coloring up to maxval
+      minval: float clamping coloring down to minval
+
+    Returns:
+    -------
+      the bokeh object if interactive else None
+
+    """
+    regdata = data.copy()
+    if minval is not None:
+        data[data<minval]=minval
+    if maxval is not None:
+        data[data > maxval] = maxval
+    data=data/data.max()
+    data = data.values
+    TOOLS = "hover,crosshair,pan,wheel_zoom,zoom_in,zoom_out,box_zoom,undo,redo,reset,save"
+    xname = []
+    yname = []
+    color = []
+    alpha = []
+    height = []
+    width = []
+    if pvals is not None:
+        print('we are assuming you want to display the pvals of your correlation with size')
+        regpvals = pvals.copy()
+        u = pvals<maxokpval
+        pvals[~u] = np.log10(1/pvals[~u])
+        pvals = pvals/pvals.max()
+        pvals[u]=1
+    if interactive:
+        xname = []
+        yname = []
+        color = []
+        for i, name1 in enumerate(regdata.index):
+            for j, name2 in enumerate(regdata.columns):
+                xname.append(name2)
+                yname.append(name1)
+                if pvals is not None:
+                    #import pdb;pdb.set_trace()
+                    height.append(max(0.1, min(0.9, pvals.loc[name1][name2])))
+                    color.append(cc.coolwarm[int((data[i, j]*128)+127)])
+                    alpha.append(min(abs(data[i, j]), 0.9))
+                elif other is not None:
+                    color.append(cc.coolwarm[int((data[i, j]*128)+127)])
+                    alpha.append(max(min(other[i, j], 0.9),0.1) if other[i, j]!=0 else 0)
+                else:
+                    alpha.append(min(abs(data[i, j]), 0.9))
+                if colors is not None:
+                    if type(colors) is list:
+                        if colors[i] == colors[j]:
+                            color.append(
+                                Category10[10][colors[i]])
+                        else:
+                            color.append('lightgrey')
+
+                elif pvals is None and other is None:
+                    color.append('grey' if data[i, j]
+                                 > 0 else Category20[3][2])
+        if pvals is not None:
+            width = height.copy()
+            data = dict(
+                xname=xname,
+                yname=yname,
+                colors=color,
+                alphas=alpha,
+                data=regdata.values.ravel(),
+                pvals=regpvals.values.ravel(),
+                width=width,
+                height=height
+            )
+        else:
+            data = dict(
+                xname=xname,
+                yname=yname,
+                colors=color,
+                alphas=alpha,
+                data=data.ravel()
+            )
+        tt = [('names', '@yname, @xname'), ('value', '@data')]
+        if pvals is not None:
+            tt.append(('pvals','@pvals'))
+        p = figure(title=title if title is not None else "Heatmap",
+                   x_axis_location="above", tools=TOOLS,
+                   x_range=list(reversed(regdata.columns.astype(str).tolist())), y_range=regdata.index.tolist(),
+                   tooltips=tt)
+
+        p.plot_width = 800
+        p.plot_height = 800
+        p.grid.grid_line_color = None
+        p.axis.axis_line_color = None
+        p.axis.major_tick_line_color = None
+        p.axis.major_label_text_font_size = "5pt"
+        p.axis.major_label_standoff = 0
+        p.xaxis.major_label_orientation = np.pi / 3
+
+        p.rect('xname', 'yname', width = 0.9 if not width else 'width',
+                height = 0.9 if not height else 'height', source=data,
+               color='colors', alpha='alphas', line_color=None,
+               hover_line_color='black', hover_color='colors')
+        save(p, folder + title.replace(' ', "_") + "_heatmap.html")
+        #export_png(p, filename=folder + title.replace(' ', "_") + "_correlation.png")
+        try:
+            show(p)
+        except:
+            show(p)
+        return p  # show the plot
+    else:
+        plt.figure(figsize=size)
+        plt.title('the correlation matrix')
+        plt.imshow(data)
+        plt.savefig(title + "_correlation.pdf")
+        plt.show()
+
+
 def venn(inp, names, title="venn", folder=''):
     """
     Plots a venn diagram using the pyvenn package
@@ -613,40 +742,3 @@ def SOMPlot(net, size, colnames, minweight=0.1, distq1=0.535, distq2=0.055, dist
         somnodes.loc[i, 'features'] = tot
     #interactive SOM with features with highest importance to the nodes, displayed when hovering
     bigScatter(somnodes, precomputed=True, features=True, binsize=1, title='Cobinding SOM cluster of '+str(size), folder=folder)
-
-
-def andrew(groups, merged, annot, enr=None, pvals=None, cols=8, precise=True, title = "sorted clustermap of cobindings clustered", folder="", rangeval=4, okpval=10**-3, size=(20,15),vmax=3, vmin=0):
-    if enr is None or pvals is None:
-        enr, pvals = chip.enrichment(merged, groups=groups)
-    rand = np.random.choice(merged.index,5000)
-    subgroups = groups[rand]
-    sorting = np.argsort(subgroups)
-    redblue = cm.get_cmap('RdBu_r',256)
-    subenr = enr.iloc[annot-cols:]
-    subenr[subenr>rangeval]=rangeval
-    subenr[subenr<-rangeval]=-rangeval
-    subenr = subenr/rangeval
-    data = []
-    #colors = []
-    impv = pvals.values
-    for i in subgroups[sorting]:
-        #colors.append(viridis(i))
-        a = redblue((128+(subenr[i]*128)).astype(int)).tolist()
-        for j in range(len(a)):
-            a[j] = [1.,1.,1.,1.] if impv[j,i] > okpval else a[j]
-        data.append(a)
-    data = pd.DataFrame(data=data,columns=list(subenr.index),index= rand[sorting])
-    #data["clusters"]  = colors
-    
-    a = np.log2(1.01+merged[merged.columns[cols:annot]].iloc[rand].iloc[sorting].T)
-    if not precise:
-        for i in set(groups):
-            e = a[a.columns[subgroups[sorting]==i]].mean(1)
-            e = pd.DataFrame([e for i in range((subgroups[sorting]==i).sum())]).T
-            a[a.columns[subgroups[sorting]==i]] = e
-    
-    fig = sns.clustermap(a, vmin=vmin, vmax=vmax, figsize=size, z_score=0, colors_ratio=0.01, col_cluster=False,col_colors=data, xticklabels=False)
-    fig.ax_col_dendrogram.set_visible(False)
-    fig.fig.suptitle(title)
-    fig.savefig(folder + str(len(set(groups))) + '_clustermap_cobinding_enrichment_andrewplot.pdf')
-    plt.show()
diff --git a/requirements.txt b/requirements.txt
index 044a359..90b55ca 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,24 +1,23 @@
-pandas==0.24.2
-gseapy==0.9.13
-scipy==1.0.0
-bokeh==1.2.0
-requests==2.22.0
-pytest==4.6.3
-setuptools==41.0.1
-statsmodels==0.9.0
-pysam==0.15.3
-gsheets==0.4
-matplotlib==2.2.0
-pybedtools==0.8.0
-seaborn==0.9.0
-venn==0.1.3
-rpy2==3.0.4
-numpy==1.17.2
+oauth2client>=4.1.3
+bokeh
+taigapy==2.12
+pandas
+scipy>=1.0.0
+firecloud_dalmatian>=0.0.17
+rpy2>=3.3.2
+statsmodels
+pytest
+venn
+gsheets==0.4.1
+pysam
+numpy
+seaborn
+pybedtools
+requests>=2.24.0
+matplotlib
+gseapy==0.9.18
+colorcet
 gspread==3.6.0
-firecloud_dalmatian==0.0.17
-Pillow==8.1.0
-colorcet==2.0.6
-google_api_python_client==1.12.8
-pyBigWig==0.3.17
-scikit_learn==0.24.1
-taigapy==3.2.1
+Pillow
+pyBigWig
+scikit_learn
diff --git a/setup.py b/setup.py
index 3b5bc3a..0c4ab00 100644
--- a/setup.py
+++ b/setup.py
@@ -1,17 +1,29 @@
 from setuptools import setup
+import sys
+import os
+if sys.version_info.major < 3 or sys.version_info.minor < 2:
+    raise ValueError("genepy is only compatible with Python 3.3 and below")
+if sys.version_info.minor < 5:
+    import warnings
+    warnings.warn("genepy may not function properly on Python < 3.5")
+
+os.system('git submodule init && git submodule sync')
 
 with open("README.md", 'r') as f:
     long_description = f.read()
 
 setup(
-    name='Broad-GenePy',
+    name='Broad-genepy',
     version='1.0',
     description='A useful module for any CompBio',
-    long_description = long_description,
+    long_description=long_description,
     author='Jeremie Kalfon',
     author_email='jkobject@gmail.com',
-    url="https://github.com/BroadInstitute/GenePy",
-    packages=['cell_line_mapping-master/python'],  # same as name
+    url="https://github.com/BroadInstitute/genepy",
+    packages=['genepy/cell_line_mapping-master/python/cell_line_mapper',
+              'genepy/epigenetics', 'genepy/mutations', 'genepy/google', 'genepy/sequencing/',
+              'genepy/terra', 'genepy/rna', 'genepy/utils'],
+    package_data={'genepy': ['data/*']},
     python_requires='>=3.5',
     install_requires=[
         'rpy2-bioconductor-extensions',
@@ -20,8 +32,7 @@
         'deeptools',
         ## from requirements.txt
         "bokeh",
-        "dalmatian",
-        "firecloud_dalmatian",
+        "firecloud-dalmatian",
         "google_api_python_client",
         "gsheets",
         "gspread",
@@ -36,7 +47,7 @@
         "pytest",
         "requests",
         "rpy2",
-        "scikit_learn",
+        "scikit-learn",
         "scipy",
         "seaborn",
         "setuptools",
@@ -44,13 +55,22 @@
         "typing",
         "venn",
         ],  # external packages as dependencies
+    classifiers=[
+        "Programming Language :: Python :: 3",
+        "Intended Audience :: Science/Research",
+        "Topic :: Scientific/Engineering :: Bio-Informatics",
+    ],
 )
 
+
 print("You might want to install Bowtie2, samtools, bwa and R to be able to use all functions of this package:\n\
   http://bowtie-bio.sourceforge.net/bowtie2/index.shtml\n\
   http://www.htslib.org/\n\
   https://github.com/lh3/bwa\n")
 
-print("once R is installed you need to have installed erccdashboard, GSEABase GSVA, DESeq2 to have access to aall the functions")
+print("trying to install R packages")
+os.system("R -e 'if(!requireNamespace(\"BiocManager\", quietly = TRUE)){install.packages(\"BiocManager\")};BiocManager::install(c(\"GSEABase\", \"erccdashboard\", \"GSVA\", \"DESeq2\"));'")
+print('if it did not work. please install R or check your R installation')
+print("once R is installed you need to install erccdashboard, GSEABase GSVA, DESeq2 to have access to all the functions")
 
 print("Finished!")


Test patch:
diff --git a/cell_line_mapping-master/celllinemapr/man/ccle.to.latest.Rd b/genepy/cell_line_mapping-master/celllinemapr/man/ccle.to.latest.Rd
similarity index 100%
rename from cell_line_mapping-master/celllinemapr/man/ccle.to.latest.Rd
rename to genepy/cell_line_mapping-master/celllinemapr/man/ccle.to.latest.Rd
diff --git a/cell_line_mapping-master/python/cell_line_mapper/test_mapper.py b/genepy/cell_line_mapping-master/python/cell_line_mapper/test_mapper.py
similarity index 100%
rename from cell_line_mapping-master/python/cell_line_mapper/test_mapper.py
rename to genepy/cell_line_mapping-master/python/cell_line_mapper/test_mapper.py
diff --git a/cell_line_mapping-master/python/test-data.csv b/genepy/cell_line_mapping-master/python/test-data.csv
similarity index 100%
rename from cell_line_mapping-master/python/test-data.csv
rename to genepy/cell_line_mapping-master/python/test-data.csv