Guided tutorial to COTAN V.2
9 June 2025
Package
COTAN 2.8.2
Contents
0.1 Preamble
library(COTAN)
library(zeallot)
library(rlang)
library(data.table)
library(Rtsne)
library(qpdf)
library(GEOquery)
library(ComplexHeatmap)
options(parallelly.fork.enable = TRUE)0.2 Introduction
This tutorial contains the same functionalities as the first release of the COTAN tutorial but done using the new and updated functions.
0.3 Get the data-set
Download the data-set for "mouse cortex E17.5".
dataDir <- tempdir()
GEO <- "GSM2861514"
fName <- "GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz"
dataSetFile <- file.path(dataDir, GEO, fName)
dir.create(file.path(dataDir, GEO), showWarnings = FALSE)
if (!file.exists(dataSetFile)) {
getGEOSuppFiles(GEO, makeDirectory = TRUE,
baseDir = dataDir, fetch_files = TRUE,
filter_regex = fName)
}
#> size
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz 1509523
#> isdir
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz FALSE
#> mode
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz 666
#> mtime
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz 2025-06-09 18:17:18
#> ctime
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz 2025-06-09 18:17:18
#> atime
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz 2025-06-09 18:17:18
#> exe
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz no
#> uname
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz biocbuild
#> udomain
#> E:\\biocbuild\\bbs-3.21-bioc\\tmpdir\\RtmpqgOHxN/GSM2861514/GSM2861514_E175_Only_Cortical_Cells_DGE.txt.gz palomino7
sample.dataset <- read.csv(dataSetFile, sep = "\t", row.names = 1L)Define a directory where the output will be stored.
outDir <- dataDir
# Log-level 2 was chosen to showcase better how the package works
# In normal usage a level of 0 or 1 is more appropriate
setLoggingLevel(2L)
#> Setting new log level to 2
# This file will contain all the logs produced by the package
# as if at the highest logging level
setLoggingFile(file.path(outDir, "vignette_v2.log"))
#> Setting log file to be: E:\biocbuild\bbs-3.21-bioc\tmpdir\RtmpqgOHxN/vignette_v2.log
message("COTAN uses the `torch` library when asked to `optimizeForSpeed`")
#> COTAN uses the `torch` library when asked to `optimizeForSpeed`
message("Run the command 'options(COTAN.UseTorch = FALSE)'",
" in your session to disable `torch` completely!")
#> Run the command 'options(COTAN.UseTorch = FALSE)' in your session to disable `torch` completely!
# this command does check whether the torch library is properly installed
c(useTorch, deviceStr) %<-% COTAN:::canUseTorch(TRUE, "cuda")
#> While trying to load the `torch` library Error in doTryCatch(return(expr), name, parentenv, handler): The `torch` library is installed but the required additional libraries are not avalable yet
#> Warning in value[[3L]](cond): The `torch` library is installed, but might
#> require further initialization
#> Warning in value[[3L]](cond): Please look at the `torch` package installation
#> guide to complete the installation
#> Warning in COTAN:::canUseTorch(TRUE, "cuda"): Falling back to legacy
#> [non-torch] code.
if (useTorch) {
message("The `torch` library is available and ready to use")
if (deviceStr == "cuda") {
message("The `torch` library can use the `CUDA` GPU")
} else {
message("The `torch` library can only use the CPU")
message("Please ensure you have the `OpenBLAS` libraries",
" installed on the system")
}
}
rm(useTorch, deviceStr)1 Analytical pipeline
Initialize the COTAN object with the row count table and
the metadata for the experiment.
cond <- "mouse_cortex_E17.5"
obj <- COTAN(raw = sample.dataset)
obj <- initializeMetaDataset(obj,
GEO = GEO,
sequencingMethod = "Drop_seq",
sampleCondition = cond)
#> Initializing `COTAN` meta-data
logThis(paste0("Condition ", getMetadataElement(obj, datasetTags()[["cond"]])),
logLevel = 1L)
#> Condition mouse_cortex_E17.5Before we proceed to the analysis, we need to clean the data. The analysis will use a matrix of raw UMI counts as the input. To obtain this matrix, we have to remove any potential cell doublets or multiplets, as well as any low quality or dying cells.
1.1 Data cleaning
We can check the library size (UMI number) with an empirical cumulative distribution function
plot(ECDPlot(obj))
plot(cellSizePlot(obj))
plot(genesSizePlot(obj))
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_point()`).
plot(scatterPlot(obj))
During the cleaning, every time we want to remove cells or genes
we can use the dropGenesCells()function.
Drop cells with too many reads as they are probably doublets or multiplets
cellsSizeThr <- 6000L
obj <- addElementToMetaDataset(obj, "Cells size threshold", cellsSizeThr)
cells_to_rem <- getCells(obj)[getCellsSize(obj) > cellsSizeThr]
obj <- dropGenesCells(obj, cells = cells_to_rem)
plot(cellSizePlot(obj))
To drop cells by gene number: high genes count might also indicate doublets…
genesSizeHighThr <- 3000L
obj <- addElementToMetaDataset(obj, "Num genes high threshold",
genesSizeHighThr)
cells_to_rem <- getCells(obj)[getNumExpressedGenes(obj) > genesSizeHighThr]
obj <- dropGenesCells(obj, cells = cells_to_rem)
plot(genesSizePlot(obj))
Drop cells with too low genes expression as they are probably dead
genesSizeLowThr <- 500L
obj <- addElementToMetaDataset(obj, "Num genes low threshold", genesSizeLowThr)
cells_to_rem <- getCells(obj)[getNumExpressedGenes(obj) < genesSizeLowThr]
obj <- dropGenesCells(obj, cells = cells_to_rem)
plot(genesSizePlot(obj))
Cells with a too high percentage of mitochondrial genes are likely dead (or at the last problematic) cells. So we drop them!
c(mitPlot, mitSizes) %<-% mitochondrialPercentagePlot(obj, genePrefix = "^Mt")
plot(mitPlot)
mitPercThr <- 1.0
obj <- addElementToMetaDataset(obj, "Mitoc. perc. threshold", mitPercThr)
cells_to_rem <- rownames(mitSizes)[mitSizes[["mit.percentage"]] > mitPercThr]
obj <- dropGenesCells(obj, cells = cells_to_rem)
c(mitPlot, mitSizes) %<-% mitochondrialPercentagePlot(obj, genePrefix = "^Mt")
plot(mitPlot)
If we do not want to consider the mitochondrial genes we can remove them before starting the analysis.
genes_to_rem <- getGenes(obj)[grep("^Mt", getGenes(obj))]
cells_to_rem <- getCells(obj)[which(getCellsSize(obj) == 0L)]
obj <- dropGenesCells(obj, genes_to_rem, cells_to_rem)We want also to log the current status.
logThis(paste("n cells", getNumCells(obj)), logLevel = 1L)
#> n cells 859The clean() function estimates all the parameters for the data. Therefore, we have to run it again every time we remove any genes or cells from the data.
obj <- addElementToMetaDataset(obj, "Num drop B group", 0L)
obj <- clean(obj)
#> Genes/cells selection done: dropped [4787] genes and [0] cells
#> Working on [12235] genes and [859] cells
c(pcaCellsPlot, pcaCellsData, genesPlot,
UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(obj)
#> PCA: START
#> PCA: DONE
#> Hierarchical clustering: START
#> Hierarchical clustering: DONE
plot(pcaCellsPlot)
plot(genesPlot)
We can observe here that the red cells are really enriched in hemoglobin genes so we prefer to remove them. They can be extracted from the pcaCellsData object and removed.
cells_to_rem <- rownames(pcaCellsData)[pcaCellsData[["groups"]] == "B"]
obj <- dropGenesCells(obj, cells = cells_to_rem)
obj <- addElementToMetaDataset(obj, "Num drop B group", 1L)
obj <- clean(obj)
#> Genes/cells selection done: dropped [5] genes and [0] cells
#> Working on [12230] genes and [855] cells
c(pcaCellsPlot, pcaCellsData, genesPlot,
UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(obj)
#> PCA: START
#> PCA: DONE
#> Hierarchical clustering: START
#> Hierarchical clustering: DONE
plot(pcaCellsPlot)
To color the PCA based on nu (so the cells’ efficiency)
plot(UDEPlot)
UDE (color) should not correlate with principal components! This is very important.
The next part is used to remove the cells with efficiency too low.
plot(nuPlot)
We can zoom on the smallest values and, if COTAN detects a clear elbow, we can decide to remove the cells.
plot(zoomedNuPlot)
We also save the defined threshold in the metadata and re-run the estimation
UDELowThr <- 0.30
obj <- addElementToMetaDataset(obj, "Low UDE cells' threshold", UDELowThr)
obj <- addElementToMetaDataset(obj, "Num drop B group", 2L)
obj <- estimateNuLinear(obj)
cells_to_rem <- getCells(obj)[getNu(obj) < UDELowThr]
obj <- dropGenesCells(obj, cells = cells_to_rem)Repeat the estimation after the cells are removed
obj <- clean(obj)
#> Genes/cells selection done: dropped [3] genes and [0] cells
#> Working on [12227] genes and [852] cells
c(pcaCellsPlot, pcaCellsData, genesPlot,
UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(obj)
#> PCA: START
#> PCA: DONE
#> Hierarchical clustering: START
#> Hierarchical clustering: DONE
plot(pcaCellsPlot)
plot(scatterPlot(obj))
logThis(paste("n cells", getNumCells(obj)), logLevel = 1L)
#> n cells 8521.2 COTAN analysis
In this part, all the contingency tables are computed
and used to get the statistics necessary to COEX evaluation and storing
obj <- proceedToCoex(obj, calcCoex = TRUE,
optimizeForSpeed = TRUE, cores = 6L, deviceStr = "cuda",
saveObj = FALSE, outDir = outDir)
#> Cotan analysis functions started
#> Genes/cells selection done: dropped [0] genes and [0] cells
#> Working on [12227] genes and [852] cells
#> Estimate `dispersion`: START
#> Warning in handleMultiCore(cores): On this system multi-core is not currently
#> supported; this can happen on some systems like 'windows'
#> Warning in handleMultiCore(cores): In case you might try
#> 'options(parallelly.fork.enable = TRUE)' to enable multi-core support
#> Warning in handleMultiCore(cores): The number of cores used will be set 1!
#> Estimate `dispersion`: DONE
#> `dispersion` | min: -0.056854248046875 | max: 372.75 | % negative: 19.5469043919195
#> Cotan genes' coex estimation started
#> While trying to load the `torch` library Error in doTryCatch(return(expr), name, parentenv, handler): The `torch` library is installed but the required additional libraries are not avalable yet
#> Warning in canUseTorch(optimizeForSpeed, deviceStr): Falling back to legacy
#> [non-torch] code.
#> Calculate genes' coex (legacy): START
#> Warning in asMethod(object): sparse->dense coercion: allocating vector of size
#> 1.1 GiB
#> Total calculations elapsed time: 184.007210016251
#> Calculate genes' coex (legacy): DONEWhen saveObj == TRUE, in the previous step, this step can be skipped
as the COTAN object has already been saved in the outDir.
# saving the structure
saveRDS(obj, file = file.path(outDir, paste0(cond, ".cotan.RDS")))1.3 Automatic run
It is also possible to run directly a single function if we don’t want to clean anything.
obj2 <- automaticCOTANObjectCreation(
raw = sample.dataset,
GEO = GEO,
sequencingMethod = "Drop_seq",
sampleCondition = cond,
calcCoex = TRUE, cores = 6L, optimizeForSpeed = TRUE,
saveObj = TRUE, outDir = outDir)2 Analysis of the elaborated data
2.1 GDI
To calculate and store the Global Differentiation Index (GDI) we can run:
gdiDF <- calculateGDI(obj)
#> Calculate GDI dataframe: START
#> Calculate GDI dataframe: DONE
head(gdiDF)
#> sum.raw.norm GDI exp.cells
#> 0610007N19Rik 3.332671 1.626640 2.230047
#> 0610007P14Rik 5.279001 1.349530 18.779343
#> 0610009B22Rik 4.733316 1.281530 11.619718
#> 0610009D07Rik 6.313149 1.364269 43.427230
#> 0610009E02Rik 2.590242 1.016320 1.173709
#> 0610009L18Rik 3.509489 1.256757 3.051643
# This will store only the $GDI column
obj <- storeGDI(obj, genesGDI = gdiDF)The next function can either plot the GDI for the dataset directly or
use the pre-computed dataframe.
It marks the given threshold 1.43 (in red) and
the 50% and 75% quantiles (in blue).
We can also specify some gene sets (three in this case) that
we want to label explicitly in the GDI plot.
genesList <- list(
"NPGs" = c("Nes", "Vim", "Sox2", "Sox1", "Notch1", "Hes1", "Hes5", "Pax6"),
"PNGs" = c("Map2", "Tubb3", "Neurod1", "Nefm", "Nefl", "Dcx", "Tbr1"),
"hk" = c("Calm1", "Cox6b1", "Ppia", "Rpl18", "Cox7c", "Erh", "H3f3a",
"Taf1", "Taf2", "Gapdh", "Actb", "Golph3", "Zfr", "Sub1",
"Tars", "Amacr")
)
GDIPlot(obj, cond = cond, genes = genesList, GDIThreshold = 1.40)
#> GDI plot
#> Removed 0 low GDI genes (such as the fully-expressed) in GDI plot
The percentage of cells expressing the gene in the third column of this data-frame is reported.
2.2 Heatmaps
To perform the Gene Pair Analysis, we can plot a heatmap of the COEX values
between two gene sets.
We have to define the different gene sets (list.genes) in a list.
Then we can choose which sets to use in the function parameter sets
(for example, from 1 to 3).
We also have to provide an array of the file name prefixes for each condition
(for example, “mouse_cortex_E17.5”).
In fact, this function can plot genes relationships across many different
conditions to get a complete overview.
plot(heatmapPlot(obj, genesLists = genesList))
#> heatmap plot: START
#> Hangling COTAN object with condition: mouse_cortex_E17.5
#> calculating PValues: START
#> Get p-values on a set of genes on columns and on a set of genes on rows
#> calculating PValues: DONE
#> min coex: -0.475089564305104 max coex: 0.456662313361152
We can also plot a general heatmap of COEX values based on some markers like
the following one.
invisible(
genesHeatmapPlot(obj, primaryMarkers = c("Satb2", "Bcl11b", "Vim", "Hes1"),
pValueThreshold = 0.001, symmetric = TRUE))
#> calculating PValues: START
#> Get p-values genome wide on columns and genome wide on rows
#> calculating PValues: DONE
invisible(
genesHeatmapPlot(obj, primaryMarkers = c("Satb2", "Bcl11b", "Fezf2"),
secondaryMarkers = c("Gabra3", "Meg3", "Cux1", "Neurod6"),
pValueThreshold = 0.001, symmetric = FALSE))2.3 Get data tables
Sometimes we can also be interested in the numbers present directly in the contingency tables for two specific genes. To get them we can use two functions:
contingencyTables() to produce the observed and expected data
print("Contingency Tables:")
#> [1] "Contingency Tables:"
contingencyTables(obj, g1 = "Satb2", g2 = "Bcl11b")
#> $observed
#> Satb2.yes Satb2.no
#> Bcl11b.yes 47 149
#> Bcl11b.no 287 369
#>
#> $expected
#> Satb2.yes Satb2.no
#> Bcl11b.yes 82.78537 113.2154
#> Bcl11b.no 251.21413 404.7851
print("Corresponding Coex")
#> [1] "Corresponding Coex"
getGenesCoex(obj)["Satb2", "Bcl11b"]
#> [1] -0.2028008Another useful function is getGenesCoex(). This can be used to extract
the whole or a partial COEX matrix from a COTAN object.
# For the whole matrix
coex <- getGenesCoex(obj, zeroDiagonal = FALSE)
coex[1000L:1005L, 1000L:1005L]
#> 6 x 6 Matrix of class "dspMatrix"
#> Ap3s1 Ap3s2 Ap4b1 Ap4e1 Ap4m1 Ap4s1
#> Ap3s1 0.91872636 -0.02997150 0.028065790 -0.02865879 0.01372465 -0.032402737
#> Ap3s2 -0.02997150 0.91460404 -0.045830195 -0.03058359 -0.05045147 0.095316128
#> Ap4b1 0.02806579 -0.04583020 0.907981121 -0.03650627 0.02944754 0.009052527
#> Ap4e1 -0.02865879 -0.03058359 -0.036506269 0.82865653 -0.04320115 -0.025651002
#> Ap4m1 0.01372465 -0.05045147 0.029447542 -0.04320115 0.91486865 0.047203683
#> Ap4s1 -0.03240274 0.09531613 0.009052527 -0.02565100 0.04720368 0.908152183# For a partial matrix
coex <- getGenesCoex(obj, genes = c("Satb2", "Bcl11b", "Fezf2"))
coex[1000L:1005L, ]
#> 6 x 3 Matrix of class "dgeMatrix"
#> Bcl11b Fezf2 Satb2
#> Ap3s1 -0.04615274 -0.003499993 0.01372725
#> Ap3s2 0.01281042 0.026728977 0.01784518
#> Ap4b1 0.04062971 0.010720802 -0.03173374
#> Ap4e1 -0.05290965 -0.015074415 -0.03450185
#> Ap4m1 0.05811928 0.020840815 -0.03963970
#> Ap4s1 -0.06858325 0.006087159 -0.015877052.4 Establishing genes’ clusters
COTAN provides a way to establish genes’ clusters given some lists of markers
layersGenes <- list(
"L1" = c("Reln", "Lhx5"),
"L2/3" = c("Satb2", "Cux1"),
"L4" = c("Rorb", "Sox5"),
"L5/6" = c("Bcl11b", "Fezf2"),
"Prog" = c("Vim", "Hes1", "Dummy")
)
c(gSpace, eigPlot, pcaGenesClDF, treePlot) %<-%
establishGenesClusters(obj, groupMarkers = layersGenes,
numGenesPerMarker = 25L, kCuts = 5L)
#> Establishing gene clusters - START
#> Calculating gene co-expression space - START
#> calculating PValues: START
#> Get p-values on a set of genes on columns and genome wide on rows
#> calculating PValues: DONE
#> Number of selected secondary markers: 184
#> Calculating gene co-expression space - DONE
#> Establishing gene clusters - DONE
plot(eigPlot)
plot(treePlot)
colSelection <- vapply(pcaGenesClDF, is.numeric, logical(1L))
genesUmapPl <- UMAPPlot(pcaGenesClDF[, colSelection, drop = FALSE],
clusters = getColumnFromDF(pcaGenesClDF, "hclust"),
elements = layersGenes,
title = "Genes' clusters UMAP Plot",
numNeighbors = 32L, minPointsDist = 0.25)
#> UMAP plot
#> 18:27:35 UMAP embedding parameters a = 1.121 b = 1.057
#> 18:27:35 Read 1223 rows and found 12 numeric columns
#> 18:27:35 Using Annoy for neighbor search, n_neighbors = 32
#> 18:27:35 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 18:27:36 Writing NN index file to temp file E:\biocbuild\bbs-3.21-bioc\tmpdir\RtmpqgOHxN\file170ac2001d52
#> 18:27:36 Searching Annoy index using 1 thread, search_k = 3200
#> 18:27:36 Annoy recall = 100%
#> 18:27:37 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 32
#> 18:27:40 Initializing from normalized Laplacian + noise (using RSpectra)
#> 18:27:40 Commencing optimization for 500 epochs, with 50894 positive edges
#> 18:27:40 Using rng type: pcg
#> 18:27:44 Optimization finished
plot(genesUmapPl)
2.5 Uniform Clustering
It is possible to obtain a cell clusterization based on the concept of
uniformity of expression of the genes across the cells.
That is the cluster satisfies the null hypothesis of the COTAN model:
the genes expression is not dependent on the cell in consideration.
There are two functions involved into obtaining a proper clusterization:
the first is cellsUniformClustering that uses standard tools clusterization
methods, but then discards and re-clusters any non-uniform cluster.
Please note that the most important parameters for the users are the
GDIThresholds inside the Uniform Transcript checkers: they define how
strict is the check. Default constructed advance check gives a pretty strong
guarantee of uniformity for the cluster.
# This code is a little too computationally heavy to be used in an example
# So we stored the result and we can load it in the next section
# default constructed checker is OK
advChecker <- new("AdvancedGDIUniformityCheck")
c(splitClusters, splitCoexDF) %<-%
cellsUniformClustering(obj, initialResolution = 0.8, checker = advChecker,
optimizeForSpeed = TRUE, deviceStr = "cuda",
cores = 6L, genesSel = "HGDI",
saveObj = TRUE, outDir = outDir)
obj <- addClusterization(obj, clName = "split",
clusters = splitClusters, coexDF = splitCoexDF)
table(splitClusters)In the case one already has an existing clusterization, it is possible to
calculate the DEA data.frame and add it to the COTAN object.
data("vignette.split.clusters", package = "COTAN")
splitClusters <- vignette.split.clusters[getCells(obj)]
splitCoexDF <- DEAOnClusters(obj, clusters = splitClusters)
#> Differential Expression Analysis - START
#> **********
#> Differential Expression Analysis - DONE
obj <- addClusterization(obj, clName = "split", clusters = splitClusters,
coexDF = splitCoexDF, override = FALSE)It is possible to have some statistics about the clusterization
c(summaryData, summaryPlot) %<-%
clustersSummaryPlot(obj, clName = "split", plotTitle = "split summary")
summaryData
#> split NoCond CellNumber CellPercentage MeanUDE MedianUDE ExpGenes25 ExpGenes
#> 1 -1 NoCond 82 9.6 1.09 1.10 1688 10626
#> 2 1 NoCond 69 8.1 1.01 0.91 1469 10620
#> 3 2 NoCond 27 3.2 1.74 1.78 2998 9594
#> 4 3 NoCond 153 18.0 0.62 0.59 865 10675
#> 5 4 NoCond 44 5.2 0.86 0.81 1287 9042
#> 6 5 NoCond 196 23.0 0.75 0.76 1133 11053
#> 7 6 NoCond 126 14.8 1.46 1.38 2248 11322
#> 8 7 NoCond 55 6.5 0.45 0.44 488 8826
#> 9 8 NoCond 52 6.1 1.91 1.94 2807 10627
#> 10 9 NoCond 48 5.6 1.20 1.18 1834 9763The ExpGenes column contains the number of genes that are expressed in any
cell of the relevant cluster, while the ExpGenes25 column contains the number of genes expressed in at the least 25% of the cells of the relevant cluster
plot(summaryPlot)
It is possible to visualize how relevant are the marker genes’ lists with
respect to the given clusterization
c(splitHeatmap, scoreDF, pValueDF) %<-%
clustersMarkersHeatmapPlot(obj, groupMarkers = layersGenes,
clName = "split", kCuts = 5L,
adjustmentMethod = "holm")
draw(splitHeatmap)
Since the algorithm that creates the clusters is not directly geared to achieve cluster uniformity, there might be some clusters that can be merged back together and still be uniform.
This is the purpose of the function mergeUniformCellsClusters that takes a
clusterization and tries to merge cluster pairs after checking that together
the pair forms a uniform cluster.
In order to avoid running the totality of the possible checks (as it can explode quickly with the number of clusters), the function relies on a related distance the find the cluster pairs that have the highest chance to be merged.
c(mergedClusters, mergedCoexDF) %<-%
mergeUniformCellsClusters(obj, clusters = splitClusters,
checkers = advChecker,
optimizeForSpeed = TRUE, deviceStr = "cuda",
cores = 6L, saveObj = TRUE, outDir = outDir)
# merges are:
# 1 <- 06 + 07
# 2 <- '-1' + 08 + 11
# 3 <- 09
# 4 <- 01
# 5 <- 02
# 6 <- 12
# 7 <- 03 + 04 + 10 + 13
# 8 <- 05
obj <- addClusterization(obj, clName = "merge", override = TRUE,
clusters = mergedClusters, coexDF = mergedCoexDF)
table(mergedClusters)Again, here, the most important parameters for the users are the GDIThresholds
inside the Uniform Transcript checkers: they define how strict is the
check. Default constructed advance check gives a pretty strong guarantee of
uniformity for the cluster. If one wants to achieve a more relaxed merge,
it is possible to define a sequence of checkers, each less stringent than the
previous, that will be applied sequentially: First all the merges with the
stricter checker are performed, than those that satisfy the second, and so on.
checkersList <- list(advChecker,
shiftCheckerThresholds(advChecker, 0.01),
shiftCheckerThresholds(advChecker, 0.03))
prevCheckRes <- data.frame()
# In this case we want to re-use the already calculated merge checks
# Se we reload them from the output files. This, along a similar facility for
# the split method, is also useful in the cases the execution is interrupted
# prematurely...
#
if (TRUE) {
# read from the last file among those named all_check_results_XX.csv
mergeDir <- file.path(outDir, cond, "leafs_merge")
resFiles <- list.files(path = mergeDir, pattern = "all_check_results_.*csv",
full.names = TRUE)
prevCheckRes <- read.csv(resFiles[length(resFiles)],
header = TRUE, row.names = 1L)
}
c(mergedClusters2, mergedCoexDF2) %<-%
mergeUniformCellsClusters(obj, clusters = splitClusters,
checkers = checkersList,
allCheckResults = prevCheckRes,
optimizeForSpeed = TRUE, deviceStr = "cuda",
cores = 6L, saveObj = TRUE, outDir = outDir)
# merges are:
# 1 <- '-1' + 06 + 09
# 2 <- 07
# 3 <- 03 + 04 + 10 + 13
# 4 <- 12
# 5 <- 01 + 08 + 11
# 6 <- 02 + 05
obj <- addClusterization(obj, clName = "merge2", override = TRUE,
clusters = mergedClusters2, coexDF = mergedCoexDF2)
table(mergedClusters2)In the case one already has existing clusterizations, it is possible to
calculate their DEA data.frame and add them to the COTAN object.
data("vignette.merge.clusters", package = "COTAN")
mergedClusters <- vignette.merge.clusters[getCells(obj)]
mergedCoexDF <- DEAOnClusters(obj, clusters = mergedClusters)
#> Differential Expression Analysis - START
#> *******
#> Differential Expression Analysis - DONE
obj <- addClusterization(obj, clName = "merge", clusters = mergedClusters,
coexDF = mergedCoexDF, override = FALSE)
data("vignette.merge2.clusters", package = "COTAN")
mergedClusters2 <- vignette.merge2.clusters[getCells(obj)]
mergedCoexDF2 <- DEAOnClusters(obj, clusters = mergedClusters2)
#> Differential Expression Analysis - START
#> ****
#> Differential Expression Analysis - DONE
obj <- addClusterization(obj, clName = "merge2", clusters = mergedClusters2,
coexDF = mergedCoexDF2, override = FALSE)
table(mergedClusters2, mergedClusters)
#> mergedClusters
#> mergedClusters2 1 2 3 4 5 6 7
#> 1 151 0 0 0 0 0 0
#> 2 0 153 44 0 196 0 0
#> 3 0 0 0 153 0 0 100
#> 4 0 0 0 0 0 55 0Here is possible to visualize how the split clusters are merged in to
merge and merge2
c(mergeHeatmap, mergeScoresDF, mergePValuesDF) %<-%
clustersMarkersHeatmapPlot(obj, clName = "merge", condNameList = "split",
conditionsList = list(splitClusters))
draw(mergeHeatmap)
c(merge2Heatmap, merge2ScoresDF, merge2PValuesDF) %<-%
clustersMarkersHeatmapPlot(obj, clName = "merge2", condNameList = "split",
conditionsList = list(splitClusters))
draw(merge2Heatmap)
In the following graph, it is possible to check that the found clusters align very well to the expression of the layers’ genes defined above…
It is possible to plot the clusterization on a UMAP plot
c(umapPlot, cellsPCA) %<-% cellsUMAPPlot(obj, clName = "merge2",
clusters = NULL,
dataMethod = "LogLikelihood",
genesSel = "HGDI",
numGenes = 200L,
colors = NULL,
numNeighbors = 30L,
minPointsDist = 0.3)
#> Running genes' selection: START
#> Running genes' selection: DONE
#> UMAP plot
#> 18:28:18 UMAP embedding parameters a = 0.9922 b = 1.112
#> 18:28:18 Read 852 rows and found 25 numeric columns
#> 18:28:18 Using Annoy for neighbor search, n_neighbors = 30
#> 18:28:18 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 18:28:18 Writing NN index file to temp file E:\biocbuild\bbs-3.21-bioc\tmpdir\RtmpqgOHxN\file170ac1e84535
#> 18:28:18 Searching Annoy index using 1 thread, search_k = 3000
#> 18:28:18 Annoy recall = 100%
#> 18:28:19 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
#> 18:28:21 Initializing from normalized Laplacian + noise (using RSpectra)
#> 18:28:21 Commencing optimization for 500 epochs, with 36310 positive edges
#> 18:28:21 Using rng type: pcg
#> 18:28:25 Optimization finished
plot(umapPlot)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
2.6 Vignette clean-up stage
The next few lines are just to clean.
if (file.exists(file.path(outDir, paste0(cond, ".cotan.RDS")))) {
#Delete file if it exists
file.remove(file.path(outDir, paste0(cond, ".cotan.RDS")))
}
if (file.exists(file.path(outDir, paste0(cond, "_times.csv")))) {
#Delete file if it exists
file.remove(file.path(outDir, paste0(cond, "_times.csv")))
}
if (dir.exists(file.path(outDir, cond))) {
unlink(file.path(outDir, cond), recursive = TRUE)
}
if (dir.exists(file.path(outDir, GEO))) {
unlink(file.path(outDir, GEO), recursive = TRUE)
}
# stop logging to file
setLoggingFile("")
#> Closing previous log file - Setting log file to be:
file.remove(file.path(outDir, "vignette_v2.log"))
#> [1] TRUE
options(parallelly.fork.enable = FALSE)Sys.time()
#> [1] "2025-06-09 18:28:28 EDT"sessionInfo()
#> R version 4.5.0 RC (2025-04-04 r88126 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows Server 2022 x64 (build 20348)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] grid stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ComplexHeatmap_2.24.0 GEOquery_2.76.0 Biobase_2.68.0
#> [4] BiocGenerics_0.54.0 generics_0.1.4 qpdf_1.3.5
#> [7] Rtsne_0.17 data.table_1.17.4 rlang_1.1.6
#> [10] zeallot_0.2.0 COTAN_2.8.2 BiocStyle_2.36.0
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.22 splines_4.5.0
#> [3] later_1.4.2 tibble_3.3.0
#> [5] polyclip_1.10-7 XML_3.99-0.18
#> [7] fastDummies_1.7.5 lifecycle_1.0.4
#> [9] doParallel_1.0.17 processx_3.8.6
#> [11] globals_0.18.0 lattice_0.22-7
#> [13] MASS_7.3-65 dendextend_1.19.0
#> [15] magrittr_2.0.3 limma_3.64.1
#> [17] plotly_4.10.4 sass_0.4.10
#> [19] rmarkdown_2.29 jquerylib_0.1.4
#> [21] yaml_2.3.10 httpuv_1.6.16
#> [23] Seurat_5.3.0 sctransform_0.4.2
#> [25] askpass_1.2.1 spam_2.11-1
#> [27] sp_2.2-0 spatstat.sparse_3.1-0
#> [29] reticulate_1.42.0 cowplot_1.1.3
#> [31] pbapply_1.7-2 RColorBrewer_1.1-3
#> [33] abind_1.4-8 GenomicRanges_1.60.0
#> [35] purrr_1.0.4 coro_1.1.0
#> [37] torch_0.14.2 circlize_0.4.16
#> [39] GenomeInfoDbData_1.2.14 IRanges_2.42.0
#> [41] S4Vectors_0.46.0 ggrepel_0.9.6
#> [43] irlba_2.3.5.1 listenv_0.9.1
#> [45] spatstat.utils_3.1-4 rentrez_1.2.3
#> [47] goftest_1.2-3 RSpectra_0.16-2
#> [49] spatstat.random_3.4-1 fitdistrplus_1.2-2
#> [51] parallelly_1.45.0 codetools_0.2-20
#> [53] DelayedArray_0.34.1 xml2_1.3.8
#> [55] tidyselect_1.2.1 shape_1.4.6.1
#> [57] UCSC.utils_1.4.0 farver_2.1.2
#> [59] ScaledMatrix_1.16.0 viridis_0.6.5
#> [61] matrixStats_1.5.0 stats4_4.5.0
#> [63] spatstat.explore_3.4-3 jsonlite_2.0.0
#> [65] GetoptLong_1.0.5 gghalves_0.1.4
#> [67] progressr_0.15.1 ggridges_0.5.6
#> [69] survival_3.8-3 iterators_1.0.14
#> [71] foreach_1.5.2 tools_4.5.0
#> [73] ica_1.0-3 Rcpp_1.0.14
#> [75] glue_1.8.0 gridExtra_2.3
#> [77] SparseArray_1.8.0 xfun_0.52
#> [79] MatrixGenerics_1.20.0 ggthemes_5.1.0
#> [81] GenomeInfoDb_1.44.0 dplyr_1.1.4
#> [83] withr_3.0.2 BiocManager_1.30.26
#> [85] fastmap_1.2.0 callr_3.7.6
#> [87] digest_0.6.37 rsvd_1.0.5
#> [89] parallelDist_0.2.6 R6_2.6.1
#> [91] mime_0.13 colorspace_2.1-1
#> [93] Cairo_1.6-2 scattermore_1.2
#> [95] tensor_1.5 dichromat_2.0-0.1
#> [97] spatstat.data_3.1-6 tidyr_1.3.1
#> [99] httr_1.4.7 htmlwidgets_1.6.4
#> [101] S4Arrays_1.8.1 uwot_0.2.3
#> [103] pkgconfig_2.0.3 gtable_0.3.6
#> [105] lmtest_0.9-40 SingleCellExperiment_1.30.1
#> [107] XVector_0.48.0 htmltools_0.5.8.1
#> [109] dotCall64_1.2 bookdown_0.43
#> [111] zigg_0.0.2 clue_0.3-66
#> [113] SeuratObject_5.1.0 scales_1.4.0
#> [115] png_0.1-8 spatstat.univar_3.1-3
#> [117] knitr_1.50 tzdb_0.5.0
#> [119] reshape2_1.4.4 rjson_0.2.23
#> [121] curl_6.3.0 nlme_3.1-168
#> [123] proxy_0.4-27 cachem_1.1.0
#> [125] zoo_1.8-14 GlobalOptions_0.1.2
#> [127] stringr_1.5.1 KernSmooth_2.23-26
#> [129] parallel_4.5.0 miniUI_0.1.2
#> [131] pillar_1.10.2 vctrs_0.6.5
#> [133] RANN_2.6.2 promises_1.3.3
#> [135] BiocSingular_1.24.0 beachmat_2.24.0
#> [137] xtable_1.8-4 cluster_2.1.8.1
#> [139] evaluate_1.0.3 magick_2.8.7
#> [141] tinytex_0.57 readr_2.1.5
#> [143] cli_3.6.5 compiler_4.5.0
#> [145] crayon_1.5.3 future.apply_1.20.0
#> [147] labeling_0.4.3 ps_1.9.1
#> [149] plyr_1.8.9 stringi_1.8.7
#> [151] viridisLite_0.4.2 deldir_2.0-4
#> [153] BiocParallel_1.42.1 assertthat_0.2.1
#> [155] lazyeval_0.2.2 spatstat.geom_3.4-1
#> [157] Matrix_1.7-3 RcppHNSW_0.6.0
#> [159] hms_1.1.3 patchwork_1.3.0
#> [161] bit64_4.6.0-1 future_1.58.0
#> [163] ggplot2_3.5.2 statmod_1.5.0
#> [165] shiny_1.10.0 SummarizedExperiment_1.38.1
#> [167] ROCR_1.0-11 Rfast_2.1.5.1
#> [169] igraph_2.1.4 RcppParallel_5.1.10
#> [171] bslib_0.9.0 bit_4.6.0