2.1 Create scGPS objects

# load mixed population 1 (loaded from day_2_cardio_cell_sample dataset, 
# named it as day2)
library(scGPS)

day2 <- day_2_cardio_cell_sample
mixedpop1 <- new_scGPS_object(ExpressionMatrix = day2$dat2_counts,
    GeneMetadata = day2$dat2geneInfo, CellMetadata = day2$dat2_clusters)

# load mixed population 2 (loaded from day_5_cardio_cell_sample dataset, 
# named it as day5)
day5 <- day_5_cardio_cell_sample
mixedpop2 <- new_scGPS_object(ExpressionMatrix = day5$dat5_counts,
    GeneMetadata = day5$dat5geneInfo, CellMetadata = day5$dat5_clusters)

2.2 Run prediction

# select a subpopulation
c_selectID <- 1
# load gene list (this can be any lists of user selected genes)
genes <- training_gene_sample
genes <- genes$Merged_unique
# load cluster information 
cluster_mixedpop1 <- colData(mixedpop1)[,1]
cluster_mixedpop2 <- colData(mixedpop2)[,1]
#run training (running nboots = 3 here, but recommend to use nboots = 50-100)
LSOLDA_dat <- bootstrap_prediction(nboots = 3, mixedpop1 = mixedpop1, 
    mixedpop2 = mixedpop2, genes = genes, c_selectID  = c_selectID,
    listData = list(), cluster_mixedpop1 = cluster_mixedpop1, 
    cluster_mixedpop2 = cluster_mixedpop2, trainset_ratio = 0.7)
names(LSOLDA_dat)
#> [1] "Accuracy"          "ElasticNetGenes"   "Deviance"         
#> [4] "ElasticNetFit"     "LDAFit"            "predictor_S1"     
#> [7] "ElasticNetPredict" "LDAPredict"        "cell_results"

2.3 Summarise results

# summary results LDA
sum_pred_lda <- summary_prediction_lda(LSOLDA_dat = LSOLDA_dat, nPredSubpop = 4)
# summary results Lasso to show the percent of cells
# classified as cells belonging 
sum_pred_lasso <- summary_prediction_lasso(LSOLDA_dat = LSOLDA_dat,
    nPredSubpop = 4)
# plot summary results 
plot_sum <-function(sum_dat){
    sum_dat_tf <- t(sum_dat)
    sum_dat_tf <- na.omit(sum_dat_tf)
    sum_dat_tf <- apply(sum_dat[, -ncol(sum_dat)],1,
        function(x){as.numeric(as.vector(x))})
    sum_dat$names <- gsub("ElasticNet for subpop","sp",  sum_dat$names )
    sum_dat$names <- gsub("in target mixedpop","in p",  sum_dat$names) 
    sum_dat$names <- gsub("LDA for subpop","sp",  sum_dat$names )
    sum_dat$names <- gsub("in target mixedpop","in p",  sum_dat$names)
    colnames(sum_dat_tf) <- sum_dat$names
    boxplot(sum_dat_tf, las=2)
}
plot_sum(sum_pred_lasso)

plot_sum(sum_pred_lda)

# summary accuracy to check the model accuracy in the leave-out test set 
summary_accuracy(object = LSOLDA_dat)
#> [1] 67.72727 62.73585 61.32075
# summary maximum deviance explained by the model 
summary_deviance(object = LSOLDA_dat)
#> $allDeviance
#> [1] "8.26" "8.25" "6.94"
#> 
#> $DeviMax
#>          dat_DE$Dfd          Deviance           DEgenes
#> 1                 0              8.26    genes_cluster1
#> 2                 1              8.26    genes_cluster1
#> 3                 2              8.26    genes_cluster1
#> 4                 3              8.26    genes_cluster1
#> 5                 4              8.26    genes_cluster1
#> 6 remaining DEgenes remaining DEgenes remaining DEgenes
#> 
#> $LassoGenesMax
#> NULL

3.1 Identify clusters in a dataset using CORE

(skip this step if clusters are known)

# find clustering information in an expresion data using CORE
day5 <- day_5_cardio_cell_sample
cellnames <- colnames(day5$dat5_counts)
cluster <-day5$dat5_clusters
cellnames <-data.frame("Cluster"=cluster, "cellBarcodes" = cellnames)
mixedpop2 <-new_scGPS_object(ExpressionMatrix = day5$dat5_counts,
                    GeneMetadata = day5$dat5geneInfo, CellMetadata = cellnames)

CORE_cluster <- CORE_clustering(mixedpop2, remove_outlier = c(0), PCA=FALSE)

# to update the clustering information, users can ...
key_height <- CORE_cluster$optimalClust$KeyStats$Height
optimal_res <- CORE_cluster$optimalClust$OptimalRes
optimal_index = which(key_height == optimal_res)

clustering_after_outlier_removal <- unname(unlist(
 CORE_cluster$Cluster[[optimal_index]]))
corresponding_cells_after_outlier_removal <- CORE_cluster$cellsForClustering
original_cells_before_removal <- colData(mixedpop2)[,2]
corresponding_index <- match(corresponding_cells_after_outlier_removal,
                            original_cells_before_removal )
