NB: This answer has two parts. The first part addresses the integral in the question and the second provides some simple lower bounds on $\Gamma(s,x)$.
Part I: Bounds on the integral in the question.
Here we provide some tight bounds on the integral
$$
\newcommand{\rd}{\,\mathrm d} I(a) := \int_a^\infty y^{-r} e^{-b(y-a)^2}\rd y
$$
and, in particular, give an asymptotic equivalence.
Claim: For $a > 0$, $b > 0$, and $r > 0$,
$$
\sqrt{\frac{\pi}{4b}}\, \Big(a+\frac{1}{\sqrt{b\pi}} \Big)^{-r} \leq I(a) \leq \sqrt{\frac{\pi}{4b}} \,a^{-r} \>.
$$
In particular, we have that $a^r I(a) \uparrow \sqrt{\pi/4b}$ as $a \to \infty$.
Proof. (Lower bound.) Use the substitution $z = \sqrt{2b}(y-a)$ to rewrite the integral in the "canonical form"
$$
\sqrt{\frac{\pi}{4 b}} \int_0^\infty (a+z/\sqrt{2b})^{-r} \frac{2 e^{-\frac{1}{2}z^2}}{\sqrt{2\pi}} \rd z = \sqrt{\frac{\pi}{4 b}} \mathbb E (a+Z/\sqrt{2b})^{-r}\>,
$$
where $Z$ is a standard half-normal random variable. The function $(a+z/c)^{-r}$ is convex in $z$ and so Jensen's inequality yields the stated lower bound by noting that $\mathbb E Z = \sqrt{2/\pi}$.
(Upper bound). Trivially $(a+z/\sqrt{2b})^{-r} \leq a^{-r}$ and so monotonicity of expectation yields the stated upper bound.
(Asymptotic result). That $I(a) \sim \sqrt{\pi/4b} \,a^{-r}$ as $a \to \infty$ follows directly from the upper and lower bounds. That $a^r I(a) \uparrow \sqrt{\pi/4b}$ follows from the fact that $(1+z/a)^{-r}$ is monotone increasing in $a$ and from monotonicity of expectation.
Part II: Lower bounds on $\Gamma(s,x)$.
Below, we assume real $s$ and $x > 0$. We use similar tricks to those in the previous part. Generally, the bounds we obtain below are asymptotically of the form $g(s) e^{-x} x^{s-1}$ as $x \to \infty$.
First of all, note that
$$
\Gamma(s,x) = \int_x^\infty t^{s-1} e^{-t} \rd t = e^{-x} \int_0^\infty (t+x)^{s-1} e^{-t} \rd t = e^{-x} \mathbb E(T+x)^{s-1} \>,
$$
where $T \sim \mathrm{Exp}(1)$.
Claim: For $s \leq 1$ and $s \geq 2$, we have $\Gamma(s,x) \geq e^{-x} (1+x)^{s-1}$.
Proof. $(t+x)^{s-1}$ is convex in $t$ for the corresponding $s$. Applying Jensen's inequality to the right-hand side of the above display equation and noting $\mathbb E T = 1$ yields the result.
Claim: For $s \geq 1$, we have $\Gamma(s,x) \geq e^{-x} \Gamma(s) (1+x/s)^{s-1}$. The inequality reverses for $0 < s < 1$.
Proof: The function $f(t) = (1+x/t)^{s-1}$ is convex for $s \geq 1$ and concave for $s < 1$. Apply Jensen's inequality to the right-hand side of
$$
\Gamma(s,x) = \int_x^\infty t^{s-1} e^{-t} \rd t = e^{-x} \int_0^\infty (1+x/t)^{s-1} t^{s-1} e^{-t} \rd t = e^{-x} \Gamma(s) \mathbb E(1+x/G)^{s-1} \>,
$$
where $G$ is a gamma random variable with shape parameter $s$ and scale parameter 1, and hence mean $\mathbb E G = s$. (For the reverse inequality when $s < 1$, the restriction to $s > 0$ is necessary to maintain the probabilistic interpretation and, hence, the application of Jensen's inequality.)
You can just use the general Leibniz rule
$$
(f_1\dots f_m)^{(n)}=\sum_{k_1+k_2+\dots+k_m=n}\binom{n}{k_1,k_2,\dots,k_n}\prod_{i=1}^m f_i^{(k_i)}
$$
and get
$$
\left.\frac{\mathrm{d}^n}{\mathrm{d}s^n}[\Gamma(s)^n]\right\rvert_{s=1} = \sum_{k_1+k_2+\dots+k_n=n}\binom{n}{k_1,k_2,\dots,k_n}\prod_{i=1}^n \left.\frac{\mathrm{d}^{k_i}}{\mathrm{d}s^{k_i}}\Gamma(s)\right\rvert_{s=1}.
$$
There is an error in your derivation. Since $\log t<0$ for $t\in(0,1)$, and you bound $0<e^{-t}<1$, your bound for $K_n$ doesn't work with the $(-1)^n$. If you put everything in modulus sign you get $\lvert K_n\rvert\leq n!$ which seems to be what you were doing at the end.
Best Answer
Equation 6.5.1 here gives a reasonably tight bound except for very small n. I suppose section 5.11 here also could be unpacked to yield some upper bounds, depending on exactly what you're looking for. (I assume you know that for integer n, $\Gamma(n) = (n-1)!$.)