This paper considers open quantum systems whose dynamic variables satisfy canonical commutation relations and are governed by Markovian Hudson-Parthasarathy quantum stochastic differential equations driven by external bosonic fields. The dependence of the Hamiltonian and the system-field coupling operators on the system variables is represented using the Weyl functional calculus. This leads to an integro-differential equation (IDE) for the evolution of the quasi-characteristic function (QCF) which encodes the dynamics of mixed moments of the system variables. Unlike quantum master equations for reduced density operators, this IDE involves only complex-valued functions on finite-dimensional Euclidean spaces and extends the Wigner-Moyal phase-space approach for quantum stochastic systems. The dynamics of the QCF and the related Wigner quasi-probability density function (QPDF) are discussed in more detail for the case when the coupling operators depend linearly on the system variables and the Hamiltonian has a nonquadratic part represented in the Weyl quantization form. For this class of quantum stochastic systems, we also consider an approximate computation of invariant states and discuss the deviation from Gaussian quantum states in terms of the $χ^2$-divergence (or the second-order Renyi relative entropy) applied to the QPDF. The results of the paper may find applications to investigating different aspects of the moment stability, relaxation dynamics and invariant states in open quantum systems.