Skip to content

Instantly share code, notes, and snippets.

@mikelove
Created September 27, 2022 21:09
Show Gist options
  • Save mikelove/31003a09103943b99ea853d8733bdd4f to your computer and use it in GitHub Desktop.
Save mikelove/31003a09103943b99ea853d8733bdd4f to your computer and use it in GitHub Desktop.
Example of using mpralm with a threshold defined from shuffled seq
library(mpra)
E <- 1e3 # Number of elements
B <- 3 # Number of barcodes
s <- 4 # Samples per tissue
nt <- 2 # Number of tissues
set.seed(434)
rna <- matrix(rpois(E*B*s*nt, lambda = 2 * 30), nrow = E*B, ncol = s*nt)
dna <- matrix(rpois(E*B*s*nt, lambda = 30), nrow = E*B, ncol = s*nt)
rn <- as.character(outer(paste0("barcode_", seq_len(B), "_"),
paste0("elem_", seq_len(E)), FUN = "paste0"))
cn <- c(paste0("liver_", seq_len(s)), paste0("kidney_", seq_len(s)))
rownames(rna) <- rn
rownames(dna) <- rn
colnames(rna) <- cn
colnames(dna) <- cn
eid <- rep(paste0("elem_", seq_len(E)), each = B)
m <- MPRASet(DNA = dna, RNA = rna, eid = eid, eseq = NULL, barcode = NULL)
# note bug in printing code from use of cat()
design <- data.frame(intcpt = rep(1, length(cn)))
fit <- mpralm(object = m, design = design, aggregate = "sum",
normalize = FALSE, model_type = "indep_groups", plot = TRUE)
tr <- treat(fit, lfc=1)
tt <- topTreat(tr, coef = 1, number = Inf)
with(tt, plot(logFC, -log10(P.Value))); abline(v=1, col="magenta", lwd=2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment