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

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

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

seuratObj <- as.Seurat(sceObj)

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')

## -----------------------------------------------------------------------------
geneSetExp <- expMat(seuratObj, acinarMarkers)
class(geneSetExp)
dim(geneSetExp)

## -----------------------------------------------------------------------------
cellSets <- percentileSets(geneSetExp, percentile=90)
vapply(cellSets, length, integer(1))

## ----out.height='90%', out.width='90%', fig.height=6, fig.width=5-------------
mat <- cellDistribution(cellSets, allCells = colnames(seuratObj))
basicHeatmap(mat, c('Gene', 'Cell', 'Presence'), 
             title='Distribution of cell sets among all cells') + NoLegend()

## -----------------------------------------------------------------------------
cso <- cellSetsOverlaps(cellSets, nCells=dim(seuratObj)[2])

## -----------------------------------------------------------------------------
head(cso)

## -----------------------------------------------------------------------------
cso2 <- generateOverlaps(geneSetExp, percentile = 90)
identical(cso, cso2)

## -----------------------------------------------------------------------------
po <- processOverlaps(cso)
head(po)

## ----out.height='100%', out.width='100%', fig.height=5, fig.width=5-----------
overlapCutoffPlot(po)

## -----------------------------------------------------------------------------
scoreDF <- scoreCells(geneSetExp, cso, setPairs=list(overlapPairs(cso)), 
                      geneSetNames='CSOA_acinar')
head(scoreDF)

## ----message=FALSE------------------------------------------------------------
seuratObj <- attachCellScores(seuratObj, scoreDF)
head(seuratObj$CSOA_acinar)

## -----------------------------------------------------------------------------
geneSet <- list(CSOA_acinar = acinarMarkers)
seuratObj <- runCSOA(seuratObj, geneSet)
identical(seuratObj[[]]$CSOA_acinar, scoreDF$CSOA_acinar)

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