Investigating Quiescent OPCs, Como 06.24.24

1 Load data and necessary packages

Follow this pipeline to download data https://cole-trapnell-lab.github.io/monocle3/docs/clustering/ # 1.1 load libraries

library(monocle3)
library(Seurat)

#1.2 read in RDS file, directly downloaded from GEO website

opcs <- readRDS('/Users/christinacomo/Downloads/GSE249268_Final_aging_OPC_cds_mm106.rds')

1.3 Change monocle dataframes to match with format of Seurat objects

#Change 
count.mat <- assay(o
                   pcs)

Error: unexpected symbol in:
"count.mat <- assay(o
                   pcs"

1.4 Create Seurat object

2 Quality Control

2.1 Mitochondria, Features, and Counts

Assign mitochondrial percentage to data set and look at quality control metrics (number of RNA molecules per cell, number of genes per cell and mitochondrial percent)

opcs.so[["percent.mt"]] <- PercentageFeatureSet(opcs.so, pattern = "^MT-")
opcs.so[["percent.mt"]] <- PercentageFeatureSet(opcs.so, pattern = "^MT-")
VlnPlot(opcs.so, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
Warning in SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents = idents,  :
  All cells have the same value of percent.mt.

3 Process data

3.1 Follow the Seurat vignette to analyze data, https://satijalab.org/seurat/articles/pbmc3k_tutorial#run-non-linear-dimensional-reduction-umaptsne

# Because this is mostly a homogenous population, we are only going to use a variable feature count of 2,000
opcs.so <- NormalizeData(opcs.so)
opcs.so <- FindVariableFeatures(opcs.so, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(opcs.so)
opcs.so <- ScaleData(opcs.so, features = all.genes)
opcs.so <- RunPCA(opcs.so, features = VariableFeatures(object = opcs.so))

ElbowPlot(opcs.so)

opcs.so <- FindNeighbors(opcs.so, dims = 1:10)
opcs.so <- FindClusters(opcs.so, resolution = 0.3)

opcs.so <- RunUMAP(opcs.so, dims = 1:10)
DimPlot(opcs.so)

3.2 Save RDS files for local loading

3.3 change identities of seurat object by accessing the metadata to look at how the data is represented

4 Subset data to only include the ages P180 and 360

opcs.so.sub <- subset(opcs.so, idents = c('2_P180', '3_P360'))

opcs.so.sub <- RunPCA(opcs.so.sub)
opcs.so.sub <- FindNeighbors(opcs.so.sub, dims = 1:10)
opcs.so.sub <- FindClusters(opcs.so.sub, resolution = 0.3)
opcs.so.sub <- RunUMAP(opcs.so.sub, dims = 1:10)
DimPlot(opcs.so.sub)

saveRDS(opcs.so.sub, '/Users/christinacomo/Library/CloudStorage/OneDrive-TheUniversityofColoradoDenver/10xVisium/Ethan Hughes/RDSfiles/opcs.so.ages.rds')

opcs.so.sub <- readRDS('/Users/christinacomo/Library/CloudStorage/OneDrive-TheUniversityofColoradoDenver/10xVisium/Ethan Hughes/RDSfiles/opcs.so.ages.rds')

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Warning message:
R graphics engine version 16 is not supported by this version of RStudio. The Plots tab will be disabled until a newer version of RStudio is installed. 

4.1 Look at identities of subsetted clusters

This is the final_identity/labels that was created by the authors of the paper

5 Subset data to only look at quiescent OPCs

opcs.so.qui <- subset(opcs.so.sub, ident='quiescent_OPC')

Error in WhichCells.Seurat(object = x, cells = cells, idents = idents,  : 
  Cannot find the following identities in the object: quiescent_OPC

5.1 Run standard workflow

You do not have to rescale or normalize your data after subsetting but you must re-run PCA as you are chaning the UMAP space At the 0.1 resolution, we get 4 clusters – this looks good to me and does not appear to be over or under clustered. For reference I ran a few different resolutions, but I like 0.1 the best. It makes sense that it is a low number, as it is a fairly homogeneous population. making it any higher can create false positives or identification of clusters that are not real. However, whatever you think is best in terms of the biology and your expertise is good with me, a case could be made for resolution 0.2.

opcs.so.qui <- RunPCA(opcs.so.qui)

PC_ 1 
Positive:  ENSMUSG00000101939.1, ENSMUSG00000090125.3, ENSMUSG00000030605.15, ENSMUSG00000064215.13, ENSMUSG00000097971.3, ENSMUSG00000044349.16, ENSMUSG00000032060.10, ENSMUSG00000017390.15, ENSMUSG00000057098.14, ENSMUSG00000069049.11 
       ENSMUSG00000048001.7, ENSMUSG00000006356.10, ENSMUSG00000032766.9, ENSMUSG00000042249.11, ENSMUSG00000064365.1, ENSMUSG00000062825.15, ENSMUSG00000033208.7, ENSMUSG00000052684.4, ENSMUSG00000020486.19, ENSMUSG00000056966.7 
       ENSMUSG00000110246.1, ENSMUSG00000069045.11, ENSMUSG00000068457.14, ENSMUSG00000046756.10, ENSMUSG00000045672.15, ENSMUSG00000049265.7, ENSMUSG00000105942.1, ENSMUSG00000029816.10, ENSMUSG00000015149.14, ENSMUSG00000024395.9 
Negative:  ENSMUSG00000064360.1, ENSMUSG00000059291.15, ENSMUSG00000004980.16, ENSMUSG00000062078.15, ENSMUSG00000025203.6, ENSMUSG00000052146.16, ENSMUSG00000069662.5, ENSMUSG00000049517.8, ENSMUSG00000039630.10, ENSMUSG00000067288.13 
       ENSMUSG00000067336.6, ENSMUSG00000040225.15, ENSMUSG00000106106.2, ENSMUSG00000037805.15, ENSMUSG00000028234.6, ENSMUSG00000037563.15, ENSMUSG00000061787.15, ENSMUSG00000041841.8, ENSMUSG00000015829.13, ENSMUSG00000017404.12 
       ENSMUSG00000090862.3, ENSMUSG00000068823.10, ENSMUSG00000042082.7, ENSMUSG00000066551.12, ENSMUSG00000061477.4, ENSMUSG00000090137.8, ENSMUSG00000058799.16, ENSMUSG00000005610.17, ENSMUSG00000022283.15, ENSMUSG00000063870.12 
PC_ 2 
Positive:  ENSMUSG00000028649.21, ENSMUSG00000026131.20, ENSMUSG00000015222.18, ENSMUSG00000015829.13, ENSMUSG00000039202.12, ENSMUSG00000068748.7, ENSMUSG00000021614.16, ENSMUSG00000010505.16, ENSMUSG00000040225.15, ENSMUSG00000061601.15 
       ENSMUSG00000062078.15, ENSMUSG00000030077.11, ENSMUSG00000056966.7, ENSMUSG00000035954.10, ENSMUSG00000022263.10, ENSMUSG00000035109.14, ENSMUSG00000023951.17, ENSMUSG00000002028.13, ENSMUSG00000026872.18, ENSMUSG00000021188.14 
       ENSMUSG00000062151.13, ENSMUSG00000002265.17, ENSMUSG00000005871.15, ENSMUSG00000004980.16, ENSMUSG00000029231.15, ENSMUSG00000009470.17, ENSMUSG00000063077.15, ENSMUSG00000050310.9, ENSMUSG00000038708.10, ENSMUSG00000022762.18 
Negative:  ENSMUSG00000057863.6, ENSMUSG00000060636.14, ENSMUSG00000037742.14, ENSMUSG00000071415.6, ENSMUSG00000039001.12, ENSMUSG00000046330.10, ENSMUSG00000057322.12, ENSMUSG00000003429.11, ENSMUSG00000049775.16, ENSMUSG00000060126.14 
       ENSMUSG00000003970.9, ENSMUSG00000093674.7, ENSMUSG00000062006.12, ENSMUSG00000057841.5, ENSMUSG00000062997.6, ENSMUSG00000079641.3, ENSMUSG00000058600.13, ENSMUSG00000008683.16, ENSMUSG00000031320.9, ENSMUSG00000007892.8 
       ENSMUSG00000016252.7, ENSMUSG00000028081.6, ENSMUSG00000060938.14, ENSMUSG00000025290.17, ENSMUSG00000040952.16, ENSMUSG00000067274.10, ENSMUSG00000079523.8, ENSMUSG00000059070.16, ENSMUSG00000006333.14, ENSMUSG00000032518.6 
PC_ 3 
Positive:  ENSMUSG00000064339.1, ENSMUSG00000064337.1, ENSMUSG00000008668.15, ENSMUSG00000003545.3, ENSMUSG00000069049.11, ENSMUSG00000024608.11, ENSMUSG00000021250.13, ENSMUSG00000069045.11, ENSMUSG00000046330.10, ENSMUSG00000060938.14 
       ENSMUSG00000063457.14, ENSMUSG00000068457.14, ENSMUSG00000040952.16, ENSMUSG00000093674.7, ENSMUSG00000025290.17, ENSMUSG00000012848.15, ENSMUSG00000003970.9, ENSMUSG00000057841.5, ENSMUSG00000007892.8, ENSMUSG00000052684.4 
       ENSMUSG00000056673.14, ENSMUSG00000006333.14, ENSMUSG00000003429.11, ENSMUSG00000025794.9, ENSMUSG00000064360.1, ENSMUSG00000060126.14, ENSMUSG00000039001.12, ENSMUSG00000071415.6, ENSMUSG00000032399.8, ENSMUSG00000030744.13 
Negative:  ENSMUSG00000015656.17, ENSMUSG00000064358.1, ENSMUSG00000064357.1, ENSMUSG00000059291.15, ENSMUSG00000047215.14, ENSMUSG00000062647.16, ENSMUSG00000060036.14, ENSMUSG00000020460.15, ENSMUSG00000064365.1, ENSMUSG00000090862.3 
       ENSMUSG00000037563.15, ENSMUSG00000098274.7, ENSMUSG00000071658.5, ENSMUSG00000074884.12, ENSMUSG00000017404.12, ENSMUSG00000064354.1, ENSMUSG00000066551.12, ENSMUSG00000012405.16, ENSMUSG00000029614.13, ENSMUSG00000038400.15 
       ENSMUSG00000052146.16, ENSMUSG00000062328.8, ENSMUSG00000090137.8, ENSMUSG00000029878.6, ENSMUSG00000037805.15, ENSMUSG00000024176.10, ENSMUSG00000038274.13, ENSMUSG00000023019.12, ENSMUSG00000058546.8, ENSMUSG00000017778.14 
PC_ 4 
Positive:  ENSMUSG00000064358.1, ENSMUSG00000064357.1, ENSMUSG00000029580.14, ENSMUSG00000064354.1, ENSMUSG00000030551.14, ENSMUSG00000062647.16, ENSMUSG00000041841.8, ENSMUSG00000017404.12, ENSMUSG00000052146.16, ENSMUSG00000037805.15 
       ENSMUSG00000060036.14, ENSMUSG00000059291.15, ENSMUSG00000037742.14, ENSMUSG00000030605.15, ENSMUSG00000045136.6, ENSMUSG00000015656.17, ENSMUSG00000090862.3, ENSMUSG00000020427.11, ENSMUSG00000031604.10, ENSMUSG00000039323.18 
       ENSMUSG00000022351.14, ENSMUSG00000047261.9, ENSMUSG00000032518.6, ENSMUSG00000045294.11, ENSMUSG00000037894.13, ENSMUSG00000021250.13, ENSMUSG00000018293.4, ENSMUSG00000067288.13, ENSMUSG00000021273.11, ENSMUSG00000053129.5 
Negative:  ENSMUSG00000064339.1, ENSMUSG00000021614.16, ENSMUSG00000023019.12, ENSMUSG00000098973.1, ENSMUSG00000024063.13, ENSMUSG00000034000.15, ENSMUSG00000033208.7, ENSMUSG00000024395.9, ENSMUSG00000035202.8, ENSMUSG00000015829.13 
       ENSMUSG00000064337.1, ENSMUSG00000091475.3, ENSMUSG00000106106.2, ENSMUSG00000060279.14, ENSMUSG00000012422.14, ENSMUSG00000029819.6, ENSMUSG00000022708.16, ENSMUSG00000038400.15, ENSMUSG00000046008.8, ENSMUSG00000105843.1 
       ENSMUSG00000064365.1, ENSMUSG00000026204.15, ENSMUSG00000036856.4, ENSMUSG00000042312.9, ENSMUSG00000022129.4, ENSMUSG00000056966.7, ENSMUSG00000024176.10, ENSMUSG00000020598.16, ENSMUSG00000015149.14, ENSMUSG00000064215.13 
PC_ 5 
Positive:  ENSMUSG00000038400.15, ENSMUSG00000030342.8, ENSMUSG00000046008.8, ENSMUSG00000001493.9, ENSMUSG00000024063.13, ENSMUSG00000029819.6, ENSMUSG00000039202.12, ENSMUSG00000053279.8, ENSMUSG00000036752.4, ENSMUSG00000070348.5 
       ENSMUSG00000023951.17, ENSMUSG00000091475.3, ENSMUSG00000000184.12, ENSMUSG00000027962.14, ENSMUSG00000024395.9, ENSMUSG00000022708.16, ENSMUSG00000029878.6, ENSMUSG00000026547.15, ENSMUSG00000021613.9, ENSMUSG00000021614.16 
       ENSMUSG00000030077.11, ENSMUSG00000012422.14, ENSMUSG00000029178.14, ENSMUSG00000027956.11, ENSMUSG00000007041.9, ENSMUSG00000031613.9, ENSMUSG00000100303.3, ENSMUSG00000020591.11, ENSMUSG00000022129.4, ENSMUSG00000063229.15 
Negative:  ENSMUSG00000030605.15, ENSMUSG00000034275.18, ENSMUSG00000064339.1, ENSMUSG00000030551.14, ENSMUSG00000045672.15, ENSMUSG00000060279.14, ENSMUSG00000064349.1, ENSMUSG00000064365.1, ENSMUSG00000020460.15, ENSMUSG00000035202.8 
       ENSMUSG00000047215.14, ENSMUSG00000087365.7, ENSMUSG00000048001.7, ENSMUSG00000098973.1, ENSMUSG00000056966.7, ENSMUSG00000064337.1, ENSMUSG00000042249.11, ENSMUSG00000046607.6, ENSMUSG00000029608.10, ENSMUSG00000064351.1 
       ENSMUSG00000049265.7, ENSMUSG00000010505.16, ENSMUSG00000105843.1, ENSMUSG00000022231.10, ENSMUSG00000029816.10, ENSMUSG00000041773.8, ENSMUSG00000003545.3, ENSMUSG00000060860.8, ENSMUSG00000037664.13, ENSMUSG00000110246.1 

opcs.so.qui <- FindNeighbors(opcs.so.qui, dims = 1:10)

Computing nearest neighbor graph
Computing SNN

opcs.so.qui <- FindClusters(opcs.so.qui, resolution = 0.1)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 11537
Number of edges: 344557

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9396
Number of communities: 3
Elapsed time: 1 seconds

opcs.so.qui <- RunUMAP(opcs.so.qui, dims = 1:10)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
12:22:18 UMAP embedding parameters a = 0.9922 b = 1.112
12:22:18 Read 11537 rows and found 10 numeric columns
12:22:18 Using Annoy for neighbor search, n_neighbors = 30
12:22:18 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:22:19 Writing NN index file to temp file /var/folders/2b/gcm1xy2n6s50jqf9f3ln87140000gn/T//RtmpDwIazj/file69e376982f9d
12:22:19 Searching Annoy index using 1 thread, search_k = 3000
12:22:23 Annoy recall = 100%
12:22:23 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
12:22:24 Initializing from normalized Laplacian + noise (using RSpectra)
12:22:24 Commencing optimization for 200 epochs, with 461856 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:22:30 Optimization finished

DimPlot(opcs.so.qui)

6 Find markers from clusters

marks.1 <- FindMarers(opcs.so.qui, ident.1 = "1")

Error in FindMarers(opcs.so.qui, ident.1 = "1") : 
  could not find function "FindMarers"

6.1 You can see that the tables made with the gene markers are missing a column header above the ensembl ID and are not the gene names we want

#6.2 Write gene tables as a csv and edit in excel (easier than doing it here)

6.3 Read back in csv files with updated column headers

7 Change ensebml IDs to gene IDs

We can see that the ensembl IDs have the transcript information after it and we need to remove it to be able to convert to gene IDs #7.1 Write a function to remove transcript infomration for our dataframe

write.csv(marks.0, '/Users/christinacomo/Library/CloudStorage/OneDrive-TheUniversityofColoradoDenver/10xVisium/Ethan Hughes/outputs/markers0.csv')
write.csv(marks.1, '/Users/christinacomo/Library/CloudStorage/OneDrive-TheUniversityofColoradoDenver/10xVisium/Ethan Hughes/outputs/markers1.csv')
write.csv(marks.2, '/Users/christinacomo/Library/CloudStorage/OneDrive-TheUniversityofColoradoDenver/10xVisium/Ethan Hughes/outputs/markers2.csv')

7.2 Convert gene IDs

I used this https://www.biotools.fr/mouse/ensembl_symbol_converter you can use tools like Biomart, but this is faster and easier

7.3 Interactive table of top genes based on increased expression, refer to excel sheets for all

8 Session Info

sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Denver
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] org.Mm.eg.db_3.18.0  AnnotationDbi_1.64.1 IRanges_2.36.0       S4Vectors_0.40.2     Biobase_2.62.0       BiocGenerics_0.48.1 
 [7] biomaRt_2.58.0       Seurat_5.0.1         SeuratObject_5.0.1   sp_2.1-2            

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22        splines_4.3.1           later_1.3.2             filelock_1.0.3          bitops_1.0-7           
  [6] tibble_3.2.1            polyclip_1.10-6         XML_3.99-0.16           fastDummies_1.7.3       lifecycle_1.0.4        
 [11] globals_0.16.3          lattice_0.22-5          MASS_7.3-60.0.1         magrittr_2.0.3          limma_3.58.1           
 [16] plotly_4.10.4           sass_0.4.9              rmarkdown_2.26          jquerylib_0.1.4         yaml_2.3.8             
 [21] httpuv_1.6.15           sctransform_0.4.1       spam_2.10-0             spatstat.sparse_3.0-3   reticulate_1.34.0      
 [26] cowplot_1.1.3           pbapply_1.7-2           DBI_1.2.1               RColorBrewer_1.1-3      abind_1.4-5            
 [31] zlibbioc_1.48.0         Rtsne_0.17              purrr_1.0.2             presto_1.0.0            RCurl_1.98-1.14        
 [36] rappdirs_0.3.3          GenomeInfoDbData_1.2.11 ggrepel_0.9.5           irlba_2.3.5.1           listenv_0.9.1          
 [41] spatstat.utils_3.0-4    goftest_1.2-3           RSpectra_0.16-1         spatstat.random_3.2-3   fitdistrplus_1.1-11    
 [46] parallelly_1.37.1       leiden_0.4.3.1          codetools_0.2-19        xml2_1.3.6              tidyselect_1.2.1       
 [51] farver_2.1.1            BiocFileCache_2.10.1    matrixStats_1.2.0       spatstat.explore_3.2-5  jsonlite_1.8.8         
 [56] ellipsis_0.3.2          progressr_0.14.0        ggridges_0.5.6          survival_3.5-7          tools_4.3.1            
 [61] progress_1.2.3          ica_1.0-3               Rcpp_1.0.12             glue_1.7.0              gridExtra_2.3          
 [66] xfun_0.43               GenomeInfoDb_1.38.5     dplyr_1.1.4             withr_3.0.0             BiocManager_1.30.22    
 [71] fastmap_1.1.1           fansi_1.0.6             digest_0.6.35           R6_2.5.1                mime_0.12              
 [76] colorspace_2.1-0        scattermore_1.2         tensor_1.5              spatstat.data_3.0-4     RSQLite_2.3.5          
 [81] utf8_1.2.4              tidyr_1.3.0             generics_0.1.3          renv_1.0.3              data.table_1.14.10     
 [86] prettyunits_1.2.0       httr_1.4.7              htmlwidgets_1.6.4       uwot_0.1.16             pkgconfig_2.0.3        
 [91] gtable_0.3.4            blob_1.2.4              lmtest_0.9-40           XVector_0.42.0          htmltools_0.5.8        
 [96] dotCall64_1.1-1         scales_1.3.0            png_0.1-8               knitr_1.45              rstudioapi_0.15.0      
[101] reshape2_1.4.4          curl_5.2.1              nlme_3.1-164            zoo_1.8-12              cachem_1.0.8           
[106] stringr_1.5.1           KernSmooth_2.23-22      parallel_4.3.1          miniUI_0.1.1.1          pillar_1.9.0           
[111] grid_4.3.1              vctrs_0.6.5             RANN_2.6.1              promises_1.2.1          dbplyr_2.4.0           
[116] xtable_1.8-4            cluster_2.1.6           evaluate_0.23           cli_3.6.2               compiler_4.3.1         
[121] rlang_1.1.3             crayon_1.5.2            future.apply_1.11.1     labeling_0.4.3          plyr_1.8.9             
[126] stringi_1.8.3           viridisLite_0.4.2       deldir_2.0-4            munsell_0.5.0           Biostrings_2.70.1      
[131] lazyeval_0.2.2          spatstat.geom_3.2-9     Matrix_1.6-5            RcppHNSW_0.5.0          hms_1.1.3              
[136] patchwork_1.2.0         bit64_4.0.5             future_1.33.1           ggplot2_3.5.0           KEGGREST_1.42.0        
[141] statmod_1.5.0           shiny_1.8.0             ROCR_1.0-11             igraph_1.6.0            memoise_2.0.1          
[146] bslib_0.6.2             bit_4.0.5