plae Analysis Examples

Welcome

This document contains some simple examples for how to use the advanced filtering options in https://plae.nei.nih.gov and a brief introduction to using the Seurat file to do custom tests.

Gene Focused

By “gene focused” we mean that:

1. You have a gene you are interested in
2. You want to learn more about where / when / what it is expressed in

In silico In Situ

For those who are most comfortable with stained in situ slices of the retina this visualization may be useful. The major cell types of the retina are laid on in rough anatomical positioning. The cell types are colored by intensity, with the brighter colors meaning that the gene is more highly expressed in that cell type. As an example we show expression of RHO (rhodopsin, a rod marker) and RPE65 (RPE marker).

UMAP - Table

If you are curious about a gene, then there are several ways you can learn about its retinal cell type expression patterning. We will use ATOH7, a transcription factor that regulates retinal ganglion development as our example gene.

The UMAP view is a two dimensional representation of the individual cells in the scEiaD. Cells that are closer together have more related gene expression profiles (and thus are likely to be similar cell types).

Let’s go the UMAP - Table viewer in plae:

Viz -> UMAP - Tables

Dark gray are cells which have no detectable ATOH7.

Yellow is the highest expression and dark purple is the lowest expression.

Show highest expressing cells

What if we want to see which cells have the highest expression? We can use the “Filter Gene Expression” slider to only show cells with expression above a log2(expression) value.

We see that the highest expressing cells are in the “center” before the branching happens.

Species Filtering

By default plae shows data for all organisms in the database (human, mouse, macaque).

If we only want to see ATOH7 expression in human data, then that is very easy with the powerful “Scatter Filter Category” and “Gene Filter On” sections.

Table Information

While the UMAP view is cool looking, it can’t show you everything….what if we want to know what kind of cells are expressing ATOH7?

We can have quantified information on where ATOH7 is expressed by Cell Type (predicted) (this is our machine learned cell type labels) and organism.

We see that 40 - 70 % of the mouse and human neurogenic cell type express ATOH7. In raw counts that is 10,117 of 14,502 total mouse neurogenic cells. They have an average expression of 2.15. You can sort or filter the table based on queries. If you wanted to see ATOH7 expression in the RGCs this is trivial to do by typing in the box below.

This shows us that ATOH7 expression seems to be dropping in the maturing/mature RGCs (and is much lower in the macaque) relative to the neurogenic population.

Study filtering

As scEiaD is constructed from publicly available datasets, you can also filter the data to only show results from a specific paper. This may be useful if you using the results from that paper and want to check or confirm a finding.

You can see information about the papers / studies in scEiaD by using the adjacent “Make Meta Table” section as follows:

We see that the Clark et al. 2019 study did Smart-seq2 and 10X across many developmental time points in mouse. They study_accession ID is SRP158081. We can use this ID to look at ATOH7 expression only in this study in both the UMAP view and the table view

Expression Plot

As we have a huge number of studies and samples, we can use this (for single cell data) unusual view: a boxplot! We can see how ATOH7 expression changes across celltype and study.

How do we get here?

Make a plot….that shows nothing?

We’ve entered ATOH7 as the gene to plot (1). We are faceting (splitting the plot into separate sub-plots) on Cell Type (predict) (2). We are coloring the data points by study_accession (each study’s average gene expression across the Cell Type (predict) is plotted separately) (3). But we see … nothing. Why?

That is because the Plot Height (400) is not high enough. The text is prioritized over the data, so they are hidden. As it is extremely difficult to “auto” pick the correct height, it was more straightforward to have the user pick it. Usually a value around 700 - 1200 will give a reasonable view.

Useable plot

So yes, now we can see the data.

Some cones have ATOH7 expression?

So each point is an independent study. We see high ATOH7 expression in the neurogenic population, across many studies. But we also see one of studies with ATOH7 expression in cones. The legend shows which colors correspond to which study.

One of the studies is a bright teal…probably SRP200599. We can confirm that by replotting the data with a filter that only shows study_accession SRP200599.

Yep, that is it. We can jump to the UMAP - Table view to pull up the metadata we have extracted about the study.

This is from Buenaventura et al.  and is a study that enriched early mouse cones. Some work suggests that loss of ATOH7 inhibits cone specification. These cones may be “early” cones or late neurogenic cells that are developing into cones.

