Skip to content

Instantly share code, notes, and snippets.

@mikelove
Last active February 9, 2024 20:00
Show Gist options
  • Save mikelove/e2d8fa7143e1a70cca20c28b98c336bc to your computer and use it in GitHub Desktop.
Save mikelove/e2d8fa7143e1a70cca20c28b98c336bc to your computer and use it in GitHub Desktop.
Element level analysis with mpralm
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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment