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$.
The formula:
$$ \Delta U = mgh $$
is an approximation that applies when the distance $h$ is small enough that changes in $g$ can be ignored. As you say, the expression for $U$ is:
$$ U= -\frac{G M m}{r} $$
So the change when moving a distance $h$ upwards is:
$$ \Delta U = \frac{GMm}{r} - \frac{GMm}{r + h} $$
We rearrange this to get:
$$\begin{align}
\Delta U &= GMm \left( \frac{1}{r} - \frac{1}{r + h} \right) \\
&= GMm \frac{h}{r^2 + rh} \\
&= \frac{GM}{r^2} m \frac{h}{1 + h/r} \\
&\approx \frac{GM}{r^2} m h
\end{align}$$
where the last approximation is because $h \ll r$ so $1 + h/r \approx 1$. And since $GM/r^2$ is just the gravitational acceleration $g$ at a distance $r$, we get:
$$ \Delta U = g m h $$
Best Answer
Your equation (2) is the change in potential energy when the object moves vertically by a distance $h$ i.e. when the object moves from $r$ to $r+h$. Let's use equation (1) to calculate this:
$$ \Delta U = GMm\left(\frac{1}{r}-\frac{1}{r+h}\right) $$
Subtracting the two fractions inside the bracket gives:
$$\begin{align} \Delta U &= GMm\left(\frac{r+h}{r(r+h)}\frac{r}{r(r+h)}\right) \\ &= GMm\frac{h}{r(r+h)} \end{align}$$
Since $h \ll r$ that means $r+h\approx r$ and our equation becomes:
$$\begin{align} \Delta U &\approx GMm\frac{h}{r^2} \\ &\approx \frac{GM}{r^2}mh \\ &\approx gmh \end{align}$$
Footnote:
I've just noticed that in your comment to ItachÃ's answer you ask if you can use a Taylor series. You can use a binomial expansion to make the aproximation more obvious. You rewrite:
$$ \Delta U = GMm\frac{h}{r(r+h)} $$
as:
$$ \Delta U = \frac{GM}{r^2}mh\left(1+\frac{h}{r}\right)^{-1} $$
then a binomial expansion gives:
$$ \Delta U = \frac{GM}{r^2}mh\left(1-\frac{h}{r} + O\left(\frac{h}{r}\right)^2 \right) $$
And as before since $h \ll r$ the term in the brackets is approximately one and we once again get:
$$ \Delta U = \frac{GM}{r^2}mh $$