numerical integration - Multidimensional NIntegral with singularity

I am dealing with an integral

$F(k)=-\int d^2p d^2q\frac1{\lvert p\rvert}\frac1{\lvert p-q\rvert}(\frac1{\lvert k-q\rvert}-\frac1{\lvert q\rvert})$

when I am evaluating it at some point, say k=30

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]](1/
Sqrt[30^2 + q^2 - 60 q Cos[\[Phi]q]]-1/q), {p, 0, 50}, {q, 0, 
50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0, 2 \[Pi]}]

MMA gives 1397.63 with warning

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one 
of the following: singularity, value of the integration is 0, highly 
oscillatory integrand, or WorkingPrecision too small.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has 

increased more than 2000 times. The global error is expected to decrease 
monotonically after a number of integrand evaluations. Suspect one of the 
following: the working precision is insufficient for the specified precision 
goal; the integrand is highly oscillatory or it is not a (piecewise) smooth 
function; or the true value of the integral is 0. Increasing the value of the 
GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical 
integration. NIntegrate obtained -3120.06+0.187845 I and 1266.8675450228204` 
for the integral and error estimates.

The integtal has singularity but should be converget. I tried to increase Maxrecursion Minrecursion or Working precision does not help.

Besides, if I split the integral in to 2 parts from bracket of the formula and integrate separately

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]]/
Sqrt[30^2 + q^2 - 60 q Cos[\[Phi]q]], {p, 0, 50}, {q, 0, 
50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0, 2 \[Pi]}]

answer is -3120.06 + 0.187845 I

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]]/
q, {p, 0, 50}, {q, 0, 50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0, 
2 \[Pi]}]

answer is -4638.95 + 1.66533*10^-16 I

with imaginary parts

I saw some other post with similar question, some answer suggest using

    Method -> "MonteCarlo", PrecisionGoal -> 1

It does work for me, gives result 1220.21 with no warning I don't know how much can I trust this result Because the answer changes slightly together with warning when I increase Precisiongoal

NIntegrate::maxp: The integral failed to converge after 50100 integrand 
evaluations. NIntegrate obtained -3233.79 and 35.28825355386109` for the 
integral and error estimates.

I need trustful answer for other part of the program, any suggestions?