Add Vignette with Multiple Groups and Specifying Null Genes

[dcgerard]

A frequent question concerns setting up multiple groups and specifying the null genes. Adding a vignette would help out a lot.

Here is the script I usually send folks, and this could easily be turned into a vignette if it's cleaned up a little.

library(seqgendiff)
set.seed(1)

## Generate some simulated data. In practice you would use real data.
n <- 100
p <- 1000
Y <- rpois(n = n*p, lambda = 100)
dim(Y) <- c(p, n)

## Subsample individuals so that simulation does not depend on quirks
## of some genes.
Ysub <- select_counts(mat = Y, nsamp = 50, ngene = 1000)

## Generate design and coefficient matrices
group <- sample(c(1:4), size = ncol(Ysub), replace = TRUE)
group <- as.factor(group)
X <- model.matrix(~group)
X <- X[, -1, drop = FALSE] ## Remove intercept
betamat <- rnorm(ncol(X) * nrow(Ysub))
dim(betamat) <- c(nrow(Ysub), ncol(X))

## Choose which genes are null by setting those rows in betamat to 0
nullvec <- sample(x = c(TRUE, FALSE),
                  size = nrow(betamat),
                  replace = TRUE,
                  prob = c(0.9, 0.1))
betamat[nullvec, ] <- 0

## X is the design matrix representing group.
X

## betamat is the coefficient matrix containing effect sizes for
## each group at each gene
head(betamat)

## Thin
thout <- thin_diff(mat = Ysub, design_perm = X, coef_perm = betamat)
Ynew <- thout$mat
Xnew <- thout$designmat
betanew <- thout$coefmat

## Compare linear regression estimates with true coefficients
coefest <- t(coef(lm(log2(t(Ynew) + 0.5) ~ Xnew))[-1, , drop = FALSE])
plot(coefest, betanew)
abline(0, 1, lty = 2, col = 2)