Hello, Bioconductor, I'm working lately with data from small-RNA seq in an organism without small RNA annotation. I would like to perform DE analysis on Differentially Expresed Regions and I would like to ask about the creation of the model matrix. Following the example of the vignette :

Build models

mod <- model.matrix(~ pheno$group + pheno$gender)
mod0 <- model.matrix(~ pheno$gender)

I only have pheno$group, does that mean I just follow the "main" instructions from the limma vignette?

getF <- function(fit, fit0, theData) {
    rss1 <- rowSums((fitted(fit) - theData)^2)
    df1 <- ncol(fit$coef)
    rss0 <- rowSums((fitted(fit0) - theData)^2)
    df0 <- ncol(fit0$coef)
    fstat <- ((rss0 - rss1) / (df1 - df0)) / (rss1 / (ncol(theData) - df1))
    f_pval <- pf(fstat, df1 - df0, ncol(theData) - df1, lower.tail = FALSE)
    fout <- cbind(fstat, df1 - 1, ncol(theData) - df1, f_pval)
    colnames(fout)[2:3] <- c("df1", "df0")
    fout <- data.frame(fout)
    return(fout)
}

If I want to perform the single base-level F-statistics, instead of using the:

models <- makeModels(sampleDepths,
    testvars = pheno$group,
    adjustvars = pheno[, c("gender")]
)

Do I put the testvars only?

Thank you in advance!

Konstantinos

Hi,

Can you provide an example that breaks? The following R code works for me:

library("derfinder")
set.seed(20210801)
pheno <- data.frame(
    group = sample(c("A", "B"), 10, replace = TRUE),
    gender = sample(c("M", "F"), 10, replace = TRUE)
)

sampleDepths <- runif(10, 1e7, 1e8)

models <- makeModels(sampleDepths,
    testvars = pheno$group,
    adjustvars = pheno[, c("gender")]
)

as well as

models_noadj <- makeModels(sampleDepths,
    testvars = pheno$group,
    adjustvars = NULL
)

which I think is the situation you are describing of not adjusting for any covariates (adjustment variables). You can see how the models change with http://bioconductor.org/packages/ExploreModelMatrix/. Here are how the models look:

> sapply(models, colnames)
$mod
[1] "(Intercept)"  "testvarsB"    "sampleDepths" "adjustVar1M"

$mod0
[1] "(Intercept)"  "sampleDepths" "adjustVar1M"

> sapply(models_noadj, colnames)
$mod
[1] "(Intercept)"  "testvarsB"    "sampleDepths"

$mod0
[1] "(Intercept)"  "sampleDepths"

Package info:

> packageVersion("derfinder")
[1] ‘1.26.0’

Best, Leo