text: #(c) 2016 by Authors
  #This file is a part of ABruijn program.
  #Released under the BSD license (see LICENSE file)
  
  """
  Runs Minimap2 and parses its output
  """
  
  from __future__ import absolute_import
  from __future__ import division
  import os
  from collections import namedtuple
  import subprocess
  import logging
  
  import flye.utils.fasta_parser as fp
  from flye.utils.utils import which
  from flye.utils.sam_parser import AlignmentException, preprocess_sam
  from flye.six import iteritems
  from flye.six.moves import range
  
  
  logger = logging.getLogger()
  MINIMAP_BIN = "flye-minimap2"
  
  
  ContigInfo = namedtuple("ContigInfo", ["id", "length", "type"])
  
  
  def check_binaries():
      if not which(MINIMAP_BIN):
          raise AlignmentException("Minimap2 is not installed")
      if not which("sort"):
          raise AlignmentException("UNIX sort utility is not available")
  
  
  def make_alignment(reference_file, reads_file, num_proc,
                     work_dir, platform, out_alignment, reference_mode,
                     sam_output):
      """
      Runs minimap2 and sorts its output
      """
      minimap_ref_mode = {False: "ava", True: "map"}
      minimap_reads_mode = {"nano": "ont", "pacbio": "pb"}
      mode = minimap_ref_mode[reference_mode] + "-" + minimap_reads_mode[platform]
  
      _run_minimap(reference_file, reads_file, num_proc, mode,
                   out_alignment, sam_output)
  
      if sam_output:
          preprocess_sam(out_alignment, work_dir)
  
  
  def get_contigs_info(contigs_file):
      contigs_info = {}
      contigs_fasta = fp.read_sequence_dict(contigs_file)
      for ctg_id, ctg_seq in iteritems(contigs_fasta):
          contig_type = ctg_id.split("_")[0]
          contigs_info[ctg_id] = ContigInfo(ctg_id, len(ctg_seq),
                                            contig_type)
  
      return contigs_info
  
  
  def shift_gaps(seq_trg, seq_qry):
      """
      Shifts all ambigious query gaps to the right
      """
      lst_trg, lst_qry = list("$" + seq_trg + "$"), list("$" + seq_qry + "$")
      is_gap = False
      gap_start = 0
      for i in range(len(lst_trg)):
          if is_gap and lst_qry[i] != "-":
              is_gap = False
              swap_left = gap_start - 1
              swap_right = i - 1
  
              while (swap_left > 0 and swap_right >= gap_start and
                     lst_qry[swap_left] == lst_trg[swap_right]):
                  lst_qry[swap_left], lst_qry[swap_right] = \
                              lst_qry[swap_right], lst_qry[swap_left]
                  swap_left -= 1
                  swap_right -= 1
  
          if not is_gap and lst_qry[i] == "-":
              is_gap = True
              gap_start = i
  
      return "".join(lst_qry[1 : -1])
  
  
  def get_uniform_alignments(alignments, seq_len):
      """
      Leaves top alignments for each position within contig
      assuming uniform coverage distribution
      """
      def _get_median(lst):
          if not lst:
              raise ValueError("_get_median() arg is an empty sequence")
          sorted_list = sorted(lst)
          if len(lst) % 2 == 1:
              return sorted_list[len(lst) // 2]
          else:
              mid1 = sorted_list[(len(lst) // 2) - 1]
              mid2 = sorted_list[(len(lst) // 2)]
              return (mid1 + mid2) / 2
  
      WINDOW = 100
      MIN_COV = 10
      COV_RATE = 1.25
  
      #split contig into windows, get median read coverage over all windows and
      #determine the quality threshold cutoffs for each window
      wnd_primary_cov = [0 for _ in range(seq_len // WINDOW + 1)]
      wnd_aln_quality = [[] for _ in range(seq_len // WINDOW + 1)]
      wnd_qual_thresholds = [1.0 for _ in range(seq_len // WINDOW + 1)]
      for aln in alignments:
          for i in range(aln.trg_start // WINDOW, aln.trg_end // WINDOW):
              if not aln.is_secondary:
                  wnd_primary_cov[i] += 1
              wnd_aln_quality[i].append(aln.err_rate)
  
      #for each window, select top X alignmetns, where X is the median read coverage
      cov_threshold = max(int(COV_RATE * _get_median(wnd_primary_cov)), MIN_COV)
      for i in range(len(wnd_aln_quality)):
          if len(wnd_aln_quality[i]) > cov_threshold:
              wnd_qual_thresholds[i] = sorted(wnd_aln_quality[i])[cov_threshold]
  
      #for each alignment, count in how many windows it passes the threshold
      filtered_alignments = []
      total_sequence = 0
      filtered_sequence = 0
      for aln in alignments:
          good_windows = 0
          total_windows = aln.trg_end // WINDOW - aln.trg_start // WINDOW
          total_sequence += aln.trg_end - aln.trg_start
          for i in range(aln.trg_start // WINDOW, aln.trg_end // WINDOW):
              if aln.err_rate <= wnd_qual_thresholds[i]:
                  good_windows += 1
  
          if good_windows > total_windows // 2:
              filtered_alignments.append(aln)
              filtered_sequence += aln.trg_end - aln.trg_start
  
      #filtered_reads_rate = 1 - float(len(filtered_alignments)) / len(alignments)
      #filtered_seq_rate = 1 - float(filtered_sequence) / total_sequence
      #logger.debug("Filtered {0:7.2f}% reads, {1:7.2f}% sequence"
      #                .format(filtered_reads_rate * 100, filtered_seq_rate * 100))
  
      return filtered_alignments
  
  
  def split_into_chunks(fasta_in, chunk_size):
      out_dict = {}
      for header, seq in iteritems(fasta_in):
          #print len(seq)
          for i in range(0, max(len(seq) // chunk_size, 1)):
              chunk_hdr = "{0}$chunk_{1}".format(header, i)
              start = i * chunk_size
              end = (i + 1) * chunk_size
              if len(seq) - end < chunk_size:
                  end = len(seq)
  
              #print(start, end)
              out_dict[chunk_hdr] = seq[start : end]
  
      return out_dict
  
  
  def merge_chunks(fasta_in, fold_function=lambda l: "".join(l)):
      """
      Merges sequence chunks. Chunk names are in format `orig_name$chunk_id`.
      Each chunk is as dictionary entry. Value type is arbitrary and
      one can supply a custom fold function
      """
      def name_split(h):
          orig_hdr, chunk_id = h.rsplit("$", 1)
          return orig_hdr, int(chunk_id.rsplit("_", 1)[1])
  
      out_dict = {}
      cur_seq = []
      cur_contig = None
      for hdr in sorted(fasta_in, key=name_split):
          orig_name, dummy_chunk_id = name_split(hdr)
          if orig_name != cur_contig:
              if cur_contig != None:
                  out_dict[cur_contig] = fold_function(cur_seq)
              cur_seq = []
              cur_contig = orig_name
          cur_seq.append(fasta_in[hdr])
  
      if cur_seq:
          out_dict[cur_contig] = fold_function(cur_seq)
  
      return out_dict
  
  
  def _run_minimap(reference_file, reads_files, num_proc, mode, out_file,
                   sam_output):
      cmdline = [MINIMAP_BIN, reference_file]
      cmdline.extend(reads_files)
      cmdline.extend(["-x", mode, "-t", str(num_proc)])
      if sam_output:
          #a = SAM output, p = min primary-to-seconday score
          #N = max secondary alignments
          cmdline.extend(["-a", "-p", "0.7", "-N", "10"])
  
      try:
          devnull = open(os.devnull, "wb")
          #logger.debug("Running: " + " ".join(cmdline))
          subprocess.check_call(cmdline, stderr=devnull,
                                stdout=open(out_file, "wb"))
      except (subprocess.CalledProcessError, OSError) as e:
          if e.returncode == -9:
              logger.error("Looks like the system ran out of memory")
          raise AlignmentException(str(e))
  