mashr

mashr.R

library(mashr)
library(corrplot)
library(remaCor)
library(ggplot2)

load("/sc/arion/projects/roussp01a/sanan/230706_NPSAD_filePrep/mashr/ckpts/mashr_ckpt2.RData")

# how many genes are non-zero in each tissue?
ngenes = apply(betaTab, 2, function(x) sum(x!=0))

df = data.frame(CellType = names(ngenes), ngenes = ngenes)

fig = ggplot(df, aes(CellType, ngenes )) +
	geom_bar(stat='identity') +
	coord_flip() 
ggsave(fig, file="~/www/test.png")

# keep tissues with > 7000 genes
# keep genes that are not all zero in the retained tissues
idx = which(ngenes > 7000)
keep = apply(betaTab[,idx], 1, function(x) max(abs(x))) > 0

# mashr
data = mash_set_data(betaTab[keep,idx], seTab[keep,idx])
U.c = cov_canonical(data) 
m.c = mash(data, U.c)


C = cor(get_pm(m.c))
corrplot(C)



gene = "CACNA1C"
m.c$result$lfsr[gene,]
beta = m.c$result$PosteriorMean[gene,]
se = m.c$result$PosteriorSD[gene,]

plotForest(beta, se) + ggtitle(gene)