Skip to content

Instantly share code, notes, and snippets.

Avatar

Mike Love mikelove

View GitHub Profile
@mikelove
mikelove / tximeta_plyranges.R
Created Jun 25, 2021
tximeta plyranges demo
View tximeta_plyranges.R
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
mikelove / data_gen_mutsig.R
Created Jun 25, 2021
Data generating model for mutational signatures
View data_gen_mutsig.R
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
mikelove / ihw_example.R
Created Jun 22, 2021
IHW on two covariates
View ihw_example.R
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
mikelove / boot_workflow.R
Created Jun 22, 2021
bootRanges workflow thoughts
View boot_workflow.R
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
mikelove / airpart.R
Last active Jun 25, 2021
airpart demo
View airpart.R
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
)
View goseq.R
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
mikelove / apeglm.R
Last active Jun 1, 2021
apeglm code for ASE using ashr / locfdr
View apeglm.R
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)
View discordant_dge_eGenes.R
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
mikelove / eqtl_ase.R
Created Apr 25, 2021
eqtl_ase thoughts
View eqtl_ase.R
# 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
mikelove / ratio.R
Created Apr 15, 2021
Allelic ratio with inferential variance
View ratio.R
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