I don't know anything about this.

One reason why the $n$th derangement number is very close to $n!/e$ is that the exponential generating function for the derangement numbers is $$\frac{e^{-z}}{1-z} = \frac{e^{-1}}{1-z} + g(z)$$ where $g(z)$ is an entire function, and therefore $g(z)$ has coefficients that are very small compared to the coefficients of $e^{-z}/(1-z)$. By the same reasoning, if $f(z)$ is an entire function, or even a function with radius of convergence greater than 1, then the coefficients of $f(z)/(1-z)$ will be very close to $f(1)$.

Yes, I mean integrating $x^{\alpha -1}(1-x)^{-A}$ times each term of the hypergeometric series.

This is what you get by integrating termwise, using the beta integral.

It seems likely that the formula might appear in Pólya's 1937 paper (on "Pólya's theorem") but I haven't checked this. It's also always a good idea to check Richard Stanley's Enumerative Combinatorics for anything on enumerative combinatorics.

It's not really an answer, just a suggestion.

See the OEIS, oeis.org/A000246, in particular the reference to Riordan's book.

You can use John Stembridge's SF package for Maple: dept.math.lsa.umich.edu/~jrs/maple.html