It seems that the integral can be expressed using McRobert's E-function, which is a generalized hypergeometric function. Let the function defined as
\begin{equation}
I(X)=\int_0^X e^{\frac{1}{x}+\frac{1}{X-x}} x^{b-a-1}(X-x)^{a-1}\,dx \end{equation}
its value at $X=1$ is to be found. It can be considered as a convolution, its Laplace transform is thus
\begin{equation}
\mathcal{L}\left[I;X\to p\right]=F_{b-a-1}(p)F_{a-1}(p)
\end{equation}
where
\begin{equation}
F_\alpha(p)=\int_0^\infty x^\alpha e^{-px-\tfrac{1}{x}}\,dx
\end{equation}
From Ederlyi TI, p.146 (4.5.29),
\begin{equation}
F_\alpha(p)=2p^{-\frac{\alpha+1}{2}}K_{\alpha+1}\left( 2\sqrt{p} \right)
\end{equation}
$K_\alpha(.)$ is a modified Bessel function. Then
\begin{equation}
\mathcal{L}\left[I;X\to p\right]=4p^{-\frac{b}{2}}K_{b-a}\left( 2\sqrt{p} \right)K_{a}\left( 2\sqrt{p} \right)
\end{equation}
Now, using the inverse Laplace transform of the product of two modified Bessel functions derived in this paper,
\begin{equation}
I(1)=\frac{2^{b-2a}}{\sqrt{\pi}}\sum_{i,-i}\frac{1}{i}\left[E\left( \left.
\begin{array}{l}
b+1,a+1,b-a+1,1\\
\tfrac{b}{2}+1,\tfrac{b+3}{2}
\end{array}
\right|4e^{i\pi} \right)\right]
\end{equation}
where the sum is taken above and below the branch cut $\Re z<0$ of the McRobert E-function $E\left( \left.
\begin{array}{l}
u_p\\
v_q
\end{array}
\right|z \right)$. This expression may probably be simplified, at least for some specific values of $a$ and $b$ and expressed in terms of classical generalized hypergeometric functions or in terms of Meijer G-function.
Remark 1: in the Ederlyi table (61) p.285, a particular case can be obtained:
$$ \mathcal{L}^{-1}\left[\sqrt{p}K_{\nu+\tfrac{1}{2}}(2\sqrt{p})K_{\nu-\tfrac{1}{2}}(2\sqrt{p})\right]=\frac{1}{4}\sqrt{\frac{\pi}{2}}e^{-2}W_{\tfrac{1}{2},\nu}(4)$$
where $W_{\mu,\nu}(.)$ is the Whittaker function. Choosing $\nu=-1/2,b=-1,a=-1$, it comes
$$I(-1,-1)=\int_0^1 e^{\frac{1}{x}+\frac{1}{1-x}} x^{-1}(1-x)^{-2}\,dx=\sqrt{\frac{\pi}{2}}e^{-2}W_{\tfrac{1}{2},-\tfrac{1}{2}}(4) $$
However, numerically this result is exact if the r.h.s. is multiplied by a factor $\sqrt 2$. Typo in the table or mistake in my derivation?
Remark 2: Changing $x\to 1-x$ in the integral gives the identity
$$I(a,b)=I(b-a,b)$$
Moreover, by the relation $(1-x)^a=(1-x)^{a-1}-x(1-x)^{a-1}$, we have the property
$$I(a+1,b)=I(a,b-1)-I(a,b)$$
This expression can be generalized, by the use of the binomial expansion of $(1-x)^{a+k}$. With these recurrence relation and the particular cases of @Nemo and from this answer, we can access to many special cases of the integral.