# check the matching
identical(as.character(original_cells_before_removal[corresponding_index]),
         corresponding_cells_after_outlier_removal)
#> [1] TRUE
# create new object with the new clustering after removing outliers
mixedpop2_post_clustering <- mixedpop2[,corresponding_index]
colData(mixedpop2_post_clustering)[,1] <- clustering_after_outlier_removal

3.2 Identify clusters in a dataset using SCORE (Stable Clustering at Optimal REsolution)

(skip this step if clusters are known)

(SCORE aims to get stable subpopulation results by introducing bagging aggregation and bootstrapping to the CORE algorithm)

# find clustering information in an expresion data using SCORE
day5 <- day_5_cardio_cell_sample
cellnames <- colnames(day5$dat5_counts)
cluster <-day5$dat5_clusters
cellnames <-data.frame("Cluster"=cluster, "cellBarcodes" = cellnames)
mixedpop2 <-new_scGPS_object(ExpressionMatrix = day5$dat5_counts,
                    GeneMetadata = day5$dat5geneInfo, CellMetadata = cellnames )

SCORE_test <- CORE_bagging(mixedpop2, remove_outlier = c(0), PCA=FALSE,
                                bagging_run = 20, subsample_proportion = .8)

3.3 Visualise all cluster results in all iterations

dev.off()
#> null device 
#>           1
##3.2.1 plot CORE clustering
p1 <- plot_CORE(CORE_cluster$tree, CORE_cluster$Cluster, 
    color_branch = c("#208eb7", "#6ce9d3", "#1c5e39", "#8fca40", "#154975",
        "#b1c8eb"))
p1
#> $mar
#> [1] 1 5 0 1
#extract optimal index identified by CORE
key_height <- CORE_cluster$optimalClust$KeyStats$Height
optimal_res <- CORE_cluster$optimalClust$OptimalRes
optimal_index = which(key_height == optimal_res)
#plot one optimal clustering bar
plot_optimal_CORE(original_tree= CORE_cluster$tree,
                 optimal_cluster = unlist(CORE_cluster$Cluster[optimal_index]),
                 shift = -2000)
#> Ordering and assigning labels...
#> 2
#> 162335NA
#> 3
#> 162335423
#> Plotting the colored dendrogram now....
#> Plotting the bar underneath now....

##3.2.2 plot SCORE clustering
#plot all clustering bars
plot_CORE(SCORE_test$tree, list_clusters = SCORE_test$Cluster)
#plot one stable optimal clustering bar
plot_optimal_CORE(original_tree= SCORE_test$tree,
                 optimal_cluster = unlist(SCORE_test$Cluster[
                    SCORE_test$optimal_index]),
                 shift = -100)
#> Ordering and assigning labels...
#> 2
#> 162335NA
#> 3
#> 162335423
#> Plotting the colored dendrogram now....
#> Plotting the bar underneath now....

3.4 Compare clustering results with other dimensional reduction methods (e.g., tSNE)

t <- tSNE(expression.mat=assay(mixedpop2))
#> Preparing PCA inputs using the top 1500 genes ...
#> Computing PCA values...
#> Running tSNE ...
p2 <-plot_reduced(t, color_fac = factor(colData(mixedpop2)[,1]),
                      palletes =1:length(unique(colData(mixedpop2)[,1])))
#> Warning: Use of `reduced_dat_toPlot$Dim1` is discouraged. Use `Dim1` instead.
#> Warning: Use of `reduced_dat_toPlot$Dim2` is discouraged. Use `Dim2` instead.
p2

3.5 Find gene markers and annotate clusters

#load gene list (this can be any lists of user-selected genes)
genes <-training_gene_sample
genes <-genes$Merged_unique

#the gene list can also be objectively identified by differential expression
#analysis cluster information is requied for find_markers. Here, we use
#CORE results.

#colData(mixedpop2)[,1] <- unlist(SCORE_test$Cluster[SCORE_test$optimal_index])

suppressMessages(library(locfit))

DEgenes <- find_markers(expression_matrix=assay(mixedpop2),
                            cluster = colData(mixedpop2)[,1],
                            selected_cluster=unique(colData(mixedpop2)[,1]))

#the output contains dataframes for each cluster.
#the data frame contains all genes, sorted by p-values
names(DEgenes)
#> [1] "baseMean"       "log2FoldChange" "lfcSE"          "stat"          
#> [5] "pvalue"         "padj"           "id"

#you can annotate the identified clusters
DEgeneList_1vsOthers <- DEgenes$DE_Subpop1vsRemaining$id

#users need to check the format of the gene input to make sure they are
#consistent to the gene names in the expression matrix

#the following command saves the file "PathwayEnrichment.xlsx" to the
#working dir
#use 500 top DE genes
suppressMessages(library(DOSE))
suppressMessages(library(ReactomePA))
suppressMessages(library(clusterProfiler))
genes500 <- as.factor(DEgeneList_1vsOthers[seq_len(500)])
enrichment_test <- annotate_clusters(genes, pvalueCutoff=0.05, gene_symbol=TRUE)