Cell Type Focused

If you want to get a sense about what is present in scEiaD, then there are several tools you can use. For our example, we will be starting with the rods.

How many rods do we have?

These are published labels

How many rods do we have after the machine learning?

How many rods (predicted) do we have across organism?

How many rods (predicted) we we have across organism and study?

What genes are differentially expressed in the Rods?

And what does ROM1 expression look like in the UMAP?

Very high expression in the rods (and expressed in other cells too)

Advanced Stuff - Analysis in R

We provide the full data as seurat (v3) or anndata (scanpy) objects you can download for boutique analysis. Here we demonstrate how you can use the Seurat object to run a quick custom diff test on a sub-population of neurogenic cells. The Seurat object is 25 GB in size and likely will cause serious memory issues on a generic laptop.

library(Seurat)

## Attaching SeuratObject

library(tidyverse)

## Registered S3 method overwritten by 'cli':
##   method     from         
##   print.boxx spatstat.geom

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(SingleCellExperiment)

## Loading required package: SummarizedExperiment

## Loading required package: MatrixGenerics

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following object is masked from 'package:dplyr':
## 
##     count

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:tidyr':
## 
##     expand

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## The following object is masked from 'package:purrr':
## 
##     reduce

## Loading required package: GenomeInfoDb

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

## 
## Attaching package: 'SummarizedExperiment'

## The following object is masked from 'package:SeuratObject':
## 
##     Assays

## The following object is masked from 'package:Seurat':
## 
##     Assays

library(scran)

## Loading required package: scuttle

library(org.Hs.eg.db)

## Loading required package: AnnotationDbi

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:dplyr':
## 
##     select

## 

#system('wget -O ~/data/scEiaD_2022_02/scEiaD_all_seurat_v3.Rdata https://hpc.nih.gov/~mcgaugheyd/scEiaD/2022_03_22/scEiaD_all_seurat_v3.Rdata)
# the -O renames the file and puts it in the ~/data/scEiaD directory
load('~/data/scEiaD_2022_02/scEiaD_all_seurat_v3.Rdata')

ID neurogenic cells in different clusters

We see that the neurogenic cells are a small minority of these clusters. But perhaps still we can find some interesting differences between them.

scEiaD@meta.data %>% 
  group_by(cluster, CellType_predict) %>% 
  summarise(Count = n()) %>% 
  mutate(Perc = (Count / sum(Count)) * 100) %>% 
  filter(Perc > 5) %>% 
  filter(CellType_predict == 'Neurogenic Cell') %>% 
  arrange(-Count)

## `summarise()` has grouped output by 'cluster'. You can override using the `.groups` argument.

## # A tibble: 3 × 4
## # Groups:   cluster [3]
##   cluster CellType_predict Count  Perc
##     <dbl> <chr>            <int> <dbl>
## 1      62 Neurogenic Cell   6453 61.5 
## 2      30 Neurogenic Cell   4307  7.58
## 3      40 Neurogenic Cell   3003 10.9

The most common is cluster 62 and 40

Compare cluster 62 - Neurogenic vs cluster 40 - Neurogenic

This will take a few minutes to run. This is why we cannot offer on-demand custom diff testing as while it runs the whole app shuts down for everyone while it computes.

# Create new column with a pasted together cluster ID and a CellType
scEiaD <- AddMetaData(scEiaD, 
                      metadata = 
                        paste0(
                          scEiaD@meta.data$cluster, 
                          '_', 
                          scEiaD@meta.data$CellType_predict
                        ), 
                      col.name = 'clusterCT' )
# tell Seurat this new column is the default identity
Idents(scEiaD) <- scEiaD@meta.data$clusterCT
diff_test <- FindMarkers(scEiaD, ident.1 = '62_Neurogenic Cell', ident.2 = '40_Neurogenic Cell')
diff_test %>% arrange(p_val) %>% filter(abs(avg_log2FC) > 2) %>% head(5)

##                 p_val avg_log2FC pct.1 pct.2 p_val_adj
## ENSG00000198886     0       -Inf 0.950 0.935         0
## ENSG00000173404     0  -243.6912 0.640 0.240         0
## ENSG00000198034     0  -522.0720 0.872 0.691         0
## ENSG00000205927     0  -162.7358 0.503 0.091         0
## ENSG00000136938     0  -386.3705 0.873 0.640         0

We see the second hit (by p value and log2FC > 2) is INSM1 (ENSG00000173404). It was reported by Anne Morris’s group (https://iovs.arvojournals.org/article.aspx?articleid=2148558) that knocking down this gene leads to differentiation issues in the eye.

Compare macula cones vs peripheral cones

Let us now demonstrate how having a large annotated set of data can be of value in asking more specific questions, like “how do macula/fovea cones differ relative to peripheral cones.”

First, we need to see whether this is even a tractable question.

scEiaD@meta.data %>% 
  as_tibble() %>% 
  group_by(retina_region, organism, study_accession, CellType_predict) %>% 
  summarise(Count = n()) %>% 
  filter(CellType_predict %in% c("Cone")) %>% 
  filter(Count > 10, !is.na(retina_region))

## `summarise()` has grouped output by 'retina_region', 'organism', 'study_accession'. You can override using the `.groups` argument.

## # A tibble: 9 × 5
## # Groups:   retina_region, organism, study_accession [9]
##   retina_region organism            study_accession CellType_predict Count
##   <chr>         <chr>               <chr>           <chr>            <int>
## 1 Macula        Homo sapiens        EGAD00001006350 Cone               253
## 2 Macula        Homo sapiens        SRP151023       Cone                88
## 3 Macula        Homo sapiens        SRP194595       Cone                50
## 4 Macula        Homo sapiens        SRP222958       Cone                13
## 5 Macula        Homo sapiens        SRP255195       Cone               141
## 6 Peripheral    Homo sapiens        SRP151023       Cone               329
## 7 Peripheral    Homo sapiens        SRP222958       Cone                18
## 8 Peripheral    Homo sapiens        SRP255195       Cone                26
## 9 Peripheral    Macaca fascicularis SRP158528       Cone               811

Yes, we have several human studies with macula / peripheral cone cells.

Let’s demonstrate how we can quickly run a conservative differential expression test the leverages the many studies we have. The “traditional” scRNA diff tests often return insane numbers of highly differentially expressed genes (the test we ran aboves returns many many genes with a padj < 0.05), which is highly annoying as then you really have to plot each result individually to confirm that it does not look wacky. We are going to build a “pseudo bulk” experiment here where we sum all the counts by gene, cell type (in this case we are doing macula cones versus peripheral cones), and study. This “pseudo” count table can then be used in a classic differential testing system like edgeR, limma, or DESeq2 (we will use the later).

# Create new column with a pasted together retina region and a CellType
scEiaD <- AddMetaData(scEiaD, 
                      metadata = 
                        paste0(
                          scEiaD@meta.data$retina_region, 
                          '_', 
                          scEiaD@meta.data$CellType_predict,
                          '_',
                          scEiaD@meta.data$organism
                        ), 
                      col.name = 'regionCT')


# filter down scEiaD to only human and in semi-supported (more than 10 cones cells) studies
Idents(scEiaD) <- scEiaD@meta.data$regionCT
scEiaD__subset <- subset(scEiaD, idents = c('Macula_Cone_Homo sapiens','Peripheral_Cone_Homo sapiens'))
scEiaD__subset <- subset(scEiaD__subset, 
                         subset = study_accession %in% 
                           (scEiaD__subset@meta.data %>% group_by(study_accession) %>% 
                              summarise(Count = n()) %>% filter(Count > 10) %>% 
                              pull(study_accession)))
# create SCE object so can create a pseudobulk matrix to run DESeq2 on
sce <- as.SingleCellExperiment(scEiaD__subset)
summed <- scater::aggregateAcrossCells(sce, 
                                       ids=colData(sce)[,c("regionCT", "study_accession")])
# pull out the counts to build the DESeq2 object
mat <- assay(summed, 'counts')
colnames(mat) <-colData(summed) %>% as_tibble() %>% 
  mutate(names = glue::glue("{study_accession}_{retina_region}")) %>% pull(names)

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = round(mat),
                              colData = colData(summed),
                              design= ~ study_accession + retina_region) # study as covariate, testing against fovea / not fovea

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

dds <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

## fitting model and testing

deseq_res <- results(dds)

Wrinkle

scEiaD is built on the Ensembl gene id (for examples ENSG00000185527). That’s not very friendly to humans so I added a bit of code to link the Ensembl ID to the gene name and show the top 16 diff genes (padj < 0.05) as a table.

# get human gene names 
library(org.Hs.eg.db)
symbols <- mapIds(org.Hs.eg.db, keys = row.names(deseq_res), keytype = "ENSEMBL", column="SYMBOL")

## 'select()' returned 1:many mapping between keys and columns

region_table <- deseq_res %>% as_tibble(rownames = 'name') %>%  
  left_join(symbols %>% enframe()) %>% 
  dplyr::rename(Ensembl = name, Gene = value) %>% 
  relocate(Gene)

## Joining, by = "name"

region_table %>% 
  filter(padj < 0.05) %>% 
  arrange(pvalue) %>% DT::datatable()

Always look at the plots

We use the “ScaleData” function from Seurat to zero center and scale the counts for each dataset

PDE6G

Top gene more highly expressed in the macula compared to the periphery. Seems pretty consistently higher across all the studies.

scEiaD__subset <- ScaleData(scEiaD__subset)

## Centering and scaling data matrix

VlnPlot(scEiaD__subset, features = c('ENSG00000185527'), log = TRUE)

VlnPlot(scEiaD__subset, c('ENSG00000185527'), split.by = 'retina_region', group.by='study_accession', log = TRUE) 

## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

All genes

for (i in region_table %>% 
  filter(padj < 0.05) %>% 
  pull(Ensembl)){
  #print(i)
  gene_name <- region_table %>% filter(Ensembl == i) %>% pull(Gene)
  print(VlnPlot(scEiaD__subset, i, split.by = 'retina_region', group.by='study_accession', log = TRUE) + ggtitle(gene_name))
}

Output and Session Info

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DESeq2_1.34.0               org.Hs.eg.db_3.14.0        
##  [3] AnnotationDbi_1.56.1        scran_1.22.1               
##  [5] scuttle_1.4.0               SingleCellExperiment_1.16.0
##  [7] SummarizedExperiment_1.24.0 Biobase_2.54.0             
##  [9] GenomicRanges_1.46.0        GenomeInfoDb_1.30.0        
## [11] IRanges_2.28.0              S4Vectors_0.32.2           
## [13] BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
## [15] matrixStats_0.61.0          forcats_0.5.1              
## [17] stringr_1.4.0               dplyr_1.0.7                
## [19] purrr_0.3.4                 readr_2.0.2                
## [21] tidyr_1.1.4                 tibble_3.1.6               
## [23] ggplot2_3.3.5               tidyverse_1.3.1            
## [25] SeuratObject_4.0.4          Seurat_4.0.6               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                reticulate_1.22          
##   [3] tidyselect_1.1.1          RSQLite_2.2.8            
##   [5] htmlwidgets_1.5.4         grid_4.1.2               
##   [7] BiocParallel_1.28.0       Rtsne_0.15               
##   [9] munsell_0.5.0             ScaledMatrix_1.2.0       
##  [11] codetools_0.2-18          ica_1.0-2                
##  [13] DT_0.19                   statmod_1.4.36           
##  [15] future_1.23.0             miniUI_0.1.1.1           
##  [17] withr_2.4.2               colorspace_2.0-2         
##  [19] highr_0.9                 knitr_1.36               
##  [21] rstudioapi_0.13           ROCR_1.0-11              
##  [23] tensor_1.5                listenv_0.8.0            
##  [25] GenomeInfoDbData_1.2.7    polyclip_1.10-0          
##  [27] farver_2.1.0              bit64_4.0.5              
##  [29] parallelly_1.30.0         vctrs_0.3.8              
##  [31] generics_0.1.1            xfun_0.28                
##  [33] R6_2.5.1                  ggbeeswarm_0.6.0         
##  [35] rsvd_1.0.5                locfit_1.5-9.4           
##  [37] cachem_1.0.6              bitops_1.0-7             
##  [39] spatstat.utils_2.3-0      DelayedArray_0.20.0      
##  [41] assertthat_0.2.1          promises_1.2.0.1         
##  [43] scales_1.1.1              beeswarm_0.4.0           
##  [45] gtable_0.3.0              beachmat_2.10.0          
##  [47] globals_0.14.0            goftest_1.2-3            
##  [49] rlang_0.4.12              genefilter_1.76.0        
##  [51] splines_4.1.2             lazyeval_0.2.2           
##  [53] spatstat.geom_2.3-1       broom_0.7.10             
##  [55] yaml_2.2.1                reshape2_1.4.4           
##  [57] abind_1.4-5               modelr_0.1.8             
##  [59] crosstalk_1.2.0           backports_1.3.0          
##  [61] httpuv_1.6.3              tools_4.1.2              
##  [63] ellipsis_0.3.2            spatstat.core_2.3-2      
##  [65] jquerylib_0.1.4           RColorBrewer_1.1-2       
##  [67] ggridges_0.5.3            Rcpp_1.0.7               
##  [69] plyr_1.8.6                sparseMatrixStats_1.6.0  
##  [71] zlibbioc_1.40.0           RCurl_1.98-1.5           
##  [73] rpart_4.1-15              deldir_1.0-6             
##  [75] viridis_0.6.2             pbapply_1.5-0            
##  [77] cowplot_1.1.1             zoo_1.8-9                
##  [79] haven_2.4.3               ggrepel_0.9.1            
##  [81] cluster_2.1.2             fs_1.5.0                 
##  [83] magrittr_2.0.1            data.table_1.14.2        
##  [85] scattermore_0.7           lmtest_0.9-39            
##  [87] reprex_2.0.1              RANN_2.6.1               
##  [89] fitdistrplus_1.1-6        hms_1.1.1                
##  [91] patchwork_1.1.1           mime_0.12                
##  [93] evaluate_0.14             xtable_1.8-4             
##  [95] XML_3.99-0.8              readxl_1.3.1             
##  [97] gridExtra_2.3             scater_1.22.0            
##  [99] compiler_4.1.2            KernSmooth_2.23-20       
## [101] crayon_1.4.2              htmltools_0.5.2          
## [103] mgcv_1.8-38               later_1.3.0              
## [105] tzdb_0.2.0                geneplotter_1.72.0       
## [107] lubridate_1.8.0           DBI_1.1.1                
## [109] dbplyr_2.1.1              MASS_7.3-54              
## [111] Matrix_1.3-4              cli_3.1.0                
## [113] parallel_4.1.2            metapod_1.2.0            
## [115] igraph_1.2.8              pkgconfig_2.0.3          
## [117] plotly_4.10.0             spatstat.sparse_2.1-0    
## [119] xml2_1.3.2                annotate_1.72.0          
## [121] vipor_0.4.5               bslib_0.3.1              
## [123] dqrng_0.3.0               XVector_0.34.0           
## [125] rvest_1.0.2               digest_0.6.28            
## [127] sctransform_0.3.3         RcppAnnoy_0.0.19         
## [129] Biostrings_2.62.0         spatstat.data_2.1-2      
## [131] rmarkdown_2.11            cellranger_1.1.0         
## [133] leiden_0.3.9              uwot_0.1.11              
## [135] edgeR_3.36.0              DelayedMatrixStats_1.16.0
## [137] shiny_1.7.1               lifecycle_1.0.1          
## [139] nlme_3.1-153              jsonlite_1.7.2           
## [141] BiocNeighbors_1.12.0      viridisLite_0.4.0        
## [143] limma_3.50.0              fansi_0.5.0              
## [145] pillar_1.6.4              lattice_0.20-45          
## [147] KEGGREST_1.34.0           fastmap_1.1.0            
## [149] httr_1.4.4                survival_3.2-13          
## [151] glue_1.5.0                png_0.1-7                
## [153] bit_4.0.4                 bluster_1.4.0            
## [155] stringi_1.7.5             sass_0.4.0               
## [157] blob_1.2.2                BiocSingular_1.10.0      
## [159] memoise_2.0.0             irlba_2.3.5              
## [161] future.apply_1.8.1

save(scEiaD__subset,region_table, deseq_res, file = 'pseudoBulk_cone_region_files.Rdata' )
write_csv(region_table, file = 'pseudoBulk_cone_region_table.csv.gz' )