Setup

By default, AUCell is installed only with the minimum dependencies. To run AUCell in parallel or run the examples in this tutorial, we recommend to install these packages:

source("http://bioconductor.org/biocLite.R")
# To support paralell execution:
biocLite(c("doMC", "doRNG"))
# For the main example:
biocLite(c("mixtools", "GEOquery", "SummarizedExperiment"))
# For the examples in the follow-up section of the tutorial:
biocLite(c("DT", "plotly", "NMF", "d3heatmap", "shiny", "rbokeh"))

Some tips…

# Explore tutorials in the web browser:
browseVignettes(package="AUCell") 

# Commnad line-based:
vignette(package="AUCell") # list
vignette("AUCell") # open

vignetteFile <- paste(file.path(system.file('doc', package='AUCell')), "AUCell.Rmd", sep="/")
# Copy to edit as markdown
file.copy(vignetteFile, ".")
# Alternative: extract R code
Stangle(vignetteFile)

0. Load scRNA-seq dataset and gene sets

The input data for AUCell are the expression matrix and the gene-sets.

dir.create("AUCell_tutorial")
setwd("AUCell_tutorial") # or in the first code chunk (kntr options), if running as Notebook

# i.e. Reading from a text file
exprMatrix <- read.table("myCountsMatrix.tsv")
exprMatrix <- as.matrix(exprMatrix)

# or using an expression set
exprMatrix <- exprs(myExpressionSet)

# (This may take a few minutes)
library(GEOquery)
attemptsLeft <- 20
while(attemptsLeft>0)
{
  geoFile <- tryCatch(getGEOSuppFiles("GSE60361", makeDirectory=FALSE), error=identity) 
  if("error" %in% class(geoFile)) 
  {
    attemptsLeft <- attemptsLeft-1
    Sys.sleep(5)
  }
  else
    attemptsLeft <- 0
}

gzFile <- grep(".txt.gz", basename(rownames(geoFile)), fixed=TRUE, value=TRUE)
message(gzFile)
txtFile <- gsub(".gz", "", gzFile, fixed=TRUE)
message(txtFile)
gunzip(filename=gzFile, destname=txtFile, remove=TRUE)

library(data.table)
geoData <- fread(txtFile, sep="\t")
geneNames <- unname(unlist(geoData[,1, with=FALSE]))
exprMatrix <- as.matrix(geoData[,-1, with=FALSE])
rm(geoData)
dim(exprMatrix)
rownames(exprMatrix) <- geneNames
exprMatrix[1:5,1:4]

# Remove file
file.remove(txtFile)

# Save for future use
mouseBrainExprMatrix <- exprMatrix
save(mouseBrainExprMatrix, file="exprMatrix_AUCellVignette_MouseBrain.RData")

# load("exprMatrix_AUCellVignette_MouseBrain.RData")
set.seed(333)
exprMatrix <- mouseBrainExprMatrix[sample(rownames(mouseBrainExprMatrix), 5000),]

dim(exprMatrix)

## [1] 5000 3005

library(GSEABase)
genes <- c("gene1", "gene2", "gene3")
geneSets <- GeneSet(genes, setName="geneSet1")
geneSets

## setName: geneSet1 
## geneIds: gene1, gene2, gene3 (total: 3)
## geneIdType: Null
## collectionType: Null 
## details: use 'details(object)'

library(AUCell)
library(GSEABase)
gmtFile <- paste(file.path(system.file('examples', package='AUCell')), "geneSignatures.gmt", sep="/")
geneSets <- getGmt(gmtFile)

geneSets <- subsetGeneSets(geneSets, rownames(exprMatrix)) 
cbind(nGenes(geneSets))

##                       [,1]
## Astrocyte_Cahoy        555
## Neuron_Cahoy           375
## Oligodendrocyte_Cahoy  417
## Astrocyte_Lein          12
## Neuron_Lein             14
## Microglia_lavin        149

geneSets <- setGeneSetNames(geneSets, newNames=paste(names(geneSets), " (", nGenes(geneSets) ,"g)", sep=""))

# Random
set.seed(321)
extraGeneSets <- c(
  GeneSet(sample(rownames(exprMatrix), 50), setName="Random (50g)"),
  GeneSet(sample(rownames(exprMatrix), 500), setName="Random (500g)"))

