cell counts using betas

cellcounts.R

library(ggplot2)
library(devtools)
library(reshape2)
#install_github("brentp/celltypes450")
library(celltypes450)
library(meffil)
options(mc.cores=16)

# Using GSE55491
# Estimate cell counts using meffil
samplesheet <- meffil.create.samplesheet(".")
qc.objects <- meffil.qc(samplesheet, cell.type.reference="blood gse35069 complete", verbose=TRUE)
cellcounts_meffil <- t(sapply(qc.objects, function(x) x$cell.counts$counts))

# Estimate cell counts using celltypes450
norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=10)
norm.beta <- meffil.normalize.samples(norm.objects)

cellcounts_meffil_b <- meffil.estimate.cell.counts.from.betas(norm.beta, "blood gse35069 complete", verbose=TRUE)


m <- min(norm.beta[!norm.beta==0])
norm.beta[norm.beta == 0] <- m
cellcounts_brent <- adjust.beta(norm.beta, est.only=TRUE)
cellcounts_meffil <- cellcounts_meffil[, match(colnames(cellcounts_brent), colnames(cellcounts_meffil))]



plot_comparisons <- function(c1, c2, nom1, nom2)
{
	c1 <- c1[,colnames(c1) %in% colnames(c2)]
	c2 <- c2[,colnames(c2) %in% colnames(c1)]
	c1 <- c1[, match(colnames(c2), colnames(c1))]
	a <- melt(c1)
	b <- melt(c2)
	ab <- merge(a, b, by=c("Var1", "Var2"))	

	p <- ggplot(ab, aes(x = value.x, y=value.y)) +
		geom_point() +
		facet_wrap(~ Var2) +
		labs(x = nom1, y=nom2)
	return(list(p=p, dat=ab))
}


res <- plot_comparisons(cellcounts_meffil, cellcounts_meffil_b, "idat", "beta")
res$p
ggsave(filename="beta_vs_idat.pdf", plot=res$p)

# Plot estimates
a <- melt(cellcounts_brent)
b <- melt(cellcounts_meffil)
a$meffil <- b$value

p <- ggplot(a, aes(x = meffil, y=value)) +
	geom_point() +
	facet_wrap(~ Var2) +
	labs(x = "meffil", y="celltypes450")
ggsave(filename="meffil_vs_celltypes450.pdf", plot=p)



# Do a bad transformation of betas

norm.beta.squared <- 1 - sqrt(norm.beta) / 2
cellcounts_brent_squared <- adjust.beta(norm.beta.squared, est.only=TRUE)

a <- melt(cellcounts_brent_squared)
b <- melt(cellcounts_meffil)
a$meffil <- b$value

p <- ggplot(a, aes(x = meffil, y=value)) +
	geom_point() +
	facet_wrap(~ Var2) +
	labs(x = "meffil", y="celltypes450 squared")
ggsave(filename="meffil_vs_celltypes450_squared.pdf", plot=p)