After doing integration by parts, he ended up with
$$
\operatorname{Ei}(z)=\frac{e^z}{z} \biggl [\sum_{k=0}^n \frac{k!}{z^k}+\epsilon _n(z) \biggr ].
$$
[...] He then showed that $$|\epsilon _n(z)| \leq \frac{(n+1)!}{|z|^{n+1} (\sin \delta)^{n+2}},$$
where $\delta$ is some constant such that it ranges $(0,\pi)$, and such that $|\arg(-z)| \leq \delta - \pi$.
[...] I think I should point out that the $\sin \delta$ above is not in absolute values, which I think is problematic, since provided the range of $\delta $, $\sin \delta $ could take on negative values and the inequality above will not be true anymore.
If $\delta \in (0,\pi)$ then $\sin \delta > 0$, so there's no need to worry about the absolute value.
Does the last result above show that $\epsilon _n (z) = O(z^{-(n+1)})$?
The number $\delta$ is constant with respect to $z$ here, so there really is a constant
$$
C = \frac{(n+1)!}{(\sin \delta)^{n+2}}
$$
such that
$$
|e_n(z)| \leq C |z|^{-(n+1)}.
$$
This is exactly what it means to write $e_n(z) = O(z^{-(n+1)})$.
Now, owing to the fact that there's still that $\frac{e^z}{z}$ factor outside the RHS of the second equation from above, then from what I understood, it should say that:
$$
Ei(z) = \frac{e^z}{z} \sum_{k=0}^n \frac{k!}{z^k} + O \left (\frac{e^z}{z^{(n+2)}} \right),
$$
which doesn't really fit the Poincare asymptotic expansion. Or is because I shouldn't see it in light of the Poincare asymptotic expansion but in terms of some other asymptotic scale?
There are two ways to think about it. One way is to agree that it shows that
$$
\sum_{k=0}^n \frac{k!}{z^k}+\epsilon _n(z) = \sum_{k=0}^n \frac{k!}{z^k}+O(z^{-(n+1)})
$$
is a Poincaré expansion, so what the formula
$$
\operatorname{Ei}(x) \sim \frac{e^x}{x} \sum_{k=0}^{\infty} \frac{k!}{z^k}
$$
means is precisely that
$$
\operatorname{Ei}(x) = \frac{e^x}{x} \left[\sum_{k=0}^{n} \frac{k!}{z^k}+O(z^{-(n+1)})\right]
$$
for every $n \in \mathbb N$.
The other way is to expand things as you've done and write
$$
\operatorname{Ei}(z) = \sum_{k=0}^n \frac{k!e^z}{z^{k+1}} + O\left( \frac{e^z}{z^{(n+2)}} \right). \tag{$*$}
$$
This is now a generalized asymptotic expansion $\sum_k g_k(z)$, where the summands
$$
g_k(z) = \frac{k!e^z}{z^{k+1}}
$$
form an asymptotic sequence. Indeed we have
$$
O\left( \frac{e^z}{z^{(n+2)}} \right) \subset o(g_n(z))
$$
as $z \to \infty$, so a consequence of $(*)$ is that
$$
\operatorname{Ei}(z) = \sum_{k=0}^n g_k(z) + o(g_n(z))
$$
for every fixed $n \in \mathbb N$ as $z \to \infty$. This is precisely the statement that
$$
\operatorname{Ei}(z) \sim \sum_{k=0}^\infty g_k(z)
$$
as $x \to \infty$, that is
$$
\operatorname{Ei}(z) \sim \sum_{k=0}^\infty \frac{k!e^z}{z^{k+1}} \tag{$**$}
$$
as $x \to \infty$. Now, the interpretation would be equivalent either way (i.e. the Big O statements one could obtain from the asymptotic series) so one could formally factor out $e^z/z$ from this "sum" and thus write the result as
$$
\operatorname{Ei}(z) \sim \frac{e^z}{z}\sum_{k=0}^\infty \frac{k!}{z^k} \tag{$***$}
$$
as $z \to \infty$.
You should convince yourself that $(**)$ and $(***)$ really are equivalent statements.
Following a more elementary route we will first break the integral in two pieces:
$$B(\frac{N+1}{N+2};N+1,p+1)=B(N+1,p+1)-I\\I=\int^{1}_{\frac{N+1}{N+2}}x^{N}(1-x)^pdx$$
Perform the following change of variables $x=1-\frac{t}{N+2}$:
$$I=\frac{1}{(N+2)^{p+1}}\int_{0}^1t^p(1-\frac{t}{N+2})^{N}dt$$ However it is possible to estimate that
$$(1-\frac{t}{N+2})^{N}=e^{-t}(1+\frac{4t-t^2}{2N}+\mathcal{O}(N^{-2}))$$
and performing the integrals we get:
$$I\sim\frac{1}{N^{p+1}}\Big[\gamma(p+1,1)+\frac{2}{N}(\gamma(p+2,1)-p-1-\frac{1}{4}\gamma(p+3,1))\Big]$$
While on the other hand we may write
$$B(N+1,p+1)\sim \frac{\Gamma(p+1)}{N^{p+1}}\Big[1-\frac{(p+1)^2+3(p+1)}{2N}\Big]$$
So all in all as a leading approximation one gets:
$$B(\frac{N+1}{N+2};N+1,p+1)\sim\frac{\Gamma(p+1)-\gamma(p+1,1)}{N^{p+1}}=\frac{\Gamma(p+1,1)}{N^{p+1}}+\mathcal{O}(N^{-p-2})$$
where $\Gamma(z,x)=\int_x^{\infty}t^{z-1}e^{-t}dt$ is the upper incomplete Gamma function. From the above one can calculate also the subleading term. More subleading terms can also be readily calculated but probably not in closed form.
Best Answer
In general no. For example, $e^{-\nu^2/x^2} = O(x^n)$ as $x \to 0$ for all $n \geq 0$ and $\nu \neq 0$ fixed, so we have
$$ \frac{1}{1-x} + e^{-\nu^2/x^2} \sim \sum_{n=0}^{\infty} x^n$$ as $x \to 0$ for all fixed $\nu \neq 0$. This asymptotic is not uniform with respect to $\nu$ when, say, $\nu = x$.