If $\Psi(q,x)$ is the larger system wave function, with $x$ the subsystem coordinates and $q$ the coordinates outside the system, then the average of an operator $\hat{A}$ acting on the $x$ coordinates of $\Psi$ is given by
\begin{align}
<\hat{A}> &= \int dx dq \Psi^*(x,q) \hat{A} \Psi(x,q) \\
&= \int dx dq \Psi^*(x',q) \hat{A} \Psi(x,q)|_{x'=x} \\
&= \int dx \hat{A} \int dq \Psi^*(x',q) \Psi(x,q)|_{x'=x} \\
&= \int dx \hat{A} \rho(x',x)|_{x'=x}
\end{align}
What's this? The density matrix...