The Bunch-Davies (BD) vacuum is a "pure state" when it comes to all the creation and annihilation operators in the patch where the BD vacuum is relevant. Inside this patch, we can say what operators annihilate the BD state which uniquely specifies it (well, for massive fields, one has the $\alpha$-vacua as generalizations).
However, there's also a thermal radiation coming from the boundary of this patch. This is the reason for the non-conservation of the particle number. Conceptually, this non-conservation is equivalent to the ordinary thermal ensemble in the Minkowski space - which is a "mixed state". So the two-point functions are essentially thermal two-point functions, and so on.
The correlation functions are computed via Feynman rules that are analogous to those in the Minkowskian thermal state. We don't compute scattering amplitudes because there is no scattering S-matrix in the de Sitter space. That's because an observer's causal diamond misses almost all of the asymptotic regions where global in- and out-states may be defined.
In particular, it is not true that the Feynman diagrams for thermal field theory require the particle numbers to be well-defined (thermal states are always mixtures of different particle number states) or conserved (interactions always change the number of particles, and interactions had to be incorporated to get the right thermal state, anyway). The number of particles is not a good quantum number in violently cosmological backgrounds and one will find nothing interesting if he or she focuses on this bad quantum number, assuming that it carries some wisdom.
At the level of field theory in a curved spacetime, the right rules are known and analogous to thermal QFT; correlation functions and not (unmeasurable) amplitudes may be calculated. When gravity is made dynamical, the calculations only work at the semiclassical level.
Beyond semiclassical QFT in curved spaces
There are subtle things and no one knows the totally rigorous way to describe the de Sitter space. In particular, we don't know whether the relevant Hilbert space is finite-dimensional (as suggested by the finite de Sitter entropy) or infinite-dimensional; whether the degrees of freedom behind the cosmic horizon are "independent" or just reparametrizations of the degrees of freedom inside the visible universe, because of the complementarity principle; whether there can be a fully predictive framework or not, and so on.
Nevertheless, at the same moment, one should realize that all these uncertainties are "very tiny" when the de Sitter radius is large. The thermal wavelength is comparable to the de Sitter radius and the expected number of thermal particles within the whole visible part of the de Sitter space is of order one.
So we're really uncertain about roughly one photon and one graviton that may or may not live somewhere in the whole visible Universe and whose energy is so tiny that the wavelength is comparable to the de Sitter radius, too. In any practical context, even during inflation, this is a very small effect.