EWCE.R
In EWCE: Expression Weighted Celltype Enrichment

## ----style, echo=FALSE, results='asis', message=FALSE--------------------
BiocStyle::markdown()
knitr::opts_chunk$set(tidy         = FALSE,
                      warning      = FALSE,
                      message      = FALSE)

## ------------------------------------------------------------------------
library(EWCE)

## ---- fig.width = 7, fig.height = 4--------------------------------------
set.seed(1234)
library(ggplot2)
library(reshape2)
data("celltype_data")
genes = c("Apoe","Gfap","Gapdh")
exp = melt(cbind(celltype_data$all_scts[genes,],genes),id.vars="genes")
colnames(exp) = c("Gene","Cell","AvgExp")
ggplot(exp)+geom_bar(aes(x=Cell,y=AvgExp),stat="identity")+facet_grid(Gene~.)+ 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

## ---- fig.width = 7, fig.height = 4--------------------------------------
exp = melt(cbind(celltype_data$subcell_dists[genes,],genes),id.vars="genes")
colnames(exp) = c("Gene","Cell","ExpProp")
ggplot(exp)+geom_bar(aes(x=Cell,y=ExpProp*100),stat="identity")+facet_grid(Gene~.)+ 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

## ---- fig.width = 7, fig.height = 4--------------------------------------
exp = melt(cbind(celltype_data$cell_dists[genes,],genes),id.vars="genes")
colnames(exp) = c("Gene","Cell","ExpProp")
ggplot(exp)+geom_bar(aes(x=Cell,y=ExpProp*100),stat="identity")+facet_grid(Gene~.)+ 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

## ------------------------------------------------------------------------
library(EWCE)
data("example_genelist")
print(example_genelist)

## ------------------------------------------------------------------------
data("mouse_to_human_homologs")
m2h = unique(mouse_to_human_homologs[,c("HGNC.symbol","MGI.symbol")])
mouse.hits = unique(m2h[m2h$HGNC.symbol %in% example_genelist,"MGI.symbol"])
mouse.bg  = unique(setdiff(m2h$MGI.symbol,mouse.hits))

## ------------------------------------------------------------------------
print(mouse.hits)

## ------------------------------------------------------------------------
reps=1000 # <- Use 1000 bootstrap lists so it runs quickly, for publishable analysis use >10000
subCellStatus=TRUE # <- Use subcell level annotations (i.e. Interneuron type 3)

## ----results='hide'------------------------------------------------------
# Bootstrap significance testing, without controlling for transcript length and GC content
full_results = bootstrap.enrichment.test(sct_data=celltype_data,mouse.hits=mouse.hits,
                    mouse.bg=mouse.bg,reps=reps,sub=subCellStatus)

## ------------------------------------------------------------------------
print(full_results$results[order(full_results$results$p),3:5][1:6,])

## ---- fig.width = 7, fig.height = 4--------------------------------------
print(ewce.plot(full_results$results,mtc_method="BH"))

## ----results='hide'------------------------------------------------------
generate.bootstrap.plots(sct_data=celltype_data,mouse.hits=mouse.hits,mouse.bg=mouse.bg,reps=100,sub=subCellStatus,full_results=full_results,listFileName="VignetteGraphs")

## ------------------------------------------------------------------------
human.hits = unique(m2h[m2h$HGNC.symbol %in% example_genelist,"HGNC.symbol"])
human.bg = unique(setdiff(m2h$HGNC.symbol,human.hits))

## ---- fig.width = 7, fig.height = 4--------------------------------------
print(ewce.plot(cont_results$results,mtc_method="BH"))

## ----results='hide', cache=FALSE, fig.width = 5, fig.height = 4----------
# Bootstrap significance testing controlling for transcript length and GC content
cont_results = bootstrap.enrichment.test(sct_data=celltype_data,human.hits=human.hits,
                    human.bg=human.bg,reps=reps,sub=FALSE,geneSizeControl=TRUE)
print(ewce.plot(cont_results$results,mtc_method="BH"))

## ----results='hide'------------------------------------------------------
gene.list.2 = mouse.bg[1:30]
second_results = bootstrap.enrichment.test(sct_data=celltype_data,mouse.hits=gene.list.2,
                    mouse.bg=mouse.bg,reps=reps,sub=subCellStatus)
full_res2 = data.frame(full_results$results,list="Alzheimers")
scnd_res2 = data.frame(second_results$results,list="Second")
merged_results = rbind(full_res2,scnd_res2)

## ----fig.width = 7, fig.height = 4---------------------------------------
print(ewce.plot(total_res=merged_results,mtc_method="BH"))

## ----results='hide'------------------------------------------------------
data(tt_alzh)
tt_results = ewce_expression_data(sct_data=celltype_data,tt=tt_alzh)

## ----fig.width = 7, fig.height = 4---------------------------------------
ewce.plot(tt_results$joint_results)

## ----results='hide'------------------------------------------------------
# Load the data
data(tt_alzh_BA36)
data(tt_alzh_BA44)

# Run EWCE analysis
tt_results_36 = ewce_expression_data(sct_data=celltype_data,tt=tt_alzh_BA36)
tt_results_44 = ewce_expression_data(sct_data=celltype_data,tt=tt_alzh_BA44)

# Fill a list with the results
results = add.res.to.merging.list(tt_results)
results = add.res.to.merging.list(tt_results_36,results)
results = add.res.to.merging.list(tt_results_44,results)

# Perform the merged analysis
merged_res = merged_ewce(results,reps=10) # <- For publication reps should be higher
print(merged_res)

## ----fig.width = 7, fig.height = 4---------------------------------------
print(ewce.plot(merged_res))