## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(
    collapse=TRUE,
    comment="#>"
)

## ----setup, eval=FALSE--------------------------------------------------------
# if (!require("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("CSOA")

## ----message=FALSE, warning=FALSE, results=FALSE------------------------------
library(CSOA)
library(ggplot2)
library(patchwork)
library(scRNAseq)
library(scuttle)
library(Seurat)

sceObj <- BaronPancreasData('human')
sceObj <- logNormCounts(sceObj)

## -----------------------------------------------------------------------------
table(colData(sceObj)$label)

## -----------------------------------------------------------------------------
acinarMarkers <- c('PRSS1', 'KLK1', 'CTRC', 'PNLIP', 'AKR1C3', 'CTRB1', 
                   'DUOXA2', 'ALDOB', 'REG3A', 'SERPINA3', 'PRSS3', 'REG1B', 
                   'CFB',	'GDF15',	'MUC1','ANPEP', 'ANGPTL4', 'OLFM4', 
                   'GSTA1', 'LGALS2', 'PDZK1IP1', 'RARRES2', 'CXCL17', 
                   'UBD', 'GSTA2', 'LYZ', 'RBPJL', 'PTF1A', 'CELA3A', 
                   'SPINK1', 'ZG16', 'CEL', 'CELA2A', 'CPB1', 'CELA1', 
                   'PNLIPRP1', 'RNASE1', 'AMY2B', 'CPA2','CPA1', 'CELA3B', 
                   'CTRB2', 'PLA2G1B', 'PRSS2', 'CLPS', 'REG1A', 'SYCN')

alphaMarkers <- c('GCG', 'TTR', 'PCSK2', 'FXYD5', 'LDB2', 'MAFB',
                  'CHGA', 'SCGB2A1', 'GLS', 'FAP', 'DPP4', 'GPR119',
                  'PAX6', 'NEUROD1', 'LOXL4', 'PLCE1', 'GC', 'KLHL41',
                  'FEV', 'PTGER3', 'RFX6', 'SMARCA1', 'PGR', 'IRX1',
                  'UCP2', 'RGS4', 'KCNK16', 'GLP1R', 'ARX', 'POU3F4',
                  'RESP18', 'PYY', 'SLC38A5', 'TM4SF4', 'CRYBA2', 'SH3GL2', 
                  'PCSK1', 'PRRG2', 'IRX2', 'ALDH1A1','PEMT', 'SMIM24', 
                  'F10', 'SCGN', 'SLC30A8')

ductalMarkers <- c('CFTR', 'SERPINA5', 'SLPI', 'TFF1', 'CFB', 'LGALS4', 
                   'CTSH',	'PERP', 'PDLIM3',	'WFDC2', 'SLC3A1', 'AQP1', 
                   'ALDH1A3', 'VTCN1',	'KRT19', 'TFF2', 'KRT7', 'CLDN4',
                   'LAMB3', 'TACSTD2', 'CCL2', 'DCDC2','CXCL2', 'CLDN10',
                   'HNF1B', 'KRT20', 'MUC1', 'ONECUT1', 'AMBP', 'HHEX',
                   'ANXA4', 'SPP1', 'PDX1', 'SERPINA3', 'GDF15', 'AKR1C3',
                   'MMP7', 'DEFB1', 'SERPING1', 'TSPAN8', 'CLDN1', 'S100A10',
                   'PIGR')

gammaMarkers <- c('PPY', 'ABCC9', 'FGB', 'ZNF503', 'MEIS1', 'LMO3', 'EGR3', 
                  'CHN2', 'PTGFR', 'ENTPD2', 'AQP3', 'THSD7A', 'CARTPT',
                  'ISL1', 'PAX6', 'NEUROD1', 'APOBEC2', 'SEMA3E', 'SLITRK6',
                  'SERTM1', 'PXK', 'PPY2P', 'ETV1', 'ARX', 'CMTM8', 'SCGB2A1', 
                  'FXYD2', 'SCGN')

geneSets <- list(acinarMarkers, alphaMarkers, ductalMarkers, gammaMarkers)
names(geneSets) <- c('CSOA_acinar', 'CSOA_alpha', 'CSOA_ductal', 'CSOA_gamma')


## ----warning=FALSE------------------------------------------------------------
seuratObj <- as.Seurat(sceObj)
seuratObj <- runCSOA(seuratObj, geneSets)

## ----warning=FALSE, out.height='100%', out.width='100%', fig.height=8, fig.width=8----
plots <- lapply(names(geneSets), function(x) {
    VlnPlot(seuratObj, x, group.by = 'label') + NoLegend() +
        theme(axis.text = element_text(size=10),
              axis.title.x = element_blank(),
              plot.title = element_text (size=11))
})
    
(plots[[1]] | plots[[2]]) / (plots[[3]] | plots[[4]])

## -----------------------------------------------------------------------------
sceObj <- runCSOA(sceObj, geneSets)

## ----message=FALSE------------------------------------------------------------
sceObj <- runCSOA(sceObj, geneSets)
vapply(names(geneSets), 
       function(x) identical(seuratObj[[]][[x]], colData(sceObj)[[x]]),
       logical(1))

## -----------------------------------------------------------------------------
geneSetExp <- expMat(seuratObj)

## ----message=FALSE------------------------------------------------------------
res <- runCSOA(geneSetExp, geneSets)
head(res[[2]])

## -----------------------------------------------------------------------------
vapply(names(geneSets), 
       function(x) identical(seuratObj[[]][[x]], res[[2]][[x]]), logical(1))

## -----------------------------------------------------------------------------
sessionInfo()