code: def clean_bam_file(bam_in, mask=None):
      
      seq_obj = defaultdict(int)
      if mask:
          mask_file = op.splitext(bam_in)[0] + "_mask.bam"
          if not file_exists(mask_file):
              pybedtools.BedTool(bam_file).intersect(b=mask, v=True).saveas(mask_file)
          bam_in = mask_file
      out_file = op.splitext(bam_in)[0] + "_rmlw.bam"
      # bam.index(bam_in, {'algorithm':{}})
      run("samtools index %s" % bam_in)
      if not file_exists(bam_in + ".bai"):
          raise IOError("Failed to created bam index of %s. Try to do it manually" % bam_in)
      bam_handle = pysam.AlignmentFile(bam_in, "rb")
      with pysam.AlignmentFile(out_file, "wb", template=bam_handle) as out_handle:
          for read in bam_handle.fetch():
              seq_name = int(read.query_name.replace('seq_', ''))
              match_size = [nts for oper, nts in read.cigartuples if oper == 0]
              subs_size = [nts for oper, nts in read.cigartuples if oper == 4]
              if match_size[0] < 17:
                  continue
              if subs_size:
                  if subs_size[0] > 3:
                      continue
              try:
                  nh = read.get_tag('NH')
              except KeyError:
                  nh = 1
              seq_obj[seq_name] = sequence(seq_name)
              seq_obj[seq_name].align = nh
              out_handle.write(read)
      return out_file, seq_obj
docs: Remove from alignment reads with low counts and highly # of hits