import os,sys
from collections import defaultdict
import numpy as np
import optparse
import logging
from pbcore.io.align.CmpH5IO import CmpH5Reader
from pbcore.io.BasH5IO import BasH5Reader
import glob
import numpy as np
import logging
import shutil
import pickle
import math
import mbin
import motif_tools
from Bio import SeqIO

def launch():
	opts, input_files, motifs_fn = __parseArgs()
	__initLog(opts)

	motifs           = np.loadtxt(motifs_fn, dtype="str", ndmin=1)
	motifs,not_found = find_motifs_in_control(opts, motifs)
	if len(not_found)>0:
		logging.warning("")
		logging.warning("  ******************** Important *********************")
		logging.warning("  Did not find %s motifs in %s:" % (len(not_found), opts.control_pkl_name))
		for nf in not_found:
			logging.warning("       %s" % nf)
		logging.warning("  These motif(s) will be removed from further analysis.")
		logging.warning("  These %s motifs will be kept:" % len(motifs))
		for m in motifs:
			logging.warning("       %s" % m)
		logging.warning("  ****************************************************")
		logging.warning("")
	else:
		logging.info("Found entries for all %s motifs in %s" % (len(motifs), opts.control_pkl_name))


	build_profiles(opts, input_files, motifs, motifs_fn)

	print >> sys.stderr, "mBin methylation profiling has finished running. See log for details."