countsPerGene <- apply(exprMatrix, 1, function(x) sum(x>0))
# Housekeeping-like
extraGeneSets <- c(extraGeneSets,
                   GeneSet(sample(names(countsPerGene)[which(countsPerGene>quantile(countsPerGene, probs=.95))], 100), setName="HK-like (100g)"))

geneSets <- GeneSetCollection(c(geneSets,extraGeneSets))
names(geneSets)

## [1] "Astrocyte_Cahoy (555g)"       "Neuron_Cahoy (375g)"         
## [3] "Oligodendrocyte_Cahoy (417g)" "Astrocyte_Lein (12g)"        
## [5] "Neuron_Lein (14g)"            "Microglia_lavin (149g)"      
## [7] "Random (50g)"                 "Random (500g)"               
## [9] "HK-like (100g)"

1. Build gene-expression rankings for each cell

The first step to calculate the enrichment of a signature is to create the “rankings”. These rankings are only an intermediate step to calculate the AUC, but they are kept as a separate step in the workflow in order to provide more flexibility (i.e. to save them for future analyses, to merge datasets, or process them by parts).

For each cell, the genes are ranked from highest to lowest value. The genes with same expression value are shuffled. Therefore, genes with expression ‘0’ are randomly sorted at the end of the ranking. It is important to check that most cells have at least the number of expressed/detected genes that are going to be used to calculate the AUC (aucMaxRank in calcAUC()). The histogram provided by AUCell_buildRankings() allows to quickly check this distribution. plotGeneCount(exprMatrix) allows to obtain only the plot before building the rankings.

Since the rankings are created individually for each cell, in principle, it is possible to merge cell-rankings from different datasets. However, the datasets should be similar in regards to their “sensitivity” (e.g. the number of genes detected in the cells of each datasets), and the genes they include (e.g. same gene IDs).

cells_rankings <- AUCell_buildRankings(exprMatrix, nCores=1, plotStats=TRUE)

## Quantiles for the number of genes detected by cell: 
## (Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).

##     min      1%      5%     10%     50%    100% 
##  185.00  255.12  356.00  432.00  902.00 2038.00

cells_rankings

## Ranking for 5000 genes (rows) and 3005 cells (columns).
## 
## Top-left corner of the ranking:
##                     cells
## genes                1772071015_C02 1772071017_G12 1772071017_A05
##   0610010B08Rik_loc4           2499           2954           2799
##   0610010K14Rik                3293            445           3485
##   0610031J06Rik                 862           1045           1102
##   0610037L13Rik                1047           2367           4809
##   0610039K10Rik                4672           4228           4199
##   0610040J01Rik                4034           1584           3729
##                     cells
## genes                1772071014_B06 1772067065_H06
##   0610010B08Rik_loc4           4932           1351
##   0610010K14Rik                3947            239
##   0610031J06Rik                4906            898
##   0610037L13Rik                3080            391
##   0610039K10Rik                3707           4924
##   0610040J01Rik                4353           4128
## 
## Quantiles for the number of genes detected by cell:
##     min      1%      5%     10%     50%    100% 
##  185.00  255.12  356.00  432.00  902.00 2038.00

The “rankings” can be seen as a new representation of the original dataset. Once they are calculated, they can be saved for future analyses.

save(cells_rankings, file="cells_rankings.RData")

2. Calculate enrichment for the gene signatures (AUC)

To determine whether the gene set is enriched at the top of the gene-ranking for each cell, AUCell uses the “Area Under the Curve” (AUC) of the recovery curve.

The function AUCell_calcAUC calculates this score, and returns a matrix with an AUC score for each gene-set in each cell.

cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings)
save(cells_AUC, file="cells_AUC.RData")

In order to calculate the AUC, by default only the top 5% of the genes in the ranking are used (i.e. checks whether the genes in the gene-set or signature are within the top 5%). This allows faster execution on bigger datasets, and reduce the effect of the noise at the bottom of the ranking (e.g. where many genes might be tied at 0 counts). The percentage to take into account can be modified with the argument aucMaxRank. For datasets where most cells express many genes (e.g. a filtered dataset), or these have high expression values, it might be good to increase this threshold. Check the histogram provided by AUCell_buildRankings to get an estimation on where this threshold lies within the dataset.

