Contents

library(COTAN)
library(data.table)
library(Matrix)
library(ggrepel)
#> Loading required package: ggplot2
library(factoextra)
#> Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(Rtsne)
library(utils)
library(plotly)
#> 
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#> 
#>     last_plot
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following object is masked from 'package:graphics':
#> 
#>     layout
library(tidyverse)
#> ── Attaching packages
#> ───────────────────────────────────────
#> tidyverse 1.3.2 ──
#> ✔ tibble  3.1.8      ✔ dplyr   1.0.10
#> ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
#> ✔ readr   2.1.3      ✔ forcats 0.5.2 
#> ✔ purrr   0.3.5      
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::between()   masks data.table::between()
#> ✖ tidyr::expand()    masks Matrix::expand()
#> ✖ dplyr::filter()    masks plotly::filter(), stats::filter()
#> ✖ dplyr::first()     masks data.table::first()
#> ✖ dplyr::lag()       masks stats::lag()
#> ✖ dplyr::last()      masks data.table::last()
#> ✖ tidyr::pack()      masks Matrix::pack()
#> ✖ purrr::transpose() masks data.table::transpose()
#> ✖ tidyr::unpack()    masks Matrix::unpack()
library(htmlwidgets)
library(MASS)
#> 
#> Attaching package: 'MASS'
#> 
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> 
#> The following object is masked from 'package:plotly':
#> 
#>     select
library(dendextend)
#> 
#> ---------------------
#> Welcome to dendextend version 1.16.0
#> Type citation('dendextend') for how to cite the package.
#> 
#> Type browseVignettes(package = 'dendextend') for the package vignette.
#> The github page is: https://github.com/talgalili/dendextend/
#> 
#> Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
#> You may ask questions at stackoverflow, use the r and dendextend tags: 
#>   https://stackoverflow.com/questions/tagged/dendextend
#> 
#>  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
#> ---------------------
#> 
#> 
#> Attaching package: 'dendextend'
#> 
#> The following object is masked from 'package:data.table':
#> 
#>     set
#> 
#> The following object is masked from 'package:stats':
#> 
#>     cutree
mycolours <- c("A" = "#8491B4B2","B"="#E64B35FF")
my_theme = theme(axis.text.x = element_text(size = 14, angle = 0, hjust = .5, vjust = .5, face = "plain", colour ="#3C5488FF" ),
                 axis.text.y = element_text( size = 14, angle = 0, hjust = 0, vjust = .5, face = "plain", colour ="#3C5488FF"),
                 axis.title.x = element_text( size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain", colour ="#3C5488FF"),
                 axis.title.y = element_text( size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain", colour ="#3C5488FF"))
data_dir = tempdir()

This is just an example on a toy dataset. If the user want to try on a real dataset it can be used the previous command to download it and the next commented line to use it as input data. There are also available online many other examples at link.

data("ERCCraw", package = "COTAN")

ERCCraw = as.data.frame(ERCCraw)
rownames(ERCCraw) = ERCCraw$V1
ERCCraw = ERCCraw[,2:ncol(ERCCraw)]
ERCCraw[1:5,1:5]
#>            AAACCGTGAAGCCT.1 AAACGCACTGCTGA.1 AAACTTGACCGAAT.1 AAAGACGAGGAAGC.1
#> ERCC-00002             1662             4036             3919             4124
#> ERCC-00003               95              283              298              290
#> ERCC-00004               25               75               70               75
#> ERCC-00009               41               57               78               98
#> ERCC-00012                0                0                0                0
#>            AAAGACGATCCGAA.1
#> ERCC-00002             4220
#> ERCC-00003              312
#> ERCC-00004               87
#> ERCC-00009               87
#> ERCC-00012                0

Define a directory where the output will be stored.

out_dir <- paste0(tempdir(),"/")

1 Analytical pipeline

Inizialise the COTAN object with the row count table and the metadata for the experiment.

obj = new("scCOTAN",raw = ERCCraw)
#obj = initRaw(obj,GEO="GSM2861514" ,sc.method="Drop_seq",cond = "mouse cortex E17.5")
obj = initRaw(obj,GEO="ERCC" ,sc.method="10X",cond = "negative ERCC dataset")
#> [1] "Initializing S4 object"

Now we can start the cleaning. Analysis requires and starts from a matrix of raw UMI counts after removing possible cell doublets or multiplets and low quality or dying cells (with too high mtRNA percentage, easily done with Seurat or other tools).

If we do not want to consider the mitochondrial genes we can remove them before starting the analysis.

genes_to_rem = get.genes(obj)[grep('^mt', get.genes(obj))] 
cells_to_rem = names(get.cell.size(obj)[which(get.cell.size(obj) == 0)])
obj = drop.genes.cells(obj,genes_to_rem,cells_to_rem )

We want also to define a prefix to identify the sample.

#t = "E17.5_cortex"
t = "ERCC"

print(paste("Condition ",t,sep = ""))
#> [1] "Condition ERCC"
#--------------------------------------
n_cells = length(get.cell.size(object = obj))
print(paste("n cells", n_cells, sep = " "))
#> [1] "n cells 1015"

n_it = 1

1.1 Data cleaning

First we create the directory to store all information regarding the data cleaning.

if(!file.exists(out_dir)){
  dir.create(file.path(out_dir))
}

if(!file.exists(paste(out_dir,"cleaning", sep = ""))){   
  dir.create(file.path(out_dir, "cleaning"))
}
ttm = clean(obj)
#> [1] "Start estimation mu with linear method"
#> [1]   88 1015
#> [1] "Start PCA"
#> [1] "starting hclust"

obj = ttm$object

ttm$pca.cell.2

Run this when B cells need to be removed.

pdf(file.path(out_dir,"cleaning",paste(t,"_",n_it,"_plots_before_cells_exlusion.pdf", sep = "")))
ttm$pca.cell.2
ggplot(ttm$D, aes(x=n,y=means)) + geom_point() +
  geom_text_repel(data=subset(ttm$D, n > (max(ttm$D$n)- 15) ), aes(n,means,label=rownames(ttm$D[ttm$D$n > (max(ttm$D$n)- 15),])),
                  nudge_y      = 0.05,
                  nudge_x      = 0.05,
                  direction    = "x",
                  angle        = 90,
                  vjust        = 0,
                  segment.size = 0.2)+
  ggtitle("B cell group genes mean expression")+my_theme +
  theme(plot.title = element_text(color = "#3C5488FF", size = 20, face = "italic",vjust = - 5,hjust = 0.02 ))
dev.off()

if (length(ttm$cl1) < length(ttm$cl2)) {
  to_rem = ttm$cl1
}else{
  to_rem = ttm$cl2
}
n_it = n_it+1
obj = drop.genes.cells(object = obj,genes = c(),cells = to_rem)

gc()

ttm = clean(obj)
#ttm = clean.sqrt(obj, cells)
obj = ttm$object

ttm$pca.cell.2

Run this only in the last iteration, instead the previous code, when B cells group has not to be removed

pdf(file.path(out_dir,"cleaning",paste(t,"_",n_it,"_plots_before_cells_exlusion.pdf", sep = "")))
ttm$pca.cell.2
ggplot(ttm$D, aes(x=n,y=means)) + geom_point() +
  geom_text_repel(data=subset(ttm$D, n > (max(ttm$D$n)- 15) ), aes(n,means,label=rownames(ttm$D[ttm$D$n > (max(ttm$D$n)- 15),])),
                  nudge_y      = 0.05,
                  nudge_x      = 0.05,
                  direction    = "x",
                  angle        = 90,
                  vjust        = 0,
                  segment.size = 0.2)+
  ggtitle(label = "B cell group genes mean expression", subtitle = " - B group NOT removed -")+my_theme +
  theme(plot.title = element_text(color = "#3C5488FF", size = 20, face = "italic",vjust = - 10,hjust = 0.02 ),
        plot.subtitle = element_text(color = "darkred",vjust = - 15,hjust = 0.01 ))

dev.off()
#> png 
#>   2

To color the pca based on nu_j (so the cells’ efficiency)

nu_est = round(get.nu(object = obj), digits = 7)

plot.nu <-ggplot(ttm$pca_cells,aes(x=PC1,y=PC2, colour = log(nu_est)))

plot.nu = plot.nu + geom_point(size = 1,alpha= 0.8)+
  scale_color_gradient2(low = "#E64B35B2",mid =  "#4DBBD5B2", high =  "#3C5488B2" ,
                        midpoint = log(mean(nu_est)),name = "ln(nu)")+
  ggtitle("Cells PCA coloured by cells efficiency") +
  my_theme +  theme(plot.title = element_text(color = "#3C5488FF", size = 20),
                    legend.title=element_text(color = "#3C5488FF", size = 14,face = "italic"),
                    legend.text = element_text(color = "#3C5488FF", size = 11),
                    legend.key.width = unit(2, "mm"),
                    legend.position="right")

pdf(file.path(out_dir,"cleaning",paste(t,"_plots_PCA_efficiency_colored.pdf", sep = "")))
plot.nu
dev.off()
#> png 
#>   2

plot.nu

The next part is use to remove the cells with efficiency too low.

nu_df = data.frame("nu"= sort(get.nu(obj)), "n"=c(1:length(get.nu(obj))))

ggplot(nu_df, aes(x = n, y=nu)) + 
    geom_point(colour = "#8491B4B2", size=1)+
    my_theme #+ ylim(0,1) + xlim(0,70)

We can zoom on the smallest values and, if we detect a clear elbow, we can decide to remove the cells.

yset = 0.4#threshold to remove low UDE cells
plot.ude <- ggplot(nu_df, aes(x = n, y=nu)) + 
    geom_point(colour = "#8491B4B2", size=1) + 
    my_theme + ylim(0,1) + xlim(0,400) +
    geom_hline(yintercept=yset, linetype="dashed", color = "darkred") +
    annotate(geom="text", x=200, y=0.25, 
             label=paste("to remove cells with nu < ",yset,sep = " "), 
             color="darkred", size=4.5)


pdf(file.path(out_dir,"cleaning",paste(t,"_plots_efficiency.pdf", sep = "")))
plot.ude
#> Warning: Removed 668 rows containing missing values (geom_point).
dev.off()
#> png 
#>   2

plot.ude
#> Warning: Removed 668 rows containing missing values (geom_point).

We also save the defined threshold in the metadata and re run the estimation

obj = add.row.to.meta(obj,c("Threshold low UDE cells:",yset)) 

to_rem = rownames(nu_df[which(nu_df$nu < yset),])

obj = drop.genes.cells(object = obj, genes = c(),cells =  to_rem)

Repeat the estimation after the cells are removed

ttm = clean(obj)
#> [1] "Start estimation mu with linear method"
#> [1]   88 1004
#> [1] "Start PCA"
#> [1] "starting hclust"
obj = ttm$object
ttm$pca.cell.2

In case we do not want to remove anything, we can run:

pdf(file.path(out_dir,"cleaning",paste(t,"_plots_efficiency.pdf", sep = "")))
ggplot(nu_df, aes(x = n, y=nu)) + geom_point(colour = "#8491B4B2", size=1) +my_theme + #xlim(0,100)+
  annotate(geom="text", x=50, y=0.25, label="nothing to remove ", color="darkred")
dev.off()
#> png 
#>   2

Just to check again, we plot the final efficiency colored PCA

nu_est = round(get.nu(object = obj), digits = 7)
plot.nu <-ggplot(ttm$pca_cells,aes(x=PC1,y=PC2, colour = log(nu_est)))
plot.nu = plot.nu + geom_point(size = 2,alpha= 0.8)+
  scale_color_gradient2(low = "#E64B35B2",mid =  "#4DBBD5B2", high =  "#3C5488B2" ,
                        midpoint = log(mean(nu_est)),name = "ln(nu)")+
  ggtitle("Cells PCA coloured by cells efficiency: last") +
  my_theme +  theme(plot.title = element_text(color = "#3C5488FF", size = 20),
                    legend.title=element_text(color = "#3C5488FF", size = 14,face = "italic"),
                    legend.text = element_text(color = "#3C5488FF", size = 11),
                    legend.key.width = unit(2, "mm"),
                    legend.position="right")

pdf(file.path(out_dir,"cleaning",paste(t,"_plots_PCA_efficiency_colored_FINAL.pdf", sep = "")))
plot.nu
dev.off()
#> png 
#>   2

plot.nu

1.2 COTAN analysis

In this part all the contingency tables are computed and used to get the statistics (S) To storage efficiency of all the observed tables only the yes_yes is stored. Note that if will be necessary re-computing the yes-yes table, the yes-yes table should be cancelled before running cotan_analysis.

obj = cotan_analysis(obj,cores = 2)
#> [1] "cotan analysis"
#> [1] "mu estimator creation"
#> [1] "save effective constitutive genes"
#> [1] "start a minimization"
#> [1] "Final trance!"
#> [1] "a min: -0.147216796875 | a max 13.5546875 | negative a %: 26.984126984127"

COEX evaluation and storing

obj = get.coex(obj)
#> [1] "coex dataframe creation"
#> [1] "creation of observed yes/yes contingency table"
#> [1] "mu estimator creation"
#> [1] "expected contingency tables creation"
#> [1] "The distance between estimated n of zeros and observed number of zero is 0.0041370202609155 over 63"
#> [1] "Done"
#> [1] "coex estimation"
#> [1] "Cleaning RAM"

# saving the structure 
saveRDS(obj,file = paste(out_dir,t,".cotan.RDS", sep = ""))

2 Analysis of the elaborated data

2.1 GDI

The next function can directly plot the GDI for the dataset with the 1.5 threshold (in red) and the two higher quantiles (in blue).

plot_GDI(obj, cond = "ERCC")
#> [1] "GDI plot "
#> [1] "function to generate GDI dataframe"
#> [1] "Using S"
#> [1] "function to generate S "

If a more complex plot is needed, or if we want to analyze more in detail the GDI data, we can get directly the GDI dataframe.

quant.p = get.GDI(obj)
#> [1] "function to generate GDI dataframe"
#> [1] "Using S"
#> [1] "function to generate S "

head(quant.p)
#>            sum.raw.norm       GDI exp.cells
#> ERCC-00012     3.175703 1.1614957  2.390438
#> ERCC-00013     5.995188 1.3512649 27.788845
#> ERCC-00014     6.077256 1.3473405 34.860558
#> ERCC-00016     3.569692 1.1518864  3.685259
#> ERCC-00017     1.609829 0.7957451  0.498008
#> ERCC-00019     7.157739 1.4853109 66.135458

In the third column of this dataframe is reported the percentage of cells expressing the gene.

Next we can define some gene sets (in this case three) that we want to specifically label in the GDI plot.

AA=c("ERCC-00012","ERCC-00013","ERCC-00014") 
BB=c("ERCC-00016","ERCC-00017","ERCC-00019")
CC=c("ERCC-00022","ERCC-00024","ERCC-00028")

text.size = 12

quant.p$highlight = with(quant.p, ifelse(rownames(quant.p) %in% AA, "AA",
                                                 ifelse(rownames(quant.p) %in% CC,"Constitutive" ,
                                                               ifelse(rownames(quant.p) %in% BB,"BB" , "normal"))))

textdf <- quant.p[rownames(quant.p) %in% c(AA,CC,BB), ]
mycolours <- c("Con" = "#00A087FF","AA"="#E64B35FF","BB"="#F39B7FFF","normal" = "#8491B4B2")
f1 = ggplot(subset(quant.p,highlight == "normal" ), aes(x=sum.raw.norm, y=GDI)) +  geom_point(alpha = 0.1, color = "#8491B4B2", size=2.5)
GDI_plot = f1 +  geom_point(data = subset(quant.p,highlight != "normal"  ), aes(x=sum.raw.norm, y=GDI, colour=highlight),size=2.5,alpha = 0.8) +
  geom_hline(yintercept=quantile(quant.p$GDI)[4], linetype="dashed", color = "darkblue") +
  geom_hline(yintercept=quantile(quant.p$GDI)[3], linetype="dashed", color = "darkblue") +
  geom_hline(yintercept=1.5, linetype="dotted", color = "red", size= 0.5) +
  scale_color_manual("Status", values = mycolours)  +
  scale_fill_manual("Status", values = mycolours)  +
  xlab("log normalized counts")+ylab("GDI")+
  geom_label_repel(data = textdf , aes(x=sum.raw.norm, y=GDI, label = rownames(textdf),fill=highlight),
                   label.size = NA, 
                   alpha = 0.5, 
                   direction ="both",
                   na.rm=TRUE,
                   seed = 1234) +
  geom_label_repel(data =textdf , aes(x=sum.raw.norm, y=GDI, label = rownames(textdf)),
                   label.size = NA, 
                   segment.color = 'black',
                   segment.size = 0.5,
                   direction = "both",
                   alpha = 0.8, 
                   na.rm=TRUE,
                   fill = NA,
                   seed = 1234) +
  theme(axis.text.x = element_text(size = text.size, angle = 0, hjust = .5, vjust = .5, face = "plain", colour ="#3C5488FF" ),
        axis.text.y = element_text( size = text.size, angle = 0, hjust = 0, vjust = .5, face = "plain", colour ="#3C5488FF"),  
        axis.title.x = element_text( size = text.size, angle = 0, hjust = .5, vjust = 0, face = "plain", colour ="#3C5488FF"),
        axis.title.y = element_text( size = text.size, angle = 90, hjust = .5, vjust = .5, face = "plain", colour ="#3C5488FF"),
        legend.title = element_blank(),
        legend.text = element_text(color = "#3C5488FF",face ="italic" ),
        legend.position = "bottom") 
legend <- cowplot::get_legend(GDI_plot)
GDI_plot =GDI_plot + theme(
        legend.position = "none") 
GDI_plot

2.2 Heatmaps

For the Gene Pair Analysis, we can plot an heatmap with the coex values between two genes sets. To do so we need to define, in a list, the different gene sets (list.genes). Then in the function parameter sets we can decide which sets need to be considered (in the example from 1 to 3). In the condition parameter we should insert an array with each file name prefix to be considered (in the example, the file names is “E17.5_cortex”).

list.genes = list("Ref.col"= BB, "AA"=AA, "Const."=CC )

plot_heatmap(df_genes = list.genes,sets = c(1:3),conditions = "ERCC",dir = out_dir)
#> [1] "plot heatmap"
#> [1] "Loading condition ERCC"
#> [1] "ERCC-00016" "ERCC-00017" "ERCC-00019"
#> [1] "Get p-values on a set of genes on columns on a set of genes on rows"
#> [1] "Using function S"
#> [1] "function to generate S "
#> [1] "Ref.col"
#> [1] "AA"
#> [1] "Const."
#> [1] "min coex: 0 max coex 0.98246092315627"

We can also plot a general heatmap of coex values based on some markers like the following one.

plot_general.heatmap(prim.markers = c("ERCC-00014","ERCC-00019"), condition = "ERCC",dir = out_dir, p_value = 0.05)
#> [1] "ploting a general heatmap"
#> NULL
#> [1] "Get p-values genome wide on columns genome wide on rows"
#> [1] "Using function S"
#> [1] "function to generate S "

plot_general.heatmap(prim.markers = c("ERCC-00014","ERCC-00019"),markers.list =c("ERCC-00084" ,"ERCC-00085" ,"ERCC-00086") ,symmetric = FALSE, condition = "ERCC",dir = out_dir, p_value = 0.05)
#> [1] "ploting a general heatmap"
#> NULL
#> [1] "Get p-values genome wide on columns genome wide on rows"
#> [1] "Using function S"
#> [1] "function to generate S "

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:

  1. get.observed.ct to produce the observed data
 get.observed.ct(object = obj, g1 = "ERCC-00014",g2 = "ERCC-00019")
#>                ERCC-00014.yes ERCC-00014.no
#> ERCC-00019.yes            241           423
#> ERCC-00019.no             109           231
  1. get.expected.ct to produce the expected values
get.expected.ct(object = obj, g1 = "ERCC-00014",g2 = "ERCC-00019")
#>                ERCC-00014.yes ERCC-00014.no
#> ERCC-00019.yes       235.7119      428.2888
#> ERCC-00019.no        114.2885      225.7108

Another useful function is extract.coex. This can be used to extract the whole or a partial coex matrix from a COTAN object.

# For the whole matrix
coex <- extract.coex(object = obj)
coex[1:5,1:5]
#>              ERCC-00012  ERCC-00013   ERCC-00014   ERCC-00016   ERCC-00017
#> ERCC-00012  0.769954876  0.03079731 -0.008469118  0.035409255 -0.004367265
#> ERCC-00013  0.030797310  0.98637914  0.020020333 -0.018534811 -0.013478967
#> ERCC-00014 -0.008469118  0.02002033  0.983490318 -0.003222142  0.094474696
#> ERCC-00016  0.035409255 -0.01853481 -0.003222142  0.982460923  0.028389361
#> ERCC-00017 -0.004367265 -0.01347897  0.094474696  0.028389361  0.185954347
# For a partial matrix
coex <- extract.coex(object = obj,genes = c("ERCC-00017", "ERCC-00019", "ERCC-00024", "ERCC-00025", "ERCC-00028", "ERCC-00031"))
head(coex)
#>              ERCC-00017   ERCC-00019   ERCC-00024   ERCC-00025    ERCC-00028
#> ERCC-00012 -0.004367265 5.272805e-02 -0.011980011  0.017574171 -0.0120737279
#> ERCC-00013 -0.013478967 4.360570e-02 -0.005497781  0.068917037  0.0310017944
#> ERCC-00014  0.094474696 2.345844e-02 -0.016243015 -0.006247842 -0.0226624019
#> ERCC-00016  0.028389361 7.232674e-05 -0.004746929  0.033126881  0.0005889499
#> ERCC-00017  0.185954347 1.880656e-02 -0.010473707  0.012328879  0.0162721054
#> ERCC-00019  0.018806561 9.695314e-01 -0.010600940 -0.032637961  0.0473851891
#>              ERCC-00031
#> ERCC-00012  0.013066001
#> ERCC-00013 -0.020289792
#> ERCC-00014 -0.054916592
#> ERCC-00016 -0.012314141
#> ERCC-00017 -0.017522836
#> ERCC-00019 -0.003396497

The next few lines are just to clean after compilation of the documentation.

if (file.exists(paste(out_dir,t,".cotan.RDS", sep = ""))) {
  #Delete file if it exists
  file.remove(paste(out_dir,t,".cotan.RDS", sep = ""))
}
#> [1] TRUE
unlink(paste(out_dir,"cleaning", sep = ""),recursive = TRUE)

print(sessionInfo())
#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [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       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] dendextend_1.16.0 MASS_7.3-58.1     htmlwidgets_1.5.4 forcats_0.5.2    
#>  [5] stringr_1.4.1     dplyr_1.0.10      purrr_0.3.5       readr_2.1.3      
#>  [9] tidyr_1.2.1       tibble_3.1.8      tidyverse_1.3.2   plotly_4.10.0    
#> [13] Rtsne_0.16        factoextra_1.0.7  ggrepel_0.9.1     ggplot2_3.3.6    
#> [17] Matrix_1.5-1      data.table_1.14.4 COTAN_1.2.0       BiocStyle_2.26.0 
#> 
#> loaded via a namespace (and not attached):
#>  [1] matrixStats_0.62.0    fs_1.5.2              lubridate_1.8.0      
#>  [4] doParallel_1.0.17     RColorBrewer_1.1-3    httr_1.4.4           
#>  [7] tools_4.2.1           backports_1.4.1       bslib_0.4.0          
#> [10] utf8_1.2.2            R6_2.5.1              irlba_2.3.5.1        
#> [13] DBI_1.1.3             lazyeval_0.2.2        BiocGenerics_0.44.0  
#> [16] colorspace_2.0-3      GetoptLong_1.0.5      withr_2.5.0          
#> [19] gridExtra_2.3         tidyselect_1.2.0      compiler_4.2.1       
#> [22] rvest_1.0.3           cli_3.4.1             Cairo_1.6-0          
#> [25] xml2_1.3.3            labeling_0.4.2        bookdown_0.29        
#> [28] sass_0.4.2            scales_1.2.1          digest_0.6.30        
#> [31] rmarkdown_2.17        pkgconfig_2.0.3       htmltools_0.5.3      
#> [34] highr_0.9             dbplyr_2.2.1          fastmap_1.1.0        
#> [37] readxl_1.4.1          rlang_1.0.6           GlobalOptions_0.1.2  
#> [40] farver_2.1.1          shape_1.4.6           jquerylib_0.1.4      
#> [43] generics_0.1.3        jsonlite_1.8.3        googlesheets4_1.0.1  
#> [46] magrittr_2.0.3        Rcpp_1.0.9            munsell_0.5.0        
#> [49] S4Vectors_0.36.0      fansi_1.0.3           viridis_0.6.2        
#> [52] lifecycle_1.0.3       RcppZiggurat_0.1.6    stringi_1.7.8        
#> [55] yaml_2.3.6            grid_4.2.1            parallel_4.2.1       
#> [58] crayon_1.5.2          lattice_0.20-45       cowplot_1.1.1        
#> [61] haven_2.5.1           circlize_0.4.15       hms_1.1.2            
#> [64] magick_2.7.3          knitr_1.40            ComplexHeatmap_2.14.0
#> [67] pillar_1.8.1          rjson_0.2.21          codetools_0.2-18     
#> [70] stats4_4.2.1          reprex_2.0.2          glue_1.6.2           
#> [73] evaluate_0.17         modelr_0.1.9          BiocManager_1.30.19  
#> [76] tzdb_0.3.0            png_0.1-7             vctrs_0.5.0          
#> [79] foreach_1.5.2         cellranger_1.1.0      gtable_0.3.1         
#> [82] clue_0.3-62           assertthat_0.2.1      cachem_1.0.6         
#> [85] xfun_0.34             Rfast_2.0.6           broom_1.0.1          
#> [88] googledrive_2.0.0     viridisLite_0.4.1     gargle_1.2.1         
#> [91] iterators_1.0.14      IRanges_2.32.0        cluster_2.1.4        
#> [94] ellipsis_0.3.2