LASSO in the data with large number of variables (p) with lower number of samples (n) I would like to fit LASSO in the following type of data where there are large number of variables (p > n). My y variable is y and would like to fit rest of variables in myd as x variables. 
#random population of 200 subjects with 1000 variables 
M <- matrix(rep(0,200*100),200,1000)
for (i in 1:200) {
set.seed(i)
  M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) <- 1:200

#random yvars 
set.seed(1234)
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5 
set.seed(234)
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

myd <- data.frame(y=y, M)

My Question: 
 What is best method for this type of data and is there any R program that implements this? 
 A: I am extending the suggestion by Stephan Kolassa above by providing codes and some explanations. 
First your data: 
#random population of 200 subjects with 1000 variables 
M <- matrix(rep(0,200*100),200,1000)
for (i in 1:200) {
set.seed(i)
  M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) <- 1:200

#random yvars 
set.seed(1234)
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5 
set.seed(234)
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

myd <- data.frame(y=y, M)

You need to install glmnet and load it:
install.packages("glmnet", repos = "http://cran.us.r-project.org")
require(glmnet)

Glmnet is a package that fits a generalized linear model via penalized maximum likelihood. 
The regularization path is computed for the lasso or elasticnet penalty at a grid of values for the regularization parameter lambda.
Please note that the generalized linear model (GLM) is a flexible generalization of ordinary linear regression that allows for response variables that have error distribution models other than a normal distribution (see wikipedia).
The GLM generalizes linear regression by allowing the linear model to be related to the response variable via a link function and by allowing the magnitude of the variance of each measurement to be a function of its predicted value. If the link function in identity generalized linear model is same as ordinary linear model. 
Please read the nice online Vignette and manual to do your own analysis. You would use Gaussian Family as you have normal distribution expectation of y. 
The function glmnet Fit a generalized linear model via penalized maximum likelihood. The regularization path is computed for the lasso or elasticnet penalty at a grid of values for the regularization parameter $\lambda$.The function can deal with all shapes of data, including very large sparse data matrices. Fits linear, logistic and multinomial, poisson, and Cox regression models from the manual. 
The argument alpha is for the elastic-net mixing parameter α, with range α∈[0,1]. α=1 is the lasso (default) and α=0 is the ridge.
The family should be "gaussian" in your case and which is default. Thus if you like to do lasso with much of arguments as default, the following would be your code:
fit1=glmnet(M,y) # or  is same as
fit1=glmnet(M,y, family="gaussian", alpha=1)

The ridge regression would be:
fit2 =glmnet(M,y, family="gaussian", alpha=0)

You can do further steps including viewing outputs and predictions: 
 print(fit1)

This displays the call that produced the object fit and a three-column matrix with columns Df (the number of nonzero coefficients), %dev (the percent deviance explained) and Lambda (the corresponding value of λ).
  coef(fit1,s=0.01) # extract coefficients at a single value of lambda
 # Let’s plot “fit1” against the log-lambda value and with each curve labeled.
plot(fit1, xvar = "lambda", label = TRUE)


 # make predictions:
  newx <-  M[1:10,]    
 predict(fit1,newx=newx,s=c(0.01,0.005)) # make predictions

Here s is  Value(s) of the penalty parameter lambda at which predictions are required. 
Default is the entire sequence used to create the model.
The another thing you would like to do is cross validation. The function  Does k-fold cross-validation for glmnet, produces a plot, and returns a value for lambda. In addition to all the glmnet parameters, cv.glmnet has its special parameters including nfolds (the number of folds), folded (user-supplied folds), type.measure(the loss used for cross-validation): “deviance” or “mse” uses squared loss while “mae” uses mean absolute error.
     cvob1 <- cv.glmnet(M, y,  nfolds=5) # 5 fold cross validation, default is 10 
     plot(cvob1)


coef(cvob1)
predict(cvob1,newx=newx, s="lambda.min")

Here s is Value(s) of the penalty parameter lambda at which predictions are required. Default is the value s="lambda.1se" stored on the CV object. Alternatively s="lambda.min" can be used. If s is numeric, it is taken as the value(s) of lambda to be used.
The following another example that uses mae with 20 fold cross validation. 
set.seed(1011)
cvob1a=cv.glmnet(M, y, nfolds=20, type.measure="mae")
plot(cvob1a)