3. Determine the cells with the given gene signatures or active gene sets

In summary: The AUC represents the proportion of expressed genes in the signature, and their relative expression value compared to the other genes within the cell. We can use this propperty to explore the population of cells that are present in the dataset according to the expression of the gene-set.

The AUC estimates the proportion of genes in the gene-set that are highly expressed in each cell. Cells expressing many genes from the gene-set will have higher AUC values than cells expressing fewer (i.e. compensating for housekeeping genes, or genes that are highly expressed in all the cells in the dataset). Since the AUC represents the proportion of expressed genes in the signature, we can use the relative AUCs across the cells to explore the population of cells that are present in the dataset according to the expression of the gene-set.

However, determining whether the signature is active (or not) in a given cell is not always trivial. The AUC is not an absolute value, but it depends on the the cell type (i.e. sell size, amount of transcripts), the specific dataset (i.e. sensitivity of the measures) and the gene-set. It is often not straight forward to obtain a pruned signature of clear marker genes that are completely “on” in the cell type of interest and off" in every other cell. In addition, at single-cell level, most genes are not expressed or detected at a constant level.

The ideal situation will be a bi-modal distribution, in which most cells in the dataset have a low “AUC” compared to a population of cells with a clearly higher value (i.e. see “Oligodendrocites” in the next figure). This is normally the case on gene sets that are active mostly in a population of cells with a good representation in the dataset (e.g. ~ 5-30% of cells in the dataset). Similar cases of “marker” gene sets but with different proportions of cells in the datasets are the “neurons” and “microglia” (see figure). When there are very few cells within the dataset, the distribution might look normal-like, but with some outliers to the higher end (e.g. microglia). While if the gene set is marker of a high percentage of cells in the dataset (i.e. neurons), the distribution might start approaching the look of a gene-set of housekeeping genes. As example, the ‘housekeeping’ gene-set in the figure includes genes that are detected in most cells in the dataset.

Note that the size of the gene-set will also affect the results. With smaller gene-genes (fewer genes), it is more likely to get cells with AUC = 0. While this is the case of the “perfect markers” it is also easier to get it by chance with smal datasets. (i.e. Random gene set with 50 genes in the figure). Bigger gene-sets (100-2k) can be more stable and easier to evaluate, as big random gene sets will approach the normal distibution.

To ease the exploration of the distributions, the function AUCell_exploreThresholds() automatically plots all the histograms and calculates several thresholds that could be used to consider a gene-set ‘active’ (returned in $aucThr). The distributions are plotted as dotted lines over the histogram and the corresponding thresholds as vertical bars in the matching color. The thicker vertical line indicates the threshold selected by default ($aucThr$selected): the highest value to reduce the false positives.

Note: This function makes use of package “mixtools” to explore the distributions. It is not essential, but we recommend to install it: source("http://bioconductor.org/biocLite.R"); biocLite("mixtools")

set.seed(123)
par(mfrow=c(3,3)) 
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, assign=TRUE) 

warningMsg <- sapply(cells_assignment, function(x) x$aucThr$comment)
warningMsg[which(warningMsg!="")]

##                                                                                         Random (50g) 
## "The AUC might follow a normal distribution (random gene-set?).  The right distribution is taller. " 
##                                                                                        Random (500g) 
## "The AUC might follow a normal distribution (random gene-set?).  The right distribution is taller. "

The thresholds calcuated for each gene set are stored in the $aucThr slot. For example, the thresholds suggested for the oligodendrocyte gene-set:

cells_assignment$Oligodendrocyte_Cahoy$aucThr$thresholds

##             threshold nCells
## Global_k1   0.2519540    747
## L_k2        0.2322519    794
## R_k3        0.1713901    955
## minimumDens 0.2470856    763

To ease for the inclusion of AUCell in workflows (e.g. SCENIC), this function can also provide the list of cells with an AUC value over this threshold ($assignment). However, keep in mind that the threshold selection in the current version is not exhaustive, so we highly recommend to check the AUC histograms and manually select the threshold if needed.

To obtain the threshold selected automatically for a given gene set(e.g. Oligodendrocytes):

cells_assignment$Oligodendrocyte_Cahoy$aucThr$selected

