def alts(args):
    """
    %prog alts HD

    Build alternative loci based on simulation data.
    """
    import pysam
    from jcvi.utils.iter import pairwise
    from jcvi.utils.grouper import Grouper

    p = OptionParser(alts.__doc__)
    p.set_outfile(outfile="TREDs.alts.csv")
    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(not p.print_help())

    treds = args
    repo = TREDsRepo()
    if "all" in treds:
        treds = repo.names

    pad_left, pad_right = 1000, 10000
    READLEN = 150
    fw = must_open(opts.outfile, "w")
    print("TRED,alts,alts.hg19", file=fw)  # Header
    for tred in treds:
        ref_regions = []

        # Simulate a depth 1000 BAM with 300 repeats
        for ref in ("hg38", "hg19"):

            # This is the region that involves the TRED locus
            repo = TREDsRepo(ref=ref)
            t = repo[tred]
            chr, start, end = t.chr, t.repeat_start, t.repeat_end
            start -= pad_left
            end += pad_right

            tred_ref = "{}_{}".format(tred, ref)
            if not op.isdir(tred_ref):
                simulate([tred_ref, "300", "300", "--depth=1000",
                          "--ref={}".format(ref), "--tred={}".format(tred)])
            bamfile = op.join(tred_ref, "300.bam")

            # Parse the BAM file, retrieve all regions
            bamfile = pysam.AlignmentFile(bamfile, "rb")
            nreads = altreads = 0
            alt_points = set()
            for read in bamfile.fetch():
                fname, fstart = bamfile.getrname(read.reference_id), read.reference_start
                rname, rstart = bamfile.getrname(read.next_reference_id), read.next_reference_start
                f_in_region = in_region(fname, fstart, chr, start, end)
                r_in_region = in_region(rname, rstart, chr, start, end)
                if (not f_in_region) and r_in_region:
                    alt_points.add((fname, fstart))
                    altreads += 1
                if (not r_in_region) and f_in_region:
                    alt_points.add((rname, rstart))
                    altreads += 1
                nreads += 1

            logging.debug("A total of {} reads ({} alts) processed".\
                        format(nreads, altreads))
            alt_points = natsorted(alt_points)

            # Chain these points together into regions
            g = Grouper()
            for a in alt_points:
                g.join(a)
            for a, b in pairwise(alt_points):
                achr, apos = a
                bchr, bpos = b
                if achr != bchr:
                    continue
                if (bpos - apos) > READLEN:
                    continue
                g.join(a, b)

            # All regions that contain ALT
            alt_sum = 0
            regions = []
            for c in g:
                chr_min, pos_min = min(c)
                chr_max, pos_max = max(c)
                assert chr_min, chr_max
                pos_min -= READLEN
                pos_max += READLEN
                regions.append((chr_min, pos_min, pos_max))
                alt_sum += pos_max - pos_min

            regions = "|".join(["{}:{}-{}".format(c, start, end) \
                            for c, start, end in natsorted(regions)])
            ref_regions.append(regions)

        line = ",".join([tred] + ref_regions)
        print(line, file=sys.stderr)
        print(line, file=fw)
        logging.debug("Alternative region sum: {} bp".format(alt_sum))

    fw.close()