fitNBthmDE {GeoDiff}R 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 = FALSE,
  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 #'

a list with parameter estimation #'

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))


[Package GeoDiff version 0.99.5 Index]