Evaluating the log gamma integral $\int_{0}^{z} \log \Gamma (x) \, \mathrm dx$ in terms of the Hurwitz zeta function
To evaluate that limit, we can expand each function in a Laurent series at $s=0$ and use the following 3 facts about the Hurwitz zeta function:
$$ \zeta(-s,a) = \zeta(-s,a+1) + a^{s} \tag{1}$$
$$ \zeta'(-s,a) = \zeta'(-s,a+1) -a^{s} \log(a) $$
$$\zeta(-n, a) = -\frac{B_{n+1}(a)}{n+1} \ , \ n \in\mathbb{N} \tag{2}$$
Doing so, we get
$$z - z \log z - \frac{\gamma z^{2}}{2} + \lim_{s \to 0^{+}} \Big[ - \Gamma(s-1) \zeta(s-1,z+1) -z \Gamma(s) \zeta(s) + \frac{z^{2}}{2} \Gamma(s+1) \zeta(s+1)$$
$$+ \Gamma(s-1) \zeta(s-1) \Big]$$
$$ = z - z \log z - \frac{\gamma z^{2}}{2}$$
$$ + \lim_{s \to 0^{+}} \Bigg[-\Big(-\frac{1}{s} + \gamma -1 + \mathcal{O}(s) \Big) \Big( -\frac{z^{2}}{2}+\frac{z}{2}-\frac{1}{12}-z + \zeta'(-1,z)s + z \log z \ s + \mathcal{O}(s^{2}) \Big)$$
$$ - z \Big( \frac{1}{s} - \gamma + \mathcal{O}(s) \Big) \Big( - \frac{1}{2} - \frac{\log (2 \pi)}{2} s + \mathcal{O}(s^{2}) \Big) + \frac{z^{2}}{2} \Big(1- \gamma s + \mathcal{O}(s^{2}) \Big) \Big( \frac{1}{s} + \gamma + \mathcal{O} (s) \Big) $$
$$ + \Big(- \frac{1}{s} + \gamma -1 + \mathcal{O} (s) \Big) \Big( - \frac{1}{12} + \zeta'(-1) s + \mathcal{O}(s^{2}) \Big) \Bigg] $$
$$ = z - z \log z - \frac{\gamma z^{2}}{2}$$
$$ + \lim_{s \to 0^{+}} \Big[\zeta'(-1,z) + z \log z + \frac{\gamma z^{2}}{2} - \frac{z^{2}}{2} + \frac{z}{2} - z + \frac{z \log(2 \pi)}{2} - \zeta(-1)+ \mathcal{O}(s) \Big] $$
$$ = \frac{z}{2} \log(2 \pi) + \frac{z(1-z)}{2} -\zeta'(-1)+ \zeta'(-1,z) $$
$ $
$(1)$ http://dlmf.nist.gov/25.11 (25.11.3)
$(2)$ http://mathworld.wolfram.com/HurwitzZetaFunction.html (9)
Another approach
Start by the following
$$\zeta(s,z) = \frac{z^{-s}}{2}+\frac{z^{1-s}}{s-1}+2\int^\infty_0 \frac{\sin(s\arctan(x/z))}{(z^2+x^2)^{s/2}(e^{2\pi x}-1)}\,dx$$
Hence we have
$$\zeta'(-1,z) = -\frac{z\log(z)}{2}+\frac{z^2\log(z)}{2}-\frac{z^2}{4}+\int^\infty_0 \frac{x\log(x^2+z^2)+2z\arctan(x/z)}{(e^{2\pi x}-1)}\,dx$$
Now use that
$$\psi(z) = \log(z)-\frac{1}{2z}-2\int^\infty_0 \frac{x }{(z^2+x^2)(e^{2\pi x}-1)}dx$$
Which implies that
$$\int^\infty_0 \frac{2zx }{(z^2+x^2)(e^{2\pi x}-1)}dx=z\log(z)-\frac{1}{2}-z\psi(z)$$
By taking the integral
$$\int^\infty_0 \frac{x\log(x^2+z^2) -x\log(x^2)}{(e^{2\pi x}-1)}dx=\int^z_0x\log(x)\,dx-\int^z_0x\psi(x)\,dx-\frac{z}{2}$$
Which simplifies to
$$\int^\infty_0 \frac{x\log(x^2+z^2) }{(e^{2\pi x}-1)}dx=\zeta'(-1)-\frac{z^2}{4}+\frac{1}{2} z^2 \log(z)-z\log\Gamma(z)+\int^z_0\log\Gamma(x)\,dx-\frac{z}{2}$$
Also we have
$$2\int^\infty_0 \frac{x}{(x^2+z^2)(e^{2\pi x}-1)}dx=\log(z)-\frac{1}{2z}-\psi(z)$$
By integration we have
$$2\int^\infty_0 \frac{\arctan(x/z)}{(e^{2\pi x}-1)}dx=z+\frac{\log(z)}{2}-z \log(z)+\log\Gamma(z)+C$$
Let $z \to 1$ to evaluate the constant
$$2\int^\infty_0 \frac{\arctan(x/z)}{(e^{2\pi x}-1)}dx=z+\frac{\log(z)}{2}-z \log(z)+\log\Gamma(z)-\frac{1}{2}\log(2\pi)$$
Multiply by $z$
$$2\int^\infty_0 \frac{z\arctan(x/z)}{(e^{2\pi x}-1)}dx=z^2+\frac{z\log(z)}{2}-z^2 \log(z)+z\log\Gamma(z)-\frac{z}{2}\log(2\pi)$$
Substitute both integrals our formula
$$\zeta'(-1,z) =-\frac{z\log(z)}{2}+\frac{z^2\log(z)}{2}-\frac{z^2}{4}+z^2+\frac{z\log(z)}{2}-z^2 \log(z)+z\log\Gamma(z)\\-\frac{z}{2}\log(2\pi)+\zeta'(-1)-\frac{z^2}{4}+\frac{1}{2} z^2 \log(z)-z\log\Gamma(z)+\int^z_0\log\Gamma(x)\,dx-\frac{z}{2}$$
Which reduces to
$$\int_{0}^{z} \log \Gamma(x) \, \mathrm dx = \frac{z}{2} \log(2 \pi) + \frac{z(1-z)}{2} - \zeta^{'}(-1) + \zeta^{'}(-1,z)$$
Here's a derivation which relies on series operations and uses a Feynman trick.
The reation to be proved is this
$$ \int_{0}^{z} \log \Gamma(x) \, \mathrm dx = \frac{z}{2} \log(2 \pi) + \frac{z(1-z)}{2} - \zeta^{'}(-1) + \zeta^{'}(-1,z)\tag{1}$$
The l.h.s. can also be written as $ \psi^{(-2)}(z)$.
Now the one hand we have by the definition of the polygamma function the following sequence of integrations
$$\psi^{(0)}(z) = \frac{\partial}{\partial z} \log \Gamma(z)\tag{2a}$$ $$\psi^{(-1)}(z) = \int_{1}^{z} \psi^{(0)}(x) \, \mathrm dx = \log \Gamma(z)\tag{2b}$$ $$\psi^{(-2)}(z) = \int_{0}^{z} \psi^{(-1)}(x) \, \mathrm dx\tag{2c}$$
The key idea now is to use the series representation
$$\psi^{(0)}(z) = H_{z-1} -\gamma = \sum_{k=1}^{\infty} \left(\frac{1}{k} - \frac{1}{k+z-1}\right) - \gamma\tag{3a}$$
where we have used the well-known series for the harmonic number.
This relation $(3a)$ allows us to use the series as a kind of "looking glass" the see what happens with the functions on integrating.
Before starting we perform a kind of Feynman trick writing
$$\psi^{(0)}(z) =\lim_{s\to 1} \psi^{(0)}(s,z) $$
where the function extended to include the parameter $s$ naturally brings zeta and Hurwitz zeta into the game (This came as a Bonus, for in the first place I only wanted to avoid divergent sums):
$$\psi^{(0)}(s,z) = \sum_{k=1}^{\infty} \left(\frac{1}{k^s} - \frac{1}{(k+z-1)^s}\right) - \gamma = \zeta(s) - \zeta(s,z) -\gamma\tag{3b}$$
Remember that
$$\zeta(s,z) = \sum _{k=1}^{\infty } \frac{1}{(k+z-1)^s}\tag{4a}$$
from which we have also
$$\int \zeta (s,z) \, \mathrm dz = \frac{\zeta (s-1,z)}{1-s}\tag{4b}$$
Now we perform the integrations of $(2)$ term-wise under the sums (without Feynman's trick this would lead to divergent sums) giving
$$\psi^{(-1)}(s,z) = \int_{x=1}^{z}\psi^{(0)}(s,x) \, \mathrm dx= \frac{\zeta (s-1,z)-\zeta (s-1)}{s-1}+z \zeta (s)-\zeta (s)-\gamma z+\gamma\tag{5b}$$
$$\psi^{(-2)}(s,z)= \int_{x=0}^{z}\psi^{(-1)}(s,x) \, \mathrm dx \\= -\frac{2 \zeta (s-2,z)+(s-2) z ((s-1) (z-2) (\gamma -\zeta (s))+2 \zeta (s-1))-2 \zeta (s-2)}{2 (s-2) (s-1)} \tag{5c}$$
Now we finalize the Feynman trick performing the limit $\psi^{(-2)}(z) =\lim_{s\to 1} \psi^{(-2)}(z,s)$ using the series expansions
$$\zeta(s \simeq -1,z) \simeq \zeta (-1,z)+(s+1) \zeta'(-1,z)\tag{6a}$$ $$\zeta(s \simeq 1) \simeq -(s-1) \gamma _1+\frac{1}{s-1}+\Gamma\tag{6b}$$ $$\zeta(s \simeq 0) \simeq -\frac{1}{2} s \log (2 \pi )-\frac{1}{2}\tag{6c}$$ $$\zeta(s \simeq -1) \simeq -\frac{1}{12}+(s+1) \zeta'(-1,0)\tag{6d}$$
we finally arrive at the relation $(1)$ to be shown.
Note that sometimes $(6d)$ is written with $\zeta'(-1,0) \to \frac{1}{12}-\log (A)$ where $A \simeq 1.28243$ is Glaisher's constant.
Remark: Actually, the justification of Feynman's trick (which is in effect an interchange of the order of limits) should be done strictly. Here we took it as an heuristic tool which is justified in hindsight by the known correctness of the result.