def write_contig_features( mbinRunner, opts ):
	"""
	Gather contig-level data generated by the pipeline and write
	to output files.
	"""
	methyl_fn = "%scontig_methyl_features.txt" % opts.prefix
	motifs_fn = "%scontig_motif_counts.txt"    % opts.prefix
	others_fn = "%scontig_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for assembled contigs" % methyl_fn)
	logging.info("   %s  \t -- Motif counts on assembled contigs" % motifs_fn)
	logging.info("   %s\t -- Coverage and length values, as well as %s-mer frequency vectors for assembled contigs" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	ref_names  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]),   dtype="str", ndmin=1)
	ref_SCp    = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_SCp"]),     dtype="float", ndmin=2)
	ref_SCp_N  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_SCp_N"]),   dtype="int", ndmin=2)
	ref_comp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_comp"]),    dtype="float", ndmin=2)
	ref_covs   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_cov"]),     dtype="float", ndmin=1)
	ref_lens   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_lengths"]), dtype="int", ndmin=1)
	
	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("contig\tlength\t%s\n" % motif_header)
	f_motifs.write("contig\tlength\t%s\n" % motif_header)
	f_others.write("contig\tlength\tcoverage\t%s\n" % kmer_header)

	for i,name in enumerate(ref_names):
		methyl_str = "\t".join(ref_SCp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, ref_lens[i], methyl_str))

		motifs_str = "\t".join(ref_SCp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, ref_lens[i], motifs_str))

		comp_str   = "\t".join(ref_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%.2f\t%s\n" % (name, ref_lens[i], ref_covs[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def write_aligned_read_features( mbinRunner, opts ):
	"""
	If --aligned_read_barcodes is used, in addition to the contig-
	level features output, copy the data from the tmp directory into 
	alignment-level output files.
	"""
	methyl_fn = "%salign_methyl_features.txt" % opts.prefix
	motifs_fn = "%salign_motif_counts.txt"    % opts.prefix
	others_fn = "%salign_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for aligned reads" % methyl_fn)
	logging.info("   %s  \t -- Motif counts on aligned reads" % motifs_fn)
	logging.info("   %s\t -- Length values, as well as %s-mer frequency vectors for aligned reads" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	align_names = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_names"]),      dtype="str",   ndmin=1)
	align_SMp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp"]),        dtype="float", ndmin=2)
	align_SMp_N = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_counts"]), dtype="int",   ndmin=2)
	align_comp  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp"]),       dtype="float", ndmin=2)
	align_lens  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_lengths"]),    dtype="int",   ndmin=1)

	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("read\tlength\t%s\n" % motif_header)
	f_motifs.write("read\tlength\t%s\n" % motif_header)
	f_others.write("read\tlength\t%s\n" % kmer_header)

	for i,name in enumerate(align_names):
		methyl_str = "\t".join(align_SMp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, align_lens[i], methyl_str))

		motifs_str = "\t".join(align_SMp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, align_lens[i], motifs_str))

		comp_str   = "\t".join(align_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%s\n" % (name, align_lens[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def write_unaligned_read_features( mbinRunner, opts ):
	"""
	Gather unaligned read-level data generated by the pipeline and 
	write to output files.
	"""
	methyl_fn = "%sread_methyl_features.txt" % opts.prefix
	motifs_fn = "%sread_motif_counts.txt"    % opts.prefix
	others_fn = "%sread_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for unaligned reads" % methyl_fn)
	logging.info("   %s  \t -- Motif counts for unaligned reads" % motifs_fn)
	logging.info("   %s\t -- Length values, as well as %s-mer frequency vectors for unaligned reads" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	read_names = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_names"]),      dtype="str",   ndmin=1)
	read_SMp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp"]),        dtype="float", ndmin=2)
	read_SMp_N = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_counts"]), dtype="int",   ndmin=2)
	read_comp  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp"]),       dtype="float", ndmin=2)
	read_lens  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_lengths"]),    dtype="int",   ndmin=1)

	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("read\tlength\t%s\n" % motif_header)
	f_motifs.write("read\tlength\t%s\n" % motif_header)
	f_others.write("read\tlength\t%s\n" % kmer_header)

	for i,name in enumerate(read_names):
		methyl_str = "\t".join(read_SMp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, read_lens[i], methyl_str))

		motifs_str = "\t".join(read_SMp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, read_lens[i], motifs_str))

		comp_str   = "\t".join(read_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%s\n" % (name, read_lens[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def find_motifs_in_control(opts, motifs):
	"""
	Make sure all specified motifs have control values in the control
	dictionary (specified by --control_pkl_name). Remove any motifs 
	that do not have control values.
	"""
	control_d    = pickle.load( open(opts.control_pkl_name,"rb" ) )
	found_motifs = []
	not_found    = []
	for m in motifs:
		if control_d.get(m):
			found_motifs.append(m)
		else:
			not_found.append(m)

	return np.array(found_motifs), not_found

def combine_subreads_for_read_level( tup ):
	input_file   = tup[0]
	contig       = tup[1]
	tmp          = tup[2]
	aln_fn_labels    = tup[3]
	j            = tup[4]
	N_contigs    = tup[5]
	logging.info("...contig %s (%s/%s)" % (contig, j, N_contigs))
	contig_fns = glob.glob( os.path.join(tmp, "%s_*.tmp" % contig) )
	for fn in contig_fns:
		if   fn.find("_labels.tmp")>-1:
			labels_fn     = fn
		elif fn.find("_lengths.tmp")>-1:
			lengths_fn    = fn
		elif fn.find("_readnames.tmp")>-1:
			readnames_fn  = fn
		elif fn.find("_ipds.tmp")>-1:
			ipds_fn       = fn
		elif fn.find("_ipdsN.tmp")>-1:
			ipds_N_fn     = fn
		elif fn.find("_compN.tmp")>-1:
			comp_N_fn     = fn
		elif fn.find("_compkmers.tmp")>-1:
			comp_kmers_fn = fn
		elif fn.find("_ipdskmers.tmp")>-1:
			ipds_kmers_fn = fn
		elif fn.find("_strand.tmp")>-1:
			strand_fn     = fn

	readname_row_idx = defaultdict(list)
	for i,line in enumerate(open(readnames_fn, "r").xreadlines()):
		readname = line.strip()
		readname_row_idx[readname].append(i)
	n_subreads = i+1

	for line in open(ipds_kmers_fn).xreadlines():
		n_motifs = len(line.strip().split("\t"))

	sublengths = np.loadtxt(lengths_fn,    dtype="int")
	ipds       = np.loadtxt(ipds_fn,       dtype="float")
	ipds_N     = np.loadtxt(ipds_N_fn,     dtype="int")
	ipds_kmers = np.loadtxt(ipds_kmers_fn, dtype="string")
	comp_N     = np.loadtxt(comp_N_fn,     dtype="int")
	comp_kmers = np.loadtxt(comp_kmers_fn, dtype="string")
	strands    = np.loadtxt(strand_fn,     dtype="int")

	if n_subreads>1:
		if n_motifs==1:
			# Only one motif survived motif-filtering
			# Case #3
			ipds       = ipds.reshape(ipds.shape[0],1)
			ipds_N     = ipds_N.reshape(ipds_N.shape[0],1)
			ipds_kmers = ipds_kmers.reshape(1)
		elif n_motifs>1:
			# The ipd data is loaded properly in matrix form
			# Case #1
			pass
		elif n_motifs==0:
			pass
	elif n_subreads==1:
		sublengths = sublengths.reshape(1)
		strands    = strands.reshape(1)
		comp_N     = comp_N.reshape(1,comp_N.shape[0])
		if n_motifs==1:
			# Case #4
			ipds       = ipds.reshape(1,1)
			ipds_N     = ipds_N.reshape(1,1)
			ipds_kmers = ipds_kmers.reshape(1)
		elif n_motifs>1:
			# Case #2
			ipds       = ipds.reshape(1,ipds.shape[0])
			ipds_N     = ipds_N.reshape(1,ipds_N.shape[0])
		elif n_motifs==0:
			pass
		
	reads_names_fn      = os.path.join(tmp, "%s_read_names.tmp"     % contig)
	reads_labels_fn     = os.path.join(tmp, "%s_read_labels.tmp"    % contig)
	reads_lengths_fn    = os.path.join(tmp, "%s_read_lengths.tmp"   % contig)
	reads_contig_fn     = os.path.join(tmp, "%s_read_contig.tmp"    % contig)
	reads_ipds_fn       = os.path.join(tmp, "%s_read_ipds.tmp"      % contig)
	reads_ipds_N_fn     = os.path.join(tmp, "%s_read_ipdsN.tmp"     % contig)
	reads_ipds_kmers_fn = os.path.join(tmp, "%s_read_ipdskmers.tmp" % contig)
	reads_comp_N_fn     = os.path.join(tmp, "%s_read_compN.tmp"     % contig)
	reads_comp_kmers_fn = os.path.join(tmp, "%s_read_compkmers.tmp" % contig)
	reads_strands_fn    = os.path.join(tmp, "%s_read_strands.tmp"   % contig)

	f_reads      = open(reads_names_fn,   "w")
	f_labels     = open(reads_labels_fn,  "w")
	f_lengths    = open(reads_lengths_fn, "w")
	f_contig     = open(reads_contig_fn,  "w")
	f_ipds       = open(reads_ipds_fn,    "w")
	f_ipds_N     = open(reads_ipds_N_fn,  "w")
	f_comp_N     = open(reads_comp_N_fn,  "w")
	f_strands    = open(reads_strands_fn, "w")
	label        = aln_fn_labels[input_file]
	for readname,row_idx in readname_row_idx.iteritems():
		f_reads.write(  "%s\n" % readname)
		f_labels.write( "%s\n" % label)
		f_contig.write( "%s\n" % contig)
		comp_N_list  = comp_N[row_idx,:].sum(axis=0)
		strands_list = strands[row_idx]
		# ipds_list   = ipds[row_idx,:].mean(axis=0)
		ipds_N_list  = ipds_N[row_idx,:].sum(axis=0)
		ipds_list    = []
		for k in range(n_motifs):
			read_motifs_N = np.sum(ipds_N[row_idx,k])
			read_ipds_sum = np.sum(ipds[row_idx,k] * ipds_N[row_idx,k])
			if read_motifs_N > 0:
				motif_mean = read_ipds_sum / read_motifs_N
			else:
				motif_mean = 0.0
			ipds_list.append(motif_mean)

		# Normalize composition kmer counts
		normed_comp_N_list = map(lambda x: math.log( float(x)/sum(comp_N_list) ), comp_N_list)
		readlength         = sublengths[row_idx].sum()
		f_lengths.write( "%s\n" % readlength)
		f_ipds.write(    "%s\n" % "\t".join(map(lambda x: str(round(x,4)), ipds_list)))
		f_ipds_N.write(  "%s\n" % "\t".join(map(lambda x: str(x),          ipds_N_list)))
		f_comp_N.write(  "%s\n" % "\t".join(map(lambda x: str(round(x,4)), normed_comp_N_list)))
		f_strands.write( "%s\n" %  ",".join(map(lambda x: str(x),          strands_list)))
	
	shutil.copy(ipds_kmers_fn, os.path.join(tmp, "%s_read_ipdskmers.tmp" % contig))
	shutil.copy(comp_kmers_fn, os.path.join(tmp, "%s_read_compkmers.tmp" % contig))
	f_reads.close()
	f_labels.close()
	f_lengths.close()
	f_contig.close()
	f_ipds.close()
	f_ipds_N.close()
	f_comp_N.close()
	f_strands.close()

	# Remove the subread-level barcodes for each contig
	for fn in contig_fns:
		os.remove(fn)

def build_profiles(opts, input_files, motifs, motifs_fn):
	"""
	Build profiles of methylation scores for the desired motifs
	across the sequences, either contigs or unaligned reads.
	"""
	if os.path.exists(opts.tmp):
		shutil.rmtree(opts.tmp)
	os.mkdir(opts.tmp)

	opts.motifs_file = motifs_fn
	opts.motifs      = motifs
	opts.bi_motifs   = None

	logging.info("Building methylation profiles using %s motifs..." % len(opts.motifs))
	to_del = glob.glob( os.path.join(opts.tmp, "*") )
	for fn in to_del:
		os.remove(fn)

	mbinRunner     = mbin.mbinRunner(opts)
	##################################################
	# Launch analysis of <N_reads> for motif filtering
	##################################################
	for i,input_file in enumerate(input_files):
		logging.info("Creating %s barcodes (%s motifs) from %s..." % (opts.N_reads, len(opts.motifs), input_file))
		mbinRunner.launch_data_loader( input_file, opts.N_reads, i, opts )

		if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
			logging.info("Combining aligned subread-level barcodes to get read-level barcodes from each contig...")
			contig_labels_fns = glob.glob( os.path.join(opts.tmp, "*_labels.tmp") )
			contigs           = map(lambda x: os.path.basename(x).split("_labels.tmp")[0], contig_labels_fns)
			args              = [ (input_file, contig, opts.tmp, opts.aln_fn_labels, i, len(contigs)) for i,contig in enumerate(contigs)]
			results           = mbin.launch_pool( opts.procs, combine_subreads_for_read_level, args )

			logging.info("Combining read-level barcodes from all contigs...")
			mbinRunner.combine_read_level_barcodes_across_contigs()
			logging.info("Done.")

			logging.info("Creating contig-level barcodes (%s motifs) from %s..." % (len(opts.motifs), input_file))
			mbinRunner.combine_subreads_for_contig_level( input_file )
			logging.info("Done.")
			n_contigs = len(np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]), dtype="str", ndmin=1))

			if opts.cross_cov_bins!=None:
				logging.info("Creating bin-level barcodes (%s motifs) using %s..." % (len(opts.motifs), opts.cross_cov_bins))
				mbinRunner.combine_contigs_for_bin_level()
				logging.info("Done.")

	if opts.input_ftype=="bas":
		# Combine subread data across multiple movies
		logging.info("Combining subread data across all movies...")
		results = mbinRunner.combine_subread_data_across_bas_movies()
		logging.info("Done.")
		# Combine movie-merged subreads to get read-level barcodes
		logging.info("Combining subreads to get read-level barcodes...")
		results = mbinRunner.bas_combine_subreads_for_read_level()
		logging.info("Done.")

		if opts.sam!=None:
			logging.info("Writing read-contig assignments based on %s..." % opts.sam)
			mbinRunner.get_read_refs_from_SAM()
			logging.info("Done.")
			for i,input_file in enumerate(input_files):
				logging.info("Creating contig-level barcodes (%s motifs) from %s..." % (len(opts.motifs), input_file))
				mbinRunner.combine_subreads_for_contig_level( input_file )
				logging.info("Done.")
			n_contigs = len(np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]), dtype="str", ndmin=1))

	logging.info("Writing output files:")

	if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
		write_contig_features(mbinRunner, opts)
		if opts.aligned_read_barcodes:
			write_aligned_read_features( mbinRunner, opts )
	elif opts.input_ftype=="bas":
		write_unaligned_read_features( mbinRunner, opts )

	logging.info("Cleaning up temp files from methylation profiling...")
	shutil.rmtree(opts.tmp)
	logging.info("Pipeline finished.")

