In this study, the Variational Bayes (VB) approach was hybridized with the bootstrap prior procedure to improve the accuracy of subset selection as well as optimizing the algorithm time in modelling high-dimensional genomic data with inherent sparse structure. The new hybrid VB approach is shown to yields a minimal sufficient statistic which under mild regularity conditions converges to the true sparse structure. Simulation and real-life high-dimensional genomic data experiments revealed comparable empirical performance with other competing frequentist and Bayesian methods. In addition, a new fast algorithm that illustrates the procedure was developed and implemented in the environment of R statistical software as package “VBbootprior”.