$\newcommand\Si\Sigma\newcommand\si\sigma\newcommand\R{\mathbb R}\newcommand{\la}{\lambda}$ The calculations can be done as follows.
Consider the spectral decomposition of $\Si$: \begin{equation} \Si=QDQ^T,\tag{0} \end{equation} where $Q$ is an orthogonal matrix $n\times n$ matrix and $D$ is the diagonal matrix with diagonal entries $\si_1^2>0,\dots,\si_k^2>0,0,\dots,0$, so that $k$ is the rank of $\Si$. Let $Y:=Q^TX$, so that $X=QY$ and $Y=(Y_1,\dots,Y_k,0,\dots,0)\sim N(0,D)$. So, for any Borel-measurable function $g\colon\mathbb R^n\to\mathbb R$ such that $Eg(X)$ exists, we have \begin{equation} \begin{aligned} Eg(X)&=Eg(Q[Y_1,\dots,Y_k,0,\dots,0]^T) \\ &=\int_{\mathbb R^k}g(Q[y_1,\dots,y_k,0,\dots,0]^T) \\ &\times f_{\si_1}(y_1)\cdots f_{\si_k}(y_k)dy_1\cdots dy_k, \end{aligned} \tag{1} \end{equation} where $f_\si$ is the pdf of $N(0,\si^2)$.
Now you can compute the probabilities for $X$ as follows: $$P(X\in B)=E1(X\in B),$$ where $B$ is any Borel subset of $\mathbb R^n$.
For an illustration, here is a Mathematica notebook with the calculation of the probability $P(X_1>X_2+1)$ given $(X_1,X_2)\sim N\left((0,0),\left( \begin{array}{cc} 1 & 2 \\ 2 & 4 \\ \end{array} \right)\right)$:
As you insist on writing the expectation $Eg(X)$ as an integral with respect to the Lebesgue measure on the support of the distribution of $X$, here it is.
Note that the support of the distribution of $X$ is the range(=column space) $\Si\R^n$ of $\Si$, which is the same as the range $q\R^k$ of $q$, where $q$ is the map \begin{equation*} \R^k\ni z=(z_1,\dots,z_k)\mapsto qz:=QJz\in q\R^k=\Si\R^n \end{equation*} and, in turn, $J$ is the map \begin{equation*} \R^k\ni z=(z_1,\dots,z_k)\mapsto Jz:=(z_1,\dots,z_k,0,\dots,0)\in\R^n. \end{equation*} The equality $\Si\R^n=q\R^k$ follows because, in view of (0), the range $\Si\R^n$ of $\Si$ is the same as that of $QD^{1/2}$, which is the same as that of $q$ (since the range of $D^{1/2}$ is the same as that of $J$). Here, we identify matrices with the corresponding linear transformations.
In view of (1), \begin{equation*} \begin{aligned} Eg(X)&=\int_{\R^k}g(qz)\nu(dz) =\int_{\Si\R^n}g(x)\mu(dx), \end{aligned} \tag{2} \end{equation*} where $\nu(dz):=\pi(z)\la_k(dz)$, $\pi(z):=f_{\si_1}(z_1)\cdots f_{\si_k}(z_k)$, $\la_k(dz)$ is the Lebesgue measure on $\R^k$, and $\mu(dx):=\nu(q^{-1}(dx))=\pi(q^{-1}(x))\la_k(q^{-1}(dx))$.
Since the map $q$ is an isometry, for the Lebesgue measure $\la_{\Si\R^n}$ over $\Si\R^n$ we have $\la_{\Si\R^n}(dx)=\la_k(q^{-1}(dx))$. So, by (2), \begin{equation*} \begin{aligned} Eg(X)&=\int_{\Si\R^n}g(x)\pi(q^{-1}(x))\la_{\Si\R^n}(dx). \end{aligned} \tag{3} \end{equation*} In particular, letting $g:=1_B$ for any Borel set $B\subseteq\Si\R^n$, we get \begin{equation*} \begin{aligned} P(X\in B)&=\int_B \pi(q^{-1}(x))\la_{\Si\R^n}(dx). \end{aligned} \tag{3a} \end{equation*}
I think formulas (3) and, in particular, (3a) are completely useless, though -- because you will actually compute $Eg(X)$ (as in the above Mathematica notebook) by formula (1).