def SH_raxml(reftree, querytree, phy_file, shout="SH_out.txt"):
    """
    SH test using RAxML

    querytree can be a single tree or a bunch of trees (eg. from bootstrapping)
    """
    assert op.isfile(reftree)
    shout = must_open(shout, "a")

    raxml_work = op.abspath(op.join(op.dirname(phy_file), "raxml_work"))
    mkdir(raxml_work)
    raxml_cl = RaxmlCommandline(cmd=RAXML_BIN("raxmlHPC"), \
    sequences=phy_file, algorithm="h", model="GTRGAMMA", \
    name="SH", starting_tree=reftree, bipartition_filename=querytree, \
    working_dir=raxml_work)

    logging.debug("Running SH test in RAxML: %s" % raxml_cl)
    o, stderr = raxml_cl()
    # hard coded
    try:
        pval = re.search('(Significantly.*:.*)', o).group(0)
    except:
        print("SH test failed.", file=sys.stderr)
    else:
        pval = pval.strip().replace("\t"," ").replace("%","\%")
        print("{0}\t{1}".format(op.basename(querytree), pval), file=shout)
        logging.debug("SH p-value appended to %s" % shout.name)

    shout.close()
    return shout.name
SH test using RAxML

    querytree can be a single tree or a bunch of trees (eg. from bootstrapping)