code: def SH_raxml(reftree, querytree, phy_file, shout="SH_out.txt"):
      
      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
docs: SH test using RAxML
  
      querytree can be a single tree or a bunch of trees (eg. from bootstrapping)