def __parseArgs():
	"""Handle command line argument parsing"""

	usage = """%prog [--help] [options] [input_seqs] [motifs_file]

	methylprofiles compiles the methylation profiles (for contigs or reads)
	across the motifs specified in the arguments. The methylation profiles 
	can be constructed using either native reads aligned to assembled contigs 
	(*.bam or *.cmp.h5) or unaligned reads (*.bas.h5), the latter of which can 
	be supplied as a single bas.h5 file or a FOFN containing multiple *.bas.h5 
	files. Unaligned BAM files are not currently supported.

	The motifs specified in the motifs file must be line-separated and
	indicate the methylated position in the motif using a 0-based index. For 
	instance, GATC-1 indicates that the adenine position is methylated, while 
	CATAG-3	says that the second adenine in the motif is methylated. This file
	would look as follows:

	GATC-1
	CATAG-3

	Two files will be output from methylprofiles (and serve as input to the 
	routine mapfeatures): 
	   (1) <prefix>_<seq>_methyl_features.txt
	   (2) <prefix>_<seq>_other_features.txt

	<prefix> is defined by --prefix and <seq> is the sequence data type: contig, 
	align, or read.


	Usage examples:

	### Using BAM of aligned reads as input and motifs specified in motifs.txt ###

	methylprofiles -i --contigs=reference.fasta aligned_reads.bam motifs.txt
	

	### Using *.cmp.h5 of aligned reads as input and motifs specified in motifs.txt ###

	methylprofiles -i --contigs=reference.fasta aligned_reads.cmp.h5 motifs.txt


	### Using unaligned reads in a bas.h5 file ###

	methylprofiles -i m12345.bas.h5 motifs.txt


	### Using a FOFN of multiple bas.h5 files ###

	methylprofiles -i bas.h5.fofn motifs.txt
	"""

	parser = optparse.OptionParser( usage=usage, description=__doc__ )

	parser.add_option( "-d", "--debug", action="store_true", help="Increase verbosity of logging" )

	parser.add_option( "-i", "--info", action="store_true", help="Add basic logging" )

	parser.add_option( "--logFile", type="str", help="Write logging to file [log.controls]" )
	
	parser.add_option( "--prefix", type="str", help="Prefix to use for output files [None]" )

	parser.add_option( "--tmp", type="str", help="Directory where numerous temporary files will be written [profiles_tmp]" )
	
	parser.add_option( "--contigs", type="str", help="Fasta file containing entries for the assembled contigs [None]" )

	parser.add_option( "--minReadScore", type="float", help="Min read score of an unaligned read [0.0]" )

	parser.add_option( "--maxPausiness", type="float", help="Max pausiness value of an unaligned read [1000]" )

	parser.add_option( "--minQV", type="float", help="If base has QV < minQV, do not include [0]" )
	
	parser.add_option( "--subreadlength_min", type="int", help="Minimum subread length to include for analysis [100]" )

	parser.add_option( "--readlength_min", type="int", help="Minimum read length to include for analysis [100]" )

	parser.add_option( "--minContigLength", type="int", help="Min length of contig to consider [10000]" )
	
	parser.add_option( "--comp_kmer", type="int", help="Kmer size to use for sequence composition measurements [5]" )
	
	parser.add_option( "--aligned_read_barcodes", action="store_true", help="Also output features for individual aligned reads, not just contigs (requires cmp.h5 input) [False]" )

	parser.add_option( "--minAcc", type="float", help="Min subread accuracy of read [0.8]" )

	parser.add_option( "--minMapQV", type="float", help="Min mapping QV of aligned read [240]" )

	parser.add_option( "--procs", type="int", help="Number of processors to use [4]" )

	parser.add_option( "--N_reads", type="int", help="Number of qualifying reads to include (from each bas.h5 if input is FOFN of bas.h5 files) in analysis [1000000000]" )
	
	parser.add_option( "--control_pkl_name", type="str", help="Filename to save control IPD data from WGA sequencing [control_ipds.pkl]" )
	
	parser.add_option( "--no_subtract_control", action="store_true", help="Subtract control IPDs in final calculations [True]" )

	parser.add_option( "--cross_cov_bins", type="str", help="Path to file containing binning results from CONCOCT. Will use to improve motif discovery. Only works with contig-level analysis (cmp.h5 input) inputs. File format should be '<contig_name>,<bin_id>' [None]" )

	parser.set_defaults( logFile="log.methylprofiles",         \
						 info=False,                           \
						 debug=False,                          \
						 prefix="",                            \
						 tmp="profiles_tmp",                   \
						 contigs=None,                         \
						 minQV=0,                              \
						 subreadlength_min=100,                \
						 readlength_min=100,                   \
						 minContigLength=10000,                \
						 minReadScore=0.0,                     \
						 maxPausiness=1000,                    \
						 comp_kmer=5,                          \
						 aligned_read_barcodes=False,          \
						 minAcc=0.8,                           \
						 minMapQV=240,                         \
						 procs=4,                              \
						 N_reads=1000000000,                   \
						 control_pkl_name="control_ipds.pkl",  \
						 no_subtract_control=False,            \
						 cross_cov_bins=None,                  \
						 )

	opts, args = parser.parse_args( )

	input_files, motifs_fn  = __check_input( opts, args, parser )

	opts.comp_only     = False
	opts.sam           = None
	opts.skip_motifs   = None
	opts.bas_whitelist = None

	if opts.no_subtract_control:
		opts.subtract_control = False
	else:
		opts.subtract_control = True

	if opts.prefix != "":
		opts.prefix+="_"
	
	opts.control_pkl_name = os.path.abspath(opts.control_pkl_name)

	if not os.path.exists(opts.control_pkl_name):
		parser.error("Can't find control IPDs pickle file: %s" % opts.control_pkl_name)

	return opts,input_files,motifs_fn

