# Path: seqcluster/libs/sam2bed.py
# def makeBED(samFields):
#     samFlag = int(samFields.flag)
#     # Only create a BED entry if the read was aligned
#     if (not (samFlag & 0x0004)) and samFields.pos:
#             name = samFields.qname
#             seq = name.split("-")[1]
#             chrom = samFields.rname
#             start = str(int(samFields.pos))
#             end = str(int(samFields.pos) + len(seq) - 1)
#             strand = getStrand(samFlag)
#             # Write out the BED entry
#             return bedaligned("%s\t%s\t%s\t%s\t.\t%s\n" %  (chrom, start, end, name, strand))
#     else:
#             return False
# 
# Path: seqcluster/libs/do.py
# def find_cmd(cmd):
#     try:
#         return subprocess.check_output(["which", cmd]).strip()
#     except subprocess.CalledProcessError:
#         return None
# 
# def run(cmd, data=None, checks=None, region=None, log_error=True,
#         log_stdout=False):
#     """Run the provided command, logging details and checking for errors.
#     """
#     try:
#         logger.debug(" ".join(str(x) for x in cmd) if not isinstance(cmd, basestring) else cmd)
#         _do_run(cmd, checks, log_stdout)
#     except:
#         if log_error:
#             logger.info("error at command")
#         raise


from collections import defaultdict
from seqcluster.libs.sam2bed import makeBED
from seqcluster.libs.do import find_cmd, run
from Bio import pairwise2
from Bio.Seq import Seq
import json, itertools
import tempfile, os, contextlib, shutil
import operator
import logging
import pysam
import pybedtools
def get_loci_fasta(loci, out_fa, ref):
    """get fasta from precursor"""
    if not find_cmd("bedtools"):
        raise ValueError("Not bedtools installed")
    with make_temp_directory() as temp:
        bed_file = os.path.join(temp, "file.bed")
        for nc, loci in loci.iteritems():
            for l in loci:
                with open(bed_file, 'w') as bed_handle:
                    logger.debug("get_fasta: loci %s" % l)
                    nc, c, s, e, st = l
                    print("{0}\t{1}\t{2}\t{3}\t{3}\t{4}".format(c, s, e, nc, st), file=bed_handle)
                get_fasta(bed_file, ref, out_fa)
    return out_fa

def read_alignment(out_sam, loci, seqs, out_file):
    """read which seqs map to which loci and
    return a tab separated file"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():

# Path: seqcluster/libs/sam2bed.py
# def makeBED(samFields):
#     samFlag = int(samFields.flag)
#     # Only create a BED entry if the read was aligned
#     if (not (samFlag & 0x0004)) and samFields.pos:
#             name = samFields.qname
#             seq = name.split("-")[1]
#             chrom = samFields.rname
#             start = str(int(samFields.pos))
#             end = str(int(samFields.pos) + len(seq) - 1)
#             strand = getStrand(samFlag)
#             # Write out the BED entry
#             return bedaligned("%s\t%s\t%s\t%s\t.\t%s\n" %  (chrom, start, end, name, strand))
#     else:
#             return False
# 
# Path: seqcluster/libs/do.py
# def find_cmd(cmd):
#     try:
#         return subprocess.check_output(["which", cmd]).strip()
#     except subprocess.CalledProcessError:
#         return None
# 
# def run(cmd, data=None, checks=None, region=None, log_error=True,
#         log_stdout=False):
#     """Run the provided command, logging details and checking for errors.
#     """
#     try:
#         logger.debug(" ".join(str(x) for x in cmd) if not isinstance(cmd, basestring) else cmd)
#         _do_run(cmd, checks, log_stdout)
#     except:
#         if log_error:
#             logger.info("error at command")
#         raise


# Path: seqcluster/libs/read.py
from collections import defaultdict
from seqcluster.libs.sam2bed import makeBED
from seqcluster.libs.do import find_cmd, run
from Bio import pairwise2
from Bio.Seq import Seq
import json, itertools
import tempfile, os, contextlib, shutil
import operator
import logging
import pysam
import pybedtools

def get_loci_fasta(loci, out_fa, ref):
    """get fasta from precursor"""
    if not find_cmd("bedtools"):
        raise ValueError("Not bedtools installed")
    with make_temp_directory() as temp:
        bed_file = os.path.join(temp, "file.bed")
        for nc, loci in loci.iteritems():
            for l in loci:
                with open(bed_file, 'w') as bed_handle:
                    logger.debug("get_fasta: loci %s" % l)
                    nc, c, s, e, st = l
                    print("{0}\t{1}\t{2}\t{3}\t{3}\t{4}".format(c, s, e, nc, st), file=bed_handle)
                get_fasta(bed_file, ref, out_fa)
    return out_fa

def read_alignment(out_sam, loci, seqs, out_file):
    """read which seqs map to which loci and
    return a tab separated file"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():

            nm = hit[0][0]