This is just about master-mind specialty topic (much more so than the coding business anyway ;), jpromvi's answer is pretty much covers it, but I thought I'd add some more details.
You basic choices are
- One big alignment (usually with different mutation models for each gene)
- Estimate a 'gene tree' for each locus and co-estimate a species tree
Number one will likely be a lot quicker, but it can run into all sorts of trouble (because a few strong signals from genes that don't represent the organismal relationships can override the species tree).
Which use will depend on your question (you might not even care about the species-level relationships) and perhaps the shape of your tree (short internal branches increase the probability of incomplete lineage sorting, one of the sources of error).
If you go down the gene-tree route then the three choices you mention, BEAST, BEST and STEM are the main options. STEM relies on trees estimated by some other method (so you reconcile your gene trees with a species designation) it doesn't work for my question because you need branch lengths in "coalsencent units" (generations/effective population size) and don't have data to reliable estimate that (so I don't know much about it).
BEAST and BEST take alignments, and do much the same thing. I've found the default priors in BEST make it almost impossible for the MCMC runs to converge (especially the population size prior) and BEAST runs a good deal faster. (The species tree thing isn't well documented in BEAST, but you just need to import species names as a "trait" in BEAUti, the gui BEAST file generator).
There, have I managed to thoroughly confuse you?