def filterloci(args):
    """
    %prog filterloci allele_freq STR-exons.wo.bed SAMPLES

    Filter subset of loci that satisfy:
    1. no redundancy (unique chr:pos)
    2. variable (n_alleles > 1)
    3. low level of missing data (>= 50% autosomal + X, > 25% for Y)

    Write meta file with the following infor:
    1. id
    2. title
    3. gene_name
    4. variant_type
    5. motif
    6. allele_frequency

    `STR-exons.wo.bed` can be generated like this:
    $ tail -n 854476 /mnt/software/lobSTR/hg38/index.tab | cut -f1-3 > all-STR.bed
    $ intersectBed -a all-STR.bed -b all-exons.bed -wo > STR-exons.wo.bed
    """
    p = OptionParser(filterloci.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    af, wobed, samples = args
    nsamples = len([x.strip() for x in open(samples)])
    nalleles = nsamples * 2
    logging.debug("Load {} samples ({} alleles) from `{}`".\
                    format(nsamples, nalleles, samples))

    logging.debug("Load gene intersections from `{}`".format(wobed))
    fp = open(wobed)
    gene_map = defaultdict(set)
    for row in fp:
        chr1, start1, end1, chr2, start2, end2, name, ov = row.split()
        gene_map[(chr1, start1)] |= set(name.split(","))
    for k, v in gene_map.items():
        non_enst = sorted(x for x in v if not x.startswith("ENST"))
        enst = sorted(x.rsplit(".", 1)[0] for x in v if x.startswith("ENST"))
        gene_map[k] = ",".join(non_enst / enst)

    logging.debug("Filtering loci from `{}`".format(af))
    fp = open(af)
    seen = set(TREDS)
    remove = []
    fw = open("meta.tsv", "w")
    header = "id title gene_name variant_type motif allele_frequency".\
                replace(" ",  "\t")
    print >> fw, header
    variant_type = "short tandem repeats"
    title = "Short tandem repeats ({})n"
    for row in fp:
        sname, counts = row.split()
        name = sname.rsplit("_", 1)[0]
        seqid, pos, motif = name.split("_")
        countst = [x for x in counts.strip("{}").split(",") if x]
        countsd = {}
        for x in countst:
            a, b = x.split(":")
            countsd[int(a)] = int(b)

        if counts_filter(countsd, nalleles, seqid):
            remove.append(sname)
            continue

        if name in seen:
            remove.append(sname)
            continue
        seen.add(name)

        gene_name = gene_map.get((seqid, pos), "")
        print >> fw, "\t".join((name, title.format(motif),
                                gene_name, variant_type, motif, counts))
    fw.close()

    removeidsfile = "remove.ids"
    fw = open(removeidsfile, "w")
    print >> fw, "\n".join(remove)
    fw.close()
    logging.debug("A total of {} filtered loci written to `{}`".\
                    format(len(remove), removeidsfile))
dataset/ETHPy150Open tanghaibao/jcvi/variation/str.py/filterloci/operator `+` -> `/`