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

fitNBthmDER Documentation

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.