I have the following setup: $p(y_i|z_i) = N(y_i|w_y^Tz_i + u_y, \sigma_y^2), p(x_i|z_i) = N(x_i|W_xz_i+u_x, \sigma_x^2I), p(z_i) = N(0,I)$.
Now $p(y_i|x_i) = \frac{p(y_i,x_i)}{p(x_i)} = \frac{\int p(y_i, x_i|z_i)p(z_i)dz_i}{p(x_i)}$. Also, I have that $x_i,y_i$ are conditionally independent given $z_i$ and have a joint Gaussian distribution. So, $p(x_i,y_i|z_i) = N(u,\Sigma), u = \begin{pmatrix}W_xz_i + u_x \\ w_y^Tz_i+u_y\end{pmatrix} = \begin{pmatrix}W_x \\ w_y^T\end{pmatrix}z_i + \begin{pmatrix}u_x \\ u_y\end{pmatrix}, \Sigma = \begin{pmatrix}\sigma_x^2I & 0\\ 0 & \sigma_y^2\end{pmatrix}$.
Now, using Linear Guassian system we can write $p(x_i,y_i) = N(\begin{pmatrix}u_x\\u_y\end{pmatrix}, \Sigma + \begin{pmatrix}W_x^TW_x & w_yW_x\\W_x^Tw_y^T & w_yw_y^T\end{pmatrix})$.
However, I don't see how can one compute $\frac{p(x_i,y_i)}{p(x_i)}$?
In multivariate Gaussian distribution, there are known formulas to obtain the conditional from the jointly distribution. For example Section 8.1.3 of https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf
In your case, I think that your jointly distribution has a mistake because $Y$ takes values in $\mathbb{R}$. Then, $p(x,y)=N\left(\left(\begin{array}{c}u_x\\u_y\end{array}\right),\Sigma+\left(\begin{array}{cc}W_xW_x^T&W_xw_y\\w_y^TW_x^T&w_y^Tw_y\end{array}\right)\right)$. Anyway, that can be checked by looking at the dimensions.
Then, the conditional can be written as $p(y|x)=N(y|\mu(x),\sigma^2)$ with $$\mu(x)=u_y+w_y^TW_x^T(\sigma_x^2I+W_xW_x^T)^{-1}(x-u_x)$$ $$\sigma^2=\sigma_y^2+w_y^Tw_y-w_y^TW_x^T(\sigma_x^2I+W_xW_x^T)^{-1}W_xw_y$$