Guided tutorial to genes’ clustering using COTAN
Preamble
Introduction
This tutorial shows the available tools in COTAN to help with genes’ clustering
Setup
Define a directory where the data and the output will be stored
dataDir <- file.path(tempdir(), "COTAN_vignette_data")
dir.create(dataDir, recursive = TRUE, showWarnings = FALSE)
print(dataDir)
#> [1] "/tmp/RtmpVAoHEq/COTAN_vignette_data"
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_genes_clustering.log"))
#> Setting log file to be: /tmp/RtmpVAoHEq/COTAN_vignette_data/vignette_genes_clustering.logRetrieving the data-set
Download the data-set for mouse cortex E17.5 from Yuzwa et al. (2017)
GEO <- "GSM2861514"
dir.create(file.path(dataDir, GEO), showWarnings = FALSE)
datasetFileName <-
file.path(dataDir, GEO, "GSM2861514_E175_All_Cells_DGE.txt.gz")
# retries up to 5 times to get the dataset
attempts <- 0L
maxAttempts <- 5L
ok <- FALSE
while (attempts < maxAttempts && !ok) {
attempts <- attempts + 1L
if (!file.exists(datasetFileName)) {
res <- try(
GEOquery::getGEOSuppFiles(
GEO = GEO,
makeDirectory = TRUE,
baseDir = dataDir,
filter_regex = base::basename(datasetFileName),
fetch_files = TRUE
),
silent = TRUE
)
}
ok <- file.exists(datasetFileName)
if (!ok && attempts < maxAttempts) {
Sys.sleep(1)
}
}
assertthat::assert_that(
ok,
msg = paste0(
"Failed to retrieve file '", datasetFileName,
"' after ", maxAttempts, " attempts."
)
)
#> [1] TRUE
rawDataset <- read.csv(datasetFileName, sep = "\t", row.names = 1L)
print(dim(rawDataset))
#> [1] 17085 2000Create the COTAN object
Initialize the COTAN object with the row count table and
the metadata from the experiment.
cond <- "mouse_cortex_E17.5"
obj <- COTAN(raw = rawDataset)
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.5Assign to each cell its origin
# mark cells origin
data(vignette.cells.origin)
head(vignette.cells.origin, 18)
#> GCGAATTGTGAA ACCGATAACTGA CCCCGGGTGCGA CTGTAGATGTTA TCATCGAAGCGC TTCTACCGAGTC
#> Other Other Other Other Other Other
#> CTGTTCCCGGCG GCGTGTTAGTTC CTCGCGCGTTTA GATGTATAACTT GCGCTATGATTT CGTTTAGTTTAC
#> Other Other Other Other Other Cortical
#> GTGGAGGCCCAT TGTCACTACATC TCTAGAACAACG ACCTTTGTTCGT TTGTCTTCTTCG TAAAATATCGCC
#> Other Other Cortical Cortical Cortical Cortical
#> Levels: Cortical Other
obj <-
addCondition(obj, condName = "origin", conditions = vignette.cells.origin)
rm(vignette.cells.origin)Drop low quality cells
# use previously established results to determine
# which cells were dropped in the cleaning stage
data(vignette.split.clusters)
cellsToDrop <-
getCells(obj)[!(getCells(obj) %in% names(vignette.split.clusters))]
obj <- dropGenesCells(obj, cells = cellsToDrop)
# Log the remaining number of cells
logThis(paste("n cells", getNumCells(obj)), logLevel = 1L)
#> n cells 1783
table(getCondition(obj, condName = "origin"))
#>
#> Cortical Other
#> 844 939
rm(vignette.split.clusters)Establish genes’ clusters
COTAN analysis
In this part, the COTAN model is calibrated: 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 = 3L,
deviceStr = "cuda"
)
#> COTAN dataset analysis: START
#> Genes/cells selection done: dropped [4330] genes and [0] cells
#> Working on [12755] genes and [1783] cells
#> Estimate `dispersion`: START
#> Total calculations elapsed time: 8.91566443443298
#> Estimate `dispersion`: DONE
#> Estimate `dispersion`: DONE
#> `dispersion` | min: -0.0349260648055935 | max: 760.574661197634 | % negative: 12.5362602900823
#> COTAN genes' COEX estimation: START
#> 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 available 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 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.2 GiB
#> Total calculations elapsed time: 102.62933754921
#> Calculate genes' COEX (legacy): DONE
#> COTAN dataset analysis: DONEGet data tables and stored COEX
Sometimes we can also be interested in the numbers present directly in the contingency tables for two specific genes.
cat("Contingency Tables:")
#> Contingency Tables:
print(contingencyTables(obj, g1 = "Satb2", g2 = "Bcl11b"))
#> $observed
#> Satb2.yes Satb2.no
#> Bcl11b.yes 62 565
#> Bcl11b.no 308 848
#>
#> $expected
#> Satb2.yes Satb2.no
#> Bcl11b.yes 138.583 488.417
#> Bcl11b.no 231.417 924.583
cat(paste("Corresponding COEX value:",
getGenesCoex(obj, zeroDiagonal = FALSE)["Satb2", "Bcl11b"]))
#> Corresponding COEX value: -0.219641253268049Another useful function is getGenesCoex(). This can be
used to extract the whole or a partial COEX matrix from a
COTAN object.
# Extract the whole matrix
coex <- getGenesCoex(obj, zeroDiagonal = FALSE)
coex[1000L:1005L, 1000L:1005L]
#> 6 x 6 Matrix of class "dspMatrix"
#> Ankrd6 Ankrd61 Anks1 Anks1b Anks3
#> Ankrd6 0.919828628 0.006022968 0.039727285 0.019217926 0.022676686
#> Ankrd61 0.006022968 0.304643397 0.012120271 0.012826424 0.002818821
#> Anks1 0.039727285 0.012120271 0.907300176 -0.011094145 -0.006198863
#> Anks1b 0.019217926 0.012826424 -0.011094145 0.919554726 -0.017552746
#> Anks3 0.022676686 0.002818821 -0.006198863 -0.017552746 0.905659149
#> Anks6 0.014713461 -0.004394969 -0.021761229 -0.007003361 0.029217368
#> Anks6
#> Ankrd6 0.014713461
#> Ankrd61 -0.004394969
#> Anks1 -0.021761229
#> Anks1b -0.007003361
#> Anks3 0.029217368
#> Anks6 0.513438755# For a partial matrix
coex_sm <- getGenesCoex(obj, genes = c("Satb2", "Bcl11b", "Fezf2"),
zeroDiagonal = FALSE)
coex_sm[2001:2010, ]
#> 10 x 3 Matrix of class "dgeMatrix"
#> Bcl11b Fezf2 Satb2
#> Ccdc91 -0.007340495 -0.015017365 0.0081678484
#> Ccdc92 0.006863524 -0.002183516 0.0121977478
#> Ccdc93 -0.009075446 0.032998076 -0.0185992406
#> Ccdc94 -0.025712653 -0.013185117 0.0153650072
#> Ccdc96 -0.028979603 -0.015681037 0.0484463876
#> Ccdc97 -0.029885047 -0.013649152 -0.0390277387
#> Cchcr1 -0.015829490 0.019683483 0.0005360145
#> Cck -0.040912312 0.049164448 0.0376753146
#> Ccl12 -0.030998730 -0.011044525 -0.0124730051
#> Ccl25 -0.037731486 -0.006071913 0.0028277194Analysis of the elaborated data: GDI
To calculate and store the Global Differentiation Index
(GDI) use the following:
gdiDF <- calculateGDI(obj, cores = 3L)
#> Calculate GDI dataframe: START
#> Calculate `GDI`: START
#> Total calculations elapsed time: 5.72920107841492
#> Calculate `GDI`: DONE
#> Total calculations elapsed time: 10.9503879547119
#> Calculate GDI dataframe: DONE
# Show the most *differentially expressed* genes
perm <- order(gdiDF[["GDI"]], decreasing = TRUE)
head(gdiDF[perm, ], n = 10)
#> sum.raw.norm GDI exp.cells
#> Mki67 2.400195 2.810177 6.393718
#> Hmgb2 2.669711 2.789769 10.095345
#> 2810417H13Rik 2.373125 2.773485 6.505889
#> Pbk 2.189592 2.716527 5.496354
#> Prc1 2.357143 2.699208 5.384184
#> Dbi 3.145560 2.687651 23.331464
#> Ccnd2 3.220094 2.687415 21.817162
#> Vim 2.711821 2.670697 14.638250
#> Mfge8 2.370640 2.664075 5.384184
#> Cdca8 2.260013 2.658388 6.169377
# This will store only the $GDI column
obj <- storeGDI(obj, genesGDI = gdiDF)The other two columns return the sum of the normalized data rows and the percentage of cells expressing the gene respectively.
Define genes lists
genesLists <- list(
# Neural Progenitor Genes
"NPGs" = c("Nes", "Vim", "Sox2", "Sox1", "Notch1", "Hes1", "Hes5", "Pax6"),
# Pan Neural Genes
"PNGs" = c("Map2", "Tubb3", "Neurod1", "Nefm", "Nefl", "Dcx", "Tbr1"),
# House Keeping
"hk" = c("Calm1", "Cox6b1", "Ppia", "Rpl18", "Cox7c", "Erh", "H3f3a",
"Taf1", "Taf2", "Gapdh", "Actb", "Golph3", "Zfr", "Sub1",
"Tars", "Amacr")
)GDI Plot
The next function can either plot the GDI from the
dataset directly or use the already stored data.
It marks the given threshold 1.40 (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.
GDIPlot(obj, genes = genesLists, GDIThreshold = 1.40)
#> GDI plot
#> Removed 0 low GDI genes (such as the fully-expressed) in GDI plot
In the picture above it is possible to see how, as expected, the
house-keeping genes (hk) have generally a low
GDI value.
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. The first set will be used to define the columns; then we can
choose which sets to use in the function parameter sets (for example,
from 1 to 3).
plot(heatmapPlot(obj, genesLists = genesLists, sets = 1L:3L, cores = 3L))
#> 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
#> Total calculations elapsed time: 2.46536231040955
#> calculating PValues: DONE
#> min COEX: -0.335360548131135 max COEX: 0.329177031378469
#> Total calculations elapsed time: 2.5075147151947
#> Heatmap plot: DONE
layersGenes <- list(
"L1" = c("Reln", "Lhx5"),
"L2/3" = c("Cux1", "Satb2"),
"L4" = c("Rorb", "Sox5"),
"L5/6" = c("Bcl11b", "Fezf2"),
"Prog" = c("Hes1", "Vim", "Dummy")
)We can also plot a general heatmap of COEX values based
on some markers like the following one.
layersMarkers <- c("Satb2", "Bcl11b", "Vim", "Hes1")
genesHeatmapPlot(
obj,
primaryMarkers = layersMarkers,
pValueThreshold = 0.0001,
symmetric = TRUE,
cores = 3L
)
#> calculating PValues: START
#> Get p-values genome wide on columns and genome wide on rows
#> Total calculations elapsed time: 24.8557078838348
#> calculating PValues: DONE
Establishing genes’ clusters
COTAN provides a way to establish genes’ clusters given
some lists of markers
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
#> Total calculations elapsed time: 3.46485352516174
#> calculating PValues: DONE
#> Number of selected secondary markers: 190
#> Calculate `GDI`: START
#> Total calculations elapsed time: 0.159004211425781
#> Calculate `GDI`: DONE
#> Total calculations elapsed time: 7.24500346183777
#> Calculating gene co-expression space - DONE
#> Total calculations elapsed time: 7.64452624320984
#> Establishing gene clusters - DONE
plot(eigPlot)
UMAP plot
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: START
#> 07:31:33 UMAP embedding parameters a = 1.121 b = 1.057
#> 07:31:33 Read 1276 rows and found 12 numeric columns
#> 07:31:33 Using Annoy for neighbor search, n_neighbors = 32
#> 07:31:33 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 07:31:33 Writing NN index file to temp file /tmp/RtmpVAoHEq/file28504e48bb28
#> 07:31:33 Searching Annoy index using 1 thread, search_k = 3200
#> 07:31:33 Annoy recall = 100%
#> 07:31:34 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 32
#> 07:31:36 Initializing from normalized Laplacian + noise (using RSpectra)
#> 07:31:36 Commencing optimization for 500 epochs, with 55474 positive edges
#> 07:31:36 Using rng type: pcg
#> 07:31:38 Optimization finished
#> Total calculations elapsed time: 4.87094330787659
#> UMAP plot: DONE
plot(genesUmapPl)
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_genes_clustering.log"))
#> [1] TRUE
options(prevOptState)Sys.time()
#> [1] "2026-10-17 07:31:38 UTC"
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] zeallot_0.2.0 COTAN_2.12.1 BiocStyle_2.40.0
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.23 splines_4.6.1
#> [3] later_1.4.8 tibble_3.3.1
#> [5] polyclip_1.10-7 XML_3.99-0.25
#> [7] fastDummies_1.7.6 httr2_1.3.0
#> [9] lifecycle_1.0.5 doParallel_1.0.17
#> [11] processx_3.9.0 globals_0.20.0
#> [13] lattice_0.23-1 MASS_7.3-66
#> [15] ggdist_3.3.3 dendextend_1.19.1
#> [17] magrittr_2.0.5 limma_3.68.5
#> [19] plotly_4.12.1 sass_0.4.10
#> [21] rmarkdown_2.32 jquerylib_0.1.4
#> [23] yaml_2.3.12 httpuv_1.6.17
#> [25] otel_0.2.0 Seurat_5.6.0
#> [27] sctransform_0.4.3 spam_2.11-4
#> [29] sp_2.2-4 spatstat.sparse_3.2-0
#> [31] reticulate_1.47.0 cowplot_1.2.0
#> [33] pbapply_1.7-5 buildtools_1.0.0
#> [35] RColorBrewer_1.1-3 abind_1.4-8
#> [37] Rtsne_0.17 GenomicRanges_1.64.0
#> [39] purrr_1.2.2 BiocGenerics_0.58.1
#> [41] coro_1.1.0 torch_0.17.0
#> [43] circlize_0.4.18 IRanges_2.46.0
#> [45] S4Vectors_0.50.3 ggrepel_0.9.8
#> [47] irlba_2.4.1 listenv_1.1.0
#> [49] spatstat.utils_3.2-5 rentrez_1.2.4
#> [51] maketools_1.3.2 goftest_1.2-3
#> [53] RSpectra_0.16-2 spatstat.random_3.5-2
#> [55] fitdistrplus_1.2-6 parallelly_1.48.0
#> [57] codetools_0.2-20 DelayedArray_0.38.2
#> [59] xml2_1.6.0 tidyselect_1.2.1
#> [61] shape_1.4.6.1 farver_2.1.2
#> [63] ScaledMatrix_1.20.0 viridis_0.6.5
#> [65] matrixStats_1.5.0 stats4_4.6.1
#> [67] spatstat.explore_3.8-3 Seqinfo_1.2.0
#> [69] jsonlite_2.0.0 GetoptLong_1.1.1
#> [71] progressr_1.0.0 ggridges_0.5.7
#> [73] survival_3.8-12 iterators_1.0.14
#> [75] systemfonts_1.3.2 foreach_1.5.2
#> [77] tools_4.6.1 ragg_1.5.2
#> [79] ica_1.0-3 Rcpp_1.1.2
#> [81] glue_1.8.1 gridExtra_2.3.1
#> [83] SparseArray_1.12.3 xfun_0.62
#> [85] distributional_0.9.0 MatrixGenerics_1.24.0
#> [87] ggthemes_7.0.0 dplyr_1.2.1
#> [89] withr_3.0.3 BiocManager_1.30.27
#> [91] fastmap_1.2.0 callr_3.8.0
#> [93] digest_0.6.39 rsvd_1.0.5
#> [95] parallelDist_0.2.7 R6_2.6.1
#> [97] mime_0.13 textshaping_1.0.5
#> [99] colorspace_2.1-4 Cairo_1.7-0
#> [101] scattermore_1.2 tensor_1.5.1
#> [103] spatstat.data_3.1-9 tidyr_1.3.2
#> [105] generics_0.1.4 data.table_1.18.6.1
#> [107] httr_1.4.9 htmlwidgets_1.6.4
#> [109] S4Arrays_1.12.1 uwot_0.2.5
#> [111] pkgconfig_2.0.3 gtable_0.3.6
#> [113] ComplexHeatmap_2.28.0 lmtest_0.9-40
#> [115] S7_0.2.2 SingleCellExperiment_1.34.0
#> [117] XVector_0.52.0 sys_3.4.3
#> [119] htmltools_0.5.9 dotCall64_1.2
#> [121] zigg_0.0.2 clue_0.3-68
#> [123] SeuratObject_5.4.0 scales_1.4.0
#> [125] Biobase_2.72.0 png_0.1-9
#> [127] spatstat.univar_3.2-0 knitr_1.52
#> [129] tzdb_0.5.0 reshape2_1.4.5
#> [131] rjson_0.2.23 curl_8.1.0
#> [133] nlme_3.1-171 proxy_0.4-29
#> [135] cachem_1.1.0 zoo_1.9-1
#> [137] GlobalOptions_0.1.4 stringr_1.6.0
#> [139] KernSmooth_2.23-27 parallel_4.6.1
#> [141] miniUI_0.1.2 GEOquery_2.80.0
#> [143] pillar_1.11.1 grid_4.6.1
#> [145] vctrs_0.7.3 RANN_2.6.3
#> [147] promises_1.5.0 BiocSingular_1.28.1
#> [149] beachmat_2.28.0 xtable_1.8-8
#> [151] cluster_2.1.8.3 evaluate_1.0.5
#> [153] readr_2.2.0 cli_3.6.6
#> [155] compiler_4.6.1 rlang_1.3.0
#> [157] crayon_1.5.3 future.apply_1.21.0
#> [159] labeling_0.4.3 ps_1.9.3
#> [161] plyr_1.8.9 stringi_1.8.9
#> [163] viridisLite_0.4.3 deldir_2.0-4
#> [165] BiocParallel_1.46.0 assertthat_0.2.1
#> [167] spatstat.geom_3.8-3 Matrix_1.7-6
#> [169] RcppHNSW_0.7.0 hms_1.1.4
#> [171] patchwork_1.3.2 bit64_4.8.6
#> [173] future_1.80.0 conflicted_1.2.0
#> [175] ggplot2_4.0.3 statmod_1.5.2
#> [177] shiny_1.14.0 SummarizedExperiment_1.42.0
#> [179] ROCR_1.0-12 Rfast_2.1.5.2
#> [181] memoise_2.0.1 igraph_2.3.4
#> [183] RcppParallel_6.2.1 bslib_0.12.0
#> [185] bit_4.6.0