Calculate $\zeta(s) = \frac{1}{\Gamma(s)} \sum_{n=0}^\infty \frac{B_n}{n!} \int_0^\infty x^{s + n - 2} dx$ By Taylor expanding 
$$\frac{x}{e^{x}-1} = \sum_{n=0}^\infty \frac{B_n}{n!}x^n$$
in the Zeta function
$$\zeta(s) = \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} (\frac{x}{e^{x}-1})dx$$
we find 
\begin{align}
\zeta(s) &= \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} (\frac{x}{e^{x}-1})dx = \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} \sum_{n=0}^\infty \frac{B_n}{n!}x^n dx \\
&= \frac{1}{\Gamma(s)} \sum_{n=0}^\infty \frac{B_n}{n!}  \int_0^\infty x^{s + n - 2} dx
\end{align}
what exactly do we do next to arrive at the formulas for $\zeta(-n)$ and $\zeta(2n)$?
 A: By Taylor expanding, for $|x| < 2 \pi$ :
$$\frac{x}{e^{x}-1} = \sum_{n=0}^\infty \frac{B_n}{n!}x^n$$
in the Zeta function, for $Re(s) > 1$ :
$$\zeta(s)\Gamma(s) =  \int_0^\infty  \frac{x^{s-1}}{e^{x}-1}dx =  \int_0^a x^{s-2} \frac{x}{e^{x}-1}dx +  \int_a^\infty  \frac{x^{s-1}}{e^{x}-1}dx$$
we find whenever $0 < a < 2 \pi$ and $Re(s)> 1$, inverting $\sum$ and $\int$ by absolute/monotone/dominated convergence :
\begin{align}
\zeta(s)\Gamma(s) &=  \int_a^\infty \frac{x^{s-1}}{e^{x}-1}dx + \int_0^a x^{s-2} \sum_{n=0}^\infty \frac{B_n}{n!}x^n dx \\
&= \int_a^\infty  \frac{x^{s-1}}{e^{x}-1}dx  + \sum_{n=0}^\infty\frac{B_n}{n!}  \int_0^a x^{s + n - 2} dx\\
&=  \int_a^\infty  \frac{x^{s-1}}{e^{x}-1}dx  +  \sum_{n=0}^\infty\frac{B_n}{n!} \frac{a^{s+n-1}}{s+n-1}\end{align}
note how by analytic continuation it stays valid for every $s \in \mathbb{C}$ except  at the poles, 
and how it tells us the residue of $\zeta(s) \Gamma(s)$ at its poles, and hence the value of $\zeta(-k)$  for $k \in \mathbb{N}$
