What cannot be calculated using density matrix? Let's say I have a compound quantum system (CQS) in an (unknown to me) pure state $\left|\Psi\right>$. If an operator $\mathbf{A}$ acts only on variables of a subsystem (S) of CQS, then I can calculate expectation value of the corresponding observable as
$$
    \bar{A} = \textrm{Tr} (\mathbf{\rho}_S \mathbf{A})
$$
where $\mathbf{\rho}_S$ is a partial density matrix of S (trace over other subsystems of CQS). I can do this for any observable (that acts only on S).
So, if I limit myself to S what is it that I cannot calculate using $\mathbf{\rho}_S$ that I'd be able to calculate if I knew $\left|\Psi\right>$? It seems such things should exist, otherwise I'd be able to construct the wave function for S (is that correct?)
 A: Below I explain how one would reconstruct a "wavefunction" from a density matrix known to come from a pure state. However, your situation is different. If $\rho$ is the density matrix of an entangled pure state of two systemsn $S$ and $S'$, then $\rho_S$ is the density matrix of a mixed state of $S$, see also this question and its answers. A mixed state does not have a unique wavefunction, so the question of whether you can reconstruct the wavefunction is meaningless. By design, if you are only interested in observables of $S$, then $\rho_S$ contains all necessary information for that - the partial trace is the reverse of the operation of "combining" quantum systems, and if you apply it to a non-entangled state of the combined system it just gives back the corresponding pure state of $S$, which obviously contains all necessary information about $S$.

You are able to construct the wave-function for a density matrix $\rho$ if you know that it comes from a pure state that has such a wavefunction. Or at least, you are able to construct the only physically meaningful content of the wavefunction:
You can compute how likely it is that the system is in a region $R\subset\mathbb{R}^n$ by computing the expectation value of the associated spectral projector $P_R$ (which acts on wavefunctions by multiplying them with the characteristic function of $R$; in non-rigorous language you may think of this as the operator that projects onto the space spanned by position eigenstates with eigenvalues in $R$), and this is what the wavefunction is - a probability density that tells you how likely it is the system is in a region $R$ by integrating the wavefunction over it. The values of the wavefunction at single points are physically irrelevant and therefore you should not expect to be able to recover them - the wavefunction is formally a representant of an equivalence class of square-integrable functions in $L^2(\mathbb{R}^n)$, not a single function whose values at individual points would carry any meaning - all that matters is what happens when you integrate it over regions of non-zero measure, and that information is recoverable from the density matrix.
A: 
So, if I limit myself to S what is it that I cannot calculate using $\mathbf{\rho}_S$ that I'd be able to calculate if I knew $\left|\Psi\right>$? 

Locally to $S$? Nothing. As you've noted, the partial trace to $S$ gives you complete information about all local measurements.

It seems such things should exist, otherwise I'd be able to construct the wave function for S (is that correct?)

The problem here is that there might not be such a thing as a "wavefunction for $S$": even if the multipartite system is in a pure state $|\Psi\rangle$, if the state is entangled then there are states (and measurements on those states) that cannot be described by any wavefunction on the $S$ Hilbert space.
For an easy example, consider two entangled qubits on the Bell triplet state
$$
\left|\Psi\right\rangle = \frac{|0\rangle|0\rangle+|1\rangle|1\rangle}{\sqrt{2}}
$$
along with the local obervables $A_i=\sigma_i$, the Pauli matrices. Then the partial trace to the first qubit gives you the reduced density matrix
$$
\rho_s = \frac{|0\rangle\langle0|+|1\rangle\langle1|}{2},
$$
which has zero expectation value for all three $A_i$s, which is impossible for a pure state. If you can describe $\rho_s$ in the form $|\psi\rangle\langle\psi|$, then by all means go for it, but in general this is not possible.
