def mconsensus(args):
    """
    %prog mconsensus *.consensus

    Call consensus along the stacks from cross-sample clustering.
    """
    p = OptionParser(mconsensus.__doc__)
    p.add_option("--allele_counts", default="allele_counts",
                 help="Directory to generate allele counts")
    add_consensus_options(p)
    opts, args = p.parse_args(args)

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

    consensusfiles = args
    prefix = opts.prefix
    acdir = opts.allele_counts
    store = ClustStores(consensusfiles)
    pctid = find_pctid(consensusfiles)
    pf = prefix + ".P{0}".format(pctid)

    clustSfile = pf + ".clustS"
    AC = makeloci(clustSfile, store, prefix,
                  minsamp=opts.minsamp, pctid=pctid)

    mkdir(acdir)
    acfile = pf + ".allele_counts"
    fw = open(acfile, "w")
    seen = DefaultOrderedDict(list)        # chr, pos => taxa
    print("# " + "\t".join(ACHEADER), file=fw)
    # Sort allele counts into separate files
    for ac in AC:
        chrpos = ac.chr, ac.pos
        seen[chrpos].append(ac)
        print(ac.tostring(taxon=True), file=fw)
    fw.close()

    logging.debug("Populate all taxa and instantiate empty vector if missing")
    all_taxa = set([op.basename(x).split(".")[0] for x in consensusfiles])
    taxon_to_ac = defaultdict(list)
    for chrpos, aclist in seen.items():
        included_taxa = set([x.taxon for x in aclist])
        missing_taxa = all_taxa - included_taxa
        template = deepcopy(aclist[0])
        template.clear()
        for ac in aclist:
            taxon_to_ac[ac.taxon].append(ac)
        for tx in missing_taxa:
            taxon_to_ac[tx].append(template)

    logging.debug("Write allele counts for all taxa")
    for tx, aclist in sorted(taxon_to_ac.items()):
        tx_acfile = op.join(acdir, tx + ".allele_counts")
        fw = open(tx_acfile, "w")
        print("# " + "\t".join(ACHEADER_NO_TAXON), file=fw)
        for ac in aclist:
            print(ac.tostring(), file=fw)
        fw.close()
        logging.debug("Written {0} sites in `{1}`".\
                format(len(aclist), tx_acfile))
%prog mconsensus *.consensus

    Call consensus along the stacks from cross-sample clustering.