def __initLog( opts ):
	"""Sets up logging based on command line arguments. Allows for three levels of logging:
	logging.error( ): always emitted
	logging.info( ) : emitted with --info or --debug
	logging.debug( ): only with --debug"""

	if os.path.exists(opts.logFile):
		os.remove(opts.logFile)

	logLevel = logging.DEBUG if opts.debug \
				else logging.INFO if opts.info \
				else logging.ERROR

	logger = logging.getLogger("")
	logger.setLevel(logLevel)
	
	# create file handler which logs even debug messages
	fh = logging.FileHandler(opts.logFile)
	fh.setLevel(logLevel)
	
	# create console handler with a higher log level
	ch = logging.StreamHandler()
	ch.setLevel(logLevel)
	
	# create formatter and add it to the handlers
	logFormat = "%(asctime)s [%(levelname)s] %(message)s"
	formatter = logging.Formatter(logFormat, "%Y-%m-%d %H:%M:%S")
	ch.setFormatter(formatter)
	fh.setFormatter(formatter)
	
	# add the handlers to logger
	logger.addHandler(ch)
	logger.addHandler(fh)

def __check_input( opts, args, parser ):
	"""
	Verify that the input format looks OK and assess the input type.
	"""
	if len(args)!=2:
		print "ERROR -- expecting two arguments:"
		print "    (1) input file (bam, cmp.h5, bas.h5, or FOFN of bas.h5 files)"
		print "    (2) file containing the motifs to analyze, separated by newlines, e.g."
		print "          GATC-1"
		print "          CATG-1"
		print "          CATAG-3"
		print ""
		sys.exit()

	input_files, opts = mbin.check_input_ftype(opts, args, parser)

	if opts.input_ftype=="bas" and opts.cross_cov_bins!=None:
		parser.error("Use of the --cross_cov_bins option is not compatible with bas.h5 inputs!")

	if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
		try:
			for entry in SeqIO.parse(opts.contigs, "fasta"):
				x = entry.seq
				y = entry.id
		except:
			parser.error("Please make sure the --contigs input is a valid fasta file.")

	motifs_fn          = os.path.abspath(args[1])

	if not os.path.exists(motifs_fn):
		parser.error("Can't find file of motifs to include in methylation profile: %s" % motifs_fn)

	return input_files, motifs_fn

if __name__ == "__main__":
	main()
import os,sys
from collections import defaultdict
import numpy as np
import optparse
import logging
from pbcore.io.align.CmpH5IO import CmpH5Reader
from pbcore.io.BasH5IO import BasH5Reader
import glob
import numpy as np
import logging
import shutil
import pickle
import math
import mbin
import motif_tools
from Bio import SeqIO

def launch():
	opts, input_files, motifs_fn = __parseArgs()
	__initLog(opts)

	motifs           = np.loadtxt(motifs_fn, dtype="str", ndmin=1)
	motifs,not_found = find_motifs_in_control(opts, motifs)
	if len(not_found)>0:
		logging.warning("")
		logging.warning("  ******************** Important *********************")
		logging.warning("  Did not find %s motifs in %s:" % (len(not_found), opts.control_pkl_name))
		for nf in not_found:
			logging.warning("       %s" % nf)
		logging.warning("  These motif(s) will be removed from further analysis.")
		logging.warning("  These %s motifs will be kept:" % len(motifs))
		for m in motifs:
			logging.warning("       %s" % m)
		logging.warning("  ****************************************************")
		logging.warning("")
	else:
		logging.info("Found entries for all %s motifs in %s" % (len(motifs), opts.control_pkl_name))


	build_profiles(opts, input_files, motifs, motifs_fn)

	print >> sys.stderr, "mBin methylation profiling has finished running. See log for details."

