As suggested by Philip Cherian, you may want to check in the Math section. However, I want to point out one thing about how to handle the limit.
I will assume that $|\vec{k}|\ne 0$. The first thing you want to do it separate each propagator into its real and imaginary parts. In general, $$\lim_{\delta \rightarrow 0^+}\frac{1}{E-E(q)+i\delta} = \lim_{\delta \rightarrow 0^+}\left[\frac{1}{E-E(q)}-i\frac{\delta}{[E-E(q)]^2 + \delta^2}\right],$$
where I multiplied top and bottom by the complex conjugate of the denominator, and disregarded the $\delta^2$ in the denominator of the real part, since it'll just go away when you take the limit. Then you use the identity $$\lim_{\delta\rightarrow 0^+}\frac{\delta}{[E-E(q)]^2+\delta^2} = \pi\delta(E-E(q)),$$
that is, the Dirac delta. What you get is $$\begin{split}&\int d^3q\,\left[\frac{1}{E-E(q)}-i\pi\delta(E-E(q)) \right]\left[\frac{1}{E-E(|\vec{q}+\vec{k}|)}-i\pi\delta(E-E(\vec{q}+\vec{k})) \right]\\ &=\int d^3q\,\frac{1}{E-E(q)}\frac{1}{E-E(|\vec{q}+\vec{k}|)} + i\pi\int d^3q\,\left[ \frac{\delta(E-E(q))}{E-E(|\vec{q}+\vec{k}|)} + \frac{\delta(E-E(|\vec{q}+\vec{k}|))}{E-E(q)} \right]. \end{split}$$
Notice that, since I'm assuming $k\ne 0$, the product of two Dirac deltas will have to vanish, because their arguments are never simultaneously zero. You can easily evaluate the integral with the Dirac deltas. Recall that, since the dispersion is quadratic in $q$, you need to follow certain rules when you change variables. That gives you the imaginary part of the integral. I have no real advice on how to evaluate the real part $$\mathrm{Re}\{I(\vec{k})\}=\int d^3q\,\frac{1}{E-E(q)}\frac{1}{E-E(|\vec{q}+\vec{k}|)},$$ except to check that you actually need it. Depending on the problem you're solving, sometimes all you need is the imaginary part. On the other hand, I also recommend you check that both denominators in $I(\vec{k})$ have $+i \epsilon$, instead of one having $+i \epsilon$ and the other $-i\epsilon$. I say this because that expression looks like the skeleton approximation to a polarization bubble, which is the product of a retarded and an advanced propagator, in which case their imaginary parts have opposite signs.