Fermi's golden rule still applies in the relativistic limit, and can be rewritten in a Lorentz invariant fashion. Starting with the transition probability
$$ W_{i\rightarrow f} = \frac{2\pi}{\hbar} |m_{if}|^2 \rho(E) \,,$$ to have $W$ Lorentz invariant we'd like both the matrix element $|m_{if}|^2$ and the density of final states $\rho(E)$ to be invariant.
This can be done by shifting a few terms around. A little bit of handwaving to motivate it: The wave function $\psi$ (which is in the matrix element) has to be normalized by $\int |\psi|^2 dV = 1$, which gives us a density (of probability to encounter a particle) of $1/V$. Now, a boosted observer experiences length contraction of $1/\gamma$, which changes the density to $\gamma/V$. To obtain the correct probability again, we should re-normalize the wave function to $\psi' = \sqrt{\gamma}\,\psi $ by pulling the Lorentz factor out.
So we intoduce a new matrix element $$|{\cal M}_{if}|^2 = |m_{if}|^2 \prod_{i=1}^n (2 \gamma_i m_i c^2) =|m_{if}|^2 \prod_{i=1}^n (2E_i)^2 $$ (this is for an $n$-body process). Now the transition probability (here in differential form) becomes:
$$
dW = \frac{2\pi}{\hbar} \frac{|{\cal M}_{if}|^2}{ (2E_1)^2 (2E_2)^2 \cdots}
\cdot \frac{1}{(2\pi\hbar)^{3n}}
\, d^3p_1 \, d^3p_2 \, \cdots
\delta({p_1}^\mu + {p_2}^\mu + \ldots - {p}^\mu )
$$
The delta function is there to ensure conservation of momentum and energy. Now we can regroup the terms:
$$
\Rightarrow \quad
dW = \frac{2\pi}{\hbar} \frac{|{\cal M}_{if}|^2}{ 2E_1 2E_2 \cdot \ldots}
\cdot d_\mathrm{LIPS}
$$
The density of states/"phase space" $d\rho$ is replaced by a relativistic version, sometimes called the Lorentz invariant phase space $d_\mathrm{LIPS}$, which is given by
$$
d_\mathrm{LIPS} = \frac{1}{(2\pi\hbar)^{3n}}
\prod_{i=1}^n \frac{d^3p_i}{ 2E_i }
\delta\left(\prod_{i=1}^n {p_i}^\mu - {p}^\mu \right) \,.
$$
The nice thing about the relativistic formula for $dW$ is that, in the case you are scattering particles off one another, it immediately shows us three important contributions: not only the matrix element and phase space, but also the flux factor $1/s$ (where $s = ({p_1}^\mu + {p_2}^\mu)^2$ is the Mandelstam variable, and in case the masses are negligible, $ s \approx 2 E $). This flux factor is responsible for the general $1/Q^2$ falling slope when you plot cross section over momentum transfer $Q = \sqrt{s}$, which comes entirely from relativistic kinematics.
Hope this answers your questions. Here is a presentation (PDF) that sums it up, with an explicit proof that it is Lorentz invariant.