def write_contig_features( mbinRunner, opts ):
	"""
	Gather contig-level data generated by the pipeline and write
	to output files.
	"""
	methyl_fn = "%scontig_methyl_features.txt" % opts.prefix
	motifs_fn = "%scontig_motif_counts.txt"    % opts.prefix
	others_fn = "%scontig_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for assembled contigs" % methyl_fn)
	logging.info("   %s  \t -- Motif counts on assembled contigs" % motifs_fn)
	logging.info("   %s\t -- Coverage and length values, as well as %s-mer frequency vectors for assembled contigs" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	ref_names  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]),   dtype="str", ndmin=1)
	ref_SCp    = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_SCp"]),     dtype="float", ndmin=2)
	ref_SCp_N  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_SCp_N"]),   dtype="int", ndmin=2)
	ref_comp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_comp"]),    dtype="float", ndmin=2)
	ref_covs   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_cov"]),     dtype="float", ndmin=1)
	ref_lens   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_lengths"]), dtype="int", ndmin=1)
	
	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("contig\tlength\t%s\n" % motif_header)
	f_motifs.write("contig\tlength\t%s\n" % motif_header)
	f_others.write("contig\tlength\tcoverage\t%s\n" % kmer_header)

	for i,name in enumerate(ref_names):
		methyl_str = "\t".join(ref_SCp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, ref_lens[i], methyl_str))

		motifs_str = "\t".join(ref_SCp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, ref_lens[i], motifs_str))

		comp_str   = "\t".join(ref_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%.2f\t%s\n" % (name, ref_lens[i], ref_covs[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def write_aligned_read_features( mbinRunner, opts ):
	"""
	If --aligned_read_barcodes is used, in addition to the contig-
	level features output, copy the data from the tmp directory into 
	alignment-level output files.
	"""
	methyl_fn = "%salign_methyl_features.txt" % opts.prefix
	motifs_fn = "%salign_motif_counts.txt"    % opts.prefix
	others_fn = "%salign_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for aligned reads" % methyl_fn)
	logging.info("   %s  \t -- Motif counts on aligned reads" % motifs_fn)
	logging.info("   %s\t -- Length values, as well as %s-mer frequency vectors for aligned reads" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	align_names = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_names"]),      dtype="str",   ndmin=1)
	align_SMp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp"]),        dtype="float", ndmin=2)
	align_SMp_N = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_counts"]), dtype="int",   ndmin=2)
	align_comp  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp"]),       dtype="float", ndmin=2)
	align_lens  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_lengths"]),    dtype="int",   ndmin=1)

	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("read\tlength\t%s\n" % motif_header)
	f_motifs.write("read\tlength\t%s\n" % motif_header)
	f_others.write("read\tlength\t%s\n" % kmer_header)

	for i,name in enumerate(align_names):
		methyl_str = "\t".join(align_SMp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, align_lens[i], methyl_str))

		motifs_str = "\t".join(align_SMp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, align_lens[i], motifs_str))

		comp_str   = "\t".join(align_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%s\n" % (name, align_lens[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def write_unaligned_read_features( mbinRunner, opts ):
	"""
	Gather unaligned read-level data generated by the pipeline and 
	write to output files.
	"""
	methyl_fn = "%sread_methyl_features.txt" % opts.prefix
	motifs_fn = "%sread_motif_counts.txt"    % opts.prefix
	others_fn = "%sread_other_features.txt"  % opts.prefix

	logging.info("   %s\t -- Methylation profiles for unaligned reads" % methyl_fn)
	logging.info("   %s  \t -- Motif counts for unaligned reads" % motifs_fn)
	logging.info("   %s\t -- Length values, as well as %s-mer frequency vectors for unaligned reads" % (others_fn, opts.comp_kmer))

	f_methyl  = open(methyl_fn, "wb")
	f_motifs  = open(motifs_fn, "wb")
	f_others  = open(others_fn, "wb")

	read_names = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_names"]),      dtype="str",   ndmin=1)
	read_SMp   = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp"]),        dtype="float", ndmin=2)
	read_SMp_N = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_counts"]), dtype="int",   ndmin=2)
	read_comp  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp"]),       dtype="float", ndmin=2)
	read_lens  = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_lengths"]),    dtype="int",   ndmin=1)

	# Write header
	motifs       = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_SMp_kmers"]),  dtype="str", ndmin=1)
	kmers        = np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["read_comp_kmers"]), dtype="str", ndmin=1)
	motif_header = "\t".join(motifs)
	kmer_header  = "\t".join(kmers)
	f_methyl.write("read\tlength\t%s\n" % motif_header)
	f_motifs.write("read\tlength\t%s\n" % motif_header)
	f_others.write("read\tlength\t%s\n" % kmer_header)

	for i,name in enumerate(read_names):
		methyl_str = "\t".join(read_SMp[i,:].astype("str"))
		f_methyl.write("%s\t%s\t%s\n" % (name, read_lens[i], methyl_str))

		motifs_str = "\t".join(read_SMp_N[i,:].astype("str"))
		f_motifs.write("%s\t%s\t%s\n" % (name, read_lens[i], motifs_str))

		comp_str   = "\t".join(read_comp[i,:].astype("str"))
		f_others.write("%s\t%s\t%s\n" % (name, read_lens[i], comp_str))

	f_methyl.close()
	f_motifs.close()
	f_others.close()

def find_motifs_in_control(opts, motifs):
	"""
	Make sure all specified motifs have control values in the control
	dictionary (specified by --control_pkl_name). Remove any motifs 
	that do not have control values.
	"""
	control_d    = pickle.load( open(opts.control_pkl_name,"rb" ) )
	found_motifs = []
	not_found    = []
	for m in motifs:
		if control_d.get(m):
			found_motifs.append(m)
		else:
			not_found.append(m)

	return np.array(found_motifs), not_found