## minimumDens 
##   0.2470856

# getThresholdSelected(cells_assignment)

Cells assigned at this threshold:

oligodencrocytesAssigned <- cells_assignment$Oligodendrocyte_Cahoy$assignment
length(oligodencrocytesAssigned)

## [1] 763

head(oligodencrocytesAssigned)

## [1] "1772067069_B02" "1772071041_F09" "1772071041_A02" "1772067073_C04"
## [5] "1772071040_G12" "1772066077_G03"

Plotting the AUC histogram of a specific gene set, and setting a new threshold:

geneSetName <- rownames(cells_AUC)[grep("Oligodendrocyte_Cahoy", rownames(cells_AUC))]
AUCell_plotHist(cells_AUC[geneSetName,], aucThr=0.25)
abline(v=0.25)

Assigning cells to this new threshold:

newSelectedCells <- names(which(getAUC(cells_AUC)[geneSetName,]>0.08))
length(newSelectedCells)

## [1] 2961

head(newSelectedCells)

## [1] "1772071015_C02" "1772071017_G12" "1772071017_A05" "1772071014_B06"
## [5] "1772067065_H06" "1772071017_E02"

Exploring the cell-assignment (table & heatmap)

The cell assignment is stored in the $assignment slot. Here we show a few ideas on how they can be extracted and visualized.

Extract these cells for all the gene-sets and transform it into a table:

cellsAssigned <- lapply(cells_assignment, function(x) x$assignment)
assignmentTable <- reshape2::melt(cellsAssigned, value.name="cell")
colnames(assignmentTable)[2] <- "geneSet"
head(assignmentTable)

##             cell                geneSet
## 1 1772066080_D10 Astrocyte_Cahoy (555g)
## 2 1772066107_E06 Astrocyte_Cahoy (555g)
## 3 1772071036_B05 Astrocyte_Cahoy (555g)
## 4 1772062109_E07 Astrocyte_Cahoy (555g)
## 5 1772058148_A03 Astrocyte_Cahoy (555g)
## 6 1772063063_B10 Astrocyte_Cahoy (555g)

Convert into an incidence matrix and plot as a histogram:

assignmentMat <- table(assignmentTable[,"geneSet"], assignmentTable[,"cell"])
assignmentMat[,1:2]

##                               
##                                1772058148_A03 1772058148_A12
##   Astrocyte_Cahoy (555g)                    1              1
##   Astrocyte_Lein (12g)                      1              0
##   Microglia_lavin (149g)                    0              0
##   Neuron_Cahoy (375g)                       0              1
##   Neuron_Lein (14g)                         0              1
##   Oligodendrocyte_Cahoy (417g)              1              0

set.seed(123)
miniAssigMat <- assignmentMat[,sample(1:ncol(assignmentMat),100)]
library(NMF)
aheatmap(miniAssigMat, scale="none", color="black", legend=FALSE)

Some interactive options (output not shown):

library(d3heatmap)
d3heatmap(matrix(miniAssigMat), scale="none", colors=c("white", "black"))

library(DT)
datatable(assignmentTable, options = list(pageLength = 10), filter="top")

Explore cells/clusters based on the signature score

The AUC score can also be used to explore the output from previous clustering analyses (or vice-versa).

In this example, we will use the AUC obtained for the diferent signatures to color a t-SNE based previously run on the whole expression matrix (i.e. not the 5000 random genes).

# Load the tSNE (included in the package)
load(paste(file.path(system.file('examples', package='AUCell')), "cellsTsne.RData", sep="/"))
cellsTsne <- cellsTsne$Y
plot(cellsTsne, pch=16, cex=.3)

This t-SNE was created with this code (it takes a while to run):

load("exprMatrix_AUCellVignette_MouseBrain.RData")
sumByGene <- apply(mouseBrainExprMatrix, 1, sum)
exprMatSubset <- mouseBrainExprMatrix[which(sumByGene>0),]
logMatrix <- log2(exprMatSubset+1)

library(Rtsne)
set.seed(123)
cellsTsne <- Rtsne(t(logMatrix))
rownames(cellsTsne$Y) <- colnames(logMatrix)
colnames(cellsTsne$Y) <- c("tsne1", "tsne2")
save(cellsTsne, file="cellsTsne.RData")

