diff --git a/R/ZlmFit-bootstrap.R b/R/ZlmFit-bootstrap.R
index e3a0972..c617dc4 100644
--- a/R/ZlmFit-bootstrap.R
+++ b/R/ZlmFit-bootstrap.R
@@ -32,6 +32,7 @@ pbootVcov1<-function (cl,zlmfit, R = 99)
 ##' Sample cells with replacement to find bootstrapped distribution of coefficients
 ##' @param zlmfit class \code{ZlmFit}
 ##' @param R number of bootstrap replicates
+##' @param boot_index list of indices to resample.  Only one of R or boot_index can be offered.
 ##' @return array of bootstrapped coefficients
 ##' @importFrom plyr raply
 ##' @examples
@@ -40,12 +41,20 @@ pbootVcov1<-function (cl,zlmfit, R = 99)
 ##' #Only run 3 boot straps, which you wouldn't ever want to do in practice...
 ##' bootVcov1(zlmVbeta, R=3)
 ##' @export
-bootVcov1 <- function(zlmfit, R=99){
+bootVcov1 <- function(zlmfit, R=99, boot_index = NULL){
     sca <- zlmfit@sca
     N <- ncol(sca)
     LMlike <- zlmfit@LMlike
-    manyvc <- raply(R, {
-        s <- sample(N, replace=TRUE)
+    if(is.numeric(R)){
+        boot_index = lapply(1:R, function(i) sample(N, replace = TRUE))
+    } else if(!is.null(boot_index)){
+        r = range(unlist(boot_index))
+        if(!is.list(boot_index) || r[1] < 1 || r[2] > N) stop("boot_index must be a list of integer vectors between 1 and ", N)
+    } else{
+        stop("Provide only of of `boot_index` or `R`.")
+    }
+    
+    manyvc <- laply(boot_index, function(s){
         newsca <- sca[,s]
         LMlike <- update(LMlike, design=colData(newsca))
         zlm(sca=newsca, LMlike=LMlike, onlyCoef=TRUE)
diff --git a/man/bootVcov1.Rd b/man/bootVcov1.Rd
index d1b2024..a84e04e 100644
--- a/man/bootVcov1.Rd
+++ b/man/bootVcov1.Rd
@@ -7,7 +7,7 @@
 \usage{
 pbootVcov1(cl, zlmfit, R = 99)
 
-bootVcov1(zlmfit, R = 99)
+bootVcov1(zlmfit, R = 99, boot_index = NULL)
 }
 \arguments{
 \item{cl}{a \code{cluster} object created by \code{makeCluster}}
@@ -15,6 +15,8 @@ bootVcov1(zlmfit, R = 99)
 \item{zlmfit}{class \code{ZlmFit}}
 
 \item{R}{number of bootstrap replicates}
+
+\item{boot_index}{list of indices to resample.  Only one of R or boot_index can be offered.}
 }
 \value{
 array of bootstrapped coefficients
diff --git a/tests/testthat/test-GSEA-by-boot.R b/tests/testthat/test-GSEA-by-boot.R
index e657f27..a28d389 100644
--- a/tests/testthat/test-GSEA-by-boot.R
+++ b/tests/testthat/test-GSEA-by-boot.R
@@ -5,7 +5,7 @@ vb1 = subset(vbetaFA[1:24,], ncells==1)
 #vb1 = vb1[,freq(vb1)>.1]
 zf = zlm(~Stim.Condition, vb1)
 set.seed(1234)
-boots = bootVcov1(zf, 36)
+boots = bootVcov1(zf, R = 36)
 ## replace NAs for each coefficient, gene and component
 bootsUncor <- apply(boots, 2:4, function(col){
     col[!is.na(col)] <- col[!is.na(col)]-mean(col[!is.na(col)], na.rm=TRUE)
diff --git a/tests/testthat/test-bootstrap.R b/tests/testthat/test-bootstrap.R
index 83a0ea4..ffcb650 100644
--- a/tests/testthat/test-bootstrap.R
+++ b/tests/testthat/test-bootstrap.R
@@ -103,9 +103,23 @@ test_that('Bootstrap recovers covariance', {
 
 parallel::stopCluster(cl)
 
-## M <- melt(boot[,,'groupB','C'])
-## ggplot(M, aes(x=value))+geom_density() + facet_wrap(~X2)
-
+# context("Nearly singular designs")
+# 
+# test_that('Bootstrap results are padded appropriately', {
+#     N <- 12
+#     m <- 20
+#     p <- 3
+#     X <- getX(p, N/p, N)
+#     beta <- rbind(2, matrix(0, nrow = p-1, ncol = m))
+#     Y <- simYs(m, X, beta, rho=0, sigma=1, p=.7)
+#     cData <- data.frame(group = attr(X, 'group'))
+#     cData$group = factor(cData$group)
+#     sca <- suppressMessages(suppressWarnings(FromMatrix(t(Y$Y), cData=cData)))
+#     zfit <- suppressWarnings(zlm(~group, sca=sca))
+#     # Only fit on groupA/groupB samples
+#     boot <- bootVcov1(zfit,R = NULL,boot_index  = list(c(rep(1, 4), 1:8)))
+#     expect_equal(colnames(boot), colnames(coef(zfit, 'D')))
+# })