code: def mergebam(args):
      
      p = OptionParser(mergebam.__doc__)
      p.set_cpus()
      opts, args = p.parse_args(args)
  
      if len(args) not in (2, 3):
          sys.exit(not p.print_help())
  
      if len(args) == 2:
          idir1, outdir = args
          dir1 = [idir1] if idir1.endswith(".bam") else iglob(idir1, "*.bam")
          logging.debug("Homozygous mode")
          dir2 = [""] * len(dir1)
      elif len(args) == 3:
          idir1, idir2, outdir = args
          dir1 = [idir1] if idir1.endswith(".bam") else iglob(idir1, "*.bam")
          dir2 = [idir2] if idir2.endswith(".bam") else iglob(idir2, "*.bam")
          assert len(dir2) == 1, "Second pile must contain a single bam"
          dir2 = [idir2] * len(dir1)
  
      assert len(dir1) == len(dir2), "Two piles must contain same number of bams"
      cmd = "samtools merge {} {} {} && samtools index {}"
      cmds = []
      mkdir(outdir)
      for a, b in zip(dir1, dir2):
          ia = op.basename(a).split(".")[0]
          ib = op.basename(b).split(".")[0] if b else ia
          outfile = op.join(outdir, "{}_{}.bam".format(ia, ib))
          cmds.append(cmd.format(outfile, a, b, outfile))
  
      p = Parallel(cmds, cpus=opts.cpus)
      p.run()
docs: %prog mergebam dir1 homo_outdir
      or
      %prog mergebam dir1 dir2/20.bam het_outdir
  
      Merge sets of BAMs to make diploid. Two modes:
      - Homozygous mode: pair-up the bams in the two folders and merge
      - Heterozygous mode: pair the bams in first folder with a particular bam