This t-SNE can be colored based on the AUC scores. To highlight the cluster of cells that are more likely of the cell type according to the signatures, we will split the cells into the cells that pass the assignment threshold (colored in shades of pink-red), and the cells that don’t (colored in black-blue).

Of course, the origin of the signatures should also be kept in mind for the interpretation (for example, these signatures were obtained from bulk RNA-seq analyses). Also, running this tutorial only on 5000 genes probably introduced some extra noise.

For this example we have used the thresholds assigned automatically:

selectedThresholds <- getThresholdSelected(cells_assignment)
par(mfrow=c(2,3)) # Splits the plot into two rows and three columns
for(geneSetName in names(selectedThresholds))
{
  nBreaks <- 5 # Number of levels in the color palettes
  # Color palette for the cells that do not pass the threshold
  colorPal_Neg <- grDevices::colorRampPalette(c("black","blue", "skyblue"))(nBreaks)
  # Color palette for the cells that pass the threshold
  colorPal_Pos <- grDevices::colorRampPalette(c("pink", "magenta", "red"))(nBreaks)
  
  # Split cells according to their AUC value for the gene set
  passThreshold <- getAUC(cells_AUC)[geneSetName,] >  selectedThresholds[geneSetName]
  if(sum(passThreshold) >0 )
  {
    aucSplit <- split(getAUC(cells_AUC)[geneSetName,], passThreshold)
    
    # Assign cell color
    cellColor <- c(setNames(colorPal_Neg[cut(aucSplit[[1]], breaks=nBreaks)], names(aucSplit[[1]])), 
    setNames(colorPal_Pos[cut(aucSplit[[2]], breaks=nBreaks)], names(aucSplit[[2]])))
    
    # Plot
    plot(cellsTsne, main=geneSetName,
    sub="Pink/red cells pass the threshold",
    col=cellColor[rownames(cellsTsne)], pch=16) 
  }
}

This kind of plots can also be used to explore the assignment thresholds. i.e.:

selectedThresholds[2] <-  0.25
par(mfrow=c(2,3))
AUCell_plotTSNE(tSNE=cellsTsne, exprMat=exprMatrix, 
cellsAUC=cells_AUC[1:2,], thresholds=selectedThresholds)

library(shiny); library(rbokeh)
exprMat <- log2(mouseBrainExprMatrix+1) # (back to the full matrix)

# Create app 
aucellApp <- AUCell_createViewerApp(auc=cells_AUC, thresholds=selectedThresholds, 
tSNE=cellsTsne, exprMat=exprMat)

# Run (the exact commands depend on the R settings, see Shiny's doc for help)
options(shiny.host="0.0.0.0") 
savedSelections <- runApp(aucellApp) 

# Other common settings:
# host = getOption("shiny.host", "127.0.0.1") 
# shinyApp(ui=aucellApp$ui, server=aucellApp$server)

# Optional: Visualize extra information about the cells (categorical variable)
cellInfoFile <- paste(file.path(system.file('examples', package='AUCell')), 
"mouseBrain_cellLabels.tsv", sep="/")
cellInfo <- read.table(cellInfoFile, row.names=1, header=TRUE, sep="\t")
head(cellInfo)

geneSetName <- "Astrocyte_Cahoy (555g)"

library(rbokeh)
tSNE.df <- data.frame(cellsTsne, cell=rownames(cellsTsne))
figure() %>%
ly_points(tsne1, tsne2, data=tSNE.df, size=1, legend = FALSE,
hover=cell, color=getAUC(cells_AUC)[geneSetName,rownames(tSNE.df)]) %>%
set_palette(continuous_color = pal_gradient(c("lightgrey", "pink", "red")))

Comparison with mean

One of the most intuitive ways to evaluate the expression of a gene signature is to calculate the mean expression. However, calculating the average of the gene signature on the raw counts will lead to cells with more reads (either for technical or biological reasons) to have higher average values for most signatures. To reduce this effect, some kind of library-size normalization would be required. Since AUCell evaluates the dectection of the signatures on each cell individually, it automatically compensates for library-size and makes its application easier and independent of the units of the data and other normalizations applied (e.g. raw counts, TPM, UMI counts).