def combine_subreads_for_read_level( tup ):
	input_file   = tup[0]
	contig       = tup[1]
	tmp          = tup[2]
	aln_fn_labels    = tup[3]
	j            = tup[4]
	N_contigs    = tup[5]
	logging.info("...contig %s (%s/%s)" % (contig, j, N_contigs))
	contig_fns = glob.glob( os.path.join(tmp, "%s_*.tmp" % contig) )
	for fn in contig_fns:
		if   fn.find("_labels.tmp")>-1:
			labels_fn     = fn
		elif fn.find("_lengths.tmp")>-1:
			lengths_fn    = fn
		elif fn.find("_readnames.tmp")>-1:
			readnames_fn  = fn
		elif fn.find("_ipds.tmp")>-1:
			ipds_fn       = fn
		elif fn.find("_ipdsN.tmp")>-1:
			ipds_N_fn     = fn
		elif fn.find("_compN.tmp")>-1:
			comp_N_fn     = fn
		elif fn.find("_compkmers.tmp")>-1:
			comp_kmers_fn = fn
		elif fn.find("_ipdskmers.tmp")>-1:
			ipds_kmers_fn = fn
		elif fn.find("_strand.tmp")>-1:
			strand_fn     = fn

	readname_row_idx = defaultdict(list)
	for i,line in enumerate(open(readnames_fn, "r").xreadlines()):
		readname = line.strip()
		readname_row_idx[readname].append(i)
	n_subreads = i+1

	for line in open(ipds_kmers_fn).xreadlines():
		n_motifs = len(line.strip().split("\t"))

	sublengths = np.loadtxt(lengths_fn,    dtype="int")
	ipds       = np.loadtxt(ipds_fn,       dtype="float")
	ipds_N     = np.loadtxt(ipds_N_fn,     dtype="int")
	ipds_kmers = np.loadtxt(ipds_kmers_fn, dtype="string")
	comp_N     = np.loadtxt(comp_N_fn,     dtype="int")
	comp_kmers = np.loadtxt(comp_kmers_fn, dtype="string")
	strands    = np.loadtxt(strand_fn,     dtype="int")

	if n_subreads>1:
		if n_motifs==1:
			# Only one motif survived motif-filtering
			# Case #3
			ipds       = ipds.reshape(ipds.shape[0],1)
			ipds_N     = ipds_N.reshape(ipds_N.shape[0],1)
			ipds_kmers = ipds_kmers.reshape(1)
		elif n_motifs>1:
			# The ipd data is loaded properly in matrix form
			# Case #1
			pass
		elif n_motifs==0:
			pass
	elif n_subreads==1:
		sublengths = sublengths.reshape(1)
		strands    = strands.reshape(1)
		comp_N     = comp_N.reshape(1,comp_N.shape[0])
		if n_motifs==1:
			# Case #4
			ipds       = ipds.reshape(1,1)
			ipds_N     = ipds_N.reshape(1,1)
			ipds_kmers = ipds_kmers.reshape(1)
		elif n_motifs>1:
			# Case #2
			ipds       = ipds.reshape(1,ipds.shape[0])
			ipds_N     = ipds_N.reshape(1,ipds_N.shape[0])
		elif n_motifs==0:
			pass
		
	reads_names_fn      = os.path.join(tmp, "%s_read_names.tmp"     % contig)
	reads_labels_fn     = os.path.join(tmp, "%s_read_labels.tmp"    % contig)
	reads_lengths_fn    = os.path.join(tmp, "%s_read_lengths.tmp"   % contig)
	reads_contig_fn     = os.path.join(tmp, "%s_read_contig.tmp"    % contig)
	reads_ipds_fn       = os.path.join(tmp, "%s_read_ipds.tmp"      % contig)
	reads_ipds_N_fn     = os.path.join(tmp, "%s_read_ipdsN.tmp"     % contig)
	reads_ipds_kmers_fn = os.path.join(tmp, "%s_read_ipdskmers.tmp" % contig)
	reads_comp_N_fn     = os.path.join(tmp, "%s_read_compN.tmp"     % contig)
	reads_comp_kmers_fn = os.path.join(tmp, "%s_read_compkmers.tmp" % contig)
	reads_strands_fn    = os.path.join(tmp, "%s_read_strands.tmp"   % contig)

	f_reads      = open(reads_names_fn,   "w")
	f_labels     = open(reads_labels_fn,  "w")
	f_lengths    = open(reads_lengths_fn, "w")
	f_contig     = open(reads_contig_fn,  "w")
	f_ipds       = open(reads_ipds_fn,    "w")
	f_ipds_N     = open(reads_ipds_N_fn,  "w")
	f_comp_N     = open(reads_comp_N_fn,  "w")
	f_strands    = open(reads_strands_fn, "w")
	label        = aln_fn_labels[input_file]
	for readname,row_idx in readname_row_idx.iteritems():
		f_reads.write(  "%s\n" % readname)
		f_labels.write( "%s\n" % label)
		f_contig.write( "%s\n" % contig)
		comp_N_list  = comp_N[row_idx,:].sum(axis=0)
		strands_list = strands[row_idx]
		# ipds_list   = ipds[row_idx,:].mean(axis=0)
		ipds_N_list  = ipds_N[row_idx,:].sum(axis=0)
		ipds_list    = []
		for k in range(n_motifs):
			read_motifs_N = np.sum(ipds_N[row_idx,k])
			read_ipds_sum = np.sum(ipds[row_idx,k] * ipds_N[row_idx,k])
			if read_motifs_N > 0:
				motif_mean = read_ipds_sum / read_motifs_N
			else:
				motif_mean = 0.0
			ipds_list.append(motif_mean)

		# Normalize composition kmer counts
		normed_comp_N_list = map(lambda x: math.log( float(x)/sum(comp_N_list) ), comp_N_list)
		readlength         = sublengths[row_idx].sum()
		f_lengths.write( "%s\n" % readlength)
		f_ipds.write(    "%s\n" % "\t".join(map(lambda x: str(round(x,4)), ipds_list)))
		f_ipds_N.write(  "%s\n" % "\t".join(map(lambda x: str(x),          ipds_N_list)))
		f_comp_N.write(  "%s\n" % "\t".join(map(lambda x: str(round(x,4)), normed_comp_N_list)))
		f_strands.write( "%s\n" %  ",".join(map(lambda x: str(x),          strands_list)))
	
	shutil.copy(ipds_kmers_fn, os.path.join(tmp, "%s_read_ipdskmers.tmp" % contig))
	shutil.copy(comp_kmers_fn, os.path.join(tmp, "%s_read_compkmers.tmp" % contig))
	f_reads.close()
	f_labels.close()
	f_lengths.close()
	f_contig.close()
	f_ipds.close()
	f_ipds_N.close()
	f_comp_N.close()
	f_strands.close()

	# Remove the subread-level barcodes for each contig
	for fn in contig_fns:
		os.remove(fn)

