OSCR

Traumatic brain injury induces cell-type-specific remodelling of the proteostasis network.

Code ↔ Paper

10 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 10 matches
  1. [1] § Methods › Data clustering ↔ Mouse Code.R, lines 291 to 380 · score 0.76 · MRtree, adjusted Rand, clustering tree, clustree, resolutions, Jaccard · Report a mismatch
  2. [2] § Methods › Data clustering ↔ Human Code.R, lines 172 to 216 · score 0.74 · MRtree, adjusted Rand, clustering tree, clustree, resolutions, human · Report a mismatch
  3. [3] § Methods › Transcriptomic data pre-processing ↔ Mouse Code.R, lines 62 to 97 · score 0.69 · biomaRt, gene mapping, human genes, fewer, filtered, mouse · Report a mismatch
  4. [4] § Results › Identifying preferential gene expression in specific cell types ↔ 1A basic_features.R, lines 110 to 156 · score 0.60 · GM2A, N4BP1, UCHL1, A2M, KLHL2, naive · Report a mismatch
  5. [5] § Results › A core sub-proteostasis network in traumatic brain injury ↔ 1A basic_features.R, lines 158 to 218 · score 0.56 · Core Naive Sub, Core CCI Sub, TRANSL, ALP, TBI, PN · Report a mismatch
  6. [6] § Methods › Transcriptomic data pre-processing ↔ Human Code.R, lines 119 to 170 · score 0.55 · NormalizeData, ScaleData, PN genes, pre, Seurat, PCA · Report a mismatch
  7. [7] § Results › A core sub-proteostasis network in traumatic brain injury ↔ Mouse Code.R, lines 1303 to 1366 · score 0.54 · Core Naive Sub, Core CCI Sub, TRANSL, ALP, mice, PN · Report a mismatch
  8. [8] § Methods › Gene set enrichment analysis ↔ Mouse Code.R, lines 769 to 810 · score 0.52 · log2FC, FindAllMarkers, fgsea, Seurat, mouse, clusters · Report a mismatch
  9. [9] § Methods › Gene set enrichment analysis ↔ Human Code.R, lines 706 to 789 · score 0.52 · log2FC, FindAllMarkers, fgsea, human, TBI, clusters · Report a mismatch
  10. [10] § Methods › Proteostasis network and Hallmark Pathway association ↔ Figure S1.R, lines 44 to 129 · score 0.50 · enrichment odds ratio, Fisher, PN · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,395 lines · 81 KB · no license · 4 matches

  1. suppressMessages(library(Seurat))
  2. suppressMessages(library(biomaRt))
  3. suppressMessages(library(presto))
  4. suppressMessages(library(ggplot2))
  5. suppressMessages(library(RColorBrewer))
  6. suppressMessages(library(clustree))
  7. suppressMessages(library(mclust))
  8. suppressMessages(library(mrtree))
  9. suppressMessages(library(igraph))
  10. suppressMessages(library(fgsea))
  11. suppressMessages(library(data.table))
  12. suppressMessages(library(ComplexHeatmap))
  13. suppressMessages(library(DESeq2))
  14. suppressMessages(library(ggVennDiagram))
  17. mem.maxVSize(vsize=Inf)
  18. options(future.globals.maxSize = 1000 * 1024^2)
  20. ####session information####
  21. # R version 4.4.1 (2024-06-14)
  22. # Platform: aarch64-apple-darwin20
  23. # Running under: macOS Sonoma 14.5
  24. #
  25. # Matrix products: default
  26. # BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
  27. # LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
  28. #
  29. # locale:
  30. # [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
  31. #
  32. # time zone: America/New_York
  33. # tzcode source: internal
  34. #
  35. # attached base packages:
  36. # [1] stats4 grid stats graphics grDevices utils datasets methods base
  37. #
  38. # other attached packages:
  39. # [1] DESeq2_1.44.0 SummarizedExperiment_1.34.0 Biobase_2.64.0 MatrixGenerics_1.16.0 matrixStats_1.5.0 GenomicRanges_1.56.2 GenomeInfoDb_1.40.1 IRanges_2.38.1 S4Vectors_0.42.1 BiocGenerics_0.50.0 ComplexHeatmap_2.20.0
  40. # [12] fgsea_1.30.0 igraph_2.2.2 mrtree_0.0.0.9000 mclust_6.1.2 clustree_0.5.1 ggraph_2.2.2 RColorBrewer_1.1-3 ggplot2_4.0.2 presto_1.0.0 data.table_1.18.2.1 Rcpp_1.1.1
  41. # [23] biomaRt_2.60.1 Seurat_5.4.0 SeuratObject_5.3.0 sp_2.2-1
  42. #
  43. # loaded via a namespace (and not attached):
  44. # [1] fs_1.6.6 spatstat.sparse_3.1-0 httr_1.4.8 doParallel_1.0.17 numDeriv_2016.8-1.1 tools_4.4.1 sctransform_0.4.3 backports_1.5.0 R6_2.6.1 lazyeval_0.2.2 uwot_0.2.4
  45. # [12] GetoptLong_1.1.0 withr_3.0.2 prettyunits_1.2.0 gridExtra_2.3 progressr_0.18.0 cli_3.6.5 spatstat.explore_3.7-0 fastDummies_1.7.5 S7_0.2.1 spatstat.data_3.1-9 ggridges_0.5.7
  46. # [23] pbapply_1.7-4 systemfonts_1.3.1 yulab.utils_0.2.4 DEoptim_2.2-8 parallelly_1.46.1 maps_3.4.3 rstudioapi_0.18.0 RSQLite_2.4.6 optimParallel_1.0-2 generics_0.1.4 gridGraphics_0.5-1
  47. # [34] shape_1.4.6.1 combinat_0.0-8 ica_1.0-3 spatstat.random_3.4-4 SymSim_0.0.0.9000 dendextend_1.19.1 dplyr_1.2.0 Matrix_1.7-4 abind_1.4-8 lifecycle_1.0.5 yaml_2.3.12
  48. # [45] scatterplot3d_0.3-45 clusterGeneration_1.3.8 SparseArray_1.4.8 BiocFileCache_2.12.0 Rtsne_0.17 blob_1.3.0 promises_1.5.0 crayon_1.5.3 miniUI_0.1.2 lattice_0.22-9 cowplot_1.2.0
  49. # [56] KEGGREST_1.44.1 magick_2.9.1 pillar_1.11.1 rjson_0.2.23 future.apply_1.20.2 codetools_0.2-20 fastmatch_1.1-8 glue_1.8.0 ggiraph_0.9.6 ggimage_0.3.5 ggfun_0.2.0
  50. # [67] spatstat.univar_3.1-6 fontLiberation_0.1.0 treeio_1.28.0 vctrs_0.7.1 png_0.1-8 spam_2.11-3 gtable_0.3.6 cachem_1.1.0 S4Arrays_1.4.1 mime_0.13 tidygraph_1.3.1
  51. # [78] coda_0.19-4.1 survival_3.8-6 SingleCellExperiment_1.26.0 iterators_1.0.14 fitdistrplus_1.2-6 ROCR_1.0-12 nlme_3.1-168 ggtree_3.12.0 bit64_4.6.0-1 fontquiver_0.2.1 progress_1.2.3
  52. # [89] filelock_1.0.3 RcppAnnoy_0.0.23 data.tree_1.2.0 irlba_2.3.7 KernSmooth_2.23-26 otel_0.2.0 colorspace_2.1-2 DBI_1.3.0 phangorn_2.12.1 mnormt_2.1.2 tidyselect_1.2.1
  53. # [100] bit_4.6.0 compiler_4.4.1 curl_7.0.0 httr2_1.2.2 expm_1.0-0 xml2_1.5.2 fontBitstreamVera_0.1.1 DelayedArray_0.30.1 plotly_4.12.0 checkmate_2.3.4 scales_1.4.0
  54. # [111] lmtest_0.9-40 quadprog_1.5-8 rappdirs_0.3.4 stringr_1.6.0 digest_0.6.39 goftest_1.2-3 spatstat.utils_3.2-1 XVector_0.44.0 htmltools_0.5.9 pkgconfig_2.0.3 dbplyr_2.5.1
  55. # [122] fastmap_1.2.0 rlang_1.1.7 GlobalOptions_0.1.3 htmlwidgets_1.6.4 UCSC.utils_1.0.0 shiny_1.13.0 farver_2.1.2 zoo_1.8-15 jsonlite_2.0.0 BiocParallel_1.38.0 magrittr_2.0.4
  56. # [133] GenomeInfoDbData_1.2.12 ggplotify_0.1.3 dotCall64_1.2 patchwork_1.3.2 ape_5.8-1 viridis_0.6.5 gdtools_0.5.0 reticulate_1.45.0 stringi_1.8.7 zlibbioc_1.50.0 MASS_7.3-65
  57. # [144] plyr_1.8.9 parallel_4.4.1 listenv_0.10.0 ggrepel_0.9.6 deldir_2.0-4 Biostrings_2.72.1 graphlayouts_1.2.3 splines_4.4.1 tensor_1.5.1 hms_1.1.4 circlize_0.4.17
  58. # [155] locfit_1.5-9.12 spatstat.geom_3.7-0 RcppHNSW_0.6.0 reshape2_1.4.5 foreach_1.5.2 tweenr_2.0.3 httpuv_1.6.16 RANN_2.6.2 tidyr_1.3.2 purrr_1.2.1 polyclip_1.10-7
  59. # [166] clue_0.3-67 future_1.69.0 scattermore_1.2 ggforce_0.5.0 xtable_1.8-8 tidytree_0.4.7 phytools_2.5-2 RSpectra_0.16-2 later_1.4.7 viridisLite_0.4.3 tibble_3.3.1
  60. # [177] aplot_0.2.9 memoise_2.0.1 AnnotationDbi_1.66.0 cluster_2.1.8.2 globals_0.19.0
  62. #####PREPARE OBJECT#####
  64. #For gene ID mapping
  65. mousegenes <- useEnsembl(biomart = "ensembl", dataset = "mmusculus_gene_ensembl", version = 105)
  66. humangenes <- useEnsembl(biomart = "ensembl", dataset = "hsapiens_gene_ensembl", version = 105)
  68. load('all_SeuratObjects.RData') #this has the data
  70. expression_matrix <- Int_all_4C.combined[['RNA']]$counts #extract the counts for all conditions. At this stage we have 19558 genes.
  72. #Do gene mapping
  73. gene_name_map <- unique(getLDS(attributes = c("external_gene_name"),
  74. filters = "external_gene_name",
  75. values = rownames(expression_matrix),
  76. mart = mousegenes,
  77. attributesL = c("hgnc_symbol"),
  78. martL = humangenes, uniqueRows=T)) #make mapping table from biomaRt
  79. mouse_count <- as.data.frame(table(gene_name_map$Gene.name)) #counting frequency to identify duplicates in mouse names (one to many mapping)
  80. human_count <- as.data.frame(table(gene_name_map$HGNC.symbol)) #counting frequency to identify duplicates in mouse names (many to one mapping)
  81. gene_name_map$mouse_counts <- mouse_count$Freq[match(unlist(gene_name_map$Gene.name), mouse_count$Var1)] #mouse name frequencies added to data frame
  82. gene_name_map$human_counts <- human_count$Freq[match(unlist(gene_name_map$HGNC.symbol), human_count$Var1)] #human name frequencies added to data frame
  83. equals_one <- gene_name_map$mouse_counts == 1 & gene_name_map$human_counts == 1 #create boolean for one-to-one mapping
  84. gene_name_map <- gene_name_map[equals_one,] #filter to keep only one-to-one mapping
  86. in_map <- rownames(expression_matrix) %in% gene_name_map$Gene.name #which gene names in the expression matrix are mappable?
  87. expression_matrix <- expression_matrix[in_map,] #use only mappable genes. At this stage we have 14348 genes.
  89. rownames(expression_matrix) <- gene_name_map$HGNC.symbol[match(unlist(rownames(expression_matrix)), gene_name_map$Gene.name)] #convert from mouse to human gene names
  91. #But we only want to include two conditions, Naive and CCI
  92. conditions_of_interest <- c('CCI','Naive')
  93. metadata_of_interest <- subset([email hidden],condition %in% conditions_of_interest)
  94. expression_matrix <- expression_matrix[,rownames(metadata_of_interest)] #at this stage we have 14348 genes
  96. #Now we want to get rid of any genes that have fewer than 10 cells with at least count 1 in each sample.
  97. origidents <- unique(metadata_of_interest$orig.ident)
  99. #read_count_function returns the number of values in a supplied row or column that have a non-zero account. If you provide a row from the expression matrix, it's the number of cells that have a count of at least 1 for the gene. If you provide a column from the expression matrix, it's the number of genes that have a non-zero count for the cell. But I don't use the latter because I use the built-in feature of CreateSeuratObject below.
  100. read_count_function <- function(row_or_column) {
  101. the_table <- table(row_or_column)
  102. zero_count <- as.numeric(the_table[names(the_table)==0])
  103. non_zero_count <- length(row_or_column)-zero_count
  104. return(non_zero_count)
  105. }
  107. #Get rid of genes that are expressed in fewer than 10 cells per sample.
  108. too_sparse_genes_list = c()
  109. for (origident in origidents) { #for each of the idents that we are including (each naive and cci replicate)
  110. print(origident)
  111. temp_meta <- subset([email hidden], orig.ident %in% origident) #identify the subset of cells that match that particular replicate
  112. exp_temp <- expression_matrix[,rownames(temp_meta)] #subset the expression matrix for only those cells
  113. non_zero_count_list <- apply(exp_temp, 1, read_count_function) #use the count_function to get the # of cells that have a non-zero value for each gene. Supplying a row.
  114. too_sparse <- names(subset(non_zero_count_list, non_zero_count_list < 10)) #get the names of the genes that have fewer than 5 cells with non-zero values
  115. too_sparse_genes_list <- union(too_sparse_genes_list, too_sparse) #add that to a list of genes that are too sparsley expressed to be handled by default SCTransform settings
  116. }
  118. keep_genes <- subset(rownames(expression_matrix), !(rownames(expression_matrix) %in% too_sparse_genes_list))
  120. expression_matrix <- expression_matrix[keep_genes,] #11371 genes, 52384 cells
  122. ##Now we want to get rid of any cell that has fewer than 200 features. We set min.cells to 0 because we've already done this filtration on a sample by sample basis.
  123. cci_seurat_pre <- CreateSeuratObject(counts=expression_matrix, min.cells=0, min.features=200, meta.data=metadata_of_interest) #11371 genes, 52368 cells
  125. #Now I want to create a separate expression_matrix only for PN genes
  127. PN = read.table('proteostasis_list.txt',header=FALSE,sep='\n')$V1 #PN gene list
  128. in_PN <- rownames(cci_seurat_pre[['RNA']]$counts) %in% PN #list of genes in the expression_matrix that are in the PN
  129. expression_matrix_PN <- cci_seurat_pre[['RNA']]$counts[in_PN,] #include only PN genes #this is 2420 genes
  131. ##We set min.cells to 0 because we've already done this filtration on a sample by sample basis. We don't get rid of additional features here. This could potentially create the issue of cells that have a very low number of PN only features.
  132. ##You can test this out.
  133. cci_seurat_pre_PN <- CreateSeuratObject(counts=expression_matrix_PN, min.cells=0, min.features=0, meta.data=[email hidden])
  134. #2420 genes, 52368 cells. By empirically testing, we found that every cell expresses at least 48 PN genes.
  136. #now we prep the cci_seurat_pre object
  137. #Add meta data with the full cell type names
  138. lookup <- c("0" = "Microglia_1", "1" = "Microglia_2", "2" = "Ependymal_1", "3" = "Microglia_3",
  139. "4" = "Ependymal_2", "5" = "BCells", "6" = "Endothelial_1", "7" = "Neutrophils_1",
  140. "8" = "Microglia_4", "9" = "Microglia_5", "10" = "Ependymal_3", "11" = "Microglia_6",
  141. "12" = "Endothelial_2", "13" = "Astrocytes", "14" = "NKcells", "15" = "Neurons",
  142. "16" = "Neutrophils_2", "17" = "Ependymal_4", "18" = "Microglia_7", "19" = "Oligodendrocytes",
  143. "20" = "TmemoryCells", "21" = "Microglia_8", "22" = "Ependymal_5") #From Jha group, these are the cell types from their clusters
  145. celltype_col <- data.frame(celltype_col=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  146. rownames(celltype_col) <- rownames([email hidden]) #add the cell row names there
  147. cci_seurat_pre <- AddMetaData(cci_seurat_pre, celltype_col) #add this to the meta data
  149. #Add meta data with abbreviated cell type names
  150. lookup <- c("0" = "M1", "1" = "M2", "2" = "Ep1", "3" = "M3",
  151. "4" = "Ep2", "5" = "B", "6" = "En1", "7" = "Nt1",
  152. "8" = "M4", "9" = "M5", "10" = "Ep3", "11" = "M6",
  153. "12" = "En2", "13" = "As", "14" = "NK", "15" = "Nr",
  154. "16" = "Nt2", "17" = "Ep4", "18" = "M7", "19" = "Olg",
  155. "20" = "Tmem", "21" = "M8", "22" = "Ep5") #From Jha group, these are the cell types from their clusters
  157. celltype_lab <- data.frame(celltype_lab=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  158. rownames(celltype_lab) <- rownames([email hidden]) #add the cell row names there
  159. cci_seurat_pre <- AddMetaData(cci_seurat_pre, celltype_lab) #add this to the meta data
  161. #Add meta data with larger cell type categories, as defined in Jha Neuron 1B
  162. lookup <- c("0" = "M", "1" = "M", "2" = "Ep", "3" = "M",
  163. "4" = "Ep", "5" = "L", "6" = "En", "7" = "Nt",
  164. "8" = "M", "9" = "M", "10" = "Ep", "11" = "M",
  165. "12" = "En", "13" = "NG", "14" = "L", "15" = "NG",
  166. "16" = "Nt", "17" = "Ep", "18" = "M", "19" = "NG",
  167. "20" = "L", "21" = "M", "22" = "Ep") #From Jha group, these are the cell types from their clusters
  169. celltype_cat <- data.frame(celltype_cat=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  170. rownames(celltype_cat) <- rownames([email hidden]) #add the cell row names there
  171. cci_seurat_pre <- AddMetaData(cci_seurat_pre, celltype_cat) #add this to the meta data
  173. saveRDS(cci_seurat_pre, file="cci_seurat_pre.rds")
  175. #now we prep the cci_seurat_pre_PN object
  176. #Add meta data with the full cell type names
  177. lookup <- c("0" = "Microglia_1", "1" = "Microglia_2", "2" = "Ependymal_1", "3" = "Microglia_3",
  178. "4" = "Ependymal_2", "5" = "BCells", "6" = "Endothelial_1", "7" = "Neutrophils_1",
  179. "8" = "Microglia_4", "9" = "Microglia_5", "10" = "Ependymal_3", "11" = "Microglia_6",
  180. "12" = "Endothelial_2", "13" = "Astrocytes", "14" = "NKcells", "15" = "Neurons",
  181. "16" = "Neutrophils_2", "17" = "Ependymal_4", "18" = "Microglia_7", "19" = "Oligodendrocytes",
  182. "20" = "TmemoryCells", "21" = "Microglia_8", "22" = "Ependymal_5") #From Jha group, these are the cell types from their clusters
  184. celltype_col <- data.frame(celltype_col=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  185. rownames(celltype_col) <- rownames([email hidden]) #add the cell row names there
  186. cci_seurat_pre_PN <- AddMetaData(cci_seurat_pre_PN, celltype_col) #add this to the meta data
  188. #Add meta data with abbreviated cell type names
  189. lookup <- c("0" = "M1", "1" = "M2", "2" = "Ep1", "3" = "M3",
  190. "4" = "Ep2", "5" = "B", "6" = "En1", "7" = "Nt1",
  191. "8" = "M4", "9" = "M5", "10" = "Ep3", "11" = "M6",
  192. "12" = "En2", "13" = "As", "14" = "NK", "15" = "Nr",
  193. "16" = "Nt2", "17" = "Ep4", "18" = "M7", "19" = "Olg",
  194. "20" = "Tmem", "21" = "M8", "22" = "Ep5") #From Jha group, these are the cell types from their clusters
  196. celltype_lab <- data.frame(celltype_lab=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  197. rownames(celltype_lab) <- rownames([email hidden]) #add the cell row names there
  198. cci_seurat_pre_PN <- AddMetaData(cci_seurat_pre_PN, celltype_lab) #add this to the meta data
  200. #Add meta data with larger cell type categories, as defined in Jha Neuron 1B
  201. lookup <- c("0" = "M", "1" = "M", "2" = "Ep", "3" = "M",
  202. "4" = "Ep", "5" = "L", "6" = "En", "7" = "Nt",
  203. "8" = "M", "9" = "M", "10" = "Ep", "11" = "M",
  204. "12" = "En", "13" = "NG", "14" = "L", "15" = "NG",
  205. "16" = "Nt", "17" = "Ep", "18" = "M", "19" = "NG",
  206. "20" = "L", "21" = "M", "22" = "Ep") #From Jha group, these are the cell types from their clusters
  208. celltype_cat <- data.frame(celltype_cat=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  209. rownames(celltype_cat) <- rownames([email hidden]) #add the cell row names there
  210. cci_seurat_pre_PN <- AddMetaData(cci_seurat_pre_PN, celltype_cat) #add this to the meta data
  212. #Add meta data with larger cell type categories, but this time break out the neuroglial into their individual cell types because these are interesting
  213. lookup <- c("0" = "M", "1" = "M", "2" = "Ep", "3" = "M",
  214. "4" = "Ep", "5" = "L", "6" = "En", "7" = "Nt",
  215. "8" = "M", "9" = "M", "10" = "Ep", "11" = "M",
  216. "12" = "En", "13" = "As", "14" = "L", "15" = "Nr",
  217. "16" = "Nt", "17" = "Ep", "18" = "M", "19" = "Olg",
  218. "20" = "L", "21" = "M", "22" = "Ep") #From Jha group, these are the cell types from their clusters
  220. celltype_altcat <- data.frame(celltype_altcat=lookup[as.character([email hidden]$celltype)]) #make a data frame converting the celltype clusters to the celltype names
  221. rownames(celltype_altcat) <- rownames([email hidden]) #add the cell row names there
  222. cci_seurat_pre_PN <- AddMetaData(cci_seurat_pre_PN, celltype_altcat) #add this to the meta data
  224. condition_ctat <- data.frame(condition_ctat=paste0([email hidden]$condition,"_",[email hidden]$celltype_altcat))
  225. rownames(condition_ctat) <- rownames([email hidden])
  226. cci_seurat_pre_PN <- AddMetaData(cci_seurat_pre_PN, condition_ctat) #add this to the meta data
  227. saveRDS(cci_seurat_pre_PN, file="cci_seurat_pre_PN.rds")
  229. #####NORMALIZATION#####
  230. umap_theme <- theme(plot.title=element_blank(),
  231. panel.background=element_blank(),
  232. axis.line=element_line(colour='black',linewidth=0.3),
  233. axis.ticks=element_line(color='black',linewidth=0.3),
  234. axis.text=element_text(size=6,colour='black'),
  235. legend.title=element_blank(),legend.text=element_text(size=6, colour='black'),
  236. legend.position='none',title=element_blank(),
  237. axis.title=element_text(size=6,colour='black'))
  239. #Now we will identify PN genes that are in our data set --> these will be used for the PCA
  241. #log transform
  242. #First we normalize on all genes
  243. #I've decided to use this pipeline instead of the sctransform pipeline mostly because of the desire to use only the proteostasis genes as the loadings for the PCs.
  244. #The problem with using sctransform is the following. sctransform recommends splitting and re-integrated the data. The standard reintegration with CCA is very slow and may introduce some biases.
  245. #Doing RPCA is much faster, but requires doing the PCA first. But what that means is that I would either need to do the PCA on all the genes (defeating my purpose) or on the proteostasis genes alone,
  246. #which would mean that the reintegration would be limited to this smaller subset of genes, which I also think is not the ideal situation.
  247. #for this reason, I am just going to use the old pipeline.
  248. #I will continue to select genes that make some appearance in all the samples, as would be the default requirement in sctransform, because I think that's probably a reasonable thing anyway.
  250. # The other thing I note here is that I start by doing all the analysis only on PN genes. This is because even when you do the PCA with only the PN genes, when you try to find markers you find that the markers
  251. # of the clusters are largely not the PN genes, making it a bit harder to find the PN genes that drive the clustering. I guess I could look in the loadings as well.
  253. #Check the PCA explanatory power using all PN genes on the PN only seurat object.
  255. cci_seurat_pre_PN <- readRDS('cci_seurat_pre_PN.rds')
  256. cci_seurat_pre_PN <- NormalizeData(cci_seurat_pre_PN)
  258. all.genes <- rownames(cci_seurat_pre_PN) #all the genes in the seurat object
  260. cci_seurat_pre_PN <- ScaleData(cci_seurat_pre_PN, features=all.genes)
  261. cci_seurat_pre_PN <- RunPCA(cci_seurat_pre_PN, npcs=100, features=all.genes) #this does PCA using the PN_genes as features
  263. p<-ElbowPlot(cci_seurat_pre_PN, ndims=100, reduction='pca')+umap_theme
  264. p$layers[[1]]$aes_params$size <- 0.5
  266. pdf(paste0('elbow.pdf'), width=3/cm(1), height=3/cm(1))
  267. p
  268. dev.off()
  270. percent_var_pn <- 100*(cci_seurat_pre_PN[["pca"]]@stdev^2) / cci_seurat_pre_PN@reductions[["pca"]]@misc$total.variance
  271. cumulative_var_pn <- cumsum(percent_var_pn)
  273. saveRDS(cci_seurat_pre_PN, 'cci_seurat_pre_PN_PCA.rds')
  274. ##16.2% explained by 100 PCs
  276. ####FIND CLUSTERS####
  277. cci_seurat_lt_pca <- readRDS('cci_seurat_pre_PN_PCA.rds')
  278. #Now I cluster at various PC numbers and dimensions.
  280. for (dim in c(20,30,50)) {
  281. set.seed(42)
  282. cci_seurat <- FindNeighbors(cci_seurat_lt_pca, dims = 1:dim)
  284. for (i in seq(0.1,2,0.1)) {
  285. cci_seurat <- FindClusters(cci_seurat, resolution=i, cluster.name=paste0('res_',as.character(i)))
  286. }
  287. saveRDS(cci_seurat, paste0('cci_seurat_lt_pca_',as.character(dim),'.rds'))
  288. gc()
  289. }
  291. #####EVALUATE THE CLUSTERS#####
  292. if (!dir.exists('mouse_initial_cluster_tests')) {
  293. dir.create('mouse_initial_cluster_tests')
  294. }
  296. for (dim in c(20,30,50)) {
  297. cci_seurat <- readRDS(paste0('cci_seurat_lt_pca_',as.character(dim),'.rds'))
  299. rand_frame_lt_pca <- data.frame(resolution=c(), rand=c())
  301. for (i in seq(0.1,2,0.1)) {
  302. ari <- adjustedRandIndex([email hidden][[paste0('res_',as.character(i))]],[email hidden]$celltype)
  303. temp_frame <- data.frame(resolution=i,rand=ari)
  304. rand_frame_lt_pca <- rbind(rand_frame_lt_pca, temp_frame)
  305. }
  307. ari_plot <- ggplot(rand_frame_lt_pca, aes(x=resolution,y=rand))+geom_point(stat='identity', colour='black',size=1,alpha=1)+scale_x_continuous('Resolution')+scale_y_continuous('Adjusted Rand Index')+umap_theme
  308. ggsave(paste0('mouse_initial_cluster_tests/mouse_ari',as.character(dim),'.pdf'), plot=ari_plot, width=3, height=3, units='cm', dpi=600)
  310. cluster_tree<-clustree(cci_seurat,prefix='res_')
  311. ggsave(paste0('mouse_initial_cluster_tests/mouse_clustree_',as.character(dim),'.pdf'),plot=cluster_tree, width=40, height=30, units='cm', dpi=600)
  313. set.seed(42)
  314. out <- mrtree([email hidden], prefix='res_', n.cores = 12, consensus=FALSE, augment.path=FALSE)
  315. saveRDS(out, paste0('mouse_initial_cluster_tests/out_',as.character(dim),'.rds'))
  317. ks.flat = apply(out$labelmat.flat, 2, FUN=function(x) length(unique(x)))
  318. ks.mrtree = apply(out$labelmat.mrtree, 2, FUN=function(x) length(unique(x)))
  319. amri.flat = sapply(1:ncol(out$labelmat.flat), function(i) AMRI(out$labelmat.flat[,i], [email hidden]$celltype_lab)$amri)
  320. amri.flat = aggregate(amri.flat, by=list(k=ks.flat), FUN=mean)
  321. amri.recon = sapply(1:ncol(out$labelmat.mrtree), function(i) AMRI(out$labelmat.mrtree[,i], [email hidden]$celltype_lab)$amri)
  323. df = rbind(data.frame(k=amri.flat$k, amri=amri.flat$x, method='Seurat flat'),
  324. data.frame(k=ks.mrtree, amri=amri.recon, method='MRtree'))
  326. amri_plot <- ggplot(data=df, aes(x=k, y=amri, color=method)) + geom_line() + theme_bw()
  327. ggsave(paste0('mouse_initial_cluster_tests/mouse_amri_',as.character(dim),'.pdf'), plot=amri_plot, width=8, height=3, units='cm', dpi=600)
  329. stab.out = stability_plot(out)
  330. ggsave(paste0('mouse_initial_cluster_tests/mouse_stab_',as.character(dim),'.pdf'), plot=stab.out$plot, width=6, height=9, units='cm', dpi=600)
  332. cluster_numbers <- apply(out$labelmat.flat, 2, max)+1
  333. print(cluster_numbers)
  334. gc()
  335. }
  337. #20
  338. #res_0.1 res_0.2 res_0.3 res_0.4 res_0.5 res_0.6 res_0.7 res_0.8 res_0.9 res_1 res_1.1 res_1.2 res_1.3 res_1.4 res_1.5 res_1.6 res_1.7 res_1.8 res_1.9 res_2
  339. #10 14 16 17 17 19 19 21 24 25 27 27 29 30 31 31 33 35 36 37
  341. #30
  342. #res_0.1 res_0.2 res_0.3 res_0.4 res_0.5 res_0.6 res_0.7 res_0.8 res_0.9 res_1 res_1.1 res_1.2 res_1.3 res_1.4 res_1.6 res_1.5 res_1.7 res_1.8 res_1.9 res_2
  343. #12 14 16 17 17 19 21 24 25 26 27 29 30 31 33 34 35 36 37 39
  345. #50
  346. # res_0.1 res_0.2 res_0.3 res_0.4 res_0.5 res_0.6 res_0.7 res_0.8 res_0.9 res_1 res_1.2 res_1.1 res_1.3 res_1.4 res_1.6 res_1.5 res_1.7 res_1.8 res_1.9 res_2
  347. # 12 13 15 17 17 18 19 21 23 25 27 28 30 31 32 33 34 34 36 37
  349. #####MAKE UMAPs FOR THE SUBSET OF CLUSTERS THAT PERFORM REASONABLY WELL#####
  350. if (!dir.exists('mouse_initial_UMAPs')) {
  351. dir.create('mouse_initial_UMAPs')
  352. }
  354. iter_frame <- data.frame(dim_options = c("20","30","30","50","50","50","50","50","50","50","50"),
  355. res_options=c("0.9","1","1.1","0.4","0.5","0.6","0.7","0.8","1.1","1.2","1.3"))
  357. for (dim in c(20,30,50)) {
  358. cci_seurat <- readRDS(paste0('cci_seurat_lt_pca_',dim,'.rds'))
  359. cci_seurat <- RunUMAP(cci_seurat, dims=1:dim)
  360. dimplot_by_celltype <- DimPlot(cci_seurat,reduction='umap', group.by="celltype_lab", label=TRUE, label.size=2,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cell type
  361. ggsave(paste0('mouse_initial_UMAPs/mouse_UMAP_',as.character(dim),'celltype_lab.pdf'), plot=dimplot_by_celltype, width=7, height=7, units='cm', dpi=600)
  363. saveRDS(cci_seurat, paste0('cci_seurat_lt_pca_',as.character(dim),'_UMAP.rds'))
  365. for (i in 1:nrow(iter_frame)) {
  366. dim_option <- iter_frame[i,]$dim_options
  367. res_option <- iter_frame[i,]$res_options
  369. if (dim==dim_option) {
  370. dimplot_by_res <- DimPlot(cci_seurat,reduction='umap', group.by=paste0('res_',res_option), label=TRUE, label.size=2,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cell type
  371. ggsave(paste0('mouse_initial_UMAPs/mouse_UMAP_',as.character(dim_option),'_',as.character(res_option),'.pdf'), plot=dimplot_by_res, width=7, height=7, units='cm', dpi=600)
  372. }
  373. }
  374. }
  376. #####IDENTIFY CLUSTERS FOR MERGER BASED ON MARKER JACCARD SIMILARITY#####
  377. ###This function generates the markers to be used in the jaccard function
  378. if (!dir.exists('mouse_clustering_evaluation')) {
  379. dir.create('mouse_clustering_evaluation')
  380. }
  382. markers_function <- function(pca_dims, res) {
  383. prefix <- paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',pca_dims,'_res_',res)
  384. seurat_object <- readRDS(paste0('mouse_initial_UMAPs/cci_seurat_lt_pca_',pca_dims,'_UMAP.rds'))
  385. Idents(object = seurat_object) <- paste0('res_',res)
  386. the_markers_wilcox <- FindAllMarkers(seurat_object, logfc.threshold=-Inf, min.pct=-Inf, genes.use=row.names(seurat_object)) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  387. write.table(the_markers_wilcox, file=paste0(prefix,'_wilcox_markers.txt'), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  388. the_markers_wilcox_subset <- subset(the_markers_wilcox, the_markers_wilcox$p_val_adj < 0.05 & the_markers_wilcox$avg_log2FC >= 0.585 & the_markers_wilcox$pct.1 >= 0.5) #we then subset for those that are significant, with at least 1.5 fold change, and that are present in at least 50% of cells in the cluster.
  389. write.table(the_markers_wilcox_subset, file=paste0(prefix,'_wilcox_markers_significant.txt'), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  390. return (list(all_markers = the_markers_wilcox, subset_markers = the_markers_wilcox_subset))
  391. }
  393. ###The jaccard function looks for similar clusters and returns the jaccard values and a reclustering in graph form
  394. jaccard_function <- function(wilcox_markers, pca_dims, res, percentage_threshold = 33.3) {
  395. prefix <- paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',pca_dims,'_res_',res)
  397. the_clusters <- unique(wilcox_markers$cluster)
  398. jaccard_matrix <- matrix(nrow=length(the_clusters), ncol=length(the_clusters))
  399. colnames(jaccard_matrix) <- the_clusters
  400. rownames(jaccard_matrix) <- the_clusters
  402. jaccard_threshold <- matrix(nrow=length(the_clusters), ncol=length(the_clusters))
  403. colnames(jaccard_threshold) <- the_clusters
  404. rownames(jaccard_threshold) <- the_clusters
  406. for (cluster1 in the_clusters) {
  407. for (cluster2 in the_clusters) {
  408. if (cluster1 != cluster2 & is.na(jaccard_matrix[cluster2,cluster1])) {
  409. subset1 <- unique(subset(wilcox_markers, wilcox_markers$cluster == cluster1)$gene)
  410. subset2 <- unique(subset(wilcox_markers, wilcox_markers$cluster == cluster2)$gene)
  411. the_intersection <- intersect(subset1, subset2)
  412. the_union <- union(subset1, subset2)
  413. jaccard <- 100*length(the_intersection)/length(the_union)
  414. jaccard_matrix[cluster1,cluster2] = jaccard
  415. if (jaccard >= 33.3) {
  416. jaccard_threshold[cluster1, cluster2] = 1
  417. }
  418. else {
  419. jaccard_threshold[cluster1, cluster2] = 0
  420. }
  421. }
  422. }
  423. }
  425. jaccard_adjacency <- matrix(nrow=length(the_clusters), ncol=length(the_clusters))
  426. colnames(jaccard_adjacency) <- the_clusters
  427. rownames(jaccard_adjacency) <- the_clusters
  429. for (cluster1 in the_clusters) {
  430. for (cluster2 in the_clusters) {
  431. subset1 <- unique(subset(wilcox_markers, wilcox_markers$cluster == cluster1)$gene)
  432. subset2 <- unique(subset(wilcox_markers, wilcox_markers$cluster == cluster2)$gene)
  433. the_intersection <- intersect(subset1, subset2)
  434. the_union <- union(subset1, subset2)
  435. jaccard <- 100*length(the_intersection)/length(the_union)
  436. if (cluster1 == cluster2) {
  437. jaccard_adjacency[cluster1, cluster2] = 0
  438. }
  439. else if (jaccard >= percentage_threshold) {
  440. jaccard_adjacency[cluster1, cluster2] = 1
  441. }
  442. else {
  443. jaccard_adjacency[cluster1, cluster2] = 0
  444. }
  445. }
  446. }
  447. set.seed(123)
  448. tg <- graph_from_adjacency_matrix(jaccard_adjacency, mode='undirected')
  449. saveRDS(tg, paste0(prefix,'_test_graph.rds'))
  450. saveRDS(jaccard_matrix, paste0(prefix,'_jaccard_matrix.rds'))
  451. write.table(jaccard_matrix, file=paste0(prefix,'_jaccard_wilcox_matrix.txt'), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  452. write.table(jaccard_threshold, file=paste0(prefix,'_jaccard_wilcox_threshold.txt'), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  454. pdf(paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',pca_dims,'_res_',res,'_test_graph.pdf'), width=4.5, height=4.5)
  455. plot.igraph(tg,vertex.color = "orange",vertex.size=15,vertex.frame.color = "black",vertex.label.cex=0.75, layout=layout_with_fr(tg, niter=10000))
  456. dev.off()
  458. return(list(mat = jaccard_matrix, adj = jaccard_adjacency, test_graph = tg))
  459. }
  461. iter_frame <- data.frame(dim_options = c("20","30","30","50","50","50","50","50","50","50","50"),
  462. res_options=c("0.9","1","1.1","0.4","0.5","0.6","0.7","0.8","1.1","1.2","1.3"))
  464. for (i in 1:nrow(iter_frame)) {
  465. dim_option <- iter_frame[i,]$dim_options
  466. res_option <- iter_frame[i,]$res_options
  468. the_markers <- markers_function(dim_option,res_option)
  469. jaccard <- jaccard_function(the_markers$subset_markers,dim_option,res_option)
  470. }
  472. ####MERGE CLUSTERS BASED ON ABOVE#####
  473. if (!dir.exists('mouse_merged_cluster_compare_fraction')) {
  474. dir.create('mouse_merged_cluster_compare_fraction')
  475. }
  477. if (!dir.exists('mouse_merged_clusters')) {
  478. dir.create('mouse_merged_clusters')
  479. }
  481. iter_frame <- data.frame(dim_options = c("20","30","30","50","50","50","50","50","50","50","50"),
  482. res_options=c("0.9","1","1.1","0.4","0.5","0.6","0.7","0.8","1.1","1.2","1.3"))
  484. for (i in 1:nrow(iter_frame)) {
  485. dim_option <- iter_frame[i,]$dim_options
  486. res_option <- iter_frame[i,]$res_options
  488. lookup_1 <- clusters(readRDS(paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',dim_option,'_res_',res_option,'_test_graph.rds')))$membership
  490. #now if we are on a new dimension seurat object, we open that file
  492. if(i==1) {
  493. cci_seurat <- readRDS(paste0('mouse_initial_UMAPs/cci_seurat_lt_pca_',dim_option,'_UMAP.rds'))
  494. }
  495. else if(!(iter_frame[i,]$dim_options == iter_frame[i-1,]$dim_options)) {
  496. cci_seurat <- readRDS(paste0('mouse_initial_UMAPs/cci_seurat_lt_pca_',dim_option,'_UMAP.rds'))
  497. }
  499. merged_clusters <- data.frame(x=lookup_1[as.character([email hidden][[paste0('res_',res_option)]])]) #make a data frame converting the celltype clusters to the celltype names
  500. rownames(merged_clusters) <- rownames([email hidden])#add the cell row names there
  502. column_name <- paste0("merged_clusters_",res_option)
  503. colnames(merged_clusters) <- column_name #rename column
  504. cci_seurat <- AddMetaData(cci_seurat, merged_clusters) #add this to the meta data
  506. dimplot_by_column_name <- DimPlot(cci_seurat,reduction='umap', group.by=column_name, label=TRUE, label.size=3,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cluster
  507. ggsave(paste0('mouse_merged_clusters/mouse_merged_clusters_',dim_option,'_',res_option,'.pdf'), plot=dimplot_by_column_name, width=7, height=7, units='cm', dpi=600)
  509. merged_cluster_compare <- table(data.frame(clusters=[email hidden][[column_name]], celltypes=[email hidden]$celltype_lab))
  510. merged_cluster_compare_fraction <- sweep(merged_cluster_compare, 2, colSums(merged_cluster_compare), FUN='/')
  511. write.table(merged_cluster_compare_fraction, paste0('mouse_merged_cluster_compare_fraction/cci_seurat_lt_pca_',dim_option,'_merged_cluster_compare_fraction_',res_option,'.txt'),quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  513. #if we have finished that dimension seurat object, we can finally save that file
  514. if(i==nrow(iter_frame)) {
  515. saveRDS(cci_seurat, paste0('mouse_merged_clusters/cci_seurat_lt_pca_',dim_option,'_mergedclusters.rds'))
  516. }
  517. else if(!(iter_frame[i,]$dim_options == iter_frame[i+1,]$dim_options)) {
  518. saveRDS(cci_seurat, paste0('mouse_merged_clusters/cci_seurat_lt_pca_',dim_option,'_mergedclusters.rds'))
  519. }
  520. }
  522. #####SECOND MERGE: BASED ON MERGING IN LOW FIDELITY CLUSTERS#####
  523. ##One thing we notice after the new clusters are made is that some of them are low fidelity.
  524. #That is, their maximum value in any given cell type cluster is quite low.
  525. #We arbitrarily set the cutoff that if the maximum a given cluster contributes to the assignment of a cell type is < 20%, it should be merged with its closest cluster based on jaccard.
  526. #This tells us which clusters will need to be recast in this way.
  527. if (!dir.exists('mouse_merged_cluster_compare_fraction/cluster_assignment_scores')) {
  528. dir.create('mouse_merged_cluster_compare_fraction/cluster_assignment_scores')
  529. }
  531. umap_theme <- theme(plot.title=element_blank(),
  532. panel.background=element_blank(),
  533. axis.line=element_line(colour='black',linewidth=0.3),
  534. axis.ticks=element_line(color='black',linewidth=0.3),
  535. axis.text=element_text(size=6,colour='black'),
  536. legend.title=element_blank(),legend.text=element_text(size=6, colour='black'),
  537. legend.position='none',title=element_blank(),
  538. axis.title=element_text(size=6,colour='black'))
  540. merged_cluster_compare_files <- list.files('mouse_merged_cluster_compare_fraction')
  542. for (the_file in merged_cluster_compare_files) {
  543. if (!the_file=='cluster_assignment_scores') {
  544. temp_table <- read.table(paste0('mouse_merged_cluster_compare_fraction/',the_file), header=TRUE)
  545. temp_table_t <- as.data.frame(t(temp_table))
  546. the_scores <- data.frame(cluster=c(), scores=c())
  547. for (the_col in colnames(temp_table_t)) {
  548. temp_frame <- data.frame(cluster=c(the_col), scores=c(max(temp_table_t[[the_col]])))
  549. the_scores <- rbind(the_scores, temp_frame)
  550. }
  551. write.table(the_scores, paste0('mouse_merged_cluster_compare_fraction/cluster_assignment_scores/cluster_assignment_scores_',the_file), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  552. }
  553. }
  555. #The algorithm for evaluating this first we identify what original cluster the merged cluster maps from.
  556. #Then we find near neighbors of that cluster based on jaccard, and consider as potential candidates those clusters that have a jaccard >= 0.2.
  557. #Then we find the average jaccard for the clusters to which that cluster belongs, and we select the one for which this cluster has the highest average jaccard.
  559. cluster_assignment_score_files <- list.files('mouse_merged_cluster_compare_fraction/cluster_assignment_scores')
  561. reassignment_function <- function(the_dim, the_res, lookup) {
  562. cluster_scores <- as.data.frame(read.table(paste0("mouse_merged_cluster_compare_fraction/cluster_assignment_scores/cluster_assignment_scores_cci_seurat_lt_pca_",the_dim,'_merged_cluster_compare_fraction_',the_res,'.txt'), header=TRUE))
  563. poor_clusters <- subset(cluster_scores, scores < 0.2) #find the clusters with poor scores
  564. jaccard_matrix <- read.table(paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',the_dim,'_res_',the_res,'_jaccard_wilcox_matrix.txt'),header=TRUE) #we will use the original jaccard matrix below
  565. jaccard_matrix_t <- as.data.frame(t(jaccard_matrix)) #the transposed matrix is also useful
  567. message("These are the poor clusters:")
  568. print(poor_clusters)
  570. cluster_assignments <- c()
  571. for (the_cluster in poor_clusters$cluster) { #for every poor scoring cluster
  572. original_cluster <- names(which(lookup == as.character(the_cluster))) #first find out what the original cluster is that it is associated with. This is not designed to work on clusters that are associated with more than one original cluster
  573. message(paste0("For the cluster ",the_cluster," the original cluster is ",original_cluster))
  575. if (length(original_cluster) > 1) {
  576. stop("This merged cluster arises from more than one original cluster. Cannot be handled by current algorithm.")
  577. }
  579. #now we reassemble the set of jaccard scores for the original cluster
  580. beginning <- na.omit(jaccard_matrix[[paste0('X',original_cluster)]])
  581. end <- na.omit(jaccard_matrix_t[[original_cluster]])
  582. jaccard_scores <- c(beginning,NA,end)
  583. print(jaccard_scores)
  585. #we only want to look at the jaccard scores >= 20 as being good potential merge candidates
  586. jaccard_scores_select <- jaccard_scores[jaccard_scores >=20 & !(jaccard_scores %in% NA)]
  587. select_clusters <- as.character(which(jaccard_scores %in% jaccard_scores_select)-1) #Find the indices associated with the high scoring jaccards and then subtract one to get the corresponding cluster
  588. if (length(select_clusters) > 0) {
  589. message(paste0("These original clusters have high associated scores with original cluster ",original_cluster))
  590. print(select_clusters)
  591. merged_cluster_candidates <- c()
  592. for (select_cluster in select_clusters) {
  593. merged_cluster_candidates <- c(merged_cluster_candidates, lookup[[select_cluster]])
  594. }
  595. merged_cluster_candidates <- unique(merged_cluster_candidates)
  596. message("These merged clusters are associated with those original clusters, and may be candidates to merge into:")
  597. print(merged_cluster_candidates)
  598. candidate_scores = c()
  599. for (merged_candidate in merged_cluster_candidates) {
  600. candidate_clusters <- as.numeric(names(which(lookup==merged_candidate)))
  601. candidate_scores[[as.character(merged_candidate)]] <- mean(jaccard_scores[candidate_clusters+1])
  602. }
  603. candidate_scores <- as.data.frame(candidate_scores)
  604. winning_merge_candidate <- merged_cluster_candidates[which(candidate_scores == max(candidate_scores))]
  605. message(paste0("Merge original cluster ", original_cluster," into merged cluster ", winning_merge_candidate))
  606. cluster_assignments[[original_cluster]] <- winning_merge_candidate
  607. }
  608. else {message('There are not other original clusters highly associated with ',original_cluster)}
  609. }
  611. lookup_2 <- lookup
  612. for (original_cluster in names(cluster_assignments)) {
  613. lookup_2[[original_cluster]] <- cluster_assignments[[original_cluster]]
  614. }
  616. print(lookup_2)
  617. lookup_lookup <- seq(1, length(unique(lookup_2)), 1)
  618. names(lookup_lookup) <- sort(unique(lookup_2))
  620. lookup_3 <- lookup_lookup[as.character(lookup_2)]
  621. names(lookup_3) <- names(lookup_2)
  622. return(lookup_3)
  623. }
  625. ##Now I re-do the assignments
  626. if (!dir.exists('mouse_recluster')) {
  627. dir.create('mouse_recluster')
  628. }
  630. iter_frame <- data.frame(dim_options = c("20","30","30","50","50","50","50","50","50","50","50"),
  631. res_options=c("0.9","1","1.1","0.4","0.5","0.6","0.7","0.8","1.1","1.2","1.3"))
  633. iter_frame <- data.frame(dim_options = c("50"),
  634. res_options=c("0.5"))
  636. for (i in 1:nrow(iter_frame)) {
  637. dim_option <- iter_frame[i,]$dim_options
  638. res_option <- iter_frame[i,]$res_options
  640. if(i==1) {
  641. cci_seurat <- readRDS(paste0('mouse_merged_clusters/cci_seurat_lt_pca_',dim_option,'_mergedclusters.rds'))
  642. }
  643. else if(!(iter_frame[i,]$dim_options == iter_frame[i-1,]$dim_options)) {
  644. cci_seurat <- readRDS(paste0('mouse_merged_clusters/cci_seurat_lt_pca_',dim_option,'_mergedclusters.rds'))
  645. }
  647. lookup_1 <- clusters(readRDS(paste0('mouse_clustering_evaluation/cci_seurat_lt_pca_',dim_option,'_res_',res_option,'_test_graph.rds')))$membership
  648. new_lookup <- reassignment_function(dim_option, res_option, lookup_1)
  650. merged_clusters <- data.frame(x=new_lookup[as.character([email hidden][[paste0('res_',res_option)]])])#make a data frame converting the celltype clusters to the celltype names
  651. rownames(merged_clusters) <- rownames([email hidden])#add the cell row names there
  653. column_name <- paste0("recluster_merged_clusters_",res_option)
  654. colnames(merged_clusters) <- column_name #rename column
  655. cci_seurat <- AddMetaData(cci_seurat, merged_clusters) #add this to the meta data
  657. dimplot_by_column_name <- DimPlot(cci_seurat,reduction='umap', group.by=column_name, label=TRUE, label.size=3,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cluster
  658. ggsave(paste0('mouse_recluster/test_mouse_recluster_merged_clusters_',dim_option,'_',res_option,'_UMAP.pdf'), plot=dimplot_by_column_name, width=7, height=7, units='cm', dpi=600)
  660. merged_cluster_compare <- table(data.frame(clusters=[email hidden][[column_name]], celltypes=[email hidden]$celltype_lab))
  661. merged_cluster_compare_fraction <- sweep(merged_cluster_compare, 2, colSums(merged_cluster_compare), FUN='/')
  662. write.table(merged_cluster_compare_fraction, paste0('mouse_recluster/test_recluster_cci_seurat_lt_pca_',dim_option,'_merged_cluster_compare_fraction_',res_option,'.txt'),quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  664. if(i==nrow(iter_frame)) {
  665. saveRDS(cci_seurat, paste0('mouse_recluster/test_cci_seurat_lt_pca_',dim_option,'_recluster.rds'))
  666. }
  667. else if(!(iter_frame[i,]$dim_options == iter_frame[i+1,]$dim_options)) {
  668. saveRDS(cci_seurat, paste0('mouse_recluster/test_cci_seurat_lt_pca_',dim_option,'_recluster.rds'))
  669. }
  670. }
  672. #####FIND MARKERS ON THE NEW CLUSTERS#####
  673. cci_seurat_lt_pca_50_merged <- readRDS('mouse_recluster/cci_seurat_lt_pca_50_recluster.rds')
  674. Idents(cci_seurat_lt_pca_50_merged) <- 'recluster_merged_clusters_0.5'
  675. Idents(cci_seurat_lt_pca_50_merged) <- factor(Idents(cci_seurat_lt_pca_50_merged), levels=order(unique(Idents(cci_seurat_lt_pca_50_merged))))
  677. the_markers_wilcox <- FindAllMarkers(cci_seurat_lt_pca_50_merged, logfc.threshold=-Inf, min.pct=-Inf, genes.use=row.names(cci_seurat_lt_pca_50_merged)) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  678. write.table(the_markers_wilcox, file='cci_seurat_lt_pca_50_merged_wilcox_markers.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  679. the_markers_wilcox_subset <- subset(the_markers_wilcox, the_markers_wilcox$p_val_adj < 0.05 & the_markers_wilcox$avg_log2FC >= 0.585 & the_markers_wilcox$pct.1 >= 0.5) #we then subset for those that are significant, with at least 1.5 fold change, and that are present in at least 50% of cells in the cluster.
  680. write.table(the_markers_wilcox_subset, file='cci_seurat_lt_pca_50_merged_wilcox_markers_significant.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  682. test_markers_wilcox_subset <- subset(the_markers_wilcox, the_markers_wilcox$p_val_adj < 0.05 & the_markers_wilcox$avg_log2FC >= 0.585 & the_markers_wilcox$pct.1 >= 0.5 & the_markers_wilcox$pct.1/the_markers_wilcox$pct.2 > 1) #we then subset for those that are significant, with at least 1.5 fold change, and that are present in at least 50% of cells in the cluster.
  683. table(test_markers_wilcox_subset$cluster)
  684. write.table(test_markers_wilcox_subset, file='test_cci_seurat_lt_pca_50_merged_wilcox_markers_significant.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  686. test_markers_wilcox_subset <- cbind(test_markers_wilcox_subset, data.frame(ratio=test_markers_wilcox_subset$pct.1/test_markers_wilcox_subset$pct.2))
  688. ##now we select the test_features to put in the dot plot. This is based on the top features for each cluster, ranked by the ratio of pct.1 to everything else (pct.2)
  689. the_clusters <- order(unique(test_markers_wilcox_subset$cluster))
  691. test_features <- c()
  692. for (the_cluster in the_clusters) {
  693. temp_subset <- subset(test_markers_wilcox_subset, cluster %in% the_cluster)
  694. temp_subset <- temp_subset[order(temp_subset$avg_log2FC, decreasing=TRUE),]
  695. test_features <- c(test_features,temp_subset[1:2,]$gene)
  696. }
  698. dotplot_legend_theme <- theme(plot.title=element_blank(),
  699. panel.background=element_blank(),
  700. axis.line=element_line(colour='black',linewidth=0.3),
  701. axis.ticks=element_line(color='black',linewidth=0.3),
  702. axis.text=element_text(size=6,colour='black'),
  703. title=element_blank(),
  704. axis.title=element_text(size=6,colour='black'), legend.text = element_text(size=6, colour='black'))
  706. dotplot <- DotPlot(cci_seurat_lt_pca_50_merged, features = test_features, dot.scale=3) + RotatedAxis()+dotplot_legend_theme
  707. ggsave('mouse_dotplot.pdf', plot=dotplot, width=9, height=7, units='cm', dpi=600)
  709. #Cluster Fraction Per Cell Type
  710. recluster_compare_0.5 <- table(data.frame(clusters=[email hidden]$recluster_merged_clusters_0.5, celltypes=[email hidden]$celltype_lab))
  711. recluster_compare_fraction_0.5 <- sweep(recluster_compare_0.5, 2, colSums(recluster_compare_0.5), FUN='/')
  713. pdf('mouse_cluster_fraction_per_celltype.pdf', width=3.9839, height=2.0617)
  714. Heatmap(recluster_compare_fraction_0.5,
  715. circlize::colorRamp2(seq(0,1,.01), colorRampPalette(rev(brewer.pal(n=7, name='RdYlBu')))(101)), na_col='#EEEEEE',
  716. column_names_gp=gpar(fontsize=6), row_names_gp = gpar(fontsize = 6),
  717. cluster_rows=TRUE, cluster_columns=TRUE,
  718. heatmap_legend_param = list(
  719. title = "Cluster Fraction Per Cell Type", at = c(0,0.5,1),
  720. title_position='leftcenter-rot',
  721. title_gp=gpar(fontsize=6, fontface='plain'),
  722. labels_gp=gpar(fontsize=6)))
  723. dev.off()
  725. #Cluster Fraction Per Condition
  726. recluster_compare_0.5_condition <- table(data.frame(clusters=[email hidden]$recluster_merged_clusters_0.5, celltypes=[email hidden]$orig.ident))
  727. recluster_compare_fraction_0.5_condition <- sweep(recluster_compare_0.5_condition, 2, colSums(recluster_compare_0.5_condition), FUN='/')
  729. pdf('mouse_cluster_fraction_per_condition.pdf', width=2.094, height=2.13)
  730. Heatmap(recluster_compare_fraction_0.5_condition,
  731. circlize::colorRamp2(seq(0,1,.01), colorRampPalette(rev(brewer.pal(n=7, name='RdYlBu')))(101)), na_col='#EEEEEE',
  732. column_names_gp=gpar(fontsize=6), row_names_gp = gpar(fontsize = 6),
  733. cluster_rows=TRUE, cluster_columns=TRUE,
  734. heatmap_legend_param = list(
  735. title = "Cluster Fraction Per Condition", at = c(0,0.5,1),
  736. title_position='leftcenter-rot',
  737. title_gp=gpar(fontsize=6, fontface='plain'),
  738. labels_gp=gpar(fontsize=6)))
  739. dev.off()
  741. umap_theme <- theme(plot.title=element_blank(),
  742. panel.background=element_blank(),
  743. axis.line=element_line(colour='black',linewidth=0.3),
  744. axis.ticks=element_line(color='black',linewidth=0.3),
  745. axis.text=element_text(size=6,colour='black'),
  746. legend.title=element_blank(),legend.text=element_text(size=6, colour='black'),
  747. legend.position='none',title=element_blank(),
  748. axis.title=element_text(size=6,colour='black'))
  750. dimplot_by_recluster_merged_clusters <- DimPlot(cci_seurat_lt_pca_50_merged,reduction='umap', group.by="recluster_merged_clusters_0.5", label=TRUE, label.size=2,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cluster
  751. ggsave('test_mouse_recluster_merged_clusters_0.5_UMAP.pdf', plot=dimplot_by_recluster_merged_clusters, width=7, height=7, units='cm', dpi=600)
  753. dimplot_by_recluster_merged_clusters <- DimPlot(cci_seurat_lt_pca_50_merged,reduction='umap', group.by="celltype_lab", label=TRUE, label.size=2,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cluster
  754. ggsave('mouse_recluster_merged_clusters_0.5_UMAP_by_celltype_lab.pdf', plot=dimplot_by_recluster_merged_clusters, width=7, height=7, units='cm', dpi=600)
  756. dimplot_by_recluster_merged_clusters <- DimPlot(cci_seurat_lt_pca_50_merged,reduction='umap', group.by="condition", label=TRUE, label.size=2,label.box=TRUE, repel=TRUE, pt.size=1,raster=TRUE)+umap_theme+labs(x='UMAP 1',y='UMAP 2') #plot by cluster
  757. ggsave('mouse_recluster_merged_clusters_0.5_UMAP_by_condition.pdf', plot=dimplot_by_recluster_merged_clusters, width=7, height=7, units='cm', dpi=600)
  759. #####PREP OBJECTS FOR GSEA#####
  760. #I now prep seurat objects that can be used for GSEA (they contain the new clustering, but on objects that have all the genes)
  761. #First I apply the new clusters onto a seurat object that hs all the genes
  762. cci_recluster <- readRDS('mouse_recluster/cci_seurat_lt_pca_50_recluster.rds')
  763. cci_pre <- readRDS('cci_seurat_pre.rds')
  765. recluster <- data.frame(recluster=[email hidden]$recluster_merged_clusters_0.5) #make a data frame converting the celltype clusters to the celltype names
  766. rownames(recluster) <- rownames([email hidden]) #add the cell row names there
  767. cci_pre <- AddMetaData(cci_pre, recluster)
  769. #now I normalize that seurat object
  770. cci_pre <- NormalizeData(cci_pre)
  771. all.genes <- rownames(cci_pre) #all the genes in the seurat object
  773. cci_pre <- ScaleData(cci_pre, features=all.genes)
  775. cci_pre <- FindVariableFeatures(cci_pre, selection.method = "vst") #I just want to check how much variance the PCA explains if you use the most variable genes
  776. cci_pre <- RunPCA(cci_pre, npcs=100) #this does PCA on variable genes only
  778. saveRDS(cci_pre, 'cci_pre_scaled.rds')
  780. #####FGSEA WITHIN CELL TYPES/CLUSTERS USING FINDALLMARKERS#####
  781. if (!dir.exists('gsea_results')) {
  782. dir.create('gsea_results')
  783. }
  785. fgsea_pathways <- gmtPathways("PN_limited_branchclass.gmt") #this is the file that will be used as gsea pathways
  787. #First use FindAllMarkers so I can have a log2FC for all the genes
  788. the_idents <- c('celltype_col', 'recluster') #I will do GSEA on groups of cell subtypes and each of the clusters
  789. for (the_ident in the_idents) {
  790. Idents(object = cci_pre) <- the_ident
  791. the_markers_wilcox <- FindAllMarkers(cci_pre, logfc.threshold=-Inf, min.pct=-Inf, genes.use=row.names(cci_pre))
  792. write.table(the_markers_wilcox, file=paste0('cci_pre_',the_ident,'_merged_wilcox_markers.txt'), quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  793. }
  795. #Now I loop through the files and run them through fgsea
  796. for (the_ident in the_idents) {
  797. groups <- unique([email hidden][[the_ident]])
  798. the_markers <- as.data.frame(read.table(paste0('cci_pre_',the_ident,'_merged_wilcox_markers.txt'), sep='\t',header=TRUE,row.names=1))
  799. for (group in groups) {
  800. key <- paste0(the_ident,"_",group)
  801. print(group)
  802. markers_subset <- subset(the_markers, cluster == group)
  803. rank_stat_matrix <- matrix(markers_subset$avg_log2FC)
  804. rownames(rank_stat_matrix) <- markers_subset$gene
  805. fgsea_stat_ranks <- setNames(markers_subset$avg_log2FC, markers_subset$gene)
  806. fgseaRes_stat <- fgsea(fgsea_pathways, fgsea_stat_ranks, minSize=10, maxSize=10000, nPermSimple=50000, scoreType="pos") #I use scoretype=pos here because I am really trying to find pathways that define these clusters, so the ones that are enriched, as opposed to ones that are downregulated.
  807. fwrite(fgseaRes_stat, file=paste0('branchclass_',key,'_stat_fgsea.tsv'), sep='\t', sep2=c(""," ",""))
  808. write.table(rank_stat_matrix,file=paste0('branchclass_rank_stat_matrix_',key,'.rnk'),sep='\t',col.names=FALSE,quote=FALSE)
  809. }
  810. }
  812. #Now this is the plotting
  813. the_idents <- c('recluster','celltype_col')
  814. size_frame <- data.frame(recluster=c(4.383,5.175),
  815. celltype_col=c(5.565,6.14))
  817. for (the_ident in the_idents) {
  818. groups <- unique([email hidden][[the_ident]])
  820. the_pathways <- c()
  821. exclude_list <- c()
  823. for (group in groups) {
  824. key <- paste0(the_ident,"_",group)
  825. gsea_table <- read.table(paste0('gsea_results/branchclass_',key,'_stat_fgsea.tsv'), sep='\t', header=TRUE)
  826. gsea_table <- subset(gsea_table, padj < 0.05) #I take the subset of significant pathways for each of the files
  827. if (length(rownames(gsea_table)) == 0) {
  828. exclude_list <- c(exclude_list, key) #if a file has no significant pathways, it will not be plotted later
  829. }
  830. the_pathways <- c(the_pathways, gsea_table$pathway)
  831. }
  832. the_pathways <- unique(the_pathways)
  833. NES_subset <- list()
  835. for (group in groups) {
  836. key <- paste0(the_ident,"_",group)
  837. if (!(key %in% exclude_list)) {
  838. gsea_table <- read.table(paste0('gsea_results/branchclass_',key,'_stat_fgsea.tsv'), sep='\t', row.names=1, header=TRUE)
  839. gsea_table <- gsea_table[the_pathways,]
  841. p_filter <- gsea_table$padj < 0.1
  843. NES_temp <- gsea_table$NES
  844. NES_temp[!p_filter] <- 0 #I am basically going to plot the NES as 0 if the padj ≥ 0.1
  845. NES_subset[[key]] <- NES_temp
  846. }
  847. }
  849. NES_subset <- as.data.frame(NES_subset)
  850. row.names(NES_subset) <- the_pathways
  851. NES_subset <- NES_subset[order(row.names(NES_subset)), ]
  852. write.table(NES_subset, paste0('mouse_',the_ident,'_NES_subset.txt'), sep='\t', quote=FALSE)
  854. pdf(paste0('mouse_',the_ident,'_GSEA_heatmap.pdf'), width=size_frame[[the_ident]][1], height=size_frame[[the_ident]][2])
  855. draw(pheatmap(as.matrix(NES_subset), fontsize=6, treeheight_row=10, treeheight_col=10, border_color=NA))
  856. dev.off()
  857. }
  859. #Here I plot the analysis just for microglia
  860. groups <- c('Microglia_1', 'Microglia_2', 'Microglia_3', 'Microglia_4', 'Microglia_5', 'Microglia_6', 'Microglia_7', 'Microglia_8')
  862. the_pathways <- c()
  863. the_ident <- 'celltype_col'
  864. for (group in groups) {
  865. key <- paste0(the_ident,"_",group)
  866. gsea_table <- read.table(paste0('gsea_results/branchclass_',key,'_stat_fgsea.tsv'), sep='\t', header=TRUE)
  867. gsea_table <- subset(gsea_table, padj < 0.05)
  868. # gsea_table <- subset(gsea_table, !(grepl('Translation_*', pathway)))
  869. the_pathways <- c(the_pathways, gsea_table$pathway)
  870. }
  871. the_pathways <- unique(the_pathways)
  872. NES_subset <- list()
  874. for (group in groups) {
  875. key <- paste0(the_ident,"_",group)
  876. gsea_table <- read.table(paste0('gsea_results/branchclass_',key,'_stat_fgsea.tsv'), sep='\t', row.names=1, header=TRUE)
  877. gsea_table <- gsea_table[the_pathways,]
  879. p_filter <- gsea_table$padj < 0.1
  881. NES_temp <- gsea_table$NES
  882. NES_temp[!p_filter] <- 0
  883. NES_subset[[key]] <- NES_temp
  884. }
  886. NES_subset <- as.data.frame(NES_subset)
  887. row.names(NES_subset) <- the_pathways
  888. NES_subset <- NES_subset[order(row.names(NES_subset)), ]
  889. write.table(NES_subset, 'NES_subset.txt', sep='\t', quote=FALSE)
  891. pdf(paste0('microglia_GSEA_heatmap.pdf'), width=4.501, height=4.96)
  892. draw(pheatmap(as.matrix(NES_subset), fontsize=6, treeheight_row=10, treeheight_col=10, border_color=NA))
  893. dev.off()
  895. #####FGSEA BETWEEN NAIVE AND CCI USING PSEUDOBULK#####
  896. if (!dir.exists('pseudobulk')) {
  897. dir.create('pseudobulk')
  898. }
  900. cci_seurat <- readRDS('cci_pre_scaled.rds')
  902. pseudobulk_process <- function(seurat_object, group_by, condition_start_string, condition_end_string) {
  903. pseudo_seurat <- AggregateExpression(seurat_object, assays="RNA", return.seurat=T, group.by = c('orig.ident', group_by))
  904. pseudo_seurat$condition <- substr(pseudo_seurat$orig.ident, condition_start_string, condition_end_string)
  905. return(pseudo_seurat)
  906. }
  908. results_output <- function(pseudobulk_object, group_by, pathways_file, p_threshold, lfc_threshold, contrast1, contrast2, select_split, split_index) {
  909. cell_type_list <- unique([email hidden][[group_by]])
  910. counts_object <- pseudobulk_object[['RNA']]$counts #we are using the raw counts matrix after pseudobulk for our analysis
  912. print(colnames(counts_object))
  913. resAsh_list <- list()
  915. for (cell_type in cell_type_list) {#loop on cell type
  917. #set up the DESeq2 sample table
  918. select_list = grep(glob2rx(paste0('*',cell_type)), colnames(counts_object), value=TRUE)
  919. print(select_list)
  920. countData <- counts_object[, select_list]
  921. cond_list = c()
  922. for (select in select_list) {
  923. cond_list = c(cond_list, strsplit(select, select_split)[[1]][split_index])
  924. }
  925. sample_table <- data.frame(Condition=factor(cond_list))
  926. row.names(sample_table) <- colnames(countData)
  927. print(sample_table)
  929. #DESeq2
  930. dds <- DESeqDataSetFromMatrix(countData=countData, colData=sample_table, design=~Condition) #DeSeq design matrix, we only use condition as a variable, and then we do the analysis on each gene separately
  932. ##filter out some genes
  933. keep <- rowSums(counts(dds) >=10) >= 3 #selected based on it being the smallest group size in any given sample
  934. dds <- dds[keep]
  935. dds1 <- DESeq(dds)
  937. ##generate results table
  938. resAsh <- lfcShrink(dds1, contrast=c('Condition',contrast1,contrast2), type='ashr', alpha=p_threshold, lfcThreshold=lfc_threshold) #perform lfcShrink, I prefer ashr because of the format of the comparisons
  939. print(length(rownames(resAsh)))
  941. resAsh_list[[cell_type]] <- resAsh #assign this to a list where we store all the differential expression data
  943. #GSEA
  944. ##based on various suggestions, including this one from Anthony Castanza at GSEA (https://groups.google.com/g/gsea-help/c/UUNzkRPbvz4), we use stat as the metric.
  945. fgsea_pathways <- gmtPathways(pathways_file) #this is the file that will be used as gsea pathways
  947. rank_stat_matrix <- matrix(resAsh$stat)
  948. rownames(rank_stat_matrix) <- row.names(resAsh)
  949. fgsea_stat_ranks <- setNames(resAsh$stat, row.names(resAsh))
  950. fgseaRes_stat <- fgsea(fgsea_pathways, fgsea_stat_ranks, minSize=10, maxSize=10000, nPermSimple=50000) #number of permutations based on recommendations from the package GOAT
  951. fwrite(fgseaRes_stat, file=paste0('pseudobulk/',cell_type,'_stat_fgsea.tsv'),sep='\t',sep2=c(""," ",""))
  953. #output files
  954. write.table(rank_stat_matrix,file=paste0('pseudobulk/branchclass_rank_stat_matrix_',cell_type,'.rnk'),sep='\t',col.names=FALSE,quote=FALSE)
  955. write.table(resAsh, file=paste0('pseudobulk/branchclass_',cell_type,'_resAsh.tsv'),sep='\t',row.names=TRUE, col.names=NA, quote=FALSE)
  956. }
  957. return(resAsh_list)
  958. }
  960. the_groups <- c('recluster','celltype_col')
  962. for (group_by_cat in the_groups) {
  963. pseudo_seurat <- pseudobulk_process(cci_seurat, group_by_cat, 1, 3)
  964. resAsh_list <- results_output(pseudobulk_object=pseudo_seurat, group_by=group_by_cat, pathways_file='PN_limited_branchclass.gmt',p_threshold=0.01, lfc_threshold=0.585, contrast1='CCI', contrast2='Naive', select_split='-',split_index=1)
  965. }
  967. size_frame <- data.frame(recluster=c(6,6),
  968. celltype_col=c(6,6))
  970. for (group_by_cat in the_groups) {
  972. exclude_list <- c()
  973. for (group in unique([email hidden][[group_by_cat]])) {
  974. temp_subset <- subset([email hidden], [email hidden][[group_by_cat]]==group)
  975. print(group)
  976. print(table(temp_subset$orig.ident))
  977. if (min(table(temp_subset$orig.ident)) < 10) {
  978. exclude_list <- c(exclude_list, group)
  979. }
  980. }
  981. print(exclude_list)
  983. if (group_by_cat == 'celltype_col') {
  984. celltypes <- c('Microglia-1','Microglia-2','Microglia-3','Microglia-4','Microglia-5','Microglia-6','Microglia-7','Microglia-8',
  985. 'Astrocytes','Neurons','Ependymal-1','Ependymal-2','Ependymal-3','Ependymal-4','Ependymal-5',
  986. 'Endothelial-1','Endothelial-2','Neutrophils-1','Neutrophils-2','BCells','TmemoryCells','NKcells')
  987. }
  988. else {
  989. celltypes = c('1','2','3','4','5','6','7','8','9','10')
  990. }
  991. the_pathways <- c()
  993. for (celltype in celltypes) {
  994. gsea_table <- read.table(paste0('pseudobulk/',celltype,'_stat_fgsea.tsv'), sep='\t', header=TRUE)
  995. gsea_table <- subset(gsea_table, padj < 0.05)
  996. if (length(rownames(gsea_table)) == 0) {
  997. exclude_list <- c(exclude_list, celltype)
  998. }
  999. the_pathways <- c(the_pathways, gsea_table$pathway)
  1000. }
  1002. exclude_list <- unique(exclude_list)
  1003. the_pathways <- unique(the_pathways)
  1005. NES_subset <- list()
  1007. for (celltype in celltypes) {
  1008. if (!(celltype %in% exclude_list)) {
  1009. gsea_table <- read.table(paste0('pseudobulk/',celltype,'_stat_fgsea.tsv'), sep='\t', row.names=1, header=TRUE)
  1010. gsea_table <- gsea_table[the_pathways,]
  1012. p_filter <- gsea_table$padj < 0.1
  1014. NES_temp <- gsea_table$NES
  1015. NES_temp[!p_filter] <- 0
  1016. NES_subset[[celltype]] <- NES_temp
  1017. }
  1018. }
  1019. NES_subset <- as.data.frame(NES_subset)
  1020. row.names(NES_subset) <- the_pathways
  1021. NES_subset <- NES_subset[order(row.names(NES_subset)), ]
  1022. write.table(NES_subset, paste0('pseudobulk/mouse_pseudobulk_',group_by_cat,'_NES_subset.txt'), sep='\t', quote=FALSE)
  1024. pdf(paste0('pseudobulk/mouse_pseudobulk_',group_by_cat,'_NES_heatmap.pdf'),width=6, height=6)
  1025. draw(pheatmap(as.matrix(NES_subset),fontsize=6, treeheight_row=10, treeheight_col=10, border_color=NA,cellwidth=8.54,cellheight=8.54))
  1026. dev.off()
  1027. }
  1029. #####BASIC FEATURES#####
  1030. if (!dir.exists('mouse_basic_features')) {
  1031. dir.create('mouse_basic_features')
  1032. }
  1034. cci_seurat_pre_PN <- readRDS("cci_seurat_pre_PN.rds")
  1036. read_count_function_1 <- function(row) {
  1037. the_table <- table(row)
  1038. low_count <- sum(as.numeric(the_table[names(the_table) ==0]))
  1039. non_low_count <- length(row)-low_count
  1040. fraction_non_low <- non_low_count/length(row)
  1041. return(fraction_non_low)
  1042. }
  1044. #What is the common PN in naive
  1045. naive_only <- subset(cci_seurat_pre_PN,[email hidden]$condition=='Naive')
  1046. naive_common_list <- c()
  1047. for (ct in unique([email hidden]$celltype_altcat)) {
  1048. subset_seurat <- subset(naive_only, [email hidden]$celltype_altcat==ct)
  1049. print(ct)
  1050. print(length(colnames(subset_seurat)))
  1051. expression_matrix <- subset_seurat[['RNA']]$counts
  1052. fraction_high <- apply(expression_matrix, 1, read_count_function_1)
  1053. subset_high <- subset(fraction_high, fraction_high > 0.5)
  1055. if (length(naive_common_list) == 0) {
  1056. naive_common_list <- names(subset_high)
  1057. }
  1058. else {
  1059. naive_common_list <- intersect(naive_common_list, names(subset_high))
  1060. }
  1061. write.table(fraction_high,paste0('mouse_basic_features/naive_fraction_',ct,'.txt'), sep='\t', quote=FALSE)
  1062. write.table(subset_high,paste0('mouse_basic_features/naive_fraction>=0.5_',ct,'.txt'), sep='\t', quote=FALSE)
  1063. }
  1065. write.table(naive_common_list, 'mouse_basic_features/coreNaive_subPN.txt',sep='\t',quote=FALSE)
  1067. #What is the common PN in CCI
  1068. cci_only <- subset(cci_seurat_pre_PN,[email hidden]$condition=='CCI')
  1069. cci_common_list <- c()
  1070. for (ct in unique([email hidden]$celltype_altcat)) {
  1071. subset_seurat <- subset(cci_only, [email hidden]$celltype_altcat==ct)
  1072. print(ct)
  1073. print(length(colnames(subset_seurat)))
  1074. expression_matrix <- subset_seurat[['RNA']]$counts
  1075. fraction_high <- apply(expression_matrix, 1, read_count_function_1)
  1076. subset_high <- subset(fraction_high, fraction_high > 0.5)
  1078. if (length(cci_common_list) == 0) {
  1079. cci_common_list <- names(subset_high)
  1080. }
  1081. else {
  1082. cci_common_list <- intersect(cci_common_list, names(subset_high))
  1083. }
  1084. write.table(fraction_high,paste0('mouse_basic_features/cci_fraction_',ct,'.txt'), sep='\t', quote=FALSE)
  1085. write.table(subset_high,paste0('mouse_basic_features/cci_fraction>=0.5_',ct,'.txt'), sep='\t', quote=FALSE)
  1086. }
  1088. write.table(cci_common_list, 'mouse_basic_features/coreCCI_subPN.txt',sep='\t',quote=FALSE)
  1090. #Unique in naive
  1091. select_markers <- function(seurat_object, name_column, greater_than, less_than, the_folder, read_file) {
  1092. unique_list <- c()
  1093. for (ct1 in unique([email hidden][[name_column]])) {
  1094. list_1 <- read.table(paste0(the_folder,'/',read_file,'_',ct1,'.txt'),sep='\t')
  1095. list_1 <- row.names(subset(list_1, list_1$x >= greater_than))
  1096. print(list_1)
  1097. common_list = c()
  1098. for (ct2 in unique([email hidden][[name_column]])) {
  1099. if (!(ct1 == ct2)) {
  1100. list_2 <- read.table(paste0(the_folder,'/',read_file,'_',ct2,'.txt'),sep='\t')
  1101. list_2 <- row.names(subset(list_2, list_2$x <= less_than))
  1102. if (length(common_list)==0) {
  1103. common_list <- list_2
  1104. }
  1105. else {
  1106. common_list <- intersect(common_list, list_2)
  1107. }
  1108. }
  1109. }
  1110. unique_list[[ct1]] <- list_1[(list_1 %in% common_list)]
  1111. }
  1112. return(unique_list)
  1113. }
  1115. cci_seurat_lt_pca_50_merged <- readRDS('mouse_recluster/cci_seurat_lt_pca_50_recluster.rds')
  1116. Idents(cci_seurat_lt_pca_50_merged) <- 'condition'
  1118. #Unique in naive
  1119. ##Here we need to pull something in from later analysis to be able to sense check that these genes are actually up in naive
  1120. naive_unique_list <- select_markers(naive_only, "celltype_altcat", 0.6, 0.4,'mouse_basic_features', 'naive_fraction')
  1121. flat_naive_unique_list <- as.character(unlist(naive_unique_list))
  1123. naive_markers <- FindMarkers(cci_seurat_lt_pca_50_merged, ident.1='Naive', ident.2='CCI', thresh.use=-Inf, min.pct=0, features=flat_naive_unique_list) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  1124. write.table(naive_markers, file='mouse_basic_features/naive_markers6040.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  1126. exclude_naive_markers_1 <- row.names(subset(naive_markers, naive_markers$avg_log2FC <= 0))
  1127. exclude_naive_markers_2 <- flat_naive_unique_list[!(flat_naive_unique_list %in% row.names(naive_markers))]
  1128. exclude_naive_markers <- c(exclude_naive_markers_1,exclude_naive_markers_2)
  1129. new_flat_naive_unique_list <- flat_naive_unique_list[!(flat_naive_unique_list %in% exclude_naive_markers)]
  1131. naive_unique_list_6040 = c()
  1132. for (ct1 in names(naive_unique_list)) {
  1133. temp <- naive_unique_list[[ct1]]
  1134. print(temp)
  1135. naive_unique_list_6040[[ct1]] <- temp[!(temp %in% exclude_naive_markers)]
  1136. }
  1138. ##Now I am doing 50/50 so I can get my freqspec for naive
  1139. naive_unique_list <- select_markers(naive_only, "celltype_altcat", 0.5, 0.5,'mouse_basic_features', 'naive_fraction')
  1140. flat_naive_unique_list <- as.character(unlist(naive_unique_list))
  1142. naive_markers <- FindMarkers(cci_seurat_lt_pca_50_merged, ident.1='Naive', ident.2='CCI', thresh.use=-Inf, min.pct=0, features=flat_naive_unique_list) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  1143. write.table(naive_markers, file='mouse_basic_features/naive_markers5050.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  1145. exclude_naive_markers_1 <- row.names(subset(naive_markers, naive_markers$avg_log2FC <= 0))
  1146. exclude_naive_markers_2 <- flat_naive_unique_list[!(flat_naive_unique_list %in% row.names(naive_markers))]
  1147. exclude_naive_markers <- c(exclude_naive_markers_1,exclude_naive_markers_2)
  1148. new_flat_naive_unique_list <- flat_naive_unique_list[!(flat_naive_unique_list %in% exclude_naive_markers)]
  1150. naive_unique_list_5050 = c()
  1151. for (ct1 in names(naive_unique_list)) {
  1152. temp <- naive_unique_list[[ct1]]
  1153. print(temp)
  1154. naive_unique_list_5050[[ct1]] <- temp[!(temp %in% exclude_naive_markers)]
  1155. }
  1157. out_data <- data.frame(freq=c(), spec=c(), celltype=c())
  1158. for (ct1 in names(naive_unique_list_5050)) {
  1159. list_1 <- read.table(paste0('mouse_basic_features/naive_fraction_',ct1,'.txt'),sep='\t')
  1160. list_1 <- subset(list_1, row.names(list_1) %in% naive_unique_list_5050[[ct1]])
  1162. check = 0
  1163. for (ct2 in unique([email hidden]$celltype_altcat)) {
  1164. if (!ct2==ct1) {
  1165. list_2 <- read.table(paste0('mouse_basic_features/naive_fraction_',ct2,'.txt'),sep='\t')
  1166. list_2 <- subset(list_2, row.names(list_2) %in% naive_unique_list_5050[[ct1]])
  1167. names(list_2) <- ct2
  1168. if (check == 0) {
  1169. running_list <- data.frame(list_2)
  1170. }
  1171. else {
  1172. running_list <- cbind(running_list, list_2)
  1173. }
  1174. check = 1
  1175. }
  1176. }
  1177. print(running_list)
  1178. max_list <- apply(running_list, 1, max)
  1179. out_data <- rbind(out_data, data.frame(freq=list_1$x, spec=max_list, celltype=rep(ct1, length(max_list))))
  1180. }
  1181. write.table(out_data, file='freqspec_naive.txt', sep='\t',quote=FALSE, col.names=NA)
  1183. fs <- read.table("mouse_basic_features/freqspec_naive_unique_table.txt", header=TRUE, sep='\t')
  1185. color_table <- c('En'="#E41A1CFF",
  1186. "Ep"="#377EB8FF",
  1187. "Nr"="#4DAF4AFF",
  1188. "OPC"="#984EA3FF",
  1189. "M"="#FECEA8FF",
  1190. "Olg"="#F781BFFF",
  1191. "As"="#A65628FF",
  1192. "L"=,"#FF7F00FF",
  1193. "Nt"="#321701")
  1195. the_theme <- theme(plot.title=element_blank(),
  1196. panel.background=element_blank(),
  1197. panel.border = element_blank(),
  1198. panel.grid=element_blank(),
  1199. axis.line=element_line(linewidth=0.3),
  1200. axis.ticks=element_line(linewidth=0.3),
  1201. axis.text=element_text(size=6),
  1202. legend.position='none',
  1203. title=element_blank(),
  1204. axis.title=element_blank())
  1206. fs_plot <- ggplot(fs, aes(x=spec, y=freq, label=label, colour=celltype_label))+geom_point(size=0.5)+geom_vline(xintercept=0.4, linetype=2, linewidth=0.3)+geom_hline(yintercept=0.6, linetype=2, linewidth=0.3)+
  1207. geom_text(size=2)+scale_x_continuous(breaks=c(0.1,0.2,0.3,0.4,0.5), limits=c(0.05,0.55))+scale_colour_manual(values=color_table, na.value='#CCCCCC')+
  1208. the_theme
  1209. ggsave('fs_plot_naive.pdf', fs_plot, width=6, height=6, units='cm', dpi=600)
  1211. ####################
  1212. cci_naive <- subset(cci_seurat_lt_pca_50_merged, [email hidden]$condition=='Naive')
  1213. Idents(cci_naive) <- "celltype_altcat"
  1214. feature_list <- c('EPN2','GJA1','HEMK1','FNBP1L') #highest log2FC for marker for each cell type that has a marker#
  1216. dev.new(width=8, height=6, noRStudioGD = TRUE, units = "cm")
  1217. RidgePlot(cci_naive, features=feature_list)
  1218. ####################
  1221. ###########Unique in cci
  1222. ##Here we need to pull something in from later analysis to be able to sense check that these genes are actually up in cci
  1223. cci_unique_list <- select_markers(cci_only, "celltype_altcat", 0.6, 0.4,'mouse_basic_features', 'cci_fraction')
  1224. flat_cci_unique_list <- as.character(unlist(cci_unique_list))
  1226. cci_markers <- FindMarkers(cci_seurat_lt_pca_50_merged, ident.1='CCI', ident.2='Naive', thresh.use=-Inf, min.pct=0, features=flat_cci_unique_list) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  1227. write.table(cci_markers, file='mouse_basic_features/cci_markers6040.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  1229. exclude_cci_markers_1 <- row.names(subset(cci_markers, cci_markers$avg_log2FC <= 0))
  1230. exclude_cci_markers_2 <- flat_cci_unique_list[!(flat_cci_unique_list %in% row.names(cci_markers))]
  1231. exclude_cci_markers <- c(exclude_cci_markers_1,exclude_cci_markers_2)
  1232. new_flat_cci_unique_list <- flat_cci_unique_list[!(flat_cci_unique_list %in% exclude_cci_markers)]
  1234. cci_unique_list_6040 = c()
  1235. for (ct1 in names(cci_unique_list)) {
  1236. temp <- cci_unique_list[[ct1]]
  1237. cci_unique_list_6040[[ct1]] <- temp[!(temp %in% exclude_cci_markers)]
  1238. }
  1240. ##Now I am doing 50/50 so I can get my my freqspec
  1241. cci_unique_list <- select_markers(cci_only, "celltype_altcat", 0.5, 0.5,'mouse_basic_features', 'cci_fraction')
  1242. flat_cci_unique_list <- as.character(unlist(cci_unique_list))
  1244. cci_markers <- FindMarkers(cci_seurat_lt_pca_50_merged, ident.1='CCI', ident.2='Naive', thresh.use=-Inf, min.pct=0, features=flat_cci_unique_list) #first calculate all the wilcox scores. Note that we use the library presto to speed up computation as suggested by Seurat.
  1245. write.table(cci_markers, file='mouse_basic_features/cci_markers5050.txt', quote=FALSE, sep='\t', row.names=TRUE, col.names=NA)
  1247. exclude_cci_markers_1 <- row.names(subset(cci_markers, cci_markers$avg_log2FC <= 0))
  1248. exclude_cci_markers_2 <- flat_cci_unique_list[!(flat_cci_unique_list %in% row.names(cci_markers))]
  1249. exclude_cci_markers <- c(exclude_cci_markers_1,exclude_cci_markers_2)
  1250. new_flat_cci_unique_list <- flat_cci_unique_list[!(flat_cci_unique_list %in% exclude_cci_markers)]
  1252. ###now 5050
  1253. cci_unique_list_5050 = c()
  1254. for (ct1 in names(cci_unique_list)) {
  1255. temp <- cci_unique_list[[ct1]]
  1256. cci_unique_list_5050[[ct1]] <- temp[!(temp %in% exclude_cci_markers)]
  1257. }
  1259. out_data <- data.frame(freq=c(), spec=c(), celltype=c())
  1260. for (ct1 in names(cci_unique_list_5050)) {
  1261. list_1 <- read.table(paste0('mouse_basic_features/cci_fraction_',ct1,'.txt'),sep='\t')
  1262. list_1 <- subset(list_1, row.names(list_1) %in% cci_unique_list_5050[[ct1]])
  1264. check = 0
  1265. for (ct2 in unique([email hidden]$celltype_altcat)) {
  1266. if (!ct2==ct1) {
  1267. list_2 <- read.table(paste0('mouse_basic_features/cci_fraction_',ct2,'.txt'),sep='\t')
  1268. list_2 <- subset(list_2, row.names(list_2) %in% cci_unique_list_5050[[ct1]])
  1269. names(list_2) <- ct2
  1270. if (check == 0) {
  1271. running_list <- data.frame(list_2)
  1272. }
  1273. else {
  1274. running_list <- cbind(running_list, list_2)
  1275. }
  1276. check = 1
  1277. }
  1278. }
  1279. print(running_list)
  1280. max_list <- apply(running_list, 1, max)
  1281. out_data <- rbind(out_data, data.frame(freq=list_1$x, spec=max_list, celltype=rep(ct1, length(max_list))))
  1282. }
  1283. write.table(out_data, file='freqspec_cci.txt', sep='\t',quote=FALSE, col.names=NA)
  1285. fs <- read.table("mouse_basic_features/freqspec_cci_unique_table.txt", header=TRUE, sep='\t')
  1287. the_theme <- theme(plot.title=element_blank(),
  1288. panel.background=element_blank(),
  1289. panel.border = element_blank(),
  1290. panel.grid=element_blank(),
  1291. axis.line=element_line(linewidth=0.3),
  1292. axis.ticks=element_line(linewidth=0.3),
  1293. axis.text=element_text(size=6),
  1294. legend.position='none',
  1295. title=element_blank(),
  1296. axis.title=element_blank())
  1298. fs_plot <- ggplot(fs, aes(x=spec, y=freq, label=label, colour=celltype_label))+geom_point(size=0.5)+geom_vline(xintercept=0.4, linetype=2, linewidth=0.3)+geom_hline(yintercept=0.6, linetype=2, linewidth=0.3)+
  1299. geom_text(size=2, position=position_nudge(x=0.02,y=0.02))+scale_x_continuous(breaks=c(0.1,0.2,0.3,0.4,0.5), limits=c(0.05,0.55))+scale_colour_manual(values=color_table, na.value='#cccccc')+
  1300. the_theme
  1301. ggsave('fs_plot_cci.pdf', fs_plot, width=6, height=6, units='cm', dpi=600)
  1303. #############
  1305. cci_cci <- subset(cci_seurat_lt_pca_50_merged, [email hidden]$condition=='CCI')
  1306. Idents(cci_cci) <- "celltype_altcat"
  1307. feature_list <- c('TRAF1','A2M','KLHL2','N4BP1','TTC9B','UCHL1') #*top for each cell type that has a marker, except I also added UCHL1*#
  1309. dev.new(width=8, height=6, noRStudioGD = TRUE, units = "cm")
  1310. RidgePlot(cci_cci, features=feature_list)
  1312. #naive
  1313. naive_subPN <- data.frame(
  1314. group = factor(c("TRANSL","ALP", "EXTRA", "ER","UPS","CYTONUC",'CYTONUC/ALP','TRANSL/ALP','UPS/ALP','UPS/PNREG'), levels=c("TRANSL", "ALP", "EXTRA", "ER","UPS","CYTONUC",'CYTONUC/ALP','TRANSL/ALP','UPS/ALP','UPS/PNREG')),
  1315. value = c(46,3,2,2,2,1,2,1,1,1),
  1316. colors = c('#377eb8','#984ea3','#ffff33','#f781bf','#a65628','#ff7f00','#fdae61','#ece7f2','#cc9966','#ffcc99'))
  1318. the_theme <- theme(plot.title=element_blank(),
  1319. panel.background=element_blank(),
  1320. panel.border = element_blank(),
  1321. panel.grid=element_blank(),
  1322. axis.line=element_blank(),
  1323. axis.ticks=element_blank(),
  1324. axis.text=element_blank(),
  1325. legend.title=element_blank(),legend.text=element_text(size=7.5, colour='black'),
  1326. legend.position='none',
  1327. title=element_blank(),
  1328. axis.title=element_blank())
  1329. dev.new(width=1.5, height=1.5, noRStudioGD = TRUE, units = "cm")
  1330. ggplot(naive_subPN, aes(x="",y=value,fill=group))+geom_bar(width=1, stat='identity')+coord_polar("y",start=0)+scale_fill_manual(values=naive_subPN$colors)+the_theme
  1332. #coreCCI
  1333. coreCCI_subPN <- data.frame(
  1334. group = factor(c("TRANSL","ALP", "EXTRA", "ER","UPS","CYTONUC",'CYTONUC/ALP','TRANSL/ALP','UPS/ALP','UPS/PNREG'), levels=c("TRANSL","ALP", "EXTRA", "ER","UPS","CYTONUC",'CYTONUC/ALP','TRANSL/ALP','UPS/ALP','UPS/PNREG')),
  1335. value = c(27,3,3,1,6,0,1,1,1,1),
  1336. colors = c('#377eb8','#984ea3','#ffff33','#f781bf','#a65628','#ff7f00','#fdae61','#ece7f2','#cc9966','#ffcc99'))
  1338. the_theme <- theme(plot.title=element_blank(),
  1339. panel.background=element_blank(),
  1340. panel.border = element_blank(),
  1341. panel.grid=element_blank(),
  1342. axis.line=element_blank(),
  1343. axis.ticks=element_blank(),
  1344. axis.text=element_blank(),
  1345. legend.title=element_blank(),legend.text=element_text(size=7.5, colour='black'),
  1346. legend.position='none',
  1347. title=element_blank(),
  1348. axis.title=element_blank())
  1349. dev.new(width=1.5, height=1.5, noRStudioGD = TRUE, units = "cm")
  1350. ggplot(coreCCI_subPN, aes(x="",y=value,fill=group))+geom_bar(width=1, stat='identity')+coord_polar("y",start=0)+scale_fill_manual(values=coreCCI_subPN$colors)+the_theme
  1352. naive_common_list <- read.table('mouse_basic_features/coreNaive_subPN.txt',sep='\t')$x
  1353. cci_common_list <- read.table('mouse_basic_features/coreCCI_subPN.txt',sep='\t')$x
  1355. dev.new(width=2.1, height=1.6, noRStudioGD = TRUE, units = "cm")
  1356. the_theme <- theme(plot.title=element_blank(),
  1357. panel.background=element_blank(),
  1358. panel.border = element_blank(),
  1359. panel.grid=element_blank(),
  1360. axis.line=element_blank(),
  1361. axis.ticks=element_blank(),
  1362. axis.text=element_blank(),
  1363. legend.title=element_blank(),legend.text=element_text(size=7.5, colour='black'),
  1364. legend.position='none',
  1365. title=element_blank(),
  1366. axis.title=element_blank())
  1368. #values=c('Naive'='#67a9cf','CCI'='#ef8a62','Naive:CCI'='#caa099') just do this in illustrator
  1369. ggVennDiagram(x=list('Naive'=naive_common_list, 'CCI'=cci_common_list),category.names=c('Naive','CCI'),fontfamily='Arial',label='count',label_size=5,label_alpha=0)+coord_flip()+the_theme
  1371. #####FIGURE 3#####
  1372. #CCI_1 CCI_2 CCI_3 Naive_1 Naive_2 Naive_3 Naive_4
  1373. #Astrocytes 130 97 71 112 213 22 152
  1374. #BCells 601 893 902 40 53 59 59
  1375. #Endothelial_1 147 187 326 213 627 187 339
  1376. #Endothelial_2 65 82 133 106 293 76 169
  1377. #Ependymal_1 467 564 960 392 726 1204 2130
  1378. #Ependymal_2 80 78 375 161 265 283 861
  1379. #Ependymal_3 80 89 131 84 180 300 328
  1380. #Ependymal_4 36 135 233 1 3 4 15
  1381. #Ependymal_5 0 15 19 4 48 18 76
  1382. #Microglia_1 1111 1638 2112 2186 2362 3494 1869
  1383. #Microglia_2 606 1351 1214 741 1132 1127 973
  1384. #Microglia_3 954 1261 1410 300 319 407 240
  1385. #Microglia_4 217 887 640 32 36 50 29
  1386. #Microglia_5 414 604 663 5 27 10 4
  1387. #Microglia_6 64 69 81 73 329 183 174
  1388. #Microglia_7 8 3 24 16 48 24 117
  1389. #Microglia_8 9 21 39 22 48 27 37
  1390. #Neurons 129 59 67 27 23 38 32
  1391. #Neutrophils_1 317 1071 606 8 9 7 12
  1392. #Neutrophils_2 28 135 60 7 11 11 12
  1393. #NKcells 162 125 114 41 51 63 99
  1394. #Oligodendrocytes 21 43 32 11 88 19 28
  1395. #TmemoryCells 107 60 91 0 2 1 1

Mouse Code.R at commit ff25fc9, no license · at the source

Overview

Authors: Jonathan Lefkowitz1, Milo Taylor1, Aditya Kumar2,3,4, Ruchira M Jha2,3,4,5, Dhivyaa Rajasundaram6, Prajwal Ciryam1,7
ORCID iDs: Jonathan Lefkowitz, Prajwal Ciryam
  1. Department of Neurology, University of Maryland School of Medicine, Baltimore, MD, USA
  2. Departments of Neurology, Barrow Neurological Institute, Phoenix, AZ, USA
  3. Department of Neurological Surgery, Barrow Neurological Institute, Phoenix, AZ, USA
  4. Department of Translational Neuroscience, Barrow Neurological Institute, Phoenix, AZ, USA
  5. Safar Center for Resuscitation Research, UPMC Children's Hospital of Pittsburgh, School of Medicine, University of Pittsburgh, Pittsburgh, PA, USA
  6. Division of Health Informatics, Department of Pediatrics, UPMC Children's Hospital of Pittsburgh, School of Medicine, University of Pittsburgh, Pittsburgh, PA, USA
  7. Shock Trauma Neurocritical Care, Program in Trauma, R Adams Cowley Shock Trauma Center, University of Maryland Medical Center, Baltimore, MD, USA
Institutions: University of Maryland, Baltimore (United States); Barrow Neurological Institute (United States); Children's Hospital of Pittsburgh (United States); University of Maryland Medical Center (United States); International Trauma Anesthesia and Critical Care Society (United States)
Journal: Philosophical transactions of the Royal Society of London. Series B, Biological sciences, volume 381, issue 1960, article 20250264
Dates: received 2 March 2026; accepted 20 June 2026; published online 1 October 2026; in print October 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1098/rstb.2025.0264 · PMID 42817638 · PMCID PMC13628056 · OpenAlex W7215098471
Open access: hybrid, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: genetics / omics (modality), human (organism), mouse (organism), traumatic brain injury (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning
Keywords: proteostasis, traumatic brain injury, inflammation
MeSH: Brain Injuries, Traumatic*, Proteostasis*, Transcriptome*, Animals, Humans, Mice, Proteotoxic Stress (* major topic)
Topic: Genetic Neurodegenerative Diseases (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: American Academy of Neurology; Passano Foundation; Barrow Neurological Foundation; Neurocritical Care Society; Burroughs-Wellcome Scholar-AHEAD Program; NIH National Institute of Neurological Disorders and Stroke (R01NS115815, R21NS131689); Chuck Noll Foundation; GEn1E Lifesciences; Henry M. Jackson Foundation for the Advancement of Military Medicine
Citations: not cited yet (Europe PMC); 50 references in the paper

Abstract

Traumatic brain injury (TBI) triggers competing damage and repair processes that can result in chronic neurodegenerative protein conformational disorders. Acute alterations in proteostasis have been reported after injury, but a comprehensive understanding of changes to the proteostasis network (PN) has been limited by the complex cellular pathology of TBI. The recent enumeration of the human PN and advances in single-cell transcriptomics make it possible to delineate PN dynamics in multifactorial diseases such as TBI. Here, we leverage archival mouse and human transcriptomic data to identify cell-type-specific Sub-PNs in the brain and their alteration by TBI. We find characteristic PN gene expression profiles in many cell types and subtypes, reflecting tremendous diversity in the enrichment of specific PN processes. TBI induces a broad transcriptional remodelling of translation, the ubiquitin–proteasome system and nuclear proteostasis machinery across cell types. In addition to this shared response, numerous cell-type-specific responses occur, some of which are conserved between mouse and human, despite differences in injury mechanism and timing. Our extensive characterization of the PN transcriptome across cell type and disease reveals potentially modifiable mechanisms of PN alteration in TBI. More broadly, our results provide a framework for elucidating changes to the PN in complex diseases.

This article is part of the Theo Murphy meeting issue ‘ProteostaSys: a systems view of proteostasis’.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 10 matches between paragraphs and lines of code.

Ciryam-Lab/Lefkowitz-Proteostasis-2026

License: none: the authors keep all their rights
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: ff25fc9aa97e585cd71ea47710d7d811441b99bf, 5 August 2026
Languages: R (4), Python (1)
Size: 6 files, 5 scripts
Software Heritage: not archived
Found in: “Data accessibility”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (4 files), circlize (3 files), ComplexHeatmap (3 files), igraph (3 files), Seurat (3 files), data.table (2 files), DESeq2 (2 files)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
5 files
At the source: github.com/Ciryam-Lab/Lefkowitz-Proteostasis-2026/

doi:10.6084/m9.figshare.c.8672767

License: none: the authors keep all their rights
State: unreachable at the last attempt, verified on 14 October 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 20 checks, the latest on 14 October 2026: unreachable at the last attempt (HTTP 202)
  • 14 October 2026: unreachable at the last attempt (HTTP 202)
  • 14 October 2026: unreachable at the last attempt (HTTP 202)
  • 14 October 2026: unreachable at the last attempt (HTTP 202)
  • 13 October 2026: unreachable at the last attempt (HTTP 202)
  • 13 October 2026: unreachable at the last attempt (HTTP 202)
  • 8 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
  • 6 October 2026: unreachable at the last attempt (HTTP 202)
At the source: doi.org/10.6084/m9.figshare.c.8672767

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 5 scripts, each with its path and the digest of its content;
  • 10 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data accessibility

Code relevant to this manuscript is available at https://github.com/Ciryam-Lab/Lefkowitz-Proteostasis-2026/. The data anlyzed in this article are available on the Gene Expression Omnibus (GEO), accession numbers GSE209552 (https://www.omicsdi.org/dataset/geo/GSE209552) and GSE269748 (https://usebiotransfer.org/tools/geo-finder/GSE269748/).

Supplementary material is available online [50].

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 6 October 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 3 keywords, 7 MeSH terms, 9 funders, 50 references.

Cite

This paper

Lefkowitz, J., Taylor, M., Kumar, A., Jha, R. M., Rajasundaram, D., & Ciryam, P. (2026). Traumatic brain injury induces cell-type-specific remodelling of the proteostasis network. Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 381(1960), 20250264. https://doi.org/10.1098/rstb.2025.0264

BibTeX

@article{lefkowitz2026traumatic,
author = {Lefkowitz, Jonathan and Taylor, Milo and Kumar, Aditya and Jha, Ruchira M and Rajasundaram, Dhivyaa and Ciryam, Prajwal},
title = {{Traumatic brain injury induces cell-type-specific remodelling of the proteostasis network}},
journal = {Philosophical transactions of the Royal Society of London. Series B, Biological sciences},
year = {2026},
month = oct,
volume = {381},
number = {1960},
pages = {20250264},
publisher = {The Royal Society},
issn = {0962-8436},
doi = {10.1098/rstb.2025.0264},
url = {https://doi.org/10.1098/rstb.2025.0264},
pmid = {42817638},
pmcid = {PMC13628056}
}

RIS

TY - JOUR
AU - Lefkowitz, Jonathan
AU - Taylor, Milo
AU - Kumar, Aditya
AU - Jha, Ruchira M
AU - Rajasundaram, Dhivyaa
AU - Ciryam, Prajwal
TI - Traumatic brain injury induces cell-type-specific remodelling of the proteostasis network
T2 - Philosophical transactions of the Royal Society of London. Series B, Biological sciences
J2 - Philos Trans R Soc Lond B Biol Sci
PY - 2026
DA - 2026/10/01
VL - 381
IS - 1960
SP - 20250264
SN - 0962-8436
PB - The Royal Society
DO - 10.1098/rstb.2025.0264
UR - https://doi.org/10.1098/rstb.2025.0264
LA - en
ER -

CSL-JSON

{
"id": "10.1098/rstb.2025.0264",
"type": "article-journal",
"title": "Traumatic brain injury induces cell-type-specific remodelling of the proteostasis network",
"container-title": "Philosophical transactions of the Royal Society of London. Series B, Biological sciences",
"author": [
{
"family": "Lefkowitz",
"given": "Jonathan"
},
{
"family": "Taylor",
"given": "Milo"
},
{
"family": "Kumar",
"given": "Aditya"
},
{
"family": "Jha",
"given": "Ruchira M"
},
{
"family": "Rajasundaram",
"given": "Dhivyaa"
},
{
"family": "Ciryam",
"given": "Prajwal"
}
],
"container-title-short": "Philos Trans R Soc Lond B Biol Sci",
"volume": "381",
"issue": "1960",
"page": "20250264",
"DOI": "10.1098/rstb.2025.0264",
"PMID": "42817638",
"PMCID": "PMC13628056",
"ISSN": "0962-8436",
"publisher": "The Royal Society",
"URL": "https://doi.org/10.1098/rstb.2025.0264",
"language": "en",
"issued": {
"date-parts": [
[
2026,
10,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/ccs3.70006 [code]
Evaluation of altered cell-cell communication between glia and neurons in the hippocampus of 3xTg-AD mice at two time points.
Journal: Journal of cell communication and signaling
In common: circlize, ComplexHeatmap, DESeq2, 4 other tools, genetics / omics, mouse, cellular / molecular, 2 references
[2] doi:10.1016/j.celrep.2023.113070 [code]
TMEM127 suppresses tumor development by promoting RET ubiquitination, positioning, and degradation.
Journal: Cell reports
In common: igraph, Seurat, data.table, 1 other tool, genetics / omics, mouse, cellular / molecular, 1 reference, author Jonathan Lefkowitz
[3] doi:10.1007/s00401-025-02913-3 [code]
Single-cell transcriptomic landscape of the neuroimmune compartment in amyotrophic lateral sclerosis brain and spinal cord.
Journal: Acta neuropathologica
In common: circlize, ComplexHeatmap, DESeq2, 4 other tools, genetics / omics, cellular / molecular, 1 reference
[4] doi:10.1084/jem.20232357 [code]
Microglia and CD8+ T cell activation precede neuronal loss in a murine model of spastic paraplegia 15.
Journal: The Journal of experimental medicine
In common: circlize, ComplexHeatmap, igraph, 3 other tools, mouse, cellular / molecular, 2 references
[5] doi:10.1126/sciadv.adf6251 [code]
Single-cell genomics reveals region-specific developmental trajectories underlying neuronal diversity in the human hypothalamus.
Journal: Science advances
In common: circlize, ComplexHeatmap, igraph, 3 other tools, genetics / omics, mouse, 1 reference
[6] doi:10.1038/s41380-026-03585-5 [code]
Multiomics analysis identifies VPA-induced changes in neural progenitor cells, ventricular-like regions, and cellular microenvironment in dorsal forebrain organoids.
Journal: Molecular psychiatry
In common: ComplexHeatmap, igraph, Seurat, 2 other tools, genetics / omics, mouse, 2 references
[7] doi:10.1016/j.ebiom.2026.106324 [code]
Refined single-cell profiling captures a CCR5&lt;sup&gt;high&lt;/sup&gt; CD4&lt;sup&gt;+&lt;/sup&gt; cytotoxic T-cell precursor in multiple sclerosis.
Journal: EBioMedicine
In common: DESeq2, igraph, Seurat, 2 other tools, genetics / omics, cellular / molecular, 2 references
[8] doi:10.1016/j.neuron.2025.09.022 [code]
Integrated omics reveals disease-associated radial glia-like cells with epigenetically dysregulated interferon response in multiple sclerosis.
Journal: Neuron
In common: circlize, ComplexHeatmap, Seurat, 2 other tools, genetics / omics, cellular / molecular, 1 reference
[9] doi:10.1093/molbev/msaf189 [code]
A General Principle of Neuronal Evolution Reveals a Human-Accelerated Neuron Type Potentially Underlying the High Prevalence of Autism in Humans.
Journal: Molecular biology and evolution
In common: DESeq2, igraph, Seurat, 2 other tools, genetics / omics, cellular / molecular, 1 reference
[10] doi:10.1038/s41559-023-02186-7 [code]
Comparative single-cell transcriptomic analysis of primate brains highlights human-specific regulatory evolution.
Journal: Nature ecology & evolution
In common: DESeq2, igraph, Seurat, 2 other tools, genetics / omics, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1098/rstb.2025.0264.svg)](https://openscicode.org/badge/doi/10.1098/rstb.2025.0264)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1098/rstb.2025.0264.svg
   :target: https://openscicode.org/badge/doi/10.1098/rstb.2025.0264
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1098/rstb.2025.0264"><img src="https://openscicode.org/badge/doi/10.1098/rstb.2025.0264.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.