Can you find an analytic expression for the (beta)hfe of a BJT?
The equation you are thinking of is from the Ebers-Moll BJT model, simplified for the forward-active operating region:

$$I_E = I_{ES}\left(\exp\left(\frac{V_{BE}}{V_T}\right) - 1 \right)$$

Note this gives the emitter current, not the collector current.

The other two equations in the Ebers-Moll model are

$$I_C = \alpha{}I_E$$

$$I_B = (1-\alpha)I_E$$

Unfortunately these equations don't help you find \$\beta\$ because \$\alpha\$ is just \$\beta\$ in disguise:

$$\beta=\frac{I_C}{I_B}=\frac{\alpha}{1-\alpha}$$


  Small related question, is the \$I_S\$ ... stated anywhere in a transistor datasheet?


It's usually not in the datasheet, but if a SPICE model is available for your device, you can get a "typical" value from there. Look for the IS parameter in the MODEL card.