#the enrichment outputs can be displayed by running
clusterProfiler::dotplot(enrichment_test, showCategory=10, font.size = 6)

4.1 Start the scGPS prediction to find relationship between clusters

#select a subpopulation, and input gene list
c_selectID <- 1
#note make sure the format for genes input here is the same to the format
#for genes in the mixedpop1 and mixedpop2
genes = DEgenes$id[1:500]

#run the test bootstrap with nboots = 2 runs

cluster_mixedpop1 <- colData(mixedpop1)[,1]
cluster_mixedpop2 <- colData(mixedpop2)[,1]

LSOLDA_dat <- bootstrap_prediction(nboots = 2, mixedpop1 = mixedpop1,
                             mixedpop2 = mixedpop2, genes = genes, 
                             c_selectID  = c_selectID,
                             listData = list(),
                             cluster_mixedpop1 = cluster_mixedpop1,
                             cluster_mixedpop2 = cluster_mixedpop2)

4.2 Display summary results for the prediction

#get the number of rows for the summary matrix
row_cluster <-length(unique(colData(mixedpop2)[,1]))

#summary results LDA to to show the percent of cells classified as cells
#belonging by LDA classifier
summary_prediction_lda(LSOLDA_dat=LSOLDA_dat, nPredSubpop = row_cluster )
#>                 V1               V2                                names
#> 1 85.0267379679144  74.331550802139 LDA for subpop 1 in target mixedpop2
#> 2               95 32.8571428571429 LDA for subpop 2 in target mixedpop2
#> 3 86.4661654135338 83.4586466165414 LDA for subpop 3 in target mixedpop2
#> 4             77.5             52.5 LDA for subpop 4 in target mixedpop2

#summary results Lasso to show the percent of cells classified as cells
#belonging by Lasso classifier
summary_prediction_lasso(LSOLDA_dat=LSOLDA_dat, nPredSubpop = row_cluster)
#>                 V1               V2                                      names
#> 1 58.2887700534759 49.7326203208556 ElasticNet for subpop1 in target mixedpop2
#> 2 99.2857142857143 93.5714285714286 ElasticNet for subpop2 in target mixedpop2
#> 3 93.2330827067669 86.4661654135338 ElasticNet for subpop3 in target mixedpop2
#> 4               95               90 ElasticNet for subpop4 in target mixedpop2

# summary maximum deviance explained by the model during the model training
summary_deviance(object = LSOLDA_dat)
#> $allDeviance
#> [1] "51.27" "49.9" 
#> 
#> $DeviMax
#>           dat_DE$Dfd          Deviance           DEgenes
#> 1                  0             51.27    genes_cluster1
#> 2                  1             51.27    genes_cluster1
#> 3                  4             51.27    genes_cluster1
#> 4                  6             51.27    genes_cluster1
#> 5                  8             51.27    genes_cluster1
#> 6                 10             51.27    genes_cluster1
#> 7                 12             51.27    genes_cluster1
#> 8                 14             51.27    genes_cluster1
#> 9                 15             51.27    genes_cluster1
#> 10                16             51.27    genes_cluster1
#> 11                17             51.27    genes_cluster1
#> 12                18             51.27    genes_cluster1
#> 13                20             51.27    genes_cluster1
#> 14                22             51.27    genes_cluster1
#> 15                27             51.27    genes_cluster1
#> 16 remaining DEgenes remaining DEgenes remaining DEgenes
#> 
#> $LassoGenesMax
#> NULL

# summary accuracy to check the model accuracy in the leave-out test set
summary_accuracy(object = LSOLDA_dat)
#> [1] 65.62500 67.41071

4.3 Plot the relationship between clusters in one sample

Here we look at one example use case to find relationship between clusters within one sample or between two sample

#run prediction for 3 clusters
cluster_mixedpop1 <- colData(mixedpop1)[,1]
cluster_mixedpop2 <- colData(mixedpop2)[,1]
#cluster_mixedpop2 <- as.numeric(as.vector(colData(mixedpop2)[,1]))

c_selectID <- 1
#top 200 gene markers distinguishing cluster 1
genes = DEgenes$id[1:200]

LSOLDA_dat1 <- bootstrap_prediction(nboots = 2, mixedpop1 = mixedpop2,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop2,
                        cluster_mixedpop2 = cluster_mixedpop2)

c_selectID <- 2
genes = DEgenes$id[1:200]

LSOLDA_dat2 <- bootstrap_prediction(nboots = 2,mixedpop1 = mixedpop2,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop2,
                        cluster_mixedpop2 = cluster_mixedpop2)

c_selectID <- 3
genes = DEgenes$id[1:200]
LSOLDA_dat3 <- bootstrap_prediction(nboots = 2,mixedpop1 = mixedpop2,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop2,
                        cluster_mixedpop2 = cluster_mixedpop2)

c_selectID <- 4
genes = DEgenes$id[1:200]
LSOLDA_dat4 <- bootstrap_prediction(nboots = 2,mixedpop1 = mixedpop2,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop2,
                        cluster_mixedpop2 = cluster_mixedpop2)


#prepare table input for sankey plot

