OSCR

A Systems-Level Transcriptomic Framework Identifies Shared Cellular Hubs in Osteoarthritis and Alzheimer's Disease.

Code ↔ Paper

19 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 19 matches
  1. [1] § Results › Single-cell transcriptomic analysis of OA cartilage highlights Fibrochondrocyte_1 as a putative sender-like hub within inferred intratissue networks ↔ singlecell code OA.R, lines 173–239 · score 0.95 · reparative chondrocytes, proliferation chondrocytes, homeostatic chondrocytes, ProC, HomC_1, HomC_2
  2. [2] § Materials and Methods › Single-cell data analysis ↔ singlecell code OA.R, lines 84–170 · score 0.94 · FindVariableFeatures, NormalizeData, RunPCA, ScaleData, novelty score, Seurat
  3. [3] § Materials and Methods › Single-cell data analysis ↔ singlecell codeAD.R, lines 115–162 · score 0.88 · FindVariableFeatures, NormalizeData, RunPCA, ScaleData, batch, resolution
  4. [4] § Results › Neurovascular- and endothelial-associated pathway activity in OA and AD ↔ singlecell codeAD.R, lines 699–705 · score 0.85 · leukocyte transendothelial migration, VEGF signaling pathway, Wnt signaling pathway, ECM receptor interaction, adherens junction, cell
  5. [5] § Results › Neurovascular- and endothelial-associated pathway activity in OA and AD ↔ singlecell code OA.R, lines 662–667 · score 0.84 · leukocyte transendothelial migration, VEGF signaling pathway, Wnt signaling pathway, ECM receptor interaction, adherens junction, OA
  6. [6] § Results › Single-cell transcriptomic analysis of OA cartilage highlights Fibrochondrocyte_1 as a putative sender-like hub within inferred intratissue networks ↔ singlecell code OA.R, lines 361–420 · score 0.84 · quantitative centrality metrics, preHTC_2, RepC_1, preInfC_1, cell communication, preHTC_1
  7. [7] § Results › Neurovascular- and endothelial-associated pathway activity in OA and AD ↔ singlecell code OA.R, lines 662–667 · score 0.80 · leukocyte transendothelial migration, ECM receptor interaction, adherens junction, preHTC_1, InfC_1, pathway
  8. [8] § Materials and Methods › Differential expression analysis ↔ DIFF code.R, lines 137–193 · score 0.79 · log2 fold change, cor.test, phyper, Spearman, concordance, hypergeometric
  9. [9] § Results › Single-cell transcriptomic analysis of AD cortex highlights Oligodendrocyte_3 as a putative receiver-like hub within inferred intratissue networks ↔ singlecell codeAD.R, lines 241–304 · score 0.77 · microglia_1, excitatory_neuron_5, oligodendrocyte_1, endothelial_cell_1, astrocytes_1, oligodendrocyte_3
  10. [10] § Results › Neurovascular- and endothelial-associated pathway activity in OA and AD ↔ singlecell codeAD.R, lines 699–705 · score 0.77 · adherens junction, VEGF signaling, Wnt signaling, endothelial_cell_1, astrocytes_1, box
  11. [11] § Results › Single-cell transcriptomic analysis of OA cartilage highlights Fibrochondrocyte_1 as a putative sender-like hub within inferred intratissue networks ↔ singlecell code OA.R, lines 173–239 · score 0.74 · ABI3BP, COL1A1, preHTC_1, InfC_1, CRTAC1, MMP2
  12. [12] § Results › Single-cell transcriptomic analysis of AD cortex highlights Oligodendrocyte_3 as a putative receiver-like hub within inferred intratissue networks ↔ singlecell codeAD.R, lines 241–304 · score 0.73 · microglia_1, excitatory_neuron_5, oligodendrocyte_1, endothelial_cell_1, astrocytes_1, oligodendrocyte_3
  13. [13] § Materials and Methods › Single-cell data analysis ↔ singlecell codeAD.R, lines 76–113 · score 0.65 · AD13, AD19, AD6, AD8, NC11, NC16
  14. [14] § Results › Single-cell transcriptomic analysis of AD cortex highlights Oligodendrocyte_3 as a putative receiver-like hub within inferred intratissue networks ↔ singlecell codeAD.R, lines 368–438 · score 0.64 · quantitative centrality metrics, oligodendrocyte_1, astrocytes_1, oligodendrocyte_3, interaction, networks
  15. [15] § Results › Single-cell transcriptomic analysis of OA cartilage highlights Fibrochondrocyte_1 as a putative sender-like hub within inferred intratissue networks ↔ singlecell code OA.R, lines 304–346 · score 0.59 · ProC_1, preInfC_1, preHTC_1, expansion, OA, RNA
  16. [16] § Results › Single-cell transcriptomic analysis of AD cortex highlights Oligodendrocyte_3 as a putative receiver-like hub within inferred intratissue networks ↔ singlecell codeAD.R, lines 306–352 · score 0.59 · donor aware, microglia_1, oligodendrocyte_3, expanded, abundance, AD
  17. [17] § Results › Single-cell transcriptomic analysis of AD cortex highlights Oligodendrocyte_3 as a putative receiver-like hub within inferred intratissue networks ↔ singlecell codeAD.R, lines 495–541 · score 0.56 · excitatory_neuron_5, oligodendrocyte_3, endothelial cells, subtypes, receiver, Astrocytes
  18. [18] § Materials and Methods › Single-cell data analysis ↔ singlecell code OA.R, lines 361–420 · score 0.51 · CellChat, cell communication, databases, interactions, networks, Seurat
  19. [19] § Materials and Methods › Single-cell data analysis ↔ singlecell codeAD.R, lines 368–438 · score 0.51 · CellChat, cell communication, databases, interactions, networks, Seurat

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 705 lines · 33 KB · no license · 10 matches

  1. ########Data Processing, Filtering, and Integration#########
  2. #load packages#
  3. library(Seurat)
  4. library(multtest)
  5. library(dplyr)
  6. library(ggplot2)
  7. library(patchwork)
  8. library(SeuratData)
  9. library(stringr)
  11. ###load_data###
  12. #load_data+creat_seurat#
  13. rm(list = ls())
  14. AD1 <- Read10X(data.dir = "GSE157827/AD1/")
  15. AD1 <- CreateSeuratObject(counts = AD1, project = "AD1", min.cells = 3, min.features = 200)
  16. AD2 <- Read10X(data.dir = "GSE157827/AD2/")
  17. AD2 <- CreateSeuratObject(counts = AD2, project = "AD2", min.cells = 3, min.features = 200)
  18. AD4 <- Read10X(data.dir = "GSE157827/AD4/")
  19. AD4 <- CreateSeuratObject(counts = AD4, project = "AD4", min.cells = 3, min.features = 200)
  20. AD5 <- Read10X(data.dir = "GSE157827/AD5/")
  21. AD5 <- CreateSeuratObject(counts = AD5, project = "AD5", min.cells = 3, min.features = 200)
  22. AD9 <- Read10X(data.dir = "GSE157827/AD9/")
  23. AD9 <- CreateSeuratObject(counts = AD9, project = "AD9", min.cells = 3, min.features = 200)
  24. AD10 <- Read10X(data.dir = "GSE157827/AD10/")
  25. AD10 <- CreateSeuratObject(counts = AD10, project = "AD10", min.cells = 3, min.features = 200)
  26. AD20 <- Read10X(data.dir = "GSE157827/AD20/")
  27. AD20 <- CreateSeuratObject(counts = AD20, project = "AD20", min.cells = 3, min.features = 200)
  28. AD21 <- Read10X(data.dir = "GSE157827/AD21/")
  29. AD21 <- CreateSeuratObject(counts = AD21, project = "AD21", min.cells = 3, min.features = 200)
  30. NC3 <- Read10X(data.dir = "GSE157827/NC3/")
  31. NC3 <- CreateSeuratObject(counts = NC3, project = "NC3", min.cells = 3, min.features = 200)
  32. NC7 <- Read10X(data.dir = "GSE157827/NC7/")
  33. NC7 <- CreateSeuratObject(counts = NC7, project = "NC7", min.cells = 3, min.features = 200)
  34. NC12 <- Read10X(data.dir = "GSE157827/NC12/")
  35. NC12 <- CreateSeuratObject(counts = NC12, project = "NC12", min.cells = 3, min.features = 200)
  36. NC14 <- Read10X(data.dir = "GSE157827/NC14/")
  37. NC14 <- CreateSeuratObject(counts = NC14, project = "NC14", min.cells = 3, min.features = 200)
  38. NC15 <- Read10X(data.dir = "GSE157827/NC15/")
  39. NC15 <- CreateSeuratObject(counts = NC15, project = "NC15", min.cells = 3, min.features = 200)
  40. NC17 <- Read10X(data.dir = "GSE157827/NC17/")
  41. NC17 <- CreateSeuratObject(counts = NC17, project = "NC17", min.cells = 3, min.features = 200)
  42. NC18 <- Read10X(data.dir = "GSE157827/NC18/")
  43. NC18 <- CreateSeuratObject(counts = NC18, project = "NC18", min.cells = 3, min.features = 200)
  44. #merge_data#
  45. merged_seurat <- merge(AD1, y = c(AD2,AD4,AD5,AD9,AD10,AD20,AD21,
  46. NC3,NC7,NC12,NC14,NC15,NC17,NC18),
  47. add.cell.ids = c("AD1","AD2","AD4","AD5","AD9","AD10","AD20","AD21",
  48. "NC3","NC7","NC12","NC14","NC15","NC17","NC18"))
  49. #count_UMI_MT#
  50. merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
  51. merged_seurat$mitoRatio <- [email hidden]$mitoRatio / 100
  52. #creat_metadata_include_base&sample&group_information#
  53. metadata <- [email hidden]
  54. metadata$cells <- rownames(metadata)
  55. metadata <- metadata %>%
  56. dplyr::rename(seq_folder = orig.ident,
  57. nUMI = nCount_RNA,
  58. nGene = nFeature_RNA)
  59. metadata$sample <- NA
  60. metadata$sample[which(str_detect(metadata$cells, "AD1"))] <- 'AD1'
  61. metadata$sample[which(str_detect(metadata$cells, "AD2"))] <- 'AD2'
  62. metadata$sample[which(str_detect(metadata$cells, "AD4"))] <- 'AD4'
  63. metadata$sample[which(str_detect(metadata$cells, "AD5"))] <- 'AD5'
  64. metadata$sample[which(str_detect(metadata$cells, "AD9"))] <- 'AD9'
  65. metadata$sample[which(str_detect(metadata$cells, "AD10"))] <- 'AD10'
  66. metadata$sample[which(str_detect(metadata$cells, "AD20"))] <- 'AD20'
  67. metadata$sample[which(str_detect(metadata$cells, "AD21"))] <- 'AD21'
  68. metadata$sample[which(str_detect(metadata$cells, "NC3"))] <- 'NC3'
  69. metadata$sample[which(str_detect(metadata$cells, "NC7"))] <- 'NC7'
  70. metadata$sample[which(str_detect(metadata$cells, "NC12"))] <- 'NC12'
  71. metadata$sample[which(str_detect(metadata$cells, "NC14"))] <- 'NC14'
  72. metadata$sample[which(str_detect(metadata$cells, "NC15"))] <- 'NC15'
  73. metadata$sample[which(str_detect(metadata$cells, "NC17"))] <- 'NC17'
  74. metadata$sample[which(str_detect(metadata$cells, "NC18"))] <- 'NC18'
  76. metadata <- metadata %>% mutate(disease = case_when(
  77. sample %in% c("AD1","AD2","AD4","AD5","AD9","AD10",
  78. "AD20","AD21") ~ 'AD',
  79. sample %in% c("NC3","NC7","NC12","NC14","NC15","NC17","NC18") ~ 'NC'))
  81. [email hidden] <- metadata
  82. ###QC###
  83. filtered_seurat <- subset(x = merged_seurat,
  84. subset= (nUMI < 20000) &
  85. (nGene > 200) &
  86. (mitoRatio < 0.20))
  87. filtered_seurat_QC <- JoinLayers(filtered_seurat)
  88. counts <- LayerData(object = filtered_seurat_QC, layer = "counts")
  89. nonzero <- counts > 0
  90. keep_genes <- Matrix::rowSums(nonzero) >= 10
  91. filtered_counts <- counts[keep_genes, ]
  92. matched_meta <- [email hidden][colnames(filtered_counts), ]
  93. filtered_seurat <- CreateSeuratObject(counts = filtered_counts,meta.data = matched_meta)
  94. #Filter_out_samples_that_contain_fewer_than_4,000 cells#
  95. #remove(AD6,AD8,AD13,AD19,NC11)#
  96. table(filtered_seurat$sample)
  97. metadata <- [email hidden]
  98. #QC_Visualize#
  99. metadata %>% ggplot(aes(x=sample, fill=sample)) + geom_bar() +
  100. theme_classic() +theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  101. theme(plot.title = element_text(hjust=0.5, face="bold")) +
  102. ggtitle("NCells")
  104. metadata %>% ggplot(aes(color=sample, x=nUMI, fill= sample)) +
  105. geom_density(alpha = 0.2) + scale_x_log10() + theme_classic() +
  106. ylab("Cell density") +geom_vline(xintercept = 500)
  107. #(because_shift_remove_NC16)⬇#
  108. metadata %>% ggplot(aes(color=sample, x=nGene, fill= sample)) + geom_density(alpha = 0.2) +
  109. theme_classic() +scale_x_log10() + geom_vline(xintercept = 300)
  111. metadata %>% ggplot(aes(x=sample, y=log10(nGene), fill=sample)) + geom_boxplot() +
  112. theme_classic() +theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  113. theme(plot.title = element_text(hjust=0.5, face="bold")) +ggtitle("NCells vs NGenes")
  115. metadata %>%
  116. ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) + geom_point() +
  117. scale_colour_gradient(low = "gray90", high = "black") +stat_smooth(method=lm) +
  118. scale_x_log10() + scale_y_log10() + theme_classic() +geom_vline(xintercept = 500) +
  119. geom_hline(yintercept = 250) +facet_wrap(~sample)
  120. #Memory_optimization##
  121. rm(AD1,AD2,AD4,AD5,AD9,AD10,AD20,AD21,NC3,NC7,NC12,NC14,NC15,NC17,NC18,
  122. filtered_seurat_QC,counts,filtered_counts,nonzero,merged_seurat)
  123. gc()
  124. #normalize data#
  125. scRNAlist_merge <- NormalizeData(filtered_seurat)
  126. scRNAlist_merge <- FindVariableFeatures(scRNAlist_merge,nfeatures = 2000)
  127. scRNAlist_merge <- ScaleData(scRNAlist_merge,vars.to.regress = c('mitoRatio'))
  128. #run_PCA#
  129. scRNAlist_merge <- RunPCA(scRNAlist_merge,npcs = 50)
  130. print(scRNAlist_merge[["pca"]], dims = 1:5, nfeatures = 5)
  131. #visualize_check_data#
  132. VizDimLoadings(scRNAlist_merge, dims = 1:2, reduction = "pca")
  133. DimPlot(scRNAlist_merge, reduction = "pca")
  134. DimHeatmap(scRNAlist_merge, dims = 1:20, cells = 500, balanced = TRUE)
  135. #remove_Batch_effect_by_harmony#
  136. library(harmony)
  137. scRNA_harmony <- RunHarmony(object = scRNAlist_merge,
  138. group.by.vars = "orig.ident",
  139. reduction = "pca",
  140. dims.use = 1:20,
  141. reduction.save = "harmony")
  142. #joinLayers#
  143. scRNA_harmony[['RNA']] <- JoinLayers(scRNA_harmony[['RNA']])
  144. # chose_0.4_resolution#
  145. scRNA_harmony <- FindNeighbors(scRNA_harmony,reduction = 'harmony',dims = 1:20)
  146. scRNA_harmony <- FindClusters(scRNA_harmony,resolution = seq(from = 0.1,to = 1.0, by = 0.1))
  147. scRNA_harmony <- RunUMAP(scRNA_harmony,dims = 1:20,reduction = 'harmony')
  148. scRNA_harmony <- RunTSNE(scRNA_harmony,dims = 1:20,reduction = 'harmony')
  149. scRNA_harmony$RNA_snn_res.0.4
  150. Idents(scRNA_harmony) <- 'RNA_snn_res.0.4'
  151. DimPlot(scRNA_harmony,reduction = 'umap')
  152. #UMAP_and_TSNE#
  153. umap_integrated_1 <- DimPlot(scRNA_harmony,reduction = 'umap',group.by = 'orig.ident',label = T)
  154. umap_integrated_2 <- DimPlot(scRNA_harmony,reduction = 'tsne', label = T)
  155. umap_integrated_3 <- DimPlot(scRNA_harmony,reduction = 'umap', label = T)
  156. umap_integrated_4 <- DimPlot(scRNA_harmony,reduction = 'umap',group.by = 'disease',label = T)
  157. umap_integrated_1
  158. umap_integrated_2
  159. umap_integrated_3
  160. umap_integrated_4
  161. # save data
  162. save(scRNA_harmony,file = 'scRNA_harmony.Rdata')
  167. ########Run analysis#########
  168. rm(list = ls())
  169. gc()
  170. #load packages#
  171. library(Seurat)
  172. library(multtest)
  173. library(dplyr)
  174. library(ggplot2)
  175. library(patchwork)
  176. library(SeuratData)
  177. library(dplyr)
  179. ###load_data###
  180. load('scRNA_harmony.Rdata')
  181. scRNA_harmony$RNA_snn_res.0.4
  182. Idents(scRNA_harmony) <- 'RNA_snn_res.0.4'
  183. DimPlot(scRNA_harmony,reduction = 'umap')
  184. ###Annotation_was_based_on_marker_genes_reported_in_the_original_publication###
  185. current_levels <- levels(scRNA_harmony)
  186. sorted_levels <- as.character(sort(as.numeric(current_levels)))
  187. Idents(scRNA_harmony) <- factor(Idents(scRNA_harmony), levels = sorted_levels)
  188. levels(scRNA_harmony)
  189. new.cluster.ids <- c("oligodendrocyte_1", # "MBP"0
  190. "astrocytes_1", # "AQP4"1
  191. "excitatory_neuron_1", # "CAMK2A"2
  192. "excitatory_neuron_2", # "CAMK2A"3
  193. "inhibitory_neuron_1", # "GAD1"
  194. "inhibitory_neuron_2", # "GAD1"8
  195. "excitatory_neuron_3", # "CAMK2A"6
  196. "microglia_1", # "C3" 7
  197. "inhibitory_neuron_3", # "GAD1"8
  198. "excitatory_neuron_4", # "CAMK2A"9
  199. "excitatory_neuron_5", # "CAMK2A"10
  200. "inhibitory_neuron_4", # "GAD1"11
  201. "excitatory_neuron_6", # "CAMK2A"12
  202. "inhibitory_neuron_5", # "GAD1"
  203. "excitatory_neuron_7", # "CAMK2A"
  204. "excitatory_neuron_8", # "CAMK2A"
  205. "excitatory_neuron_9", # "CAMK2A"
  206. "endothelial_cell_1", # "CLDN5"
  207. "astrocytes_2", # "AQP4"
  208. "excitatory_neuron_10", # "CAMK2A"
  209. "excitatory_neuron_11", # "CAMK2A"
  210. "oligodendrocyte_2", # "MBP"
  211. "excitatory_neuron_12", # "CAMK2A"
  212. "excitatory_neuron_13", # "CAMK2A"
  213. "inhibitory_neuron_6", # "GAD1"
  214. "oligodendrocyte_3", # "MBP"
  215. "astrocytes_3") # "AQP4"
  216. names(new.cluster.ids) <- levels(scRNA_harmony)
  217. scRNA_harmony <- RenameIdents(scRNA_harmony, new.cluster.ids)
  218. scRNA_harmony<- AddMetaData(object = scRNA_harmony,
  219. metadata = [email hidden],
  220. col.name = "celltype")
  221. #visualize_Annotation_marker_genes#
  222. cell_mark <- c("ADGRV1","GPC5","RYR3", "AQP4", # astrocytes
  223. "ABCB1","EBF1", "CLDN5", # endothelial_cell
  224. "CBLN2","LDB2", "CAMK2A", # excitatory_neuron
  225. "LHFPL3","PCDH15","GAD1", # inhibitory_neuron
  226. "LRMDA","DOCK8","C3", # microglia
  227. "PLP1","ST18","MBP") # oligodendrocyte
  228. DotPlot(scRNA_harmony,
  229. features = cell_mark,
  230. assay = "RNA",
  231. cluster.idents = TRUE,
  232. scale.by = "size",
  233. scale = TRUE,
  234. col.min = -2,
  235. col.max = 2) +
  236. coord_flip() +theme_bw() + labs(x = "Genes", y = "Cell Types") +
  237. theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1, size=8),
  238. axis.text.y = element_text(size=8)) +
  239. scale_color_gradient2(low = "#2166AC",mid = "white",high = "#B2182B")
  241. ###AUC###
  242. seurat_integrated<-scRNA_harmony
  243. DefaultAssay(seurat_integrated) <- "RNA"
  244. library(readxl)
  245. library("Matrix")
  246. library("AUCell")
  247. terms<- read_excel("list.xlsx")
  248. print(terms)
  249. geneset<- terms[['comorbidity_associated_upregulated_genes']]
  250. geneset <- geneset[!is.na(geneset)]
  251. valid_genes <- intersect(geneset, rownames(seurat_integrated))
  252. if (length(valid_genes) < 5) stop("Too few valid genes).")
  253. gene_sets <- list("comorbidity_associated_upregulated_genes" = valid_genes)
  254. seurat_integrated <- JoinLayers(seurat_integrated)
  255. expr_matrix <- LayerData(seurat_integrated, layer = "data")
  256. expr_matrix <- as(expr_matrix, "dgCMatrix")
  257. class(expr_matrix)
  258. cells_rankings <- AUCell_buildRankings(expr_matrix,plotStats = FALSE, splitByBlocks = TRUE)
  259. cells_AUC <- AUCell_calcAUC(gene_sets,cells_rankings,
  260. aucMaxRank = ceiling(0.05 * nrow(cells_rankings)),nCores = 4)
  261. seurat_integrated$upregulated <- as.numeric(getAUC(cells_AUC)["comorbidity_associated_upregulated_genes", ])
  262. ##visualize_AUC_umap##
  263. umap_df <- FetchData(seurat_integrated, vars = c("umap_1", "umap_2", "upregulated"))
  264. ggplot(umap_df, aes(x = umap_1, y = umap_2, color = upregulated)) +
  265. geom_point(size = 0.8, alpha = 0.8) +
  266. scale_color_viridis_c(option = "viridis",name = "AUC Score" ) +
  267. labs(x = "UMAP1", y = "UMAP2",
  268. title = "" ) +
  269. theme_classic()
  270. ##visualize_AUC_Bubble_Plot##
  271. auc_df <- [email hidden] %>%group_by(celltype) %>%
  272. summarise(auc_mean = mean(upregulated),n_cells = n())
  273. ggplot(auc_df, aes(x = celltype, y = 1,size = n_cells,color = auc_mean)) +
  274. geom_point(alpha = 0.9) +
  275. scale_size(range = c(3, 15)) +
  276. scale_color_gradientn(colors = c("#2166AC", "white", "#B2182B")) +
  277. theme_bw() +
  278. theme(axis.title.y = element_blank(),axis.text.y = element_blank(),
  279. axis.ticks.y = element_blank(),axis.text.x = element_text(angle =90, hjust = 1, size = 10)) +
  280. labs(x = "Cell Type",color = "AUC Score",size = "Number of Cells",
  281. title = "AUC Activity Across Cell Types")
  282. ##visualize_AUC_BOX_plot##
  283. #chose_cell_type#
  284. target_cells <- c('astrocytes_1',"microglia_1","oligodendrocyte_1","oligodendrocyte_3",
  285. "excitatory_neuron_5","endothelial_cell_1")
  286. seurat_subset <- subset(seurat_integrated, subset = celltype %in% target_cells)
  287. plot_data <- FetchData(seurat_subset, vars = c("upregulated","disease","celltype"))
  288. #disease_group_in_order#
  289. plot_data$disease <- factor(plot_data$disease, levels = c("AD", "NC"))
  290. #cell_type_in_order#
  291. plot_data$cell_type <- factor(plot_data$celltype,
  292. levels = c('astrocytes_1',"microglia_1","oligodendrocyte_1",
  293. "oligodendrocyte_3","excitatory_neuron_5","endothelial_cell_1"))
  294. #visualize#
  295. ggplot(plot_data, aes(x = cell_type,y = upregulated,fill = disease)) +
  296. geom_boxplot(position = position_dodge(0.8),width = 0.7,outlier.size = 0.5)+
  297. scale_fill_manual(values = c("NC" = "#1f77b4", "AD" = "#ff7f0e"),name = " ")+
  298. labs(x = "Cell Type",y = "AUC Score",title = " ",subtitle = " ") +
  299. theme_classic() +
  300. theme(legend.position = "top",legend.title = element_text(size = 16, face = "bold"),
  301. legend.text = element_text(size = 16), axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
  302. axis.text.y = element_text(size = 16), axis.title.x = element_text(size = 16, face = "bold"),
  303. axis.title.y = element_text(size = 16, face = "bold"),plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
  304. plot.subtitle = element_text(size = 16, hjust = 0.5))
  306. ###visualize_proportions_Bar_plot###
  307. library(forcats)
  308. library(scales)
  309. metadata <- [email hidden]
  310. #chose_cell_type#
  311. selected_celltypes <- c('astrocytes_1',"microglia_1","oligodendrocyte_1", "oligodendrocyte_3",
  312. "excitatory_neuron_5","endothelial_cell_1")
  313. #Calculate_cell-type_proportions#
  314. total_per_disease <- metadata %>%group_by(disease) %>%summarise(total_cells = n(), .groups = "drop")
  315. celltype_percent <- metadata %>%filter(celltype %in% selected_celltypes) %>%
  316. group_by(disease, celltype) %>%summarise(n = n(), .groups = "drop") %>%
  317. left_join(total_per_disease, by = "disease") %>%mutate(percent = 100 * n / total_cells)
  318. #cell_type_in_order#
  319. celltype_percent <- celltype_percent %>%mutate(celltype = factor(celltype, levels = selected_celltypes))
  320. #disease_group_in_order#
  321. celltype_percent <- celltype_percent %>%mutate(disease = factor(disease, levels = c("NC", "AD")))
  322. #visualize#
  323. p <- ggplot(celltype_percent, aes(x = celltype, y = percent, fill = disease)) +
  324. geom_bar(stat = "identity",position = position_dodge(width = 0.8),color = "black",width = 0.7)+
  325. geom_text(aes(label = sprintf("%.1f", percent)),position = position_dodge(width = 0.8),
  326. vjust = -0.4, size = 4.5, color = "black")+
  327. scale_fill_manual(name = " ",values = c("NC" = "#1f77b4", "AD" ="#ff7f0e")) +
  328. labs(x = "",y = "Proportion of cells (%)",title = " "
  329. ) +theme_classic() +
  330. theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
  331. axis.text.y = element_text(size = 12),
  332. axis.title = element_text(size = 12, face = "bold"),
  333. legend.position = "top",legend.text = element_text(size = 14),
  334. legend.title = element_text(size = 15),legend.key.size = unit(0.9, "cm"),
  335. panel.grid.major.y = element_line(color = "gray90", linewidth = 0.2))+
  336. scale_y_continuous(limits = c(0, max(celltype_percent$percent, na.rm = TRUE) * 1.15),
  337. expand = expansion(mult = c(0, 0.05)))
  338. p
  340. ###donor-aware differential abundance testing###
  341. #chose_cell_type#
  342. target_celltypes <- c("oligodendrocyte_1", "astrocytes_1","excitatory_neuron_1","excitatory_neuron_2",
  343. "inhibitory_neuron_1", "inhibitory_neuron_2","excitatory_neuron_3",
  344. "microglia_1","inhibitory_neuron_3","excitatory_neuron_4", "excitatory_neuron_5",
  345. "inhibitory_neuron_4", "excitatory_neuron_6", "inhibitory_neuron_5",
  346. "excitatory_neuron_7", "excitatory_neuron_8", "excitatory_neuron_9","endothelial_cell_1",
  347. "astrocytes_2", "excitatory_neuron_10", "excitatory_neuron_11", "oligodendrocyte_2",
  348. "excitatory_neuron_12", "excitatory_neuron_13","inhibitory_neuron_6",
  349. "oligodendrocyte_3","astrocytes_3")
  350. #check_data#
  351. md <- [email hidden]
  352. donor_group <- tapply(md$disease, md$sample, function(x) unique(x)[1])
  353. #run#
  354. one_ct <- function(ct) { n_ct<- tapply(md$celltype == ct, md$sample, sum)
  355. n_total <- tapply(rep(1, nrow(md)), md$sample, sum)
  356. prop_ct <- n_ct / n_total
  357. g <- donor_group[names(prop_ct)]
  358. grp_levels <- unique(md$disease)
  359. if (length(grp_levels) < 2) return(NULL)
  360. g1 <- prop_ct[g == grp_levels[1]]
  361. g2 <- prop_ct[g == grp_levels[2]]
  362. if (sum(!is.na(g1)) < 2 || sum(!is.na(g2)) < 2) {p <- NA_real_}
  363. else {p <- suppressWarnings(wilcox.test(g1, g2, exact=FALSE)$p.value)}
  364. data.frame(celltype = ct,n_donors_group1 = sum(!is.na(g1)),
  365. n_donors_group2 = sum(!is.na(g2)),median_group1=median(g1, na.rm = TRUE),
  366. median_group2= median(g2, na.rm = TRUE),p.value = p,
  367. stringsAsFactors = FALSE)}
  368. res_list <- lapply(target_celltypes, one_ct)
  369. res_DA <- do.call(rbind, res_list)
  370. res_DA$FDR <- p.adjust(res_DA$p.value, method = "BH")
  371. write.csv(res_DA, "DA_donor_wilcoxon_SELECTED_celltypes.csv", row.names = FALSE)
  372. print(res_DA)
  374. ###CellChat ####
  375. library(CellChat)
  376. library(ggalluvial)
  377. library(patchwork)
  378. library(RColorBrewer)
  379. options(stringsAsFactors = FALSE)
  380. #Prepare_the_data_for_CellChat#
  381. seurat_integrated[["RNA"]] <- JoinLayers(seurat_integrated[["RNA"]])
  382. data.input <- LayerData(seurat_integrated, assay = "RNA", layer = "data")
  383. metadata <- data.frame(cell_type = Idents(seurat_integrated),
  384. disease = seurat_integrated$disease,
  385. row.names = colnames(seurat_integrated))
  386. print(unique(metadata$disease))
  387. #Cell–cell_communication_is_shown_for_AD_only#
  388. selected_diseases <- c('AD')
  389. disease_cells <- rownames(metadata)[metadata$disease %in% selected_diseases]
  390. #chose_cell_type#
  391. selected_celltypes <- c('astrocytes_1',"microglia_1","oligodendrocyte_1", "oligodendrocyte_3",
  392. "excitatory_neuron_5","endothelial_cell_1")
  394. celltype_cells <- rownames(metadata)[metadata$cell_type %in% selected_celltypes]
  395. selected_cells <- intersect(disease_cells, celltype_cells)
  396. data.input.subset <- data.input[, selected_cells]
  397. metadata.subset <- metadata[selected_cells, , drop = FALSE]
  398. cellchat <- createCellChat(object = data.input.subset)
  399. cell_type_df <- data.frame(cell_type = metadata.subset$cell_type,
  400. row.names = rownames(metadata.subset))
  401. full_meta_df <- data.frame(cell_type = metadata.subset$cell_type,
  402. disease = metadata.subset$disease,
  403. row.names = rownames(metadata.subset))
  404. cellchat <- addMeta(cellchat, meta = full_meta_df)
  405. cellchat <- setIdent(cellchat, ident.use = "cell_type")
  406. cellchat@idents <- factor(as.character(cellchat@meta$cell_type),
  407. levels = selected_celltypes)
  408. print(levels(cellchat@idents))
  409. groupSize <- as.numeric(table(cellchat@idents))
  410. print(groupSize)
  411. print(table(cellchat@meta$disease))
  412. #Set_database#
  413. CellChatDB <- CellChatDB.human
  414. cellchat@DB <- CellChatDB
  415. cellchat <- subsetData(cellchat)
  416. #run_analysis#
  417. cellchat <- identifyOverExpressedGenes(cellchat)
  418. cellchat <- identifyOverExpressedInteractions(cellchat)
  419. cellchat <- projectData(cellchat, PPI.human)
  420. cellchat <- computeCommunProb(cellchat)
  421. cellchat <- computeCommunProbPathway(cellchat)
  422. cellchat <- aggregateNet(cellchat)
  423. #quantitative_centrality_metrics#
  424. library(igraph)
  425. W <- cellchat@net$weight
  426. if (is.null(W) || all(W == 0)) stop("CellChat no network")
  427. g <- graph_from_adjacency_matrix(W, mode = "directed", weighted = TRUE, diag = FALSE)
  428. outgoing <- strength(g, mode = "out", weights = E(g)$weight)
  429. incoming <- strength(g, mode = "in", weights = E(g)$weight)
  430. bet <- betweenness(g, directed = TRUE, weights = 1/E(g)$weight)
  431. clo <- closeness(g, mode = "all", weights = 1/E(g)$weight)
  432. centrality <- data.frame(celltype = rownames(W),
  433. outgoing = outgoing,
  434. incoming = incoming,
  435. betweenness = bet,
  436. closeness = clo)
  437. write.csv(centrality, "CellChat_global_centrality.csv", row.names = FALSE)
  438. centrality
  440. #visualize#
  441. cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
  442. print("pathways")
  443. pathways <- cellchat@netP$pathways
  444. print(pathways)
  445. #visualize#
  446. cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
  447. colors <- RColorBrewer::brewer.pal(length(selected_celltypes), "Set1")
  448. names(colors) <- selected_celltypes
  449. #visualize_communication_weight _matrix_plot#
  450. par(mar = c(0, 0, 2.5, 0)) #
  451. netVisual_circle(cellchat@net$weight,
  452. color.use=colors,
  453. vertex.weight=groupSize,
  454. weight.scale=TRUE,
  455. edge.weight.max=max(cellchat@net$weight),
  456. title.name="Communication Weight Matrix",
  457. vertex.label.cex = 1.4)
  458. #visualize_communication_weight _heatmap_plot#
  459. netVisual_heatmap(cellchat,
  460. signaling = NULL,
  461. measure = "weight",
  462. color.heatmap = "Reds",
  463. title.name = "Communication Weight Heatmap",
  464. font.size = 14,
  465. font.size.title = 16)
  466. #Prioritization_of_sender_and_receiver_signals_of_interest#
  467. sender_cell <- "astrocytes_1"
  468. receiver_cell <- "oligodendrocyte_3"
  469. pathway_comm <- subsetCommunication(cellchat, slot.name = "netP",
  470. sources.use = sender_cell,
  471. targets.use = receiver_cell)
  472. pathway_comm_sorted <- pathway_comm[order(-pathway_comm$prob), ]
  473. print("Top 10 signaling pathways:")
  474. print(head(pathway_comm_sorted, 10))
  475. #Visualize_specific_pathways#
  476. colors <- brewer.pal(length(selected_celltypes), "Dark2")
  477. names(colors) <- selected_celltypes
  478. pathways.show<-'NRXN'
  479. netVisual_aggregate(cellchat,
  480. signaling = pathways.show,
  481. layout = "chord",
  482. vertex.size = groupSize,
  483. color.use = colors,
  484. title.name = paste(pathway, "Signaling"),
  485. arrow.size = 0.02,
  486. edge.width.max = 15)
  488. ###overlap_between_SUGS_and_ligand–receptor###
  489. caup_genes <- unique(valid_genes)
  490. caup_genes <- caup_genes[!is.na(caup_genes)]
  491. UP <- function(x) unique(toupper(na.omit(x)))
  492. split_genes <- function(x) {if (is.null(x)) return(character(0))
  493. x <- na.omit(x)
  494. unlist(strsplit(x, "\\s*[+|_,/; ]\\s*"))}
  495. #check_overlap_in_Global#
  496. comm_all <- subsetCommunication(cellchat, slot.name = "netP")
  497. lr_all <- unique(c(split_genes(comm_all$ligand), split_genes(comm_all$receptor)))
  498. overlap_global <- intersect(UP(caup_genes), UP(lr_all))
  499. cat("Global overlap size =", length(overlap_global), "\n")
  500. print(sort(overlap_global))
  501. #check_overlap_in_one_sender#
  502. sender_only <- "astrocytes_1"
  503. comm_sender <- subsetCommunication(cellchat, slot.name = "netP", sources.use = sender_only)
  504. lr_sender <- unique(c(split_genes(comm_sender$ligand), split_genes(comm_sender$receptor)))
  505. overlap_sender <- intersect(UP(caup_genes), UP(lr_sender))
  506. cat(sprintf("Overlap (sender=%s) size = %d\n", sender_only, length(overlap_sender)))
  507. print(sort(overlap_sender))
  508. #check_overlap_in_one_receiver#
  509. receiver_only <- "oligodendrocyte_1"
  510. comm_receiver <- subsetCommunication(cellchat, slot.name = "netP", targets.use = receiver_only)
  511. lr_receiver <- unique(c(split_genes(comm_receiver$ligand), split_genes(comm_receiver$receptor)))
  512. overlap_receiver <- intersect(UP(caup_genes), UP(lr_receiver))
  513. cat(sprintf("Overlap (receiver=%s) size = %d\n", receiver_only, length(overlap_receiver)))
  514. print(sort(overlap_receiver))
  515. #check_overlap_in_one_receiver#
  516. receiver_only <- "oligodendrocyte_3"
  517. comm_receiver <- subsetCommunication(cellchat, slot.name = "netP", targets.use = receiver_only)
  518. lr_receiver <- unique(c(split_genes(comm_receiver$ligand), split_genes(comm_receiver$receptor)))
  519. overlap_receiver <- intersect(UP(caup_genes), UP(lr_receiver))
  520. cat(sprintf("Overlap (receiver=%s) size = %d\n", receiver_only, length(overlap_receiver)))
  521. print(sort(overlap_receiver))
  523. ###Find_top10_Markers_for_every_cell-subtype###
  524. pbmc.markers <- FindAllMarkers(scRNA_harmony, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
  525. topmaker<-pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
  526. write.csv(topmaker, file = "topmaker.csv", row.names = FALSE)
  528. ###pseudobulk###
  529. target_celltypes <-c("oligodendrocyte_1", "astrocytes_1","excitatory_neuron_1","excitatory_neuron_2",
  530. "inhibitory_neuron_1", "inhibitory_neuron_2","excitatory_neuron_3",
  531. "microglia_1","inhibitory_neuron_3","excitatory_neuron_4", "excitatory_neuron_5",
  532. "inhibitory_neuron_4", "excitatory_neuron_6", "inhibitory_neuron_5",
  533. "excitatory_neuron_7", "excitatory_neuron_8", "excitatory_neuron_9","endothelial_cell_1",
  534. "astrocytes_2", "excitatory_neuron_10", "excitatory_neuron_11", "oligodendrocyte_2",
  535. "excitatory_neuron_12", "excitatory_neuron_13","inhibitory_neuron_6",
  536. "oligodendrocyte_3","astrocytes_3")
  538. suppressWarnings(suppressMessages({
  539. edgeR_ok <- requireNamespace("edgeR", quietly = TRUE)
  540. library(Matrix)}))
  541. .assay <- DefaultAssay(scRNA_harmony)
  542. get_counts <- function(obj, assay) {
  543. lyr_names <- tryCatch(SeuratObject::Layers(obj[[assay]]), error = function(e) character(0))
  544. if ("counts" %in% lyr_names) {
  545. return(SeuratObject::GetAssayData(obj, assay = assay, layer = "counts"))}
  546. if ("data" %in% lyr_names) {
  547. m <- SeuratObject::GetAssayData(obj, assay = assay, layer = "data")
  548. m <- as(m, "dgCMatrix")
  549. return(m)}
  550. m <- tryCatch(SeuratObject::GetAssayData(obj, assay = assay), error = function(e) NULL)
  551. if (!is.null(m)) return(as(m, "dgCMatrix"))
  552. stop("Seurat v5 no counts/data")}
  553. cm <- get_counts(scRNA_harmony, .assay)
  554. md <- [email hidden]
  555. out_dir <- "PB_validation_SELECTED"
  556. dir.create(out_dir, showWarnings = FALSE)
  557. run_one_ct <- function(ct) {
  558. cells_ct <- rownames(md)[md$celltype == ct]
  559. if (length(cells_ct) == 0) { message("jump over：", ct, "nocells"); return(invisible(NULL)) }
  560. donors <- unique(md$sample[md$celltype == ct])
  561. grp_by_donor <- tapply(md$disease, md$sample, function(x) unique(x)[1])[donors]
  562. grp_levels <- unique(md$disease)
  563. if (length(grp_levels) < 2) { message("jump over：", ct, "groupproblem"); return(invisible(NULL)) }
  564. donors_g1 <- donors[grp_by_donor == grp_levels[1]]
  565. donors_g2 <- donors[grp_by_donor == grp_levels[2]]
  566. if (length(donors_g1) < 2 || length(donors_g2) < 2) {
  567. message("jump over：", ct, "sample<2 per group"); return(invisible(NULL))}
  568. pb_cols <- lapply(donors, function(dn) {
  569. cs <- rownames(md)[md$celltype == ct & md$sample == dn]
  570. if (length(cs) == 0) return(Matrix(0, nrow = nrow(cm), ncol = 1, sparse = TRUE))
  571. Matrix::rowSums(cm[, cs, drop = FALSE])
  572. })
  573. pb <- do.call(cbind, pb_cols)
  574. colnames(pb) <- donors
  575. if (edgeR_ok) {y <- edgeR::DGEList(counts = pb,
  576. samples = data.frame(sample = donors,
  577. group = grp_by_donor[donors],
  578. row.names = donors))
  579. keep <- rowSums(y$counts) > 0
  580. if (sum(keep) < 10) { message("jump over", ct, "genes < 10"); return(invisible(NULL)) }
  581. y <- y[keep, , keep.lib.sizes = FALSE]
  582. y <- edgeR::calcNormFactors(y)
  583. design <- model.matrix(~ y$samples$group)
  584. y <- edgeR::estimateDisp(y, design)
  585. fit <- edgeR::glmQLFit(y, design)
  586. qlf <- edgeR::glmQLFTest(fit, coef = 2)
  587. tt <- edgeR::topTags(qlf, n = Inf)$table
  588. tt$gene <- rownames(tt)
  589. out <- file.path(out_dir, paste0("PB_edgeR_DEGs_", make.names(ct), ".csv"))
  590. write.csv(tt, out, row.names = FALSE)
  591. message("[edgeR]：", out)}
  592. invisible(TRUE)}
  593. invisible(lapply(target_celltypes, run_one_ct))
  595. ###BBB_related_pathway_check###
  596. library(readxl)
  597. library(dplyr)
  598. library(Matrix)
  599. library(AUCell)
  600. library(ggplot2)
  601. library(tidyr)
  602. seurat_integrated <- scRNA_harmony
  603. terms <- read_excel("Blood–Brain Barrier (BBB)-related signaling pathways.xlsx")
  604. seurat_integrated <- JoinLayers(seurat_integrated)
  605. expr_matrix <- LayerData(seurat_integrated, layer = "data")
  606. expr_matrix <- as(expr_matrix, "dgCMatrix")
  607. cells_rankings <- AUCell_buildRankings(
  608. expr_matrix,
  609. plotStats = FALSE,
  610. splitByBlocks = TRUE
  611. )
  612. bubble_list <- list()
  613. for (col_name in colnames(terms)) {
  614. geneset <- terms[[col_name]]
  615. geneset <- geneset[!is.na(geneset)]
  616. valid_genes <- intersect(geneset, rownames(expr_matrix))
  617. if (length(valid_genes) < 5) {
  618. cat("Skipping:", col_name, "(too few genes)\n")
  619. next
  620. }
  621. gene_sets <- list(current_set = valid_genes)
  622. cells_AUC <- AUCell_calcAUC(
  623. gene_sets,
  624. cells_rankings,
  625. aucMaxRank = ceiling(0.05 * nrow(cells_rankings))
  626. )
  627. auc_scores <- as.numeric(getAUC(cells_AUC)["current_set", ])
  628. auc_col_name <- paste0("AUC_", col_name)
  629. [email hidden][[auc_col_name]] <- auc_scores
  630. tmp <- [email hidden] %>%
  631. group_by(celltype) %>%
  632. summarise(
  633. auc_mean = mean(.data[[auc_col_name]], na.rm = TRUE),
  634. n_cells = n(),
  635. .groups = "drop"
  636. ) %>%
  637. mutate(pathway = col_name)
  638. bubble_list[[col_name]] <- tmp
  639. }
  640. bubble_df <- bind_rows(bubble_list)
  641. p <- ggplot(bubble_df,aes(x = celltype, y = pathway,
  642. size = n_cells, color = auc_mean)) +
  643. geom_point(alpha = 0.9) +
  644. scale_size(range = c(2, 12)) +
  645. scale_color_gradientn(colors = c("#2166AC", "white", "#B2182B")) +
  646. theme_bw() +
  647. theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
  648. axis.text.y = element_text(size = 10)) +
  649. labs(x = "Cell Types",y = "",
  650. color = "AUC Score",
  651. size = "Number of Cells",
  652. title = "")
  653. print(p)
  654. #AUC_box#
  655. plot_AUC_box <- function(
  656. seurat_obj,
  657. pathways,
  658. celltypes,
  659. disease_col = "disease",
  660. celltype_col = "celltype"
  661. ){
  662. library(ggplot2)
  663. library(dplyr)
  664. df <- [email hidden] %>%
  665. filter(.data[[celltype_col]] %in% celltypes)
  666. auc_cols <- paste0("AUC_", pathways)
  667. plot_df <- df %>%
  668. select(all_of(c(celltype_col, disease_col, auc_cols))) %>%
  669. pivot_longer(
  670. cols = all_of(auc_cols),
  671. names_to = "pathway",
  672. values_to = "AUC"
  673. ) %>%
  674. mutate(
  675. pathway = gsub("AUC_", "", pathway),
  676. celltype = factor(.data[[celltype_col]], levels = celltypes),
  677. disease = factor(.data[[disease_col]], levels = c("AD","NC"))
  678. )
  679. p <- ggplot(plot_df, aes(x = celltype, y = AUC, fill = disease)) +
  680. geom_boxplot(position = position_dodge(0.8), width = 0.7, outlier.size = 0.5) +
  681. facet_wrap(~ pathway, scales = "free_y") +
  682. scale_fill_manual(values = c("NC" = "#1f77b4", "AD" = "#ff7f0e"),name = " ") +
  683. theme_classic() +
  684. labs(
  685. x = "Cell Type",
  686. y = "AUC Score",
  687. title = " "
  688. ) +
  689. theme(
  690. legend.position = "top",
  691. legend.text = element_text(size = 14),
  692. axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
  693. axis.title.x = element_text(size = 14, face = "bold"),
  694. axis.text.y = element_text(size = 12),
  695. axis.title.y = element_text(size = 14, face = "bold"),
  696. strip.text = element_text(size = 12, face = "bold")
  697. )
  698. print(p)}
  699. plot_AUC_box(
  700. seurat_obj = seurat_integrated,
  701. pathways = c("Tight junction", "Adherens junction", "Leukocyte transendothelial migration",
  702. "VEGF signaling pathway", "Wnt signaling pathway", "ECM-receptor interaction"),
  703. celltypes = c('astrocytes_1',"microglia_1","oligodendrocyte_1", "oligodendrocyte_3",
  704. "excitatory_neuron_5","endothelial_cell_1")
  705. )