def build_profiles(opts, input_files, motifs, motifs_fn):
	"""
	Build profiles of methylation scores for the desired motifs
	across the sequences, either contigs or unaligned reads.
	"""
	if os.path.exists(opts.tmp):
		shutil.rmtree(opts.tmp)
	os.mkdir(opts.tmp)

	opts.motifs_file = motifs_fn
	opts.motifs      = motifs
	opts.bi_motifs   = None

	logging.info("Building methylation profiles using %s motifs..." % len(opts.motifs))
	to_del = glob.glob( os.path.join(opts.tmp, "*") )
	for fn in to_del:
		os.remove(fn)

	mbinRunner     = mbin.mbinRunner(opts)
	##################################################
	# Launch analysis of <N_reads> for motif filtering
	##################################################
	for i,input_file in enumerate(input_files):
		logging.info("Creating %s barcodes (%s motifs) from %s..." % (opts.N_reads, len(opts.motifs), input_file))
		mbinRunner.launch_data_loader( input_file, opts.N_reads, i, opts )

		if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
			logging.info("Combining aligned subread-level barcodes to get read-level barcodes from each contig...")
			contig_labels_fns = glob.glob( os.path.join(opts.tmp, "*_labels.tmp") )
			contigs           = map(lambda x: os.path.basename(x).split("_labels.tmp")[0], contig_labels_fns)
			args              = [ (input_file, contig, opts.tmp, opts.aln_fn_labels, i, len(contigs)) for i,contig in enumerate(contigs)]
			results           = mbin.launch_pool( opts.procs, combine_subreads_for_read_level, args )

			logging.info("Combining read-level barcodes from all contigs...")
			mbinRunner.combine_read_level_barcodes_across_contigs()
			logging.info("Done.")

			logging.info("Creating contig-level barcodes (%s motifs) from %s..." % (len(opts.motifs), input_file))
			mbinRunner.combine_subreads_for_contig_level( input_file )
			logging.info("Done.")
			n_contigs = len(np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]), dtype="str", ndmin=1))

			if opts.cross_cov_bins!=None:
				logging.info("Creating bin-level barcodes (%s motifs) using %s..." % (len(opts.motifs), opts.cross_cov_bins))
				mbinRunner.combine_contigs_for_bin_level()
				logging.info("Done.")

	if opts.input_ftype=="bas":
		# Combine subread data across multiple movies
		logging.info("Combining subread data across all movies...")
		results = mbinRunner.combine_subread_data_across_bas_movies()
		logging.info("Done.")
		# Combine movie-merged subreads to get read-level barcodes
		logging.info("Combining subreads to get read-level barcodes...")
		results = mbinRunner.bas_combine_subreads_for_read_level()
		logging.info("Done.")

		if opts.sam!=None:
			logging.info("Writing read-contig assignments based on %s..." % opts.sam)
			mbinRunner.get_read_refs_from_SAM()
			logging.info("Done.")
			for i,input_file in enumerate(input_files):
				logging.info("Creating contig-level barcodes (%s motifs) from %s..." % (len(opts.motifs), input_file))
				mbinRunner.combine_subreads_for_contig_level( input_file )
				logging.info("Done.")
			n_contigs = len(np.loadtxt(os.path.join(opts.tmp, mbinRunner.fns["contig_names"]), dtype="str", ndmin=1))

	logging.info("Writing output files:")

	if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
		write_contig_features(mbinRunner, opts)
		if opts.aligned_read_barcodes:
			write_aligned_read_features( mbinRunner, opts )
	elif opts.input_ftype=="bas":
		write_unaligned_read_features( mbinRunner, opts )

	logging.info("Cleaning up temp files from methylation profiling...")
	shutil.rmtree(opts.tmp)
	logging.info("Pipeline finished.")

def __parseArgs():
	"""Handle command line argument parsing"""

	usage = """%prog [--help] [options] [input_seqs] [motifs_file]

	methylprofiles compiles the methylation profiles (for contigs or reads)
	across the motifs specified in the arguments. The methylation profiles 
	can be constructed using either native reads aligned to assembled contigs 
	(*.bam or *.cmp.h5) or unaligned reads (*.bas.h5), the latter of which can 
	be supplied as a single bas.h5 file or a FOFN containing multiple *.bas.h5 
	files. Unaligned BAM files are not currently supported.

	The motifs specified in the motifs file must be line-separated and
	indicate the methylated position in the motif using a 0-based index. For 
	instance, GATC-1 indicates that the adenine position is methylated, while 
	CATAG-3	says that the second adenine in the motif is methylated. This file
	would look as follows:

	GATC-1
	CATAG-3

	Two files will be output from methylprofiles (and serve as input to the 
	routine mapfeatures): 
	   (1) <prefix>_<seq>_methyl_features.txt
	   (2) <prefix>_<seq>_other_features.txt

	<prefix> is defined by --prefix and <seq> is the sequence data type: contig, 
	align, or read.


	Usage examples:

	### Using BAM of aligned reads as input and motifs specified in motifs.txt ###

	methylprofiles -i --contigs=reference.fasta aligned_reads.bam motifs.txt
	

	### Using *.cmp.h5 of aligned reads as input and motifs specified in motifs.txt ###

	methylprofiles -i --contigs=reference.fasta aligned_reads.cmp.h5 motifs.txt


	### Using unaligned reads in a bas.h5 file ###

	methylprofiles -i m12345.bas.h5 motifs.txt


	### Using a FOFN of multiple bas.h5 files ###

	methylprofiles -i bas.h5.fofn motifs.txt
	"""

	parser = optparse.OptionParser( usage=usage, description=__doc__ )

	parser.add_option( "-d", "--debug", action="store_true", help="Increase verbosity of logging" )

	parser.add_option( "-i", "--info", action="store_true", help="Add basic logging" )

	parser.add_option( "--logFile", type="str", help="Write logging to file [log.controls]" )
	
	parser.add_option( "--prefix", type="str", help="Prefix to use for output files [None]" )

	parser.add_option( "--tmp", type="str", help="Directory where numerous temporary files will be written [profiles_tmp]" )
	
	parser.add_option( "--contigs", type="str", help="Fasta file containing entries for the assembled contigs [None]" )

	parser.add_option( "--minReadScore", type="float", help="Min read score of an unaligned read [0.0]" )

	parser.add_option( "--maxPausiness", type="float", help="Max pausiness value of an unaligned read [1000]" )

	parser.add_option( "--minQV", type="float", help="If base has QV < minQV, do not include [0]" )
	
	parser.add_option( "--subreadlength_min", type="int", help="Minimum subread length to include for analysis [100]" )

	parser.add_option( "--readlength_min", type="int", help="Minimum read length to include for analysis [100]" )

	parser.add_option( "--minContigLength", type="int", help="Min length of contig to consider [10000]" )
	
	parser.add_option( "--comp_kmer", type="int", help="Kmer size to use for sequence composition measurements [5]" )
	
	parser.add_option( "--aligned_read_barcodes", action="store_true", help="Also output features for individual aligned reads, not just contigs (requires cmp.h5 input) [False]" )

	parser.add_option( "--minAcc", type="float", help="Min subread accuracy of read [0.8]" )

	parser.add_option( "--minMapQV", type="float", help="Min mapping QV of aligned read [240]" )

	parser.add_option( "--procs", type="int", help="Number of processors to use [4]" )

	parser.add_option( "--N_reads", type="int", help="Number of qualifying reads to include (from each bas.h5 if input is FOFN of bas.h5 files) in analysis [1000000000]" )
	
	parser.add_option( "--control_pkl_name", type="str", help="Filename to save control IPD data from WGA sequencing [control_ipds.pkl]" )
	
	parser.add_option( "--no_subtract_control", action="store_true", help="Subtract control IPDs in final calculations [True]" )

	parser.add_option( "--cross_cov_bins", type="str", help="Path to file containing binning results from CONCOCT. Will use to improve motif discovery. Only works with contig-level analysis (cmp.h5 input) inputs. File format should be '<contig_name>,<bin_id>' [None]" )

	parser.set_defaults( logFile="log.methylprofiles",         \
						 info=False,                           \
						 debug=False,                          \
						 prefix="",                            \
						 tmp="profiles_tmp",                   \
						 contigs=None,                         \
						 minQV=0,                              \
						 subreadlength_min=100,                \
						 readlength_min=100,                   \
						 minContigLength=10000,                \
						 minReadScore=0.0,                     \
						 maxPausiness=1000,                    \
						 comp_kmer=5,                          \
						 aligned_read_barcodes=False,          \
						 minAcc=0.8,                           \
						 minMapQV=240,                         \
						 procs=4,                              \
						 N_reads=1000000000,                   \
						 control_pkl_name="control_ipds.pkl",  \
						 no_subtract_control=False,            \
						 cross_cov_bins=None,                  \
						 )

	opts, args = parser.parse_args( )

	input_files, motifs_fn  = __check_input( opts, args, parser )

	opts.comp_only     = False
	opts.sam           = None
	opts.skip_motifs   = None
	opts.bas_whitelist = None

	if opts.no_subtract_control:
		opts.subtract_control = False
	else:
		opts.subtract_control = True

	if opts.prefix != "":
		opts.prefix+="_"
	
	opts.control_pkl_name = os.path.abspath(opts.control_pkl_name)

	if not os.path.exists(opts.control_pkl_name):
		parser.error("Can't find control IPDs pickle file: %s" % opts.control_pkl_name)

	return opts,input_files,motifs_fn

def __initLog( opts ):
	"""Sets up logging based on command line arguments. Allows for three levels of logging:
	logging.error( ): always emitted
	logging.info( ) : emitted with --info or --debug
	logging.debug( ): only with --debug"""

	if os.path.exists(opts.logFile):
		os.remove(opts.logFile)

	logLevel = logging.DEBUG if opts.debug \
				else logging.INFO if opts.info \
				else logging.ERROR

	logger = logging.getLogger("")
	logger.setLevel(logLevel)
	
	# create file handler which logs even debug messages
	fh = logging.FileHandler(opts.logFile)
	fh.setLevel(logLevel)
	
	# create console handler with a higher log level
	ch = logging.StreamHandler()
	ch.setLevel(logLevel)
	
	# create formatter and add it to the handlers
	logFormat = "%(asctime)s [%(levelname)s] %(message)s"
	formatter = logging.Formatter(logFormat, "%Y-%m-%d %H:%M:%S")
	ch.setFormatter(formatter)
	fh.setFormatter(formatter)
	
	# add the handlers to logger
	logger.addHandler(ch)
	logger.addHandler(fh)

def __check_input( opts, args, parser ):
	"""
	Verify that the input format looks OK and assess the input type.
	"""
	if len(args)!=2:
		print "ERROR -- expecting two arguments:"
		print "    (1) input file (bam, cmp.h5, bas.h5, or FOFN of bas.h5 files)"
		print "    (2) file containing the motifs to analyze, separated by newlines, e.g."
		print "          GATC-1"
		print "          CATG-1"
		print "          CATAG-3"
		print ""
		sys.exit()

	input_files, opts = mbin.check_input_ftype(opts, args, parser)

	if opts.input_ftype=="bas" and opts.cross_cov_bins!=None:
		parser.error("Use of the --cross_cov_bins option is not compatible with bas.h5 inputs!")

	if opts.input_ftype=="cmp" or opts.input_ftype=="bam":
		try:
			for entry in SeqIO.parse(opts.contigs, "fasta"):
				x = entry.seq
				y = entry.id
		except:
			parser.error("Please make sure the --contigs input is a valid fasta file.")

	motifs_fn          = os.path.abspath(args[1])

	if not os.path.exists(motifs_fn):
		parser.error("Can't find file of motifs to include in methylation profile: %s" % motifs_fn)

	return input_files, motifs_fn

if __name__ == "__main__":
	main()
Gather contig-level data generated by the pipeline and write 	to output files. # Write header If --aligned_read_barcodes is used, in addition to the contig- 	level features output, copy the data from the tmp directory into  	alignment-level output files. # Write header Gather unaligned read-level data generated by the pipeline and  	write to output files. # Write header Make sure all specified motifs have control values in the control 	dictionary (specified by --control_pkl_name). Remove any motifs  	that do not have control values. # Only one motif survived motif-filtering # Case #3 # The ipd data is loaded properly in matrix form # Case #1 # Case #4 # Case #2 # ipds_list   = ipds[row_idx,:].mean(axis=0) # Normalize composition kmer counts # Remove the subread-level barcodes for each contig Build profiles of methylation scores for the desired motifs 	across the sequences, either contigs or unaligned reads. ################################################## # Launch analysis of <N_reads> for motif filtering ################################################## # Combine subread data across multiple movies # Combine movie-merged subreads to get read-level barcodes Handle command line argument parsing %prog [--help] [options] [input_seqs] [motifs_file]  	methylprofiles compiles the methylation profiles (for contigs or reads) 	across the motifs specified in the arguments. The methylation profiles  	can be constructed using either native reads aligned to assembled contigs  	(*.bam or *.cmp.h5) or unaligned reads (*.bas.h5), the latter of which can  	be supplied as a single bas.h5 file or a FOFN containing multiple *.bas.h5  	files. Unaligned BAM files are not currently supported.  	The motifs specified in the motifs file must be line-separated and 	indicate the methylated position in the motif using a 0-based index. For  	instance, GATC-1 indicates that the adenine position is methylated, while  	CATAG-3	says that the second adenine in the motif is methylated. This file 	would look as follows:  	GATC-1 	CATAG-3  	Two files will be output from methylprofiles (and serve as input to the  	routine mapfeatures):  	   (1) <prefix>_<seq>_methyl_features.txt 	   (2) <prefix>_<seq>_other_features.txt  	<prefix> is defined by --prefix and <seq> is the sequence data type: contig,  	align, or read.   	Usage examples:  	### Using BAM of aligned reads as input and motifs specified in motifs.txt ###  	methylprofiles -i --contigs=reference.fasta aligned_reads.bam motifs.txt 	  	### Using *.cmp.h5 of aligned reads as input and motifs specified in motifs.txt ###  	methylprofiles -i --contigs=reference.fasta aligned_reads.cmp.h5 motifs.txt   	### Using unaligned reads in a bas.h5 file ###  	methylprofiles -i m12345.bas.h5 motifs.txt   	### Using a FOFN of multiple bas.h5 files ###  	methylprofiles -i bas.h5.fofn motifs.txt Sets up logging based on command line arguments. Allows for three levels of logging: 	logging.error( ): always emitted 	logging.info( ) : emitted with --info or --debug 	logging.debug( ): only with --debug # create file handler which logs even debug messages # create console handler with a higher log level # create formatter and add it to the handlers # add the handlers to logger Verify that the input format looks OK and assess the input type.