mikelove

library(macrophage)
dir <- system.file("extdata", package="macrophage")
coldata <- read.csv(file.path(dir, "coldata.csv"))
coldata <- coldata[,c(1,2,3,5)]
names(coldata) <- c("names","id","line","condition")
coldata$files <- file.path(dir, "quants", coldata$names, "quant.sf.gz")
coldata$line <- factor(coldata$line)
lvls <- c("naive","IFNg")
coldata <- coldata[coldata$condition %in% lvls,]
coldata$condition <- factor(coldata$condition, lvls)

mikelove

library(rafalib)

J <- 100 # number of samples
# X - covariates
X <- cbind(rep(c(0,1,0,1),each=J/4),rep(0:1,each=J/2))
imagemat(X)
K <- 3 # number of signatures
# beta - associations with K = 3 signatures
beta <- cbind(c(1, -1), c(-1, 0), c(0, 1))
sigma <- 0.1

mikelove

library(ggplot2)
library(dplyr)

set.seed(1)
n <- 1e5
cov1 <- runif(n, -1, 1)
cov2 <- runif(n, -1, 1)
prop_alt <- 1/(1 + exp(-1 * (3 * cov1 + 2 * cov2 - 5)))
pvalue <- ifelse(rbinom(n, size=1, prop_alt),
                 rbeta(n, 0.25, 1),

mikelove

library(plyranges)
#library(nullranges)
load_all("~/proj/nullranges/nullranges")

genes <- GRanges("chr1", IRanges(1:10 * 100 + 1, width=10), strand="+")
pks <- GRanges("chr1", IRanges(1:20 + 90, width=10), id=paste0("pk",1:20))
seqlengths(pks) <- seqlengths(genes) <- c("chr1"=1100)

genes %>%
  anchor_5p %>%

mikelove

library(airpart)
suppressPackageStartupMessages(library(SingleCellExperiment))
suppressPackageStartupMessages(library(SummarizedExperiment))
p.vec <- c(rep(c(0.5, 0.8), each = 2), 0.65)
set.seed(1)
sce <- makeSimulatedData(
  mu1 = 2, mu2 = 10, nct = 5, n = 20,
  ngenecl = 20, theta = 20, ncl = 1,
  p.vec = p.vec
)

mikelove

library(goseq)
all.genes <- rownames(res)[res$baseMean > XYZ]
de.genes <- rownames(res)[res$padj < 0.1]
genes <- as.integer(all.genes %in% de.genes)
table(genes)
names(genes) <- sub("\\..*","",all.genes)
pwf <- nullp(genes, "hg19", "ensGene")
res <- goseq(pwf, "hg19", "ensGene")
head(res)

mikelove

library(apeglm)
suppressPackageStartupMessages(library(SummarizedExperiment))
load("se.oracle.rda")
n <- 10
cts <- assay(wide)[,1:n] + assay(wide)[,(n+1):(2*n)]
dim(cts)
ase_cts <- assay(wide)[,1:n]
mode(ase_cts) <- "integer"

theta_hat <- matrix(nrow=nrow(cts), ncol=niter+1)

mikelove

load("invivo_invitro_overlap_mvalue.RDat")
rownames(DGE.lm.vivo) <- substr(rownames(DGE.lm.vivo),1,15)
all.equal(rownames(DGE.lm.vivo),rownames(DGE.lm.vitro))
idx <- DGE.lm.vivo$adj.P.Val < 0.01 & DGE.lm.vitro$adj.P.Val < 0.01 &
  sign(DGE.lm.vivo$logFC) != sign(DGE.lm.vitro$logFC)
plot(sign(DGE.lm.vitro$logFC)*-log10(DGE.lm.vitro$adj.P.Val),
     sign(DGE.lm.vivo$logFC)*-log10(DGE.lm.vivo$adj.P.Val),
     col=idx+1,xlab="vitro",ylab="vivo",cex=.2)
abline(v=0,h=0,col="dodgerblue")
discordant <- rownames(DGE.lm.vivo)[idx]

mikelove

# simulation of a haploid population
n <- 1000
# genotype at causal variant 1
causal1 <- rbinom(n, size=1, prob=.5)
# how many alleles to flip to create variant 2
flip_c <- .33
causal2 <- abs(causal1 - rbinom(n, size=1, prob=flip_c))
# correlation btwn genotypes at variants 1 and 2
sprintf("1-2 cor: %.2f", cor(causal1, causal2))
# how many alleles to flip to create exonic variant

mikelove

library(fishpond)
library(SingleCellExperiment)
se <- makeSimSwishData(m=20, n=10)
sce <- as(se, "SingleCellExperiment")
assayNames(sce)
# make up the structure we will be having:
rownames(sce) <- paste0("txp",1:10,"_",rep(c("M","P"),each=10))

# about here it would be good to filter out any transcripts with very low expression
# in either species... i don't have code for that here, but we can work on that when