I'm using R and I wonder whether I should first remove highly correlated variables via
vif()and apply then lasso, forward selection etc.
Or doesn't it matter? Do the packages handle multicollinearity automatically?
At first, please remember, that multicollinearity is not a binary thing, therefore it can not be handled or not-handled as 0 and 1. Every method used to handle multicollinearity will handle it to some extent - better or worse - often risking removing of important variables, when the filter is stronger or just by bad luck. One may think about it as of a trade-off and every cure at some point might be worse than disease.
This is the reason of such articles as Performance of some variable selection methods when multicollinearity is present, which give more complex, but better-than-general answers. The main conclusion from their existence is that different methods may perform differently in different contexts.
There are are also many methods and their variants. You should look at every method separately. You mentioned two methods - LASSO and forward selection:
LASSO is a special case of Penalized Regresion methods. These methods were designed to handle multicollinearity. The question is if what they do is what you want. The basic property of LASSO is that from a bunch of correlated variables it tends to pick one of them and discard the others. You do not have to be satisfied with LASSO picks (see comments). You can control them to some extent with trade-off parameter $\lambda$, which describes the strength of penalty for the coefficient values. It controls how severe levels of multicollinearity LASSO will address, but do not control what variables will be picked, when multicollinearity is met.
For a quick informations about LASSO in R, I would suggest Glmnet Vignette by Hastie & Qian.
This method is part of group of methods called Stepwise Regression. They differ not only by step procedure (forward, backwards, all possibilities and others), but also by criterion - they use for example p-values, R$^2$, MSE, AIC, BIC. Then they will perform differently when challenged by multicollinearity.
To show the last very shortly - forward selection adding variables by R$^2$ is said to be good choice for cases with multicollinearity, while backward elimination removing by p-value will never drop two collinear variables in situation when they both are significant, especially when highly biased. As in the case of LASSO, you might not be satisfied with the chosen variables for many reasons.
PRACTICAL NOTE: Sometimes is much easier just to run some things and then look how they performed. This might be the case.
Multicollinearity cannot, logically, be "solved". Rather, it indicates that you cannot readily tease apart the effect of two or more predictors. Regularisation (lasso, ridge, elastic net) try to reduce the effect of collinarity on parameter estimation and thereby improve prediction.
Also, when the correlation structure changes (in time, in space, in resolution), your model will succumb to not having been able to identify the TRUE underlying correlation (if I may blushingly refer to our paper on the subject: https://doi.org/10.1111/j.1600-0587.2012.07348.x). No model type is immune to this.
Finally, the "inflation" of variance inflation factor has been criticised as suggesting something is wrong with the estimates. It isn't: their error bars are so large because high correlation among predictors prevents the model from identifying the TRUE relationship, and this uncertainty is reflected in the standard errors (https://doi.org/10.3998/ptpbio.16039257.0010.003).