fitNBthmDE-methods: Negative Binomial threshold mixed model for differential...

In Nanostring-Biostats/GeoDiff: Count model based differential expression and normalization on GeoMx RNA data

Negative Binomial threshold mixed model for differential expression analysis

Description

Negative Binomial threshold mixed model for differential expression analysis

Negative Binomial threshold mixed model for differential expression analysis

Usage

fitNBthmDE(object, ...)

## S4 method for signature 'NanoStringGeoMxSet'
fitNBthmDE(
  object,
  form,
  split,
  ROIs_high = NULL,
  features_all = NULL,
  sizefact = NULL,
  sizefact_BG = NULL,
  preci1,
  threshold_mean = NULL,
  preci2 = 10000,
  sizescalebythreshold = TRUE,
  controlRandom = list()
)

## S4 method for signature 'matrix'
fitNBthmDE(
  form,
  annot,
  object,
  probenum = rep(1, NROW(object)),
  features_all,
  sizefact,
  sizefact_BG,
  preci1,
  threshold_mean = NULL,
  preci2 = 10000,
  sizescalebythreshold = TRUE,
  controlRandom = list()
)

Arguments

object	count matrix with features in rows and samples in columns
...	additional argument list that might be used
form	model formula
split	indicator variable on whether it is for multiple slides (Yes, TRUE; No, FALSE)
ROIs_high	ROIs with high expressions defined based on featfact and featfact
features_all	vector of all features to be run
sizefact	size factor
sizefact_BG	size factor for background
preci1	precision matrix for regression coefficients
threshold_mean	average background level
preci2	precision for the background, default=10000
sizescalebythreshold	whether to scale the size factor, default=TRUE
controlRandom	list of random effect control parameters
annot	annotations files with variables in the formula
probenum	a vector of numbers of probes in each gene, default = rep(1, NROW(object))

Value

a list with parameter estimation #'

• X, design matrix for fixed effect

• Z, design matrix for random effect

• rt, random effect terms

• para0, =NA

• para, estimated parameters, including regression coefficients, r and threshold in rows and features in columns

• sizefact, same as input sizefact

• sizefact0, NA

• preci1, input precision matrix for regression coefficients

• Im0, NA

• Im, Information matrix of parameters

• conv0, NA

• conv, vector of convergence, 0 converged, 1 not converged

• features_high, NA

• features_all, same as the input features_all

• theta, list of estimated random effect parameters

• MAP random effect

a list with parameter estimation #'

• X, design matrix for fixed effect

• Z, design matrix for random effect

• rt, random effect terms

• para0, =NA

• para, estimated parameters, including regression coefficients, r and threshold in rows and features in columns

• sizefact, same as input sizefact

• sizefact0, NA

• preci1, input precision matrix for regression coefficients

• Im0, NA

• Im, Information matrix of parameters

• conv0, NA

• conv, vector of convergence, 0 converged, 1 not converged

• features_high, NA

• features_all, same as the input features_all

• theta, list of estimated random effect parameters(for relative covariance matrix)

• varcov, list of estimated variance covariance parameter estimation

• MAP random effect

Examples

library(Biobase)
library(dplyr)
data(demoData)
demoData <- demoData[, c(1:5, 33:37)]
demoData <- fitPoisBG(demoData, size_scale = "sum")
demoData <- aggreprobe(demoData, use = "cor")
demoData <- BGScoreTest(demoData)
demoData$slidename <- substr(demoData[["slide name"]], 12, 17)
thmean <- 1 * mean(fData(demoData)$featfact, na.rm = TRUE)
demo_pos <- demoData[which(!fData(demoData)$CodeClass == "Negative"), ]
demo_neg <- demoData[which(fData(demoData)$CodeClass == "Negative"), ]
sc1_scores <- fData(demo_pos)[, "scores"]
names(sc1_scores) <- fData(demo_pos)[, "TargetName"]
features_high <- ((sc1_scores > quantile(sc1_scores, probs = 0.4)) &
   (sc1_scores < quantile(sc1_scores, probs = 0.95))) |>
    which() |>
    names()
set.seed(123)
demoData <- fitNBth(demoData,
                    features_high = features_high,
                    sizefact_BG = demo_neg$sizefact,
                    threshold_start = thmean,
                    iterations = 5,
                    start_para = c(200, 1),
                    lower_sizefact = 0,
                    lower_threshold = 100,
                    tol = 1e-8)
ROIs_high <- sampleNames(demoData)[which(demoData$sizefact_fitNBth * thmean > 2)]
features_all <- rownames(demo_pos)

pData(demoData)$group <- c(rep(1, 5), rep(2, 5))


NBthDEmod2 <- fitNBthDE(form = ~group,
                     split = FALSE,
                     object = demoData,
                     ROIs_high = ROIs_high,
                     features_high = features_high,
                     features_all = features_all,
                     sizefact_start = demoData[, ROIs_high][['sizefact_fitNBth']],
                     sizefact_BG = demoData[, ROIs_high][['sizefact']],
                     threshold_mean = notes(demoData)[["threshold"]],
                     preci2=10000,
                     prior_type="contrast",
                     covrob=FALSE,
                     preci1con=1/25,
                     sizescalebythreshold=TRUE)

set.seed(123)
NBthmDEmod1 <- fitNBthmDE(
    form = ~ group + (1 | `slide name`),
    split = FALSE,
    object = demoData,
    ROIs_high = ROIs_high,
    features_all = features_all[1:5],
    sizefact = demoData[, ROIs_high][["sizefact_fitNBth"]],
    sizefact_BG = demoData[, ROIs_high][["sizefact"]],
    preci1=NBthDEmod2$preci1,
    threshold_mean = thmean,
    preci2=10000,
    sizescale = TRUE,
    controlRandom=list(nu=12, nmh_e=400, thin_e=60))

Nanostring-Biostats/GeoDiff documentation built on April 11, 2024, 5:31 a.m.