logMat <- log2(exprMatrix+2)
meanByGs <- t(sapply(geneSets, function(gs) colMeans(logMat[geneIds(gs),])))
rownames(meanByGs) <- names(geneSets)

colorPal <- grDevices::colorRampPalette(c("black", "red"))(nBreaks)
par(mfrow=c(3,3))
for(geneSetName in names(geneSets))
{
  cellColor <- setNames(colorPal[cut(meanByGs[geneSetName,], breaks=nBreaks)], names(meanByGs[geneSetName,]))
  plot(cellsTsne, main=geneSetName, axes=FALSE, xlab="", ylab="",
  sub="Expression mean",
  col=cellColor[rownames(cellsTsne)], pch=16)
}

AUCell_plotTSNE(tSNE=cellsTsne, exprMat=exprMatrix, plots = "AUC",
cellsAUC=cells_AUC[1:9,], thresholds=selectedThresholds)

nGenesPerCell <- apply(exprMatrix, 2, function(x) sum(x>0))
colorPal <- grDevices::colorRampPalette(c("darkgreen", "yellow","red"))
cellColorNgenes <- setNames(adjustcolor(colorPal(10), alpha=.8)[as.numeric(cut(nGenesPerCell,breaks=10, right=FALSE,include.lowest=TRUE))], names(nGenesPerCell))
plot(cellsTsne, axes=FALSE, xlab="", ylab="",
sub="Number of detected genes",
col=cellColorNgenes[rownames(cellsTsne)], pch=16)

How reliable is AUCell? (Confusion matrix)

When using AUCell to identify cell types, the reliability of the cells identified will depend on the signatures used, and the differences within the cell types in the dataset. In any case, AUCell is quite robust to noise, both in the data and in the gene-set.

This is the comparison of the cells identified with AUCell (this example) and the cell-types assigned in the initial publication:

# "Real" cell type (e.g. provided in the publication)
cellLabels <- paste(file.path(system.file('examples', package='AUCell')), "mouseBrain_cellLabels.tsv", sep="/")
cellLabels <- read.table(cellLabels, row.names=1, header=TRUE, sep="\t")

# Confusion matrix:
cellTypeNames <- unique(cellLabels[,"level1class"])
confMatrix <- t(sapply(cells_assignment[c(1,4,2,5,3,6:9)], 
function(x) table(cellLabels[x$assignment,])[cellTypeNames]))
colnames(confMatrix) <- cellTypeNames
confMatrix[which(is.na(confMatrix), arr.ind=TRUE)] <- 0
confMatrix

##                              interneurons pyramidal SS pyramidal CA1
## Astrocyte_Cahoy (555g)                  0            0             3
## Astrocyte_Lein (12g)                    3            7            35
## Neuron_Cahoy (375g)                   290          397           937
## Neuron_Lein (14g)                     285          398           934
## Oligodendrocyte_Cahoy (417g)            0            5             4
## Microglia_lavin (149g)                  0            0             0
## Random (50g)                            0            0             0
## Random (500g)                           0            0             0
## HK-like (100g)                          0            0             0
##                              oligodendrocytes microglia endothelial-mural
## Astrocyte_Cahoy (555g)                     12         8                14
## Astrocyte_Lein (12g)                       55         4                 6
## Neuron_Cahoy (375g)                        57        11                45
## Neuron_Lein (14g)                          85        14                44
## Oligodendrocyte_Cahoy (417g)              748         1                 1
## Microglia_lavin (149g)                      0        40                 0
## Random (50g)                                0         0                 0
## Random (500g)                               0         0                 0
## HK-like (100g)                              0         0                 0
##                              astrocytes_ependymal
## Astrocyte_Cahoy (555g)                        188
## Astrocyte_Lein (12g)                          201
## Neuron_Cahoy (375g)                            29
## Neuron_Lein (14g)                              36
## Oligodendrocyte_Cahoy (417g)                    4
## Microglia_lavin (149g)                          0
## Random (50g)                                    0
## Random (500g)                                   0
## HK-like (100g)                                  0

Notes on interpretation:

• Interneurons & pyramidal are “neurons”.

• According to other publications that have used this dataset to evaluate new methods, a relevant group of cells marked as oligodendrocytes are likely doublets (i.e. two cells, for example a mix of an oligodendrocyte + another cell type).