I wasn't clear exactly where you were stuck in the proof, so I posted a (more-or-less) full proof.
We assume one easily proved lemma:
lemma 1: For $|z| < 1/2$, $\left|e^z -1\right| \geq |z|/2$.
Next, we define a version of the Hankel contour. For small positive $\delta$ select two points in the complex plane:
\begin{align*}
p_1&=\delta e^{i\delta} = \delta_x + i\delta_y\quad\text{where}\quad \delta_x = \delta\cos(\delta),
\delta_y = \delta\sin(\delta)\\
p_2&=\delta e^{- i\delta} = \overline{\delta e^{i\delta}} = \delta_x - i\delta_y
\end{align*}
Our contour $C$ consists of three piecewise smooth curves $C_1+C_2+C_3$ defined as follows:
$C_1$ follows a path parallel to (and just above) the real axis from $\infty$ to $p_1$.
$C_2$ follows a counterclockwise circular path from $p_1$ to $p_2$, tracing the boundary of a circle that is centered at 0 with radius $\delta$. (Making just barely less than a full circle).
$C_3$ follows a path parallel to (and just below) the real axis from $p_2$ to $\infty$.
Lemma 2
Let $ z\in\mathbb{C} \backslash [0,\infty)$ and $Re(s) > 1$. Then
$\displaystyle \lim_{\delta \to 0^+}\int_{|z|=\delta}\frac{-z^{s-1}}{e^z -1}\, dz =0$.
Proof - Lemma 2
Fix $s$. Let $A = | Re(s) - 1|$ and let $B = Im(s)$. We use the ML Inequality to obtain an upper bound for our integral. We also assume that $\delta < 1/2$ and apply lemma 1:
\begin{equation*}
\int_{|z|=\delta}\frac{-z^{s-1}}{e^z -1}\, dz \leq 2\pi\delta \cdot \max_{|z| = \delta}
\left| \frac{-z^{s-1}}{e^z -1} \right| \leq
2\pi\delta \cdot \max_{|z| = \delta} \left| \frac{-z^{s-1}}{\delta /2} \right|
= 4\pi \cdot \max_{|z| = \delta} \left|-z^{s-1}\right|.
\end{equation*}
So it only remains to show that
\begin{equation*}
\lim_{\delta \to 0^+} \max_{|z| = \delta} \left|-z^{s-1}\right| = 0.
\end{equation*}
Using the principal value of the complex logarithm (for argument $\theta$: $-\pi < \theta \leq\pi$), we rewrite $-z = \delta e^{i\theta} = e^{\log(\delta) +i\theta}$, so that
\begin{align*}
\left|-z^{s-1}\right| = \left| e^{ \left( \log(\delta) +i\theta \right) [ A + i B] } \right|
= \left| e^{A\log(\delta) - B\theta} e^{ i \left( B\log(\delta) + A\theta \right) } \right|
= \left| e^{A \log(\delta) - B\theta } \right|
= \delta^A e^{- B\theta} \leq \delta^A e^{\pi \left| B \right| }.
\end{align*}
But $e^{\pi \left| B \right| }$ is fixed. Because $A > 0$, we have $\lim_{\delta \to 0^+} \delta^A = 0$, completing the proof.
End of Proof - Lemma 2
Main Theorem
Let $s \in \mathbb{C}$. We recall the traditional definition of $\zeta(s)$:
\begin{align*}
\zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^s} = \prod_{p}\frac{1}{ \left( 1-p^{-s} \right)}
\quad\text{(where $p$ ranges over an ordered list of all the primes)},
\end{align*}
convergent for $Re(s) > 1$. For the contour of integration $C$ (a Hankel contour oriented in the positive direction defined above), we have that:
\begin{align}
\zeta(s)=\frac{\Gamma(1-s)}{ -2 \pi i} \int_C \frac{-z^{s-1}}{e^z - 1} \, dz
\end{align}
extends $\zeta(s)$ to a meromorphic function on $\mathbb{C}$, holomorphic except for a simple pole at
$s = 1$.
Proof - Main Theorem.
Unless otherwise stated, we assume $Re(s) > 1$. We start with the Gamma function:
$$ \Gamma(s) = \int_0^\infty e^{-x} x^{(s-1)} \, dx. $$
Next, we replace $x$ by $nt$ in the integral (so that $dx=ndt$):
\begin{align*}
\Gamma(s) &= \int_0^\infty e^{-nt} (nt)^{(s-1)} n\, dt = \int_0^\infty e^{-nt}n^s t^{(s-1)} \, dt = n^s \int_0^\infty e^{-nt} t^{(s-1)} \, dt \\
n^{-s}\Gamma(s) &=\int_0^\infty e^{-nt} t^{(s-1)} \, dt \\
\frac{1}{n^s} &= \frac{1}{\Gamma(s)}\int_0^\infty e^{-nt} t^{(s-1)} \, dt.
\end{align*}
Now sum both sides over $n$:
\begin{align}
\sum_{n=1}^\infty \frac{1}{n^s} &= \frac{1}{\Gamma(s)}\sum_{n=1}^\infty \int_0^\infty e^{-nt} t^{(s-1)} \, dt \\
\zeta(s) &= \frac{1}{\Gamma(s)} \int_0^\infty \sum_{n=1}^\infty e^{-nt} t^{(s-1)} \, dt \\
\zeta(s) &= \frac{1}{\Gamma(s)} \int_0^\infty t^{(s-1)}\sum_{n=1}^\infty e^{-nt} \, dt.\tag{1}
\end{align}
Because we are assuming $Re(s) > 1$, the interchange of the sum and integral is justified by the uniform convergence of the sum. Notice that the sum over $n$ in the last integral is just a geometric series with the $n=0$ term missing. (We are justified in using a geometric series, because $0 < e^{-t} < 1$). So, we have:
$$
\sum_{n=1}^\infty e^{-nt} = \sum_{n=0}^\infty e^{-nt} - e^{-0t} = \sum_{n=0}^\infty (e^{-t})^n - 1= \frac{1}{1-e^t} -1 = \frac{e^t}{1-e^t}=\frac{1}{e^t-1}.
$$
But that means we can restate equation (1) (and while we are at it replace $dt$ with the slightly more standard $dx$), and have:
$$
\zeta(s) = \frac{1}{\Gamma(s)} \int_0^\infty x^{(s-1)}\sum_{n=1}^\infty e^{-nx} \, dx = \frac{1}{\Gamma(s)} \int_0^\infty \frac{x^{s-1}}{e^x - 1} \, dx.
$$
If we set $\displaystyle I(s)= \int_0^\infty \frac{x^{s-1}}{e^x - 1} \, dx$, then we have: $\displaystyle \zeta(s) = \frac{1}{\Gamma(s)} I(s) $.
Our next goal is to evaluate $I(s)$ in the complex plane and see where that takes us. Riemann starts with a related but not identical integral, which we call $I_0(s)$:
$$
I_0(s) = \int_C \frac{-z^{s-1}}{e^z - 1} \, dz,
$$
using the contour of integration $C$ described above.
This contour stays well away from the singularities of the integrand at $2\pi i \mathbb{Z}$. We have:
$$
I_0(s) = \int_C \frac{-z^{s-1}}{e^z - 1} \, dz = \int_\infty^{p_1} \frac{-z^{s-1}}{e^z - 1} \, dz + \int_{|z|=\delta} \frac{-z^{s-1}}{e^z - 1} \, dz + \int_{p_2}^\infty \frac{-z^{s-1}}{e^z - 1} \, dz.
$$
We now study $I_0(s)$ as $\delta\rightarrow 0$. By Lemma 2, the second integral can be disregarded so the contour includes only the horizontal lines $\infty\rightarrow p_1$ and $p_2\rightarrow\infty$.
Note that $(-z)^{(s-1)}=e^{(s-1)log(-z)}$ and that $log(-z)=log|z|+iarg(-z)$. We will use $Arg(z)$, the principal value of $arg(z)$, defined as the value $\theta$ satisfying $-\pi < \theta \leq \pi$. Recall that $z \mapsto Arg(z)$ is discontinuous at each point on the nonpositive real axis. Let $z=x_0+iy$ for some fixed $x_0 < 0$. If $y\downarrow 0$ then $Arg(z) \rightarrow \pi$, whereas, if $y \uparrow 0$ then $Arg(z) \rightarrow -\pi$.
In the first integral, $-z=-\delta_x - i\delta_y$ approaches the nonpositive real axis from below, so (as $\delta\rightarrow 0$) we have
$Arg(-z) \rightarrow -\pi$. In the third integral, $-z=-\delta_x + i\delta_y$ approaches the nonpositive real axis from above, so we have $Arg(-z) \rightarrow \pi$.
\begin{align*}
I_0(s) &= %\int_C \frac{-z^{s-1}}{e^z - 1} \, dz =
\int_{\infty}^{p_1} \frac{e^{(s-1)(\log |z| -\pi i)}}{e^z - 1} \, dz
+ \int_{p_2}^{\infty} \frac{e^{(s-1)(\log |z| +\pi i)}}{e^z - 1} \, dz \\
&= -\int_{p_1}^{\infty} \frac{ e^{(s-1) \log |z| } e^{(s-1)(-\pi i)} }{e^z - 1} \, dz
+ \int_{p_2}^{\infty} \frac{ e^{(s-1) \log |z| } e^{(s-1)(\pi i)} }{e^z - 1} \, dz \\
&= -e^{(s-1)(-\pi i)} \int_{p_1}^{\infty} \frac{ e^{(s-1) \log |z| } }{e^z - 1} \, dz
+ e^{(s-1)(\pi i)} \int_{p_2}^{\infty} \frac{ e^{(s-1) \log |z| } }{e^z - 1} \, dz.
\end{align*}
Taken to the limit, our integrals are over a horizontal line with $dy = 0$, so we can replace $dz$ with $dx$ (and $z$ with $x$), replace $p_1$ and $p_2$ with $0$,
and consider $e^{(s-1) \log |z|} = x^{(s-1)}$:
\begin{align*}
\lim_{\delta \to 0+} I_0(s)
&=-e^{(s-1)(-\pi i)} \int_{0}^{\infty} \frac{ e^{(s-1) \log |z| } }{e^z - 1} \, dz
+ e^{(s-1)(\pi i)} \int_{0}^{\infty} \frac{ e^{(s-1) \log |z| } }{e^z - 1} \, dz \\
&= \left[ e^{i\pi(s-1)} - e^{-i\pi(s-1)} \right] \int_0^\infty \frac{x^{(s-1)}}{e^x - 1} \, dx.
\end{align*}
But we have: $e^{iw}-e^{-iw} = \cos(w) + i \sin(w) - \left[ \cos(w)-i\sin(w) \right] = 2i\sin(w)$.
So, the value in brackets is: $2i\sin(\pi (s-1)) = 2i\sin(\pi s - \pi) = -2i\sin(\pi s)$.
This gives our value for $I_0(s)$:
\begin{equation*}
\lim_{\delta \to 0+} I_0(s) =
-2i\sin(\pi s)\int_0^\infty \frac{x^{(s-1)}}{e^x - 1} \, dx = -2i\sin(\pi s) I(s).
\end{equation*}
Reviewing, we have:
\begin{equation*}
\zeta(s)= \frac{1}{\Gamma(s)} I(s)\quad\text{and}\quad I_0(s) = -2i\sin(\pi s) I(s)
\quad\text{so that}\quad \zeta(s)=\frac{1}{\Gamma(s)}\cdot\frac{1}{ -2i\sin(\pi s) } I_0(s).
\end{equation*}
Now use a rearrangement of Euler's Reflection Formula: $\displaystyle \frac{1}{\Gamma(s)} = \frac{\Gamma(1-s)\sin(\pi s)}{\pi}$, giving
\begin{equation*}
\zeta(s)= \left[ \frac{1}{\Gamma(s)} \right] \frac{1}{ -2i\sin(\pi s) } I_0(s)
= \left[ \frac{\Gamma(1-s)\sin(\pi s)}{\pi} \right] \frac{1}{ -2i\sin(\pi s) } I_0(s)
= \frac{\Gamma(1-s)}{ -2 \pi i} I_0(s).
\end{equation*}
We therefore have:
\begin{align}
\zeta(s)=\frac{\Gamma(1-s)}{ -2 \pi i} \int_C \frac{-z^{s-1}}{e^z - 1} \, dz.
\end{align}
Up to this point, we have assumed $Re(s) > 1$. But the integral is an entire function (it is uniformly convergent in any compact subset of $\mathbb{C}$ because $e^z$ grows faster than any power of $z$). That means that $\zeta(s)$ is analytic, except possibly at the positive integers where $\Gamma(1-s)$ has simple poles. But we already know that $\zeta(s)$ is analytic for $Re(s) > 1$, so the possible poles at $2,3,4...$ must be removable and cancel against zeros of the integral. Thus, the only pole of $\zeta(s)$ is at $s= 1$.
End of Proof - Main Theorem.
Best Answer
Please allow me to elaborate a little on Riemann's and Jeb's elegant exposition.
So imagine the contour $\gamma$ as coming just "above" the real axis from $+\infty$ to $0$ and then swinging counterclockwise around the origin, and then going back just "below" the real axis back to $+\infty$.
As Riemann suggests, we can write $(-x)^{s-1} = e^{(s-1)\log(-x)}$. But what should $\log(-x)$ be equal to for positive $x$?
As Jeb suggested, $-1 = e^{i\pi} = e^{-i\pi}$, so we can write $\log(-x)$ as $\log(e^{i\pi}x) = i\pi+\log{x}$, or as $\log(e^{-i\pi}x) = -i\pi+\log{x}$. However, we must make choices that preserve the continuity of $\log(-x)$ as $x$ moves along the contour in the complex plane. Also, as Riemann specifies, we want $\log(-x)$ to be real when $x$ is negative.
For the piece of the contour from $+\infty$ to $0$, we must choose $$\log(-x) = -i\pi+\log{x}.$$
On the way back from $0$ to $+\infty$, in order to preserve continuity we must have $$\log(-x) = i\pi+\log{x}.$$
Thus, the integral for the part of the contour from $+\infty$ to $0$ is equal to $$\int_{+\infty}^0 \frac{e^{-(s-1)i\pi}e^{(s-1)\log{x}}}{e^x-1}\,dx = e^{-\pi si} \int_0^{+\infty} \frac{e^{(s-1)\log{x}}}{e^x-1}\,dx.$$ And the integral for the part of the contour from $0$ to $+\infty$ is equal to $$\int_0^{+\infty} \frac{e^{(s-1)i\pi}e^{(s-1)\log{x}}}{e^x-1}\,dx = -e^{\pi si} \int_0^{+\infty} \frac{e^{(s-1)\log{x}}}{e^x-1}\,dx.$$
Now add those two pieces together to get Riemann's result.