text: import logging
  
  logger = logging.getLogger(__name__)
  
  import operator
  import os
  import copy
  import bisect
  import pybedtools
  import vcf
  import json
  import base64
  
  svs_of_interest = ["DEL", "INS", "DUP", "DUP:TANDEM", "INV", "ITX", "CTX"]
  sv_sources = ["Pindel", "BreakSeq", "HaplotypeCaller", "BreakDancer", "CNVnator",
                "Manta", "Lumpy", "WHAM", "CNVkit"]  # order is important!
  precise_sv_sources = ["Pindel", "BreakSeq", "HaplotypeCaller", "SoftClip"]
  sv_sources_to_type = {"Pindel": ["SR"], "BreakSeq": ["JM"], "BreakDancer": ["RP"],
                        "CNVnator": ["RD"], "HaplotypeCaller": ["AS"],
                        "Manta": ["SR", "RP"], "Lumpy": ["SR", "RP"], "CNVkit": ["RD"],
                        "WHAM": ["SR", "RP"], "SoftClip": ["SC"]}
  
  mydir = os.path.dirname(os.path.realpath(__file__))
  gaps_b37 = os.path.join(mydir, "resources/b37.gaps.bed")
  gaps_hg19 = os.path.join(mydir, "resources/hg19.gaps.bed")
  
  
  def get_gaps_file(contig_names):
      hg19_major = set(["chr%d" % i for i in xrange(1, 23)] + ["chrX", "chrY", "chrM"])
      b37_major = set([str(i) for i in xrange(1, 23)] + ["X", "Y", "MT"])
      if set(contig_names) & hg19_major: return gaps_hg19
      if set(contig_names) & b37_major: return gaps_b37
  
      logger.warn("Could not guess gaps file for reference. No gap filtering will be done.")
      return None
  
  
  class SVInterval:
      def __init__(self, chrom=None,  start=0, end=0, name=None, sv_type=None, length=0, sources=set(), gt="./1", wiggle=0,
                   info=None, cipos=[], ciend=[], native_sv=None, chrom2=None):
          self.chrom = chrom
          self.chrom2 = chrom2 if chrom2 else chrom
          self.start = start
          self.end = end
          self.name = name
          self.sv_type = sv_type
          self.length = length
          self.info = info
          self.sources = sources
          self.gt = gt
          self.wiggle = wiggle
          self.sub_intervals = []
          self.is_precise = False
          self.is_validated = False
          self.validating_interval = None
          self.cipos = cipos
          self.ciend = ciend
          self.native_sv = native_sv
          
      def merge(self, interval):
          self.start = min(self.start, interval.start - interval.wiggle)
          self.end = max(self.end, interval.end + interval.wiggle)
          self.length = max(self.length, interval.length)
          self.name = self.name + "," + interval.name
          self.sub_intervals.append(interval)
          self.sources |= interval.sources
  
      def set_merged(self, interval1, interval2):
          self.chrom = interval1.chrom
          self.start = min(interval1.start - interval1.wiggle, interval2.start - interval2.wiggle)
          self.end = max(interval1.end + interval1.wiggle, interval2.end + interval2.wiggle)
          self.length = max(interval1.length, interval2.length)
          self.name = interval1.name + "," + interval2.name
          self.sv_type = interval1.sv_type
          self.info = None
          self.sub_intervals = [interval1, interval2]
          self.sources = interval1.sources | interval2.sources
          self.gt = interval1.gt
  
      def __lt__(self, other):
          if self.chrom != other.chrom: return self.chrom < other.chrom
          if self.start != other.start: return self.start < other.start
          return self.end < other.end
  
      def __str__(self):
          if self.sub_intervals:
              return ",".join([str(interval) for interval in self.sub_intervals])
          return "%s-%d-%s-%d-%d-%s" % (self.chrom, self.start, self.chrom2, self.end, self.length, ",".join(list(self.sources)))
  
      def __repr__(self):
          return "<" + self.__class__.__name__ + " " + str(self.__dict__) + ">"
  
      def overlaps(self, other, min_fraction_self=1e-9, min_fraction_other=1e-9, min_overlap_length_self=1,
                   min_overlap_length_other=1):
          if self.sv_type not in ["ITX","CTX"] :
              if self.chrom != other.chrom:
                  return False
  
              if max(self.start - self.wiggle, other.start - other.wiggle) \
                      >= min(self.end + self.wiggle, other.end + other.wiggle):
                  return False
  
              self_length = float(self.end - self.start + 2 * self.wiggle)
              other_length = float(other.end - other.start + 2 * other.wiggle)
              overlap_length = min(self.end + self.wiggle, other.end + other.wiggle) - max(self.start - self.wiggle,
                                                                                           other.start - other.wiggle)
              return float(overlap_length) >= max(min_fraction_self * self_length,
                                                  min_fraction_other * other_length) and overlap_length >= max(
                  min_overlap_length_self, min_overlap_length_other)
          else:
              if self.sv_type == "CTX" and not (self.chrom == other.chrom and self.chrom2 == other.chrom2):
                  return False
              return ( abs(self.end - other.end) < self.wiggle ) and ( abs(self.start - other.start) < self.wiggle )
          
  
      def is_adjacent(self, other, gap=0):
          if self.sv_type in ["ITX","CTX"] :
              return False
      
          if self.chrom != other.chrom:
              return False
          return (other.start <= gap + self.end + gap < other.end) or (
              self.start <= gap + other.end + gap < self.end)
  
      def get_start(self):
          if not self.sub_intervals:
              return self.start
  
          return min([interval.get_start() for interval in self.sub_intervals])
  
      def get_end(self):
          if not self.sub_intervals:
              return self.end
  
          return max([interval.get_end() for interval in self.sub_intervals])
  
      def do_validation(self, overlap_ratio=0.5):
          self.start = self.get_start()
          self.end = self.get_end()
  
          if self.sv_type not in svs_of_interest:
              return
  
          if not self.sub_intervals:
              # This interval did not overlap anything
              self.is_precise = list(self.sources)[0] in precise_sv_sources
              # if self.sv_type == "INV" and self.length >= 1000: self.is_validated = True
              return
  
          lists = {source: [] for source in sv_sources}
  
          precise_merged = merge_intervals(
              [interval for interval in self.sub_intervals if list(interval.sources)[0] in precise_sv_sources])
  
          for interval in self.sub_intervals:
              # logger.debug("Examining subinterval " + str(interval))
              lists[list(interval.sources)[0]].append(interval)
  
          for source in sv_sources:
              # If the SV is called multiple times by the same tool, don't validate
              if len(lists[source]) == 1:
                  if self.overlaps(lists[source][0], overlap_ratio, overlap_ratio) \
                          or interval_overlaps_interval_list(lists[source][0], precise_merged, overlap_ratio,
                                                             overlap_ratio):
                      self.is_validated = len(self.sources) > 1
                      self.validating_interval = lists[source][0]
                      self.is_precise = source in precise_sv_sources
                      break  # This break makes the order of sv_sources important
  
          if not self.is_validated:
              return
  
          if self.is_precise and self.validating_interval is not None:
              self.start = self.validating_interval.start
              self.end = self.validating_interval.end
              self.length = self.validating_interval.length
              self.gt = self.validating_interval.gt
              self.info = self.validating_interval.info
  
              if self.length == 0:
                  for source in sv_sources:
                      if source in self.sources:
                          if lists[source] and lists[source][0].length > 0:
                              self.length = lists[source][0].length
                              break
  
      def fix_pos(self):
          if self.sv_type == "INS":
              mid = (self.start + self.end) / 2
              if not self.is_precise:
                  self.end = self.start
                  self.cipos = ["0", "%d" % (self.end - self.start)]
                  # self.cipos = ["%d" % (mid - self.start), "%d" % (mid - self.end)]
              else:
                  self.start = mid
                  self.end = mid
  
      def get_info(self):
          temp_info = {}
          if not self.sub_intervals:
              temp_info.update(self.info)
          else:
              if self.info:
                  temp_info.update(self.info)
              for interval in self.sub_intervals:
                  # TODO: this will just overwrite the other dict entries... this should be ok for pass variants
                  # TODO: kind of strange
                  if interval.info:
                      temp_info.update(interval.info)
          temp_info.update({"SOURCES": str(self),
                            "NUM_SVMETHODS": len(self._get_svmethods()),
                            "NUM_SVTOOLS": len(self._get_tools())})
          return temp_info
  
      def _get_tools(self):
          """Retrieve tools contributing to this structural variant call.
  
          These can come from both our current sources as well as merged sub_intervals.
          """
          tools = set(self.sources).union(*[set(x.sources) for x in self.sub_intervals])
          return sorted(list(tools))
  
      def _get_svmethods(self):
          """Retrieve structural variant methods used in all of the callers input sources.
          """
          svmethods = reduce(operator.add, [sv_sources_to_type[tool] for tool in self._get_tools()])
          return sorted(list(set(svmethods)))
  
      def to_vcf_record(self, fasta_handle=None, sample=""):
          if self.start <= 0:
              return None
          if self.sv_type not in svs_of_interest:
              return None
  
          # ignore private haplotype caller calls
          if ((not self.sub_intervals) or len(self.sources) == 1) and list(self.sources)[0] == "HaplotypeCaller":
              return None
  
          # formulate the INFO field
          info = self.get_info()
          sv_len = -self.length if self.sv_type == "DEL" else self.length
          info.update({"SVLEN": sv_len,
                       "SVTYPE": self.sv_type,
                       "SVMETHOD": ",".join(self._get_svmethods())})
          if self.sv_type in ["DEL","DUP", "ITX", "CTX"]:
              info["END"] = self.end
  
          if self.sv_type in ["ITX", "CTX"]:
              info["CHR2"] = self.chrom2
  
          if not self.is_precise:
              info.update({"IMPRECISE": True})
  
          info.update({"VT": "SV"})
          info.update({"SVTOOL": "MetaSV"})
          if self.cipos:
              info.update({"CIPOS": (",".join([str(a) for a in self.cipos]))})
  
          vcf_record = vcf.model._Record(self.chrom,
                                         self.start,
                                         ".",
                                         fasta_handle.fetch(self.chrom, max(0, self.start - 1),
                                                            max(1, self.start)) if fasta_handle else "N",
                                         [vcf.model._SV(self.sv_type)],
                                         ".",
                                         "PASS" if self.is_validated else "LowQual",
                                         info,
                                         "GT",
                                         [0],
                                         [vcf.model._Call(None, sample, vcf.model.make_calldata_tuple("GT")(GT="1/1"))])
          return vcf_record
  
      def to_bed_interval(self, sample_name):
          if self.start <= 0:
              return None
          if self.sv_type not in ["DEL", "INS", "INV", "ITX", "CTX", "DUP"]: return None
          # if not self.sub_intervals and list(self.sources)[0] == "HaplotypeCaller": return None
          # if len(self.sources) == 1 and list(self.sources)[0] == "HaplotypeCaller": return None
          if self.sv_type == "INS":
              end = self.end + 1
          elif self.sv_type in ["ITX","CTX"] :
              end = self.start + 1
          else:
              end = self.end
          info = self.get_info()
          
          if self.sv_type in ["ITX", "CTX"]:
              info["POS2"] = self.end
              info["END"] = end
              info["CHR2"] = self.chrom2
              if not self.is_precise:
                  info.update({"IMPRECISE": True})
  
          return pybedtools.Interval(self.chrom, self.start, end, name="%s,%s,%d,%s" % (
              base64.b64encode(json.dumps(info)), self.sv_type, self.length,
              ";".join(self._get_svmethods())),
              score=str(len(self.sources)))
  
      def to_svp_record(self, sample_name, id_num):
          if self.start <= 0:
              return None
          if self.sv_type not in ["DEL", "INS", "INV"]: return None
          if not self.sub_intervals and list(self.sources)[0] == "HaplotypeCaller": return None
          if len(self.sources) == 1 and list(self.sources)[0] == "HaplotypeCaller": return None
  
          start_outer = self.start
          start_inner = self.start
          end_inner = self.end
          end_outer = self.end
          type_of_computational_approach = ",".join(self._get_svmethods())
  
          return "%s\t%d\t%d\t%d\t%d\t%s\t%d\t%s\t%s\t%s\t%s\tMetaSV\t%d" % (
              self.chrom, start_outer, start_inner, end_inner, end_outer, self.sv_type, self.length, "BWA", "Illumina",
              sample_name, type_of_computational_approach, id_num)
  
  
  def interval_overlaps_interval_list(interval, interval_list, min_fraction_self=1e-9, min_fraction_other=1e-9):
      index = bisect.bisect_left(interval_list, interval)
      if index > 0 and interval.overlaps(interval_list[index - 1], min_fraction_self, min_fraction_other,
                                         min_overlap_length_self=1, min_overlap_length_other=1):
          return True
      if index < len(interval_list) and interval.overlaps(interval_list[index], min_fraction_self, min_fraction_other,
                                                          min_overlap_length_self=1, min_overlap_length_other=1):
          return True
      return False
  
  
  def merge_intervals(interval_list):
      interval_list.sort()
      merged_intervals = []
      if not interval_list:
          return []
  
      current_merged_interval = copy.deepcopy(interval_list[0])
      for i in xrange(len(interval_list) - 1):
          next_interval = interval_list[i + 1]
  
          if current_merged_interval.overlaps(next_interval) or current_merged_interval.is_adjacent(next_interval):
              if current_merged_interval.sub_intervals:
                  current_merged_interval.merge(next_interval)
              else:
                  new_merged_interval = SVInterval()
                  # logger.debug("Merging %s with %s" % (repr(current_merged_interval), repr(next_interval)))
                  new_merged_interval.set_merged(current_merged_interval, next_interval)
                  current_merged_interval = new_merged_interval
          else:
              merged_intervals.append(current_merged_interval)
              current_merged_interval = copy.deepcopy(next_interval)
  
      merged_intervals.append(current_merged_interval)
      merged_intervals.sort()
      return merged_intervals
  
  
  def merge_intervals_recursively(interval_list,overlap_ratio):
      merged_intervals = merge_intervals(interval_list)
  
      # Intervals which overlap well with merged_intervals
      intervals1 = []
      # Intervals which do not overlap well with merged_intervals.
      # Used to filter out small intervals which got merged with large intervals
      intervals2 = []
      for interval in interval_list:
          if interval_overlaps_interval_list(interval, merged_intervals, overlap_ratio,overlap_ratio ):
              intervals2.append(interval)
          else:
              intervals1.append(interval)
      if len(intervals1)==0:
          return merged_intervals
      if len(intervals2)==0:
          intervals1.sort()
          return intervals1
      final_merged = merge_intervals_recursively(intervals1,overlap_ratio) + merge_intervals_recursively(intervals2,overlap_ratio)
      final_merged.sort()
      return final_merged
  