Gaussian process models are generally fine with high dimensional datasets (I have used them with microarray data etc). They key is in choosing good values for the hyper-parameters (which effectively control the complexity of the model in a similar manner that regularisation does).
Sparse methods and pseudo-input methods are more for datasets with a large number of samples (> approx 4000 for my computer) rather than a large number of features. If you have a powerful enough computer to perform a Cholesky decomposition of the covariance matrix (n by n where n is the number of samples), then you probably don't need these methods.
If you are a MATLAB user, then I'd strongly recommend the GPML toolbox and the book by Rasmussen and Williams as good places to start.
HOWEVER, if you are interested in feature selection, then I would avoid GPs. The standard approach to feature selection with GPs is to use an Automatic Relevance Determination kernel (e.g. covSEard in GPML), and then achieve feature selection by tuning the kernel parameters to maximise the marginal likelihood. Unfortunately that is very likely to end up over-fitting the marginal likelihood and ending up with a model that performs (possibly much) worse than a model with a simple spherical radial basis function (covSEiso in GPML) covariance.
My current research focus lies on over-fitting in model selection at the moment and I have found that this is as much a problem for evidence maximisation in GPs as it is for cross-validation based optimisation of hyper-paraneters in kernel models, for details see this paper, and this one.
Feature selection for non-linear models is very tricky. Often you get better performance by sticking to a linear model and using L1 regularisation type approaches (Lasso/LARS/Elastic net etc.) to achieve sparsity or random forest methods.