For a function $u\in W^{1,p}$, we always use a sequence of $C^1$ function to approach it and derive the consequences. However, can we just claim for a.e. x, y:
$$u(x)-u(y)=\int_0^1 Du(y+t(x-y))\cdot (x-y) dt$$
I don’t think it’s right but I just can not find a counter example.
Yes this is true, see for example Theorem 6.9 in the book Analysis by Lieb and Loss. It tells more precisely that if $u\in W^{1,1}_{\rm loc}$, then for every $z\in\Bbb R^d$ and a.e. $x\in\Bbb R^d$, $$ u(x+z) - u(x)= \int_0^1 z\cdot\nabla u(x+t\,z)\,\mathrm dt. $$ The theorem also tells a condition to have the equality in the sense of distributions for when $u$ is not regular (which is what you usually really want to use as the a.e. is usually harder to manipulate).