LASSO_C1S2  <- reformat_LASSO(c_selectID=1, mp_selectID = 2,
                             LSOLDA_dat=LSOLDA_dat1,
                             nPredSubpop = length(unique(colData(mixedpop2)
                                [,1])),
                             Nodes_group ="#7570b3")

LASSO_C2S2  <- reformat_LASSO(c_selectID=2, mp_selectID =2,
                             LSOLDA_dat=LSOLDA_dat2,
                             nPredSubpop = length(unique(colData(mixedpop2)
                                [,1])),
                             Nodes_group ="#1b9e77")

LASSO_C3S2  <- reformat_LASSO(c_selectID=3, mp_selectID =2,
                             LSOLDA_dat=LSOLDA_dat3,
                             nPredSubpop = length(unique(colData(mixedpop2)
                                [,1])),
                             Nodes_group ="#e7298a")

LASSO_C4S2  <- reformat_LASSO(c_selectID=4, mp_selectID =2,
                             LSOLDA_dat=LSOLDA_dat4,
                             nPredSubpop = length(unique(colData(mixedpop2)
                                [,1])),
                             Nodes_group ="#00FFFF")

combined <- rbind(LASSO_C1S2,LASSO_C2S2,LASSO_C3S2, LASSO_C4S2 )
combined <- combined[is.na(combined$Value) != TRUE,]

nboots = 2
#links: source, target, value
#source: node, nodegroup
combined_D3obj <-list(Nodes=combined[,(nboots+3):(nboots+4)],
                     Links=combined[,c((nboots+2):(nboots+1),ncol(combined))])

library(networkD3)

Node_source <- as.vector(sort(unique(combined_D3obj$Links$Source)))
Node_target <- as.vector(sort(unique(combined_D3obj$Links$Target)))
Node_all <-unique(c(Node_source, Node_target))

#assign IDs for Source (start from 0)
Source <-combined_D3obj$Links$Source
Target <- combined_D3obj$Links$Target

for(i in 1:length(Node_all)){
   Source[Source==Node_all[i]] <-i-1
   Target[Target==Node_all[i]] <-i-1
}
# 
combined_D3obj$Links$Source <- as.numeric(Source)
combined_D3obj$Links$Target <- as.numeric(Target)
combined_D3obj$Links$LinkColor <- combined$NodeGroup

#prepare node info
node_df <-data.frame(Node=Node_all)
node_df$id <-as.numeric(c(0, 1:(length(Node_all)-1)))

suppressMessages(library(dplyr))
Color <- combined %>% count(Node, color=NodeGroup) %>% select(2)
node_df$color <- Color$color

suppressMessages(library(networkD3))
p1<-sankeyNetwork(Links =combined_D3obj$Links, Nodes = node_df,
                 Value = "Value", NodeGroup ="color", LinkGroup = "LinkColor", 
                 NodeID="Node", Source="Source", Target="Target", fontSize = 22)
p1

C2_MP2 → C2_MP2
93 

C4_MP2 → C4_MP2
88 

C3_MP2 → C3_MP2
86 

C1_MP2 → C1_MP2
73 

C2_MP2 → C1_MP2
59 

C4_MP2 → C3_MP2
56 

C4_MP2 → C2_MP2
55 

C4_MP2 → C1_MP2
37 

C3_MP2 → C1_MP2
33 

C2_MP2 → C4_MP2
26 

C3_MP2 → C4_MP2
26 

C1_MP2 → C4_MP2
10 

C1_MP2 → C3_MP2
4 

C3_MP2 → C2_MP2
3 

C1_MP2 → C2_MP2
0 

C2_MP2 → C3_MP2
0 

C1_MP2
202 

C2_MP2
179 

C3_MP2
147 

C4_MP2
235 

4.3 Plot the relationship between clusters in two samples

Here we look at one example use case to find relationship between clusters within one sample or between two sample

#run prediction for 3 clusters
cluster_mixedpop1 <- colData(mixedpop1)[,1]
cluster_mixedpop2 <- colData(mixedpop2)[,1]
row_cluster <-length(unique(colData(mixedpop2)[,1]))

c_selectID <- 1
#top 200 gene markers distinguishing cluster 1
genes = DEgenes$id[1:200]
LSOLDA_dat1 <- bootstrap_prediction(nboots = 2, mixedpop1 = mixedpop1,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop1,
                        cluster_mixedpop2 = cluster_mixedpop2)


c_selectID <- 2
genes = DEgenes$id[1:200]
LSOLDA_dat2 <- bootstrap_prediction(nboots = 2,mixedpop1 = mixedpop1,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop1,
                        cluster_mixedpop2 = cluster_mixedpop2)

c_selectID <- 3
genes = DEgenes$id[1:200]
LSOLDA_dat3 <- bootstrap_prediction(nboots = 2,mixedpop1 = mixedpop1,
                        mixedpop2 = mixedpop2, genes=genes, c_selectID, 
                        listData =list(),
                        cluster_mixedpop1 = cluster_mixedpop1,
                        cluster_mixedpop2 = cluster_mixedpop2)

#prepare table input for sankey plot

LASSO_C1S1  <- reformat_LASSO(c_selectID=1, mp_selectID = 1,
                             LSOLDA_dat=LSOLDA_dat1, nPredSubpop = row_cluster, 
                             Nodes_group = "#7570b3")

LASSO_C2S1  <- reformat_LASSO(c_selectID=2, mp_selectID = 1,
                             LSOLDA_dat=LSOLDA_dat2, nPredSubpop = row_cluster, 
                             Nodes_group = "#1b9e77")

LASSO_C3S1  <- reformat_LASSO(c_selectID=3, mp_selectID = 1,
                             LSOLDA_dat=LSOLDA_dat3, nPredSubpop = row_cluster, 
                             Nodes_group = "#e7298a")


combined <- rbind(LASSO_C1S1,LASSO_C2S1,LASSO_C3S1)

nboots = 2
#links: source, target, value
#source: node, nodegroup
combined_D3obj <-list(Nodes=combined[,(nboots+3):(nboots+4)],
                     Links=combined[,c((nboots+2):(nboots+1),ncol(combined))])
combined <- combined[is.na(combined$Value) != TRUE,]


library(networkD3)

Node_source <- as.vector(sort(unique(combined_D3obj$Links$Source)))
Node_target <- as.vector(sort(unique(combined_D3obj$Links$Target)))
Node_all <-unique(c(Node_source, Node_target))

#assign IDs for Source (start from 0)
Source <-combined_D3obj$Links$Source
Target <- combined_D3obj$Links$Target

for(i in 1:length(Node_all)){
   Source[Source==Node_all[i]] <-i-1
   Target[Target==Node_all[i]] <-i-1
}

combined_D3obj$Links$Source <- as.numeric(Source)
combined_D3obj$Links$Target <- as.numeric(Target)
combined_D3obj$Links$LinkColor <- combined$NodeGroup

#prepare node info
node_df <-data.frame(Node=Node_all)
node_df$id <-as.numeric(c(0, 1:(length(Node_all)-1)))

suppressMessages(library(dplyr))
n <- length(unique(node_df$Node))
getPalette = colorRampPalette(RColorBrewer::brewer.pal(9, "Set1"))
Color = getPalette(n)
node_df$color <- Color
suppressMessages(library(networkD3))
p1<-sankeyNetwork(Links =combined_D3obj$Links, Nodes = node_df,
                 Value = "Value", NodeGroup ="color", LinkGroup = "LinkColor",
                 NodeID="Node", Source="Source", Target="Target", fontSize = 22)
p1

C2_MP1 → C3_MP2
79 

C3_MP1 → C3_MP2
60 

C3_MP1 → C1_MP2
59 

C2_MP1 → C2_MP2
51 

C3_MP1 → C4_MP2
50 

C2_MP1 → C4_MP2
49 

C1_MP1 → C2_MP2
36 

C3_MP1 → C2_MP2
34 

C1_MP1 → C4_MP2
24 

C2_MP1 → C1_MP2
23 

C1_MP1 → C1_MP2
17 

C1_MP1 → C3_MP2
6 

C1_MP1
83 

C2_MP1
201 

C3_MP1
203 

C1_MP2
99 

C2_MP2
121 

C3_MP2
144 

C4_MP2
123 

devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.1 (2021-08-10)
#>  os       Ubuntu 20.04.3 LTS          
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  C                           
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2021-10-26                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version    date       lib source        
#>  annotate               1.72.0     2021-10-26 [2] Bioconductor  
#>  AnnotationDbi        * 1.56.0     2021-10-26 [2] Bioconductor  
#>  ape                    5.5        2021-04-25 [2] CRAN (R 4.1.1)
#>  aplot                  0.1.1      2021-09-22 [2] CRAN (R 4.1.1)
#>  assertthat             0.2.1      2019-03-21 [2] CRAN (R 4.1.1)
#>  backports              1.2.1      2020-12-09 [2] CRAN (R 4.1.1)
#>  Biobase              * 2.54.0     2021-10-26 [2] Bioconductor  
#>  BiocGenerics         * 0.40.0     2021-10-26 [2] Bioconductor  
#>  BiocParallel           1.28.0     2021-10-26 [2] Bioconductor  
#>  Biostrings             2.62.0     2021-10-26 [2] Bioconductor  
#>  bit                    4.0.4      2020-08-04 [2] CRAN (R 4.1.1)
#>  bit64                  4.0.5      2020-08-30 [2] CRAN (R 4.1.1)
#>  bitops                 1.0-7      2021-04-24 [2] CRAN (R 4.1.1)
#>  blob                   1.2.2      2021-07-23 [2] CRAN (R 4.1.1)
#>  bslib                  0.3.1      2021-10-06 [2] CRAN (R 4.1.1)
#>  cachem                 1.0.6      2021-08-19 [2] CRAN (R 4.1.1)
#>  callr                  3.7.0      2021-04-20 [2] CRAN (R 4.1.1)
#>  caret                * 6.0-90     2021-10-09 [2] CRAN (R 4.1.1)
#>  checkmate              2.0.0      2020-02-06 [2] CRAN (R 4.1.1)
#>  class                  7.3-19     2021-05-03 [2] CRAN (R 4.1.1)
#>  cli                    3.0.1      2021-07-17 [2] CRAN (R 4.1.1)
#>  clusterProfiler      * 4.2.0      2021-10-26 [2] Bioconductor  
#>  codetools              0.2-18     2020-11-04 [2] CRAN (R 4.1.1)
#>  colorspace             2.0-2      2021-06-24 [2] CRAN (R 4.1.1)
#>  cowplot                1.1.1      2020-12-30 [2] CRAN (R 4.1.1)
#>  crayon                 1.4.1      2021-02-08 [2] CRAN (R 4.1.1)
#>  data.table             1.14.2     2021-09-27 [2] CRAN (R 4.1.1)
#>  DBI                    1.1.1      2021-01-15 [2] CRAN (R 4.1.1)
#>  DelayedArray           0.20.0     2021-10-26 [2] Bioconductor  
#>  dendextend             1.15.1     2021-05-08 [2] CRAN (R 4.1.1)
#>  desc                   1.4.0      2021-09-28 [2] CRAN (R 4.1.1)
#>  DESeq2                 1.34.0     2021-10-26 [2] Bioconductor  
#>  devtools               2.4.2      2021-06-07 [2] CRAN (R 4.1.1)
#>  digest                 0.6.28     2021-09-23 [2] CRAN (R 4.1.1)
#>  DO.db                  2.9        2021-08-23 [2] Bioconductor  
#>  DOSE                 * 3.20.0     2021-10-26 [2] Bioconductor  
#>  downloader             0.4        2015-07-09 [2] CRAN (R 4.1.1)
#>  dplyr                * 1.0.7      2021-06-18 [2] CRAN (R 4.1.1)
#>  dynamicTreeCut       * 1.63-1     2016-03-11 [2] CRAN (R 4.1.1)
#>  e1071                  1.7-9      2021-09-16 [2] CRAN (R 4.1.1)
#>  ellipsis               0.3.2      2021-04-29 [2] CRAN (R 4.1.1)
#>  enrichplot             1.14.0     2021-10-26 [2] Bioconductor  
#>  evaluate               0.14       2019-05-28 [2] CRAN (R 4.1.1)
#>  fansi                  0.5.0      2021-05-25 [2] CRAN (R 4.1.1)
#>  farver                 2.1.0      2021-02-28 [2] CRAN (R 4.1.1)
#>  fastcluster            1.2.3      2021-05-24 [2] CRAN (R 4.1.1)
#>  fastmap                1.1.0      2021-01-25 [2] CRAN (R 4.1.1)
#>  fastmatch              1.1-3      2021-07-23 [2] CRAN (R 4.1.1)
#>  fgsea                  1.20.0     2021-10-26 [2] Bioconductor  
#>  foreach                1.5.1      2020-10-15 [2] CRAN (R 4.1.1)
#>  fs                     1.5.0      2020-07-31 [2] CRAN (R 4.1.1)
#>  future                 1.22.1     2021-08-25 [2] CRAN (R 4.1.1)
#>  future.apply           1.8.1      2021-08-10 [2] CRAN (R 4.1.1)
#>  genefilter             1.76.0     2021-10-26 [2] Bioconductor  
#>  geneplotter            1.72.0     2021-10-26 [2] Bioconductor  
#>  generics               0.1.1      2021-10-25 [2] CRAN (R 4.1.1)
#>  GenomeInfoDb         * 1.30.0     2021-10-26 [2] Bioconductor  
#>  GenomeInfoDbData       1.2.7      2021-09-23 [2] Bioconductor  
#>  GenomicRanges        * 1.46.0     2021-10-26 [2] Bioconductor  
#>  ggforce                0.3.3      2021-03-05 [2] CRAN (R 4.1.1)
#>  ggfun                  0.0.4      2021-09-17 [2] CRAN (R 4.1.1)
#>  ggplot2              * 3.3.5      2021-06-25 [2] CRAN (R 4.1.1)
#>  ggplotify              0.1.0      2021-09-02 [2] CRAN (R 4.1.1)
#>  ggraph                 2.0.5      2021-02-23 [2] CRAN (R 4.1.1)
#>  ggrepel                0.9.1      2021-01-15 [2] CRAN (R 4.1.1)
#>  ggtree                 3.2.0      2021-10-26 [2] Bioconductor  
#>  glmnet                 4.1-2      2021-06-24 [2] CRAN (R 4.1.1)
#>  globals                0.14.0     2020-11-22 [2] CRAN (R 4.1.1)
#>  glue                   1.4.2      2020-08-27 [2] CRAN (R 4.1.1)
#>  GO.db                  3.14.0     2021-09-23 [2] Bioconductor  
#>  GOSemSim               2.20.0     2021-10-26 [2] Bioconductor  
#>  gower                  0.2.2      2020-06-23 [2] CRAN (R 4.1.1)
#>  graph                  1.72.0     2021-10-26 [2] Bioconductor  
#>  graphite               1.40.0     2021-10-26 [2] Bioconductor  
#>  graphlayouts           0.7.1      2020-10-26 [2] CRAN (R 4.1.1)
#>  gridExtra              2.3        2017-09-09 [2] CRAN (R 4.1.1)
#>  gridGraphics           0.5-1      2020-12-13 [2] CRAN (R 4.1.1)
#>  gtable                 0.3.0      2019-03-25 [2] CRAN (R 4.1.1)
#>  highr                  0.9        2021-04-16 [2] CRAN (R 4.1.1)
#>  htmltools              0.5.2      2021-08-25 [2] CRAN (R 4.1.1)
#>  htmlwidgets            1.5.4      2021-09-08 [2] CRAN (R 4.1.1)
#>  httr                   1.4.2      2020-07-20 [2] CRAN (R 4.1.1)
#>  igraph                 1.2.7      2021-10-15 [2] CRAN (R 4.1.1)
#>  ipred                  0.9-12     2021-09-15 [2] CRAN (R 4.1.1)
#>  IRanges              * 2.28.0     2021-10-26 [2] Bioconductor  
#>  iterators              1.0.13     2020-10-15 [2] CRAN (R 4.1.1)
#>  jquerylib              0.1.4      2021-04-26 [2] CRAN (R 4.1.1)
#>  jsonlite               1.7.2      2020-12-09 [2] CRAN (R 4.1.1)
#>  KEGGREST               1.34.0     2021-10-26 [2] Bioconductor  
#>  knitr                  1.36       2021-09-29 [2] CRAN (R 4.1.1)
#>  labeling               0.4.2      2020-10-20 [2] CRAN (R 4.1.1)
#>  lattice              * 0.20-45    2021-09-22 [2] CRAN (R 4.1.1)
#>  lava                   1.6.10     2021-09-02 [2] CRAN (R 4.1.1)
#>  lazyeval               0.2.2      2019-03-15 [2] CRAN (R 4.1.1)
#>  lifecycle              1.0.1      2021-09-24 [2] CRAN (R 4.1.1)
#>  listenv                0.8.0      2019-12-05 [2] CRAN (R 4.1.1)
#>  locfit               * 1.5-9.4    2020-03-25 [2] CRAN (R 4.1.1)
#>  lubridate              1.8.0      2021-10-07 [2] CRAN (R 4.1.1)
#>  magrittr               2.0.1      2020-11-17 [2] CRAN (R 4.1.1)
#>  MASS                   7.3-54     2021-05-03 [2] CRAN (R 4.1.1)
#>  Matrix                 1.3-4      2021-06-01 [2] CRAN (R 4.1.1)
#>  MatrixGenerics       * 1.6.0      2021-10-26 [2] Bioconductor  
#>  matrixStats          * 0.61.0     2021-09-17 [2] CRAN (R 4.1.1)
#>  memoise                2.0.0      2021-01-26 [2] CRAN (R 4.1.1)
#>  ModelMetrics           1.2.2.2    2020-03-17 [2] CRAN (R 4.1.1)
#>  munsell                0.5.0      2018-06-12 [2] CRAN (R 4.1.1)
#>  networkD3            * 0.4        2017-03-18 [2] CRAN (R 4.1.1)
#>  nlme                   3.1-153    2021-09-07 [2] CRAN (R 4.1.1)
#>  nnet                   7.3-16     2021-05-03 [2] CRAN (R 4.1.1)
#>  org.Hs.eg.db         * 3.14.0     2021-09-23 [2] Bioconductor  
#>  parallelly             1.28.1     2021-09-09 [2] CRAN (R 4.1.1)
#>  patchwork              1.1.1      2020-12-17 [2] CRAN (R 4.1.1)
#>  pillar                 1.6.4      2021-10-18 [2] CRAN (R 4.1.1)
#>  pkgbuild               1.2.0      2020-12-15 [2] CRAN (R 4.1.1)
#>  pkgconfig              2.0.3      2019-09-22 [2] CRAN (R 4.1.1)
#>  pkgload                1.2.3      2021-10-13 [2] CRAN (R 4.1.1)
#>  plyr                   1.8.6      2020-03-03 [2] CRAN (R 4.1.1)
#>  png                    0.1-7      2013-12-03 [2] CRAN (R 4.1.1)
#>  polyclip               1.10-0     2019-03-14 [2] CRAN (R 4.1.1)
#>  prettyunits            1.1.1      2020-01-24 [2] CRAN (R 4.1.1)
#>  pROC                   1.18.0     2021-09-03 [2] CRAN (R 4.1.1)
#>  processx               3.5.2      2021-04-30 [2] CRAN (R 4.1.1)
#>  prodlim                2019.11.13 2019-11-17 [2] CRAN (R 4.1.1)
#>  proxy                  0.4-26     2021-06-07 [2] CRAN (R 4.1.1)
#>  ps                     1.6.0      2021-02-28 [2] CRAN (R 4.1.1)
#>  purrr                  0.3.4      2020-04-17 [2] CRAN (R 4.1.1)
#>  qvalue                 2.26.0     2021-10-26 [2] Bioconductor  
#>  R6                     2.5.1      2021-08-19 [2] CRAN (R 4.1.1)
#>  rappdirs               0.3.3      2021-01-31 [2] CRAN (R 4.1.1)
#>  RColorBrewer           1.1-2      2014-12-07 [2] CRAN (R 4.1.1)
#>  Rcpp                   1.0.7      2021-07-07 [2] CRAN (R 4.1.1)
#>  RcppArmadillo          0.10.7.0.0 2021-09-30 [2] CRAN (R 4.1.1)
#>  RcppParallel           5.1.4      2021-05-04 [2] CRAN (R 4.1.1)
#>  RCurl                  1.98-1.5   2021-09-17 [2] CRAN (R 4.1.1)
#>  reactome.db            1.77.0     2021-10-07 [2] Bioconductor  
#>  ReactomePA           * 1.38.0     2021-10-26 [2] Bioconductor  
#>  recipes                0.1.17     2021-09-27 [2] CRAN (R 4.1.1)
#>  remotes                2.4.1      2021-09-29 [2] CRAN (R 4.1.1)
#>  reshape2               1.4.4      2020-04-09 [2] CRAN (R 4.1.1)
#>  rlang                  0.4.12     2021-10-18 [2] CRAN (R 4.1.1)
#>  rmarkdown              2.11       2021-09-14 [2] CRAN (R 4.1.1)
#>  rpart                  4.1-15     2019-04-12 [2] CRAN (R 4.1.1)
#>  rprojroot              2.0.2      2020-11-15 [2] CRAN (R 4.1.1)
#>  RSQLite                2.2.8      2021-08-21 [2] CRAN (R 4.1.1)
#>  Rtsne                  0.15       2018-11-10 [2] CRAN (R 4.1.1)
#>  S4Vectors            * 0.32.0     2021-10-26 [2] Bioconductor  
#>  sass                   0.4.0      2021-05-12 [2] CRAN (R 4.1.1)
#>  scales                 1.1.1      2020-05-11 [2] CRAN (R 4.1.1)
#>  scatterpie             0.1.7      2021-08-20 [2] CRAN (R 4.1.1)
#>  scGPS                * 1.8.0      2021-10-26 [1] Bioconductor  
#>  sessioninfo            1.1.1      2018-11-05 [2] CRAN (R 4.1.1)
#>  shadowtext             0.0.9      2021-09-19 [2] CRAN (R 4.1.1)
#>  shape                  1.4.6      2021-05-19 [2] CRAN (R 4.1.1)
#>  SingleCellExperiment * 1.16.0     2021-10-26 [2] Bioconductor  
#>  stringi                1.7.5      2021-10-04 [2] CRAN (R 4.1.1)
#>  stringr                1.4.0      2019-02-10 [2] CRAN (R 4.1.1)
#>  SummarizedExperiment * 1.24.0     2021-10-26 [2] Bioconductor  
#>  survival               3.2-13     2021-08-24 [2] CRAN (R 4.1.1)
#>  testthat               3.1.0      2021-10-04 [2] CRAN (R 4.1.1)
#>  tibble                 3.1.5      2021-09-30 [2] CRAN (R 4.1.1)
#>  tidygraph              1.2.0      2020-05-12 [2] CRAN (R 4.1.1)
#>  tidyr                  1.1.4      2021-09-27 [2] CRAN (R 4.1.1)
#>  tidyselect             1.1.1      2021-04-30 [2] CRAN (R 4.1.1)
#>  tidytree               0.3.5      2021-09-08 [2] CRAN (R 4.1.1)
#>  timeDate               3043.102   2018-02-21 [2] CRAN (R 4.1.1)
#>  treeio                 1.18.0     2021-10-26 [2] Bioconductor  
#>  tweenr                 1.0.2      2021-03-23 [2] CRAN (R 4.1.1)
#>  usethis                2.1.2      2021-10-25 [2] CRAN (R 4.1.1)
#>  utf8                   1.2.2      2021-07-24 [2] CRAN (R 4.1.1)
#>  vctrs                  0.3.8      2021-04-29 [2] CRAN (R 4.1.1)
#>  viridis                0.6.2      2021-10-13 [2] CRAN (R 4.1.1)
#>  viridisLite            0.4.0      2021-04-13 [2] CRAN (R 4.1.1)
#>  withr                  2.4.2      2021-04-18 [2] CRAN (R 4.1.1)
#>  xfun                   0.27       2021-10-18 [2] CRAN (R 4.1.1)
#>  XML                    3.99-0.8   2021-09-17 [2] CRAN (R 4.1.1)
#>  xtable                 1.8-4      2019-04-21 [2] CRAN (R 4.1.1)
#>  XVector                0.34.0     2021-10-26 [2] Bioconductor  
#>  yaml                   2.2.1      2020-02-01 [2] CRAN (R 4.1.1)
#>  yulab.utils            0.0.4      2021-10-09 [2] CRAN (R 4.1.1)
#>  zlibbioc               1.40.0     2021-10-26 [2] Bioconductor  
#> 
#> [1] /tmp/RtmpAskcKd/Rinst1e6eab3779946f
#> [2] /home/biocbuild/bbs-3.14-bioc/R/library