Asymptotic for $y'' + \frac{\epsilon y'}{y^2} - y' = 0$, $y(-\infty) = 1$, $y(+\infty) = \epsilon$. 
Asymptotic for $y'' + \frac{\epsilon y'}{y^2} - y' = 0$, $y(-\infty) = 1$, $y(+\infty) = \epsilon$.

I started with a regular expansion for 
$$y^2y'' + \epsilon y' - y^2 y' = 0$$
and 
$$y = y_0 + \epsilon y_1 + O(\epsilon^2)$$ 
with $y_0(-\infty) = 1$, $y_0(+\infty) = 0$ and $y_1(-\infty) = 0$, $y_1 (+\infty) = 1$.
The zero order ODE is 
$$ y_0'' -  y_0' = 0, $$
this gives $y_0 = Ae^{x} + B$, but I can this can not satisfy the boundary condition. How should I approach differently? So this problem is a singular perturbation problem. We can not use WKB because this is not linear, and I am not familar with explicit boundary layer calculation on an unbounded domain.

Let us integrate the differential equation once to get rid of the second derivative: $y'-\epsilon/y-y = c$, where $c$ is a constant. To satisfy the boundary condition $y(-\infty) = 1$, $y'(-\infty)=0$, there is no other choice than $c=-(1+\epsilon)$. This Chini equation rewrites as $y' = (y-1)(y-\epsilon)/y$, so that $y$ is monotonously decreasing if $0<\epsilon<1$. The implicit analytical solution is defined up to a constant $k$:
$$
k + x = \frac{1}{2}\ln((y-\epsilon)(y-1)) - \frac{1+\epsilon}{1-\epsilon} \tanh^{-1}\left(\frac{2y - (1+\epsilon)}{1-\epsilon}\right)
.
$$
After some algebra,
we get
$$
k + x = \ln\left(-\left(\frac{1-y}{(y - \epsilon)^\epsilon}\right)^{1/(1-\epsilon)}\right) 
$$
for $\epsilon <y<1$.
To have $x$ real, the constant can be chosen as $k = \text{i}\pi$. For instance, in the case $\epsilon=1/2$, we get the solution $y(x) = \frac{1}{2} \left(e^x + 2 - \sqrt{e^x (e^x + 2)}\right)$,
which matches both boundary conditions:

This may give an insight about how to use perturbation methods in this case.
Let us introduce the method of matched asymptotic expansions. The idea is the following. As specified in OP, a perturbed solution $y_O = y_0 + \epsilon y_1 + \dots$ satisfies
\begin{aligned}
    y''_0-y'_0 &= 0, \\
    y''_1-y'_1 &= -y'_0/{y_0}^2, 
\end{aligned}
etc., so that
$$
y_O = a + be^x + \epsilon\left[c + de^x + \frac{a}{b^2} \left(\ln(a+be^x) - x\right) e^x\right] + \dots \, ,
$$
where the constants $a$, $b$, $c$, $d$ need to be determined. It is impossible that such an expansion matches both boundary conditions at $\pm\infty$. By rescaling the abscissas as $\xi = x\epsilon^\alpha$, the variable $Y(\xi) = y(x)$ satisfies $\epsilon^{2\alpha}Y''+\epsilon^{1+\alpha}Y'/Y^2 - \epsilon^\alpha Y' = 0$ where derivatives are computed w.r.t. $\xi$. The distinguished limit is $\epsilon\left( Y''+ Y'/Y^2\right) - Y' = 0$, corresponding to $\alpha=1$. A perturbed solution $Y_I = Y_0 + \epsilon Y_1 + \dots$ satisfies
\begin{aligned}
    Y'_0 &= 0,\\
    Y'_1 &= Y''_0-Y'_0/{Y_0}^2,
\end{aligned}
etc., so that $Y_I = A + \epsilon B + \dots,$ where the constants $A$, $B$ need to be determined. Thus, one can construct a piecewise asymptotic expansion of the solution, with an expansion $(Y_I)^-$ in the vicinity of $-\infty$, an expansion $Y_O$ in the vicinity of $x=0$ and an expansion $(Y_I)^+$ in the vicinity of $+\infty$. Then, the solutions are matched at the transition from one domain to another.
