Simple answer
No, you cannot. Since in your problem setting, $\Sigma | \vec{\mathbf{y}}$ and $\theta | \vec{\mathbf{y}}$ are by no means independent. Thus we cannot simply marginalize it.
Your try
Here the problem setting is follows:
$$
\begin{aligned}
\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt]
\Sigma &\sim \mathcal{W}^{-1} \left( \nu_0, \mathbf{S}_0 \right)
\end{aligned}
$$
where $\mathcal{W}^{-1}(\nu, \mathcal{\Psi})$ denotes the inverse-Wishart distribution, which is a generalization of inverse gamma$(\alpha, \beta)$ distribution when dimension $p=1$, with $\nu = 2\alpha, \mathcal{\Psi} = 2\beta$.
Note that we have assumed a given (nonrandom) $\theta$. This is why we can write
$$
\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma \sim \mathcal{N}\left( \theta, \Sigma \right)
$$
in the first place, rather than $\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta \sim \mathcal{N}\left( \theta, \Sigma \right)$.
Posterior when both $\theta, \Sigma$ : unknown
Alternatively, let us suppose $\theta$ is random, so that we can assume a prior on $\theta$. This is exactly where your idea may kick in. Note that when we start assuming
$$
\begin{aligned}
\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt]
\Sigma &\sim \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0)
\end{aligned}
$$
we also simply get the posterior of $\Sigma | \vec{\mathbf{y}}, \theta$ , since by the Bayes rule,
$$
\begin{aligned}
p(\Sigma | \vec{\mathbf{y}}, \theta) &= \frac{p(\vec{\mathbf{y}}|\Sigma, \theta) p(\Sigma | \theta)}{\int_{\tilde{\Sigma}} p(\vec{\mathbf{y}}|\tilde{\Sigma}, \theta) p(\tilde{\Sigma}| \theta)} \\[8pt]
&\propto \mathcal{N} \left( \nu_0 + n, [\mathbf{S} + \mathbf{S}_\theta]^{-1} \right)
\end{aligned}
$$
as in the equation (7.9). Note that unlike the above, we cannot say $\Sigma|\theta \equiv \Sigma$. This is exactly why your idea might need caveats. Rather, if we allow an additional information to be
$$
\theta | \Sigma \sim \mathcal{N} (\theta_0, \Sigma/m)
$$
for a nonrandom hyperparameters $\theta_0$,$m$, things get simpler, i.e. we can get a closed-form solution.
Posterior when $\theta | \Sigma$ : normal
The joint density of $\theta$, $\Sigma$ is commonly referred to as Normal-inverse-Wishart distribution, which is simply formulated as
$$
\theta, \Sigma \sim \mathcal{N} (\theta_0, \Sigma/m) \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \equiv NIW(\theta, m, \mathbf{S}_0, \nu_0)
$$
It can be shown that the posterior distribution is also the Normal-inverse-Wishart, as
$$
\begin{aligned}
\theta | \Sigma , \vec{\mathbf{y}} &\sim \mathcal{N} (\theta_n, \Sigma/m_n) \\[8pt]
\Sigma | \vec{\mathbf{y}} &\sim \mathcal{W}^{-1} (\nu_n, \mathbf{S}_n)
\end{aligned}
$$
where
$$
\begin{aligned}
\theta_n &= \frac{m \theta_0 n + \bar{\mathbf{y}} }{ m_n } \\
m_n &= m + n\\
\nu_n &= \nu_0 + n \\
\mathbf{S}_n &= \mathbf{S}_0 + \mathcal{S} + \frac{mn}{m+n} (\bar{\mathbf{y}} - \theta_0) (\bar{\mathbf{y}} - \theta_0)^T
\end{aligned}
$$
where $\mathcal{S} := (\vec{\mathbf{y}}-\bar{\mathbf{y}}\mathbb{1})(\vec{\mathbf{y}}-\bar{\mathbf{y}}\mathbb{1})^T$.
On your comment
So you are now assuming the model
$$
\begin{aligned}
\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt]
\Sigma &\sim \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \\[10pt]
\theta &\propto 1 \\[10pt]
\end{aligned}
$$
where definitely the improper prior $\theta$ on $p$-dimensional real space $\mathbb{R}^p$ assumed. Most Bayesians accept improper priors if the resulting posterior is proper(see here).
So let us see if this setting leads to the proper posterior. We have
$$
\begin{aligned}
p(\theta, \Sigma) &\propto \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \prod_{i=1}^n \mathcal{N}_\mathbf{y}\left( \theta, \Sigma \right) \\[10pt]
\end{aligned}
$$
thus we have
$$
\Sigma , \theta | \vec{\mathbf{y}} \sim \mathcal{W}^{-1} (\nu_0 + n, [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1})
$$
Note that the difference between here and in (7.9) is the LHS in (7.9) is $\Sigma | \theta, \vec{\mathbf{y}}$, whereas we have joint here. Thus, we can integrate out $\theta$, (kind of) as what you said,
$$
\begin{aligned}
p(\Sigma | \vec{\mathbf{y}}) &= \int_\Theta \mathcal{W}^{-1} (\nu_0 + n, [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1}) d\theta \\[10pt]
&= \int_{\mathbb{R}^p} \frac{\left|\mathbf{S}_0 + \mathbf{S}_\theta \right|^{-(\nu_0 + n)/2}}{2^{(\nu_0 + n) p/2}\Gamma_p(\frac{\nu_0 + n}{2})} \left| \Sigma \right|^{-(\nu_0 + n+p+1)/2} e^{-\frac{1}{2}\operatorname{tr}( [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1} \Sigma^{-1})} d\theta \\[10pt]
&= \frac{\left| \Sigma \right|^{-(\nu_0 + n+p+1)/2}}{2^{(\nu_0 + n) p/2}\Gamma_p(\frac{\nu_0 + n}{2})} \int_{\mathbb{R}^p} \left|\mathbf{S}_0 + \mathbf{S}_\theta \right|^{-(\nu_0 + n)/2} e^{-\frac{1}{2}\operatorname{tr}( [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1} \Sigma^{-1})} d\theta
\end{aligned}
$$
This integral does not converge, which means it is an improper posterior. Thus, most Bayesians may not accept this as the posterior.
To see why the integral does not converge in a little bit heuristic way, say for the one dimensional case($p=1$), the integral is the form
$$
p(\sigma^2 | \vec{\mathbf{y}}) = C \cdot \int_\mathbf{R} \theta^{-(\nu_0 + n)} exp(-\theta^2) d\theta
$$
for some constant $C$. It is known that the above integral converges iff
$$
\nu_0 + n < 1
$$
which is not attainable for $n \ge 1$ cases.
PS. Interestingly, we can show that the marginal prior density of the $\theta$ is multivariate $t$ distribution (reference: Bernardo and Smith[1]), where
$$
\theta \sim T(\theta_0, \kappa, q)
$$
where $\kappa=\nu_0 \mathbf{S}_0/m, q=2\nu_0 - p + 1$.
Reference
[1] Bernardo and Smith, Bayesian theory(1994), p435
