The procedure would be the same as the procedure for finding the uncertainty principle for pure states. We define $\Delta O$ as $\sqrt{\langle O^2\rangle- \langle O \rangle^2 }$. Since the two expectation values that go into calculating the uncertainty can be calculated from the density matrix via the usual expression for the expectation value of an operator $\langle O\rangle = \text{Tr}(O\rho)$.
The expression for $\Delta O$ as $\sqrt{\langle O^2\rangle- \langle O \rangle^2 }$ is so hostile to linearity that we cannot hope for a simple linear expression of the kind you propose for the uncertainty of an observable. This has been explicitly demonstrated in the answer by @user1723984. However, the Robertson-Schrodinger uncertainty relation for the product of the uncertainties of two observables has a nicer expression for a pure state $\psi$, in particular,
\begin{align}\sigma^2_A\sigma^2_B \ge {1\over 4}\Bigl|\langle\psi|[A,B]|\psi\rangle\Bigr|^2\end{align}
as you notice. Let's see if we can recover an analog of this relation for the case of the density matrix. We write
\begin{align}
\sigma^2_A&=\Bigr(\text{Tr}\big(\rho A^2\big)-\text{Tr}\big(\rho A\big)\text{Tr}\big(\rho A\big)\Bigr)\\
&=\sum_i\bigg\langle i\bigg|\rho A^2-\sum_j\Big\langle j\Big|\rho A\Big|j\Big\rangle \rho A\bigg|i\bigg\rangle\\
&=\sum_i\bigg\langle i\bigg|\rho A\Big(A-\sum_j\Big\langle j\Big|\rho A\Big|j\Big\rangle \Big) \bigg|i\bigg\rangle\\
&=\text{Tr}\Big(\rho A\big(A-\langle A\rangle\big)\Big)\\
&=\text{Tr}\Big(\rho\big(A-\langle A\rangle\big)^2\Big)\\
\end{align}
Thus,
\begin{align}
\sigma^2_A\sigma^2_B&=\text{Tr}\Big(\rho\big(A-\langle A\rangle\big)^2\Big)\text{Tr}\Big(\rho\big(B-\langle B\rangle\big)^2\Big)\\
\end{align}
Now, the Cauchy-Schwarz inequalities hold true for all inner products, and it can be shown that for Hermitian operators $X$ and $Y$, it takes the following form
\begin{align}
\big|\text{Tr}\big(SX Y\big)\big|^2 \leq \text{Tr}\big(SX^2\big)\text{Tr}\big(SY^2\big)
\end{align}
where $S$ is a semi-positive definite matrix. Since $\rho$ is a semi-positive definite matrix and $A-\langle A\rangle$, and $B-\langle B\rangle$ are Hermitian operators, an application of this inequality yields
\begin{align}
\sigma_A^2\sigma_B^2&\geq \bigg|\text{Tr}\Big(\rho \big(A-\langle A\rangle\big)\big(B-\langle B\rangle\big)\Big)\bigg|^2\\
&= \Big|\text{Tr}( \rho AB ) - \text{Tr} (\rho A)\text{Tr} (\rho B)\Big|^2
\end{align}
Exploiting the facts that
$1.$ modulus square should be bigger than the square of its imaginary part,
$2.$ $\overline{\text{Tr}(\rho A B)}=\text{Tr}(\rho BA)$, and
$3.$ $\overline{\text{Tr}(\rho A)}=\text{Tr}(\rho A)$, we can write
\begin{align}
\sigma_A^2\sigma_B^2&\geq\frac{1}{4}\Big|{\text{Tr}(\rho AB)-\text{Tr}(\rho BA)}\Big|^2\\
\implies\sigma_A\sigma_B&\geq\frac{1}{2}\Big|\text{Tr}(\rho [A,B])\Big|
\end{align}
which is the exact equivalent of the Robertson-Schrodinger uncertainty relation. Now, expanding the density matrix as $\rho = \sum_n p_n|\psi_n\rangle\langle \psi_n |$, we can easily see that $$\text{Tr}(\rho [A,B])=\sum_np_n\langle\psi_n|[A,B]|\psi_n\rangle$$
Thus, even tho not for individual uncertainties, for the uncertainty principle, we can write down
$$\Delta_\rho A\Delta_\rho B \geq \sum_n p_n \Delta_{\psi_n} A \Delta_{\psi_n} B$$
where I've switched the notation in a self-explanatory manner.
I might have made some sign/modulus mistakes, so kindly verify.
