Signs of the lasso solution Let us consider the following lasso estimator:
$$
\hat{\beta}_{L} = \arg\min \, \frac{1}{n}\sum_{i}^{n}||y_{i} - \textbf{x}_{i}\beta||_{2}^{2} + \frac{\lambda_{n}}{n}\sum_{j=1}^{p}|\beta_{j}|
$$
Assume, that $\beta \in \mathbb{R}^{p}$ and $rank(\mathbf{X}) = p$.
Then, the MLE solution $\hat{\beta}_{ML}$ is unique and assume that it is  such that all the components of the vector $\hat{\beta}_{ML}$ are non-negative. 
Does the same hold for LASSO solution (the sign will be the same or they are zero)?
EDIT: It seems that for general $\textbf{X}$ it is possible. 
Would it still be true for orthogonal design? 
Would it be true in the case when $p=n$ and $\textbf{X} = \textbf{I}$?
PS. It seems so from the geometry of the problem in 2d and 3d cases. Though I can not to show it for a general dimension.
 A: There are two questions here.  Let's look at the general one first: Is it possible for the Lasso to reverse the signs of estimated coefficients?
The answer is yes.  It can occur when two variables are strongly negatively correlated with one another but positively correlated with a third.  Here is an example of such a design matrix $X:$
$$X = \pmatrix{0.4 & 0.6 & -0.8 \\ -1.0 & 0.9 & -0.9 \\ 1.0 & -1.0 & 0.9 \\ -0.6 & -0.1 & 0.9}$$
Taking the response to be 
$$y = (0.1, -1.0, 1.0, 0.3)^\prime$$
the OLS estimates for the model (with intercept, which is going to be part of the Lasso) are approximately 
$$\hat\beta = (0.04, 0.89, 0.80, 0.97)^\prime,$$
all of which are positive. 
Using the glmnet package for R in its default (Lasso) mode produces sequences of estimates for a large range of possible $\lambda.$  The Lasso trace plots these sequences as line graphs.  Given a design matrix x and response vector y, here is the code to produce the Lasso trace:
fit <- glmnet(x,y)
plot(fit)

The horizontal axis reflects the size of the penalty ("L1 Norm") while the vertical axis represents the estimated coefficients $\hat\beta(\lambda).$  This plot is best read right to left: at the right the penalty is large because $\lambda$ is small (but nonzero, so the coefficients aren't quite the same as the OLS estimates).  Scanning to the left, the values of $\lambda$ increase, putting more weight on the penalty and pushing the estimates to zero.
Typically, once an estimate hits zero it stays there.  But this one is different:

The estimated coefficients of the columns of $X$ are (from left to right in the matrix) shown in gold, light blue, and gray. The light blue estimate initially hits zero (where the $L^1$ Norm is between $0.8$ and $1.1$), but then the gray estimate  drops while the light blue estimate becomes negative.
In effect, the Lasso is trading the second and third variables off for one another.  This can be anticipated by examining the correlation matrix of $X:$
$$\operatorname{Cor}(X) = \pmatrix{1 & -0.68 & 0.34 \\ -0.68 & 1 & -0.89 \\ 0.34 & -0.89 & 1},$$
showing how strongly anticorrelated the last two variables are.

The second question asks whether the same phenomenon can occur when $X$ is orthogonal.  The answer is no.
In this case, varying one parameter estimate in OLS does not affect the other parameter estimates. Thus, if you suppose there is a case where the Lasso solution involves a negative parameter estimate, simply switch its sign and observe (1) the penalty remains the same; (2) all other OLS estimates remain the same; but (3) because the original OLS solution had a positive coefficient, the OLS term decreases. This shows the supposed solution did not minimize the objective. Thus a negative Lasso parameter estimate is impossible in this circumstance. 
This analysis supports the claim in the first part of the question that the sign reversal must be due to the correlation patterns among the explanatory variables.
