First I presume that
If I know the posterior distribution, and from the Gibbs sampling, I
got the sampled parameter
simply means that you can compute $\pi(\theta|x)$ up to a constant and derive a Markov chain $\theta¹,\theta²,\ldots,\theta^T$ that converges to $\pi(\theta|x)$.
About your questions:
- How to draw the histogram with y axis as marginal density in R ?
is unclear. You can draw the histogram of your simulations, one component at a time, but you do not know the marginal density $\pi_1(\theta_1|x)$ in most cases. You can always produce and display an estimate of this marginal density, e.g., by Rao-Blackwellisation:$$\hat{\pi}_1(\theta_1|x)=\frac{1}{T}\sum_{t=1}^T \pi(\theta_1|\theta_{-1}^t,x)$$where the conditional density on the right represents the full conditional of $\theta_1$ given all other components that you used in Gibbs sampling. But this display cannot be used as a proof of convergence as the histogram _and_ the alternative density estimate are based on the same simulation.
In Example 7.2 of our book Introducing Monte Carlo methods with R, we
consider the pair of distributions $$ X\vert\theta \sim
\mathcal{B}\text{in}(n,\theta)\,,\quad \theta \sim {\mathcal{B}}e(a,b)
$$ (i.e., the joint distribution) as our target (this is a toy example where there is no observation, no parameter, no posterior, despite the possibly confusing notations of $X$ and $\theta$: both are random variables to be simulated there) and this leads to the joint distribution $$
f(x,\theta) = {n \choose x} \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\,
\theta^{x+a-1} (1-\theta)^{n-x+b-1} $$ that we can simulate by Gibbs
sampling since the corresponding conditional distribution of $X \vert
\theta$ is given above, while $\theta \vert x \sim {\cal
B}e(x+a,n-x+b)$. The histograms in the associated figure (Fig. 7.1)
are constructed by Gibbs sampling, while the overlaid curves are the
true marginal densities, because this is a toy problem where the exact marginals are known to be a Beta-Binomial for $X$ and a Beta for $\theta$. This hardly ever happens in practical situations.

[Reproduced from Introducing Monte Carlo methods with R]
And
- How to further apply Monte Carlo integration on the resulting Markov Chain, which are sampled from this posterior. It seems I
couldn't get a connection between these two.
is also unclear. If you mean producing an approximation of a posterior expectation $\mathbb{E}[\mathfrak{h}(\theta)|x]$, you simply apply the regular Monte Carlo formula:
$$\mathbb{E}[\mathfrak{h}(\theta)|x]\approx\frac{1}{T}\sum_{t=1}^T \mathfrak{h}(\theta^t)$$