Element level analysis with mpralm

element_level.R

set.seed(5)

n <- 1000
reps <- 10
rna <- matrix(
  rnbinom(n * reps, mu = 10, size = 100),
  ncol=reps
)
dna <- matrix(
  rnbinom(n * reps, mu = 10, size = 100),
  ncol=reps
)

# spike in a few ones above the threshold
rna[1:3,] <- rnbinom(3 * reps, mu=100, size=100)

rownames(rna) <- rownames(dna) <- paste0("variant", 1:n)

library(mpra)
mpraset <- MPRASet(DNA = dna, 
                   RNA = rna, 
                   eid = rownames(dna), 
                   eseq = NULL, 
                   barcode = NULL)

thr <- 1.234 # threshold determined from upper quantile of shuffled

design <- model.matrix(~1, data=data.frame(sample=1:reps))

fit <- mpralm(object = mpraset, design = design, aggregate = "none",
              normalize = FALSE, model_type = "indep_groups", plot = TRUE)

tr <- treat(fit, lfc=thr)

mpra_result <- topTreat(tr, coef = 1, number = Inf)

library(tibble)
mpra_result |> as_tibble()

mpra_result[1:3,1]
log2(rna[1:3,]/dna[1:3,])

 mpra_result[1:3,1]
 w <- mpra:::.fit_standard(mpraset, design, aggregate="none", normalize=FALSE, plot=FALSE, return_weights=TRUE)
 rowSums(log2((rna[1:3,]+1)/(dna[1:3,]+1)) * w[1:3,]) / rowSums(w[1:3,])
 rowMeans(log2((rna[1:3,]+1)/(dna[1:3,]+1)))