About negative energies: they set no problem:
On this context, only energy differences have significance. Negative energy appears because when you've made the integration, you've set one point where you set your energy to 0. In this case, you have chosen that $PE_1 = 0$ for $r = \infty$. If you've set $PE_1 = 1000$ at $r = \infty$, the energy was positive for some r.
However, the minus sign is important, as it is telling you that the test particle is losing potential energy when moving to $r = 0$, this is true because it is accelerating, causing an increase in $KE$:
let's calculate the $\Delta PE_1$ for a particle moving in direction of $r = 0$: $r_i = 10$ and $r_f = 1$:
$\Delta PE_1 = PE_f - PE_i = Gm(-1 - (-0.1)) = -Gm\times0.9 < 0$
as expected: we lose $PE$ and win $KE$.
Second bullet: yes, you are right. However, it is only true IF they are point particles: has they normally have a definite radius, they collide when $r = r_1 + r_2$, causing an elastic or inelastic collision.
Third bullet: you are right with $PE_2 = mgh$, however, again, you are choosing a given referential: you are assuming $PE_2 = 0$ for $y = 0$, which, on the previous notation, means that you were setting $PE_1 = 0$ for $ r = r_{earth}$.
The most important difference now is that you are saying that an increase in h is moving farther in r (if you are higher, you are farther from the Earth center).
By making the analogy to the previous problem, imagine you want to obtain the $\Delta PE_2$. In this case, you begin at $h_i = 10$ and you want move to $h_f = 1$ (moving in direction to Earth center, like $\Delta PE_1$:
$\Delta PE_2 = PE_{f} - PE_{i} = 1mg - 10mg = -9mg < 0$.
As expected, because we are falling, we are losing $PE$ and winning $KE$, the same result has $PE_1$
Fourth bullet: they both represent the same thing. The difference is that $gh$ is the first term in the Taylor series of the expansion of $PE_1$ near $r = r_{Earth}$. As exercise, try to expand $PE_1(r)$ in a taylor series, and show that the linear term is:
$PE_1 = a + \frac{Gm(r-r_{earth})}{r_{earth}^2}$.
Them numerically calculate $Gm/r_{earth}^2$ (remember that $m=m_{earth}$). If you haven't made this already, I guess you will be surprised.
So, from what I understood, your logic is totally correct, apart from two key points:
energy is defined apart of a constant value.
in the $PE_1$, increase r means decrease $1/r$, which means increase $PE_2 = -Gm/r$. In $PE_2$, increase h means increase $PE_2=mgh$.
Best Answer
There are two problems with the manipulations you've done.
First, the variables in equation (2) are ambiguously named. Equation (2) calculates the potential energy between a single mass $m$ and a mass distribution with total mass $M$. Then the equation should actually read $$\Omega = - Gm \int \frac{dM}{R}.$$ If we instead write the differential as $dm$, it looks like $m$ is being integrated as well. This results in a meaningless extra factor of $1/2$ when the integration is performed.
Next, the integral over $dM$ shouldn't be naively performed as if $R$ is constant, $$\int \frac{dM}{R} \neq \frac{M}{R}$$ in general. The issue is that every piece of mass $dM$ has its own radius $R$, so $R$ should be thought of as a function of $M$. If this doesn't make sense, just think about the discrete case, $$\sum_i \frac{m_i}{R_i}$$ where a radius $R_i$ is associated with every bit of mass $m_i$.
In your particular case, where we're thinking about two point masses separated by a distance $R$, the quantity $R$ in the integrand really is constant, so we can pull it out for $$\Omega = - \frac{Gm}{R} \int dM = -\frac{GMm}{R}$$ as expected. For a more general configuration, we would parametrize the masses and radii somehow to get a concrete integral, e.g. we could use the chain rule for $$\int \frac{dM}{R} = \int \frac{dM}{dR} \frac{dR}{R}$$ where $dM/dR$ tells us the amount of mass in thin spherical shells of radius $R$. I explain how to do this kind of integral a bit more in this answer.