<gh_stars>1-10
from __future__ import print_function

from .degenerate_tools import fix_seq_object

import errno
import os
import subprocess
import io

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

SAMTOOLS_MAP = {"FLAG": 1,
                "RNAME": 2,
                "POS": 3,
                "SEQ": 9}

BOWTIE_SEED_LENGTH = 10
BOWTIE_SEED_MISMATCH = 1


def bowtie_make_index(dir_path, genome):
    """
    Makes a bowtie index in a specified directory from a FASTA file

    :param dir_path: str
        Path to target directory
    :param genome: str
        Path to sequence FASTA file

    :return out_db_path:
        Path to the bowtie index that can be passed to bowtie
    """

    if not os.path.isfile(genome):
        raise OSError(errno.ENOENT)

    if not os.path.isdir(dir_path):
        raise OSError(errno.ENOTDIR)

    _, gen_name = os.path.split(genome)

    out_db_path = os.path.join(dir_path, gen_name.split(".")[0])
    build_cmd = ["bowtie2-build", "-f", genome, out_db_path]
    proc = subprocess.run(build_cmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
    proc.check_returncode()

    return out_db_path


def bowtie_get_hits(query, db_path, cores=1):
    """
    Align a query sequence to a target reference genome database and parse the resulting SAM file.

    :param query: [SeqRecord]
        A list of SeqRecord objects which should be aligned to the reference genome
    :param db_path: str
        A path to the reference genome database. Passed to bowtie2.
    :param cores: int
        Number of CPU cores to run. Passed to bowtie2.

    :return: {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}
        Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the
        chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches
        to the query sequence
    """

    return _bowtie_sam_parse(_bowtie_align(query, db_path, cores=cores))


def bowtie_check_index(db_path):
    """
    Check to make sure that the bowtie index exists. Raises a FileNotFoundError if there is any problem.

    :param db_path: str
        A path to the reference genome database. Passed to bowtie-inspect.
    :return: bool
        Returns TRUE if no error
    """
    try:
        inspect_command = ["bowtie2-inspect", "-s", db_path]
        proc = subprocess.run(inspect_command, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL)
        proc.check_returncode()
    except subprocess.CalledProcessError as err:
        raise OSError(errno.ENOENT, "Unable to inspect bowtie index")

    return True


def _bowtie_align(query, db_path, sam_path=None, cores=1):
    """
    Runs bowtie2 with a list of sequencs as queries against a target reference genome database.

    :param query: [SeqRecord]
        A list of SeqRecord objects
    :param db_path: str
        A path to the reference genome database. Passed to bowtie2.
    :param sam_path: str
        A path to the output SAM file. If None (or not set), a tempfile will be created
    :param cores: int
        Number of CPU cores to run. Passed to bowtie2.

    :return: str
        Returns a path to the output SAM file (unchanged, unless the param sam_path was not set)
    """

    # Write the queries into a string as a FASTA file
    temp_str = io.StringIO()
    SeqIO.write(query, temp_str, format="fasta")

    # Pipe the FASTA string into bowtie 2 and pipe the output.
    bowtie_cmd = map(str, ["bowtie2", "-f", "/dev/stdin", "-x", db_path, "-a", "-L", BOWTIE_SEED_LENGTH,
                  "-N", BOWTIE_SEED_MISMATCH, "--mp", "1,1", "--threads", cores])
    proc = subprocess.run(bowtie_cmd, input=temp_str.getvalue(), stdout=subprocess.PIPE, stderr=subprocess.DEVNULL,
                          universal_newlines=True)
    proc.check_returncode()

    # Stash the piped output (SAM file) in a string
    sam = io.StringIO(proc.stdout.strip())

    return sam


def _bowtie_sam_parse(sam, remove=False):
    """
    Takes a path to a sam file produced by bowtie and parses it into a dict keyed by query of hits

    :param sam_path: str
        Path to SAM file generated by bowtie
    :param remove: bool
        Remove the SAM file after parsing

    :return {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}
        Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the
        chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches
        to the query sequence
    """
    sequences = {}

    for line in sam:
        try:
            query_seq, map_seq, rname, pos, strand, nm = _parse_sam_line(line)
        except ValueError:
            continue

        try:
            sequences[str(query_seq)].append((map_seq, rname, pos, strand, nm))
        except KeyError:
            sequences[str(query_seq)] = [(map_seq, rname, pos, strand, nm)]

    return sequences


def _parse_sam_line(line_str):
    if line_str[0] is "@":
        raise ValueError

    li_arr = line_str.strip().split()

    opt_flags = {}

    for val in li_arr:

        try:
            if val[2] is ":":
                _v_arr = val.split(":")
                opt_flags[_v_arr[0]] = _v_arr[-1]
        except IndexError:
            pass

    query_seq = Seq(li_arr[SAMTOOLS_MAP["SEQ"]])
    rname = li_arr[SAMTOOLS_MAP["RNAME"]]
    pos = int(li_arr[SAMTOOLS_MAP["POS"]])
    _flag = int(li_arr[SAMTOOLS_MAP["FLAG"]])

    try:
        nm = int(opt_flags["NM"])
    except KeyError:
        nm = 0

    try:
        map_seq = Seq(_rewrite_to_reference(query_seq, opt_flags["MD"]), alphabet=generic_dna)
    except KeyError:
        map_seq = query_seq

    if (_flag % 256) == 16:
        query_seq = query_seq.reverse_complement()
        map_seq = map_seq.reverse_complement()
        strand = "-"
    else:
        strand = "+"

    return query_seq, map_seq, rname, pos, strand, nm


def _rewrite_to_reference(seq_str, mm):
    seq_str = fix_seq_object(seq_str, make_type="str")

    _tmp_mm = ""
    pointy = 0
    srr = list(seq_str)

    for i in range(len(mm)):
        if mm[i].isalpha():
            if _tmp_mm is not "":
                pointy += int(_tmp_mm)
                _tmp_mm = ""
            srr[pointy] = mm[i]
            pointy += 1
        else:
            _tmp_mm += mm[i]

    return "".join(srr)

from __future__ import print_function

from .degenerate_tools import fix_seq_object

import errno
import os
import subprocess
import io

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

SAMTOOLS_MAP = {"FLAG": 1,
                "RNAME": 2,
                "POS": 3,
                "SEQ": 9}

BOWTIE_SEED_LENGTH = 10
BOWTIE_SEED_MISMATCH = 1


def bowtie_make_index(dir_path, genome):
    """
    Makes a bowtie index in a specified directory from a FASTA file

    :param dir_path: str
        Path to target directory
    :param genome: str
        Path to sequence FASTA file

    :return out_db_path:
        Path to the bowtie index that can be passed to bowtie
    """

    if not os.path.isfile(genome):
        raise OSError(errno.ENOENT)

    if not os.path.isdir(dir_path):
        raise OSError(errno.ENOTDIR)

    _, gen_name = os.path.split(genome)

    out_db_path = os.path.join(dir_path, gen_name.split(".")[0])
    build_cmd = ["bowtie2-build", "-f", genome, out_db_path]
    proc = subprocess.run(build_cmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
    proc.check_returncode()

    return out_db_path


def bowtie_get_hits(query, db_path, cores=1):
    """
    Align a query sequence to a target reference genome database and parse the resulting SAM file.

    :param query: [SeqRecord]
        A list of SeqRecord objects which should be aligned to the reference genome
    :param db_path: str
        A path to the reference genome database. Passed to bowtie2.
    :param cores: int
        Number of CPU cores to run. Passed to bowtie2.

    :return: {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}
        Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the
        chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches
        to the query sequence
    """

    return _bowtie_sam_parse(_bowtie_align(query, db_path, cores=cores))


def bowtie_check_index(db_path):
    """
    Check to make sure that the bowtie index exists. Raises a FileNotFoundError if there is any problem.

    :param db_path: str
        A path to the reference genome database. Passed to bowtie-inspect.
    :return: bool
        Returns TRUE if no error
    """
    try:
        inspect_command = ["bowtie2-inspect", "-s", db_path]
        proc = subprocess.run(inspect_command, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL)
        proc.check_returncode()
    except subprocess.CalledProcessError as err:
        raise OSError(errno.ENOENT, "Unable to inspect bowtie index")

    return True


def _bowtie_align(query, db_path, sam_path=None, cores=1):
    """
    Runs bowtie2 with a list of sequencs as queries against a target reference genome database.

    :param query: [SeqRecord]
        A list of SeqRecord objects
    :param db_path: str
        A path to the reference genome database. Passed to bowtie2.
    :param sam_path: str
        A path to the output SAM file. If None (or not set), a tempfile will be created
    :param cores: int
        Number of CPU cores to run. Passed to bowtie2.

    :return: str
        Returns a path to the output SAM file (unchanged, unless the param sam_path was not set)
    """

    # Write the queries into a string as a FASTA file
    temp_str = io.StringIO()
    SeqIO.write(query, temp_str, format="fasta")

    # Pipe the FASTA string into bowtie 2 and pipe the output.
    bowtie_cmd = map(str, ["bowtie2", "-f", "/dev/stdin", "-x", db_path, "-a", "-L", BOWTIE_SEED_LENGTH,
                  "-N", BOWTIE_SEED_MISMATCH, "--mp", "1,1", "--threads", cores])
    proc = subprocess.run(bowtie_cmd, input=temp_str.getvalue(), stdout=subprocess.PIPE, stderr=subprocess.DEVNULL,
                          universal_newlines=True)
    proc.check_returncode()

    # Stash the piped output (SAM file) in a string
    sam = io.StringIO(proc.stdout.strip())

    return sam


def _bowtie_sam_parse(sam, remove=False):
    """
    Takes a path to a sam file produced by bowtie and parses it into a dict keyed by query of hits

    :param sam_path: str
        Path to SAM file generated by bowtie
    :param remove: bool
        Remove the SAM file after parsing

    :return {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}
        Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the
        chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches
        to the query sequence
    """
    sequences = {}

    for line in sam:
        try:
            query_seq, map_seq, rname, pos, strand, nm = _parse_sam_line(line)
        except ValueError:
            continue

        try:
            sequences[str(query_seq)].append((map_seq, rname, pos, strand, nm))
        except KeyError:
            sequences[str(query_seq)] = [(map_seq, rname, pos, strand, nm)]

    return sequences


def _parse_sam_line(line_str):
    if line_str[0] is "@":
        raise ValueError

    li_arr = line_str.strip().split()

    opt_flags = {}

    for val in li_arr:

        try:
            if val[2] is ":":
                _v_arr = val.split(":")
                opt_flags[_v_arr[0]] = _v_arr[-1]
        except IndexError:
            pass

    query_seq = Seq(li_arr[SAMTOOLS_MAP["SEQ"]])
    rname = li_arr[SAMTOOLS_MAP["RNAME"]]
    pos = int(li_arr[SAMTOOLS_MAP["POS"]])
    _flag = int(li_arr[SAMTOOLS_MAP["FLAG"]])

    try:
        nm = int(opt_flags["NM"])
    except KeyError:
        nm = 0

    try:
        map_seq = Seq(_rewrite_to_reference(query_seq, opt_flags["MD"]), alphabet=generic_dna)
    except KeyError:
        map_seq = query_seq

    if (_flag % 256) == 16:
        query_seq = query_seq.reverse_complement()
        map_seq = map_seq.reverse_complement()
        strand = "-"
    else:
        strand = "+"

    return query_seq, map_seq, rname, pos, strand, nm


def _rewrite_to_reference(seq_str, mm):
    seq_str = fix_seq_object(seq_str, make_type="str")

    _tmp_mm = ""
    pointy = 0
    srr = list(seq_str)

    for i in range(len(mm)):
        if mm[i].isalpha():
            if _tmp_mm is not "":
                pointy += int(_tmp_mm)
                _tmp_mm = ""
            srr[pointy] = mm[i]
            pointy += 1
        else:
            _tmp_mm += mm[i]

    return "".join(srr)
Makes a bowtie index in a specified directory from a FASTA file      :param dir_path: str         Path to target directory     :param genome: str         Path to sequence FASTA file      :return out_db_path:         Path to the bowtie index that can be passed to bowtie Align a query sequence to a target reference genome database and parse the resulting SAM file.      :param query: [SeqRecord]         A list of SeqRecord objects which should be aligned to the reference genome     :param db_path: str         A path to the reference genome database. Passed to bowtie2.     :param cores: int         Number of CPU cores to run. Passed to bowtie2.      :return: {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}         Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the         chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches         to the query sequence Check to make sure that the bowtie index exists. Raises a FileNotFoundError if there is any problem.      :param db_path: str         A path to the reference genome database. Passed to bowtie-inspect.     :return: bool         Returns TRUE if no error Runs bowtie2 with a list of sequencs as queries against a target reference genome database.      :param query: [SeqRecord]         A list of SeqRecord objects     :param db_path: str         A path to the reference genome database. Passed to bowtie2.     :param sam_path: str         A path to the output SAM file. If None (or not set), a tempfile will be created     :param cores: int         Number of CPU cores to run. Passed to bowtie2.      :return: str         Returns a path to the output SAM file (unchanged, unless the param sam_path was not set) # Write the queries into a string as a FASTA file # Pipe the FASTA string into bowtie 2 and pipe the output. # Stash the piped output (SAM file) in a string Takes a path to a sam file produced by bowtie and parses it into a dict keyed by query of hits      :param sam_path: str         Path to SAM file generated by bowtie     :param remove: bool         Remove the SAM file after parsing      :return {query_seq: [(mapped_seq, chr_name, position, strand, mismatch#)] : dict{str:[(Bio.Seq, str, int, str, int)]}         Returns a dict keyed by query sequence of a list of hit tuples that consist of the mapped sequence, the         chromosome of the mapped location, the position of the mapped location, the strand, and the number of mismatches         to the query sequence