singlecell codeAD.R at commit c3be27e, no license · at the source

Overview

Authors: Zhangzheng Wang1, Krisztián Juhász Zoltán1, Csaba Matta1,2,3, Luca Paluska1, Ahmed Al-Mnaseer1, Roland Takács1, László Ducza1
ORCID iDs: Zhangzheng Wang, Krisztián Juhász Zoltán, Csaba Matta, Ahmed Al-Mnaseer, Roland Takács, László Ducza
  1. Department of Anatomy, Histology and Embryology, Faculty of Medicine, University of Debrecen, H-4032 Debrecen, Hungary
  2. Department of Biomedical Materials Science, Graduate School of JABA, Wonkwang University, Iksan, Republic of Korea
  3. Integrated Omics Institute, Wonkwang University, Iksan, Republic of Korea
Institutions: University of Debrecen (Hungary); Wonkwang University (South Korea)
Journal: Computational and structural biotechnology journal, volume 35, issue 1, article 0085
Dates: received 5 December 2025; accepted 10 April 2026; published online 13 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.34133/csbj.0085 · PMID 42137879 · PMCID PMC13168759 · OpenAlex W7153734934
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, Machine learning
Topic: Bone Metabolism and Diseases (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Nemzeti Kutatási Fejlesztési és Innovációs Hivatal (grant no. 2025-1.2.3-TÉT-IPARI-KR-2025-00002, 152166); Nemzeti Kutatási, Fejlesztési és Innovaciós Alap (EKÖP-25-4-II-DE-249); Általános Orvostudományi Kar, Debreceni Egyetem (8XEABK00PUBT/320)
Citations: not cited yet (Europe PMC); 69 references in the paper

Abstract

Osteoarthritis (OA) and Alzheimer’s disease (AD) are prevalent age-associated disorders that frequently co-occur, yet the molecular basis of their comorbidity remains incompletely understood. To explore potential shared cellular programs, we performed an integrative analysis of publicly available bulk and single-cell transcriptomic datasets derived from human OA cartilage and AD cortex. Cross-disease comparison identified 60 overlapping differentially expressed genes, including 18 consistently up-regulated genes, which we defined as a shared up-regulated gene set (SUGS). Functional enrichment analyses indicated convergence on extracellular matrix remodeling, inflammatory signaling, metabolic stress responses, and immune regulation. Single-cell analysis of OA cartilage revealed expansion of a fibrochondrocyte subpopulation enriched for SUGS activity and extracellular-matrix-associated ligands. In AD cortex, a disease-associated oligodendrocyte subcluster displayed elevated SUGS activity and stress-response gene expression and occupied a prominent-receiver-like position within inferred neuronal–glial communication networks. Ligand–receptor analysis performed independently within each tissue identified collagen-related signaling in OA and neurexin-associated signaling in AD as dominant intratissue pathways. Because the OA and AD datasets are cross-sectional and were derived from independent cohorts and distinct tissues, these analyses do not establish direct inter-organ communication, temporal sequence, or causal directionality. In addition, CellChat-based ligand–receptor inference was performed independently within each tissue and does not itself infer cross-organ communication. The identification of sender-like and receiver-like cellular hubs should therefore be interpreted as a hypothesis-generating, systems-level conceptual framework that may help organize future experimental studies investigating potential links between joint inflammation and neurodegeneration in aging.

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

Repository

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

znGer-cel/Sender-and-Receiver-Cellular-Hubs-as-Potential-Therapeutic-Targets-in-OA-AD-Comorbidity

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: c3be27eff489e84ed63238cd68276760e43cf07b, 3 April 2026
Languages: R (3)
Size: 5 files, 3 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: edgeR (3 files), ggplot2 (3 files), tidyverse (3 files), Harmony (2 files), igraph (2 files), patchwork (2 files), Seurat (2 files), data.table (1 file), DESeq2 (1 file), limma (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
4 files
At the source: github.com/znGer-cel/Sender-and-Receiver-Cellular-Hubs-as-Potential-Therapeutic-Targets-in-OA-AD-Comorbidity

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 3 scripts, each with its path and the digest of its content;
  • 19 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 Availability

All raw and processed datasets utilized in this study are publicly accessible through the Gene Expression Omnibus under the following accession numbers: OA bulk RNA-seq (GSE114007 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114007)), AD bulk microarray (GSE122063 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE122063)), OA scRNA-seq (GSE255460 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255460)), and AD scRNA-seq (GSE157827 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE157827)). Processed data matrices and related analysis scripts will be made available from the corresponding author upon reasonable request. The full analysis pipeline and R scripts are available on GitHub (https://github.com/znGer-cel/Sender-and-Receiver-Cellular-Hubs-as-Potential-Therapeutic-Targets-in-OA-AD-Comorbidity).

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 2, 28 September 2026

  • Funding: added Nemzeti Kutatási Fejlesztési és Innovációs Hivatal: grant no. 2025-1.2.3-TÉT-IPARI-KR-2025-00002, 152166; Nemzeti Kutatási, Fejlesztési és Innovaciós Alap: EKÖP-25-4-II-DE-249; Általános Orvostudományi Kar, Debreceni Egyetem: 8XEABK00PUBT/320

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 69 references.

Cite

This paper

Wang, Z., Zoltán, K. J., Matta, C., Paluska, L., Al-Mnaseer, A., Takács, R., & Ducza, L. (2026). A Systems-Level Transcriptomic Framework Identifies Shared Cellular Hubs in Osteoarthritis and Alzheimer's Disease. Computational and structural biotechnology journal, 35(1), 0085. https://doi.org/10.34133/csbj.0085

BibTeX

@article{wang2026systems,
author = {Wang, Zhangzheng and Zoltán, Krisztián Juhász and Matta, Csaba and Paluska, Luca and Al-Mnaseer, Ahmed and Takács, Roland and Ducza, László},
title = {{A Systems-Level Transcriptomic Framework Identifies Shared Cellular Hubs in Osteoarthritis and Alzheimer's Disease}},
journal = {Computational and structural biotechnology journal},
year = {2026},
month = may,
volume = {35},
number = {1},
pages = {0085},
publisher = {AAAS Science Partner Journal Program},
issn = {2001-0370},
doi = {10.34133/csbj.0085},
url = {https://doi.org/10.34133/csbj.0085},
pmid = {42137879},
pmcid = {PMC13168759}
}

RIS

TY - JOUR
AU - Wang, Zhangzheng
AU - Zoltán, Krisztián Juhász
AU - Matta, Csaba
AU - Paluska, Luca
AU - Al-Mnaseer, Ahmed
AU - Takács, Roland
AU - Ducza, László
TI - A Systems-Level Transcriptomic Framework Identifies Shared Cellular Hubs in Osteoarthritis and Alzheimer's Disease
T2 - Computational and structural biotechnology journal
J2 - Comput Struct Biotechnol J
PY - 2026
DA - 2026/05/13
VL - 35
IS - 1
SP - 0085
SN - 2001-0370
PB - AAAS Science Partner Journal Program
DO - 10.34133/csbj.0085
UR - https://doi.org/10.34133/csbj.0085
LA - en
ER -

CSL-JSON

{
"id": "10.34133/csbj.0085",
"type": "article-journal",
"title": "A Systems-Level Transcriptomic Framework Identifies Shared Cellular Hubs in Osteoarthritis and Alzheimer's Disease",
"container-title": "Computational and structural biotechnology journal",
"author": [
{
"family": "Wang",
"given": "Zhangzheng"
},
{
"family": "Zoltán",
"given": "Krisztián Juhász"
},
{
"family": "Matta",
"given": "Csaba"
},
{
"family": "Paluska",
"given": "Luca"
},
{
"family": "Al-Mnaseer",
"given": "Ahmed"
},
{
"family": "Takács",
"given": "Roland"
},
{
"family": "Ducza",
"given": "László"
}
],
"container-title-short": "Comput Struct Biotechnol J",
"volume": "35",
"issue": "1",
"page": "0085",
"DOI": "10.34133/csbj.0085",
"PMID": "42137879",
"PMCID": "PMC13168759",
"ISSN": "2001-0370",
"publisher": "AAAS Science Partner Journal Program",
"URL": "https://doi.org/10.34133/csbj.0085",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
13
]
]
}
}

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.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, edgeR, limma, 7 other tools, genetics / omics, cellular / molecular, 5 references
[2] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, edgeR, limma, 7 other tools, Alzheimer's / dementia, cellular / molecular, 2 references
[3] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, edgeR, limma, 7 other tools, genetics / omics, cellular / molecular
[4] doi:10.1016/j.isci.2026.115657 [code]
Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: Harmony, edgeR, limma, 7 other tools
[5] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: edgeR, limma, igraph, 5 other tools, genetics / omics, cellular / molecular, 3 references
[6] doi:10.1186/s12974-026-03838-8 [code]
Acarbose modulates microglial Pkm2 acetylation to reshape immunometabolism and preserve retinal neurons after ischemia-reperfusion.
Journal: Journal of neuroinflammation
In common: Harmony, limma, DESeq2, 6 other tools, genetics / omics, cellular / molecular, 1 reference
[7] doi:10.1038/s41398-026-04200-5 [code]
Postmortem brain single-nucleus and bulk gene expression analyses identify shared and distinct abnormalities in bipolar disorder and major depressive disorder.
Journal: Translational psychiatry
In common: Harmony, edgeR, limma, 5 other tools, genetics / omics, cellular / molecular, 2 references
[8] doi:10.1186/s13059-026-04152-5 [code]
CroCoNet: a framework for the quantitative comparison of gene regulatory networks across species.
Journal: Genome biology
In common: Harmony, limma, DESeq2, 6 other tools, 1 reference
[9] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: Harmony, edgeR, DESeq2, 5 other tools, genetics / omics, cellular / molecular, 2 references
[10] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: edgeR, limma, DESeq2, 6 other tools, genetics / omics, 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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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