OSCR

Spatiotemporal brain transcriptomics reveal risk gene hot-spots in major neuropsychiatric disorders.

Code ↔ Paper

15 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 15 matches
  1. [1] § Results › Spatial enrichment patterns of neuropsychiatric risk gene sets implicate disease-relevant brain regions and circuits ↔ R code for data analysis.R, lines 113–194 · score 0.92 · parietal cortex, ganglionic eminence, occipital cortex, frontal cortex, brain stem, temporal cortex
  2. [2] § Methods › Function analysis of WGCNA module genes ↔ R code for data analysis.R, lines 458–545 · score 0.83 · enriched GO terms, clusterProfiler, compareCluster, KEGG pathway, cutoff, Module genes
  3. [3] § Methods › Weighted gene co-expression network (WGCNA) analysis ↔ R code for data analysis.R, lines 188–248 · score 0.81 · soft threshold power, blockwiseModules, cut height, Pearson, TOM, network
  4. [4] § Methods › Weighted gene co-expression network (WGCNA) analysis ↔ R code for data analysis.R, lines 741–801 · score 0.81 · soft threshold power, blockwiseModules, cut height, Pearson, TOM, network
  5. [5] § Methods › Statistics and reproducibility ↔ R code for data analysis.R, lines 2–62 · score 0.80 · Mann Whitney, ComplexHeatmap, remaining genes, selected gene, expression matrix, Rank Biserial
  6. [6] § Methods › Statistics and reproducibility ↔ R code for data analysis.R, lines 2–62 · score 0.80 · Mann Whitney, ComplexHeatmap, remaining genes, selected gene, expression matrix, Rank Biserial
  7. [7] § Results › Spatial enrichment patterns of neuropsychiatric risk gene sets implicate disease-relevant brain regions and circuits ↔ R code for data analysis.R, lines 64–111 · score 0.80 · parietal cortex, ganglionic eminence, occipital cortex, frontal cortex, brain stem, temporal cortex
  8. [8] § Methods › Function analysis of WGCNA module genes ↔ R code for data analysis.R, lines 1044–1105 · score 0.75 · clusterProfiler, compareCluster, KEGG pathway, cutoff, Module genes, enriched
  9. [9] § Results › Developmental expression enrichment patterns of risk gene sets align with the peak age of onset across major neuropsychiatric disorders ↔ R code for data analysis.R, lines 64–111 · score 0.70 · peak ages, early infancy, middle adulthood, early prenatal, enrichment patterns, onset
  10. [10] § Results › Developmental expression enrichment patterns of risk gene sets align with the peak age of onset across major neuropsychiatric disorders ↔ R code for data analysis.R, lines 64–111 · score 0.70 · peak ages, early infancy, middle adulthood, early prenatal, enrichment patterns, onset
  11. [11] § Methods › Risk gene sets expression enrichment analysis in the brain spatiotemporal transcriptome › Curation of genome-wide risk genes for the 15 traits ↔ R code for data analysis.R, lines 113–186 · score 0.68 · SZes, panic disorder, SZgw, date, ANO, TS
  12. [12] § Methods › Risk gene sets expression enrichment analysis in the brain spatiotemporal transcriptome › Curation of genome-wide risk genes for the 15 traits ↔ R code for data analysis.R, lines 113–194 · score 0.66 · SZes, panic disorder, SZgw, ANO, TS, neuroticism
  13. [13] § Results › Assessing the biological relevance of risk gene expression patterns using ENIGMA brain morphometry ↔ R code for data analysis.R, lines 683–739 · score 0.57 · d_icv, amyg, caud, hippo, pal, thal
  14. [14] § Results › Spatiotemporal enrichment patterns of risk gene sets differentiate neuropsychiatric disorder subgroups ↔ R code for data analysis.R, lines 197–235 · score 0.57 · variance explained, Principal component, Scree, rank biserial, Bar, PC
  15. [15] § Results › Spatiotemporal enrichment patterns of risk gene sets differentiate neuropsychiatric disorder subgroups ↔ R code for data analysis.R, lines 237–262 · score 0.51 · panic disorder, SZgw, PC4, PC1, epilepsy, PD

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 · 1,145 lines · 41 KB · CC-BY-4.0 · 9 matches

  2. ##Mann-Whitney U test
  3. # install.packages("ggcorrplot")
  4. # BiocManager::install("ComplexHeatmap")
  5. library(gmodels)
  6. library("data.table")
  7. library("dplyr")
  8. library(ggplot2)
  9. library(ggpubr)
  10. library(effectsize)
  11. library("ComplexHeatmap")
  13. # Mann-Whitney U test based on RPKM values ------------------------------------
  14. riskgene<-read.table("riskgenes.CSV",header=TRUE,sep=",",row.names=1)
  15. str(riskgene)
  16. braingene<-read.table("expression_matrix-ave.CSV",header=TRUE,sep=",",row.names=1)
  17. str(braingene)
  18. dim(braingene)
  19. gene <- rownames(braingene)
  20. length(gene)
  21. P_values <- data.frame()
  22. W<-data.frame()
  23. r_rank_biserial<-data.frame()
  24. for (i in 1:length(riskgene))
  25. {
  26. interesting_genes <- c(riskgene[, i])
  27. interesting_genes <- interesting_genes[interesting_genes[] != ""]
  28. selected_genes <- braingene[interesting_genes,]
  29. dim(selected_genes)
  30. selected_genes <- na.omit(selected_genes)
  31. dim(selected_genes)
  32. remaining_genes <-
  33. braingene[!rownames(braingene) %in% interesting_genes,]
  34. dim(remaining_genes)
  35. remaining_genes
  36. for (j in 1:524) {
  37. mwu_test_result <-
  38. wilcox.test(selected_genes[, j],
  39. remaining_genes[, j],
  40. alternative = "greater",
  41. verbose = FALSE)
  42. EFS_result <-
  43. cliffs_delta(
  44. selected_genes[, j],
  45. remaining_genes[, j],
  46. mu = 0,
  47. alternative = "greater",
  48. verbose = F
  49. )
  50. P_values[j, i] <- mwu_test_result$p.value
  51. r_rank_biserial[j, i] <- EFS_result$r_rank_biserial
  52. }
  53. }
  54. rownames(P_values)<-c(colnames(braingene))
  55. rownames(P_values)
  56. colnames(P_values)<-c(colnames(riskgene))
  57. colnames(P_values)
  58. rownames(r_rank_biserial)<-c(colnames(braingene))
  59. colnames(r_rank_biserial)<-c(colnames(riskgene))
  60. write.csv(P_values, file="P_Value.CSV", quote = F,row.names = T)
  61. write.csv(r_rank_biserial,file="r_rank_biserial.CSV",
  62. quote=F,row.names = T)
  64. ##complexheatmap of risk gene enrichment patterns
  65. brain<-read.table("brian regions.CSV",header=TRUE,sep=",")
  66. head(brain)
  67. brain_order<-brain[order(brain$column_num),]
  68. head(brain_order)
  69. data1<-read.table("r_rank_biserial.CSV",header=TRUE,sep=",",row.names=1)
  70. data1<-data1[brain_order$sample,]
  71. head(data1)
  72. data2<- t(data1)
  73. dim(data2)
  74. data2 <- as.matrix(data2)
  75. str(data2)
  76. row_anno <-read.table("disease_Annotation_peak age onset.CSV",header=TRUE,sep=",",row.names=1,encoding = "UTF-8")
  77. dim(row_anno)
  78. str(row_anno)
  79. row_anno$peak_age_onset<-factor(row_anno$peak_age_onset,
  80. levels = c("Y4.5","Y5.5","Y8","Y9.5","Y15.5","Y19.5","Y12.8_Y24.9","Y20.5","Childhood_Y60","Y85","Y87"))
  81. col_anno <-read.table("brain_region_Annotation.CSV",header=TRUE,sep=",",row.names=1)
  82. str(col_anno)
  83. col_anno<-col_anno[,-c(2:3)]
  84. col_anno<-col_anno[brain_order$sample,]
  85. head(col_anno)
  86. col_anno$structure<-factor(col_anno$structure,
  87. levels = c("frontal_cortex","parietal_cortex","temporal_cortex","occipital_cortex","ganglionic_eminence","striatum","thalamus","limbic_system","brain_stem","cerebellum"))
  88. col_anno$dev_stage<-factor(col_anno$dev_stage,
  89. levels = c("Early_prenatal","Early_mid_prenatal","Late_mid_prenatal","Late_prenatal","Early_infancy","Late_infancy","Early_childhood","Late_childhood","Adolescence","Young_adulthood","Middle_adulthood" ))
  91. ann_colors = list(peak_age_onset = c(Y4.5 = "#CCFFCC",Y5.5="#CCFF00", Y8="#00CC99", Y9.5="#00CC00", Y15.5="#FF99FF", Y19.5="#FF33FF",
  92. Y12.8_Y24.9="#00FFFF" , Y20.5="#33CCCC", Childhood_Y60="#3366CC", Y85="#333399", Y87="#111111"),
  93. dev_stage=c(Early_prenatal="#000000",Early_mid_prenatal="#333333",Late_mid_prenatal="#666666",Late_prenatal="#888888",Early_infancy="#DDDDDD", Late_infancy="#EEEEEE", Early_childhood="#00FFCC",Late_childhood="#00CC99",Adolescence="#FFCCFF",Young_adulthood="#33CCCC",Middle_adulthood="#333399"),
  94. structure= c(frontal_cortex="#FFFF00",parietal_cortex="pink",temporal_cortex="#CCFF00",occipital_cortex="orange" , limbic_system = "#CC0000", striatum = "#00CC66", thalamus = "#00FF00",brain_stem="#CC00FF", ganglionic_eminence="blue", cerebellum = "#330066")
  95. )
  96. pdf(file = "complexheatmap-r_rank_biserial.pdf", width = 12, height = 8)
  98. pheatmap(data2, scale = "row",show_rownames = T,
  99. clustering_method="ward.D2",
  100. show_colnames = F,
  101. cluster_cols=F,cluster_rows=F,drop_levels=T,
  102. col = colorRampPalette(c("blue","white","red"))(100),
  103. annotation_col = col_anno,
  104. annotation_row = row_anno,
  105. annotation_colors = ann_colors,
  106. heatmap_legend_param = list(title = "normalized r") ,
  107. annotation_names_row = T,
  108. annotation_names_col = T ,
  109. column_title = NULL
  110. )
  111. dev.off()
  113. # # FC-other structures(without BS)，pairwise t-test -----------------------------
  115. library(dplyr)
  116. library(ggplot2)
  117. library(rstatix)
  118. library(ggpubr)
  120. setwd("Major 15 disorder merged gene-RPKM-ave")
  121. data1<-read.csv("combine_r_rank_biserial-11structure-scaled-line.CSV",header=TRUE,sep=",")
  122. str(data1)
  124. data1$structure <- factor(data1$structure,
  125. levels = c("frontal_cortex", "parietal_cortex", "temporal_cortex", "occipital_cortex",
  126. "ganglionic_eminence","striatum","thalamus","amygdala","brain_stem","hippocampus","cerebellum"))
  127. data1$disorder <- factor(data1$disorder,
  128. levels = c("Neuroticism","BIP","TS", "SZes","IQ", "ASD","ADHD","ANO","PanicDisorder", "OCD","AD", "Epilepsy", "PD","MDD", "SZgw"))
  129. data1 <- data1[data1$structure != "brain_stem", ]
  130. df_summary <- data1 %>%
  131. group_by(disorder, structure) %>%
  132. summarise(
  133. mean_value = mean(r, na.rm = TRUE),
  134. sd_value = sd(r, na.rm = TRUE),
  135. se_value = sd(r, na.rm = TRUE) / sqrt(n()),
  136. n = n(),
  137. .groups =NULL
  138. )
  139. df_summary
  140. stat_res <- data1 %>%
  141. group_by(disorder) %>%
  142. pairwise_t_test(
  143. r ~ structure,
  144. ref.group = "frontal_cortex",
  145. p.adjust.method = "fdr"
  146. )
  147. stat_res
  148. stat_frontal <- stat_res %>%
  149. filter(group1 == "frontal_cortex")
  150. stat_frontal <- stat_frontal %>%
  151. mutate(stars = case_when(
  152. p.adj < 0.001 ~ "***",
  153. p.adj < 0.01 ~ "**",
  154. p.adj < 0.05 ~ "*",
  155. TRUE ~ ""
  156. ))
  157. write.csv(stat_frontal,file="stat_frontal_noBS.CSV", quote = F,row.names = F)
  159. plot_data <- df_summary %>%
  160. left_join(stat_frontal,
  161. by = c("structure" = "group2", "disorder")) %>%
  162. mutate(y_pos = mean_value + se_value + 0.02)
  165. ggplot() +
  166. # dot
  167. geom_jitter(data = data1,
  168. aes(x = structure, y = r),
  169. width = 0.15, height = 0.3,
  170. size = 0.1, alpha = 0.3, color = "black") +
  171. # bar
  172. geom_col(data = df_summary,
  173. aes(x = structure, y = mean_value, fill = structure),
  174. width = 0.7, alpha = 0.9) +
  175. # SE
  176. geom_errorbar(data = df_summary,
  177. aes(x = structure,
  178. ymin = mean_value - se_value,
  179. ymax = mean_value + se_value),
  180. width = 0.25, linewidth = 0.2) +
  181. # sig
  182. geom_text(data = plot_data,
  183. aes(x = structure, y = y_pos, label = stars),
  184. size = 3, color = "red", vjust = 0.3) +
  186. facet_grid(rows = vars(disorder), scales = "free_y") +
  187. scale_y_continuous(n.breaks =3) +
  188. theme_bw() +
  189. theme(
  190. strip.text.y = element_text(angle = 0, size =7),
  191. axis.text.x = element_text(angle = 30, size = 8,hjust = 0.9),
  192. axis.text.y = element_text(size = 5),
  193. ) +
  194. labs(y = "Mean ± SE")
  197. # ##PCA analysis -----------------------------------------------------------------
  198. library(ggplot2)
  199. library(ggrepel)
  201. data1<-read.table("r_rank_biserial_rpkm-ave.CSV",header=TRUE,sep=",",row.names=1)
  202. head(data1)
  203. data2<- scale(data1)
  204. data2<- t(data2)
  205. disor<-rownames(data2)
  206. disor
  207. mean(data2[,3])
  208. pca_result <- prcomp(data2,center=T,scale=T)
  209. pca_result
  210. summary(pca_result)
  212. pca_result$x
  213. screeplot(pca_result)
  214. plot(pca_result, type = "barplot")
  215. pca_result$x
  216. rownames(pca_result$x)
  218. var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)
  219. df_var <- data.frame(
  220. PC = paste0("PC", 1:length(var_explained)),
  221. Variance = var_explained
  222. )
  223. var_explained
  224. df_var$PC <- factor(df_var$PC, levels = df_var$PC[order(df_var$Variance, decreasing = TRUE)])
  225. # Scree plot
  226. ggplot(df_var, aes(x = PC, y = Variance)) +
  227. geom_bar(stat = "identity") +
  228. geom_line(aes(group = 1)) +
  229. geom_point(size = 2) +
  230. ylab("Proportion of Variance Explained") +
  231. xlab("Principal Components") +
  232. theme_bw(base_size = 14) +
  233. theme(
  234. axis.text.x = element_text(angle = 45, hjust = 1)
  235. )
  237. ##disorder order
  238. disor<-rownames(pca_result$x)
  239. disor
  240. ##PCA-15 disorder dot plot
  241. pca_df <- data.frame(
  242. PC1 = pca_result$x[,1],
  243. PC4 = pca_result$x[,4]
  244. )
  245. pca_df
  246. pca_result$x
  247. write.csv(pca_result$x,file="PCA result.csv")
  248. ggplot(pca_df, aes(PC1, PC4, color =disor)) +
  249. geom_point(size = 3) +
  250. geom_text_repel(
  251. aes(label = disor),
  252. size = 3,
  253. max.overlaps = 20,
  254. box.padding = 0.5,
  255. point.padding = 0.3,
  256. min.segment.length = 0.2,
  257. segment.color = "grey50"
  258. )+
  259. scale_color_manual(values = c("IQ"="darkgrey", "ASD"="red","TS" ="firebrick", "Neuroticism"="darkgreen", "ADHD"="#FFdd22", "OCD" ="#0099FF", "ANO" = "#00CC00","PanicDisorder" ="#FF99FF", "BIP"="#FF33FF","MDD"= "orange", "SZes"="#33CCCC", "SZgw"="blue",
  260. "Epilepsy"= "#3366CC", "AD"="purple", "PD"="#111111"))+
  262. theme_bw()
  265. # ##corr::PCA component and spatiotemporal variables --------------------------------------------------
  267. library(dplyr)
  268. library(ggplot2)
  270. scaled<-read.csv("combine_r_rank_biserial-11structure-scaled.CSV",header=TRUE,sep=",")
  271. str(scaled)
  272. head(scaled)
  273. ##disorder order
  274. disor<-rownames(pca_result$x)
  275. disor
  277. scaled$structure <- factor(scaled$structure,
  278. levels = c("frontal_cortex", "parietal_cortex", "temporal_cortex", "occipital_cortex",
  279. "ganglionic_eminence","striatum","thalamus","amygdala","hippocampus","brain_stem","cerebellum"))
  281. scaled_summary <- scaled %>%
  282. group_by(structure) %>%
  283. summarise(across(disor, mean, na.rm = TRUE))
  284. scaled_summary_t<-t(scaled_summary)
  285. scaled_summary_t
  286. colnames(scaled_summary_t)
  287. colnames(scaled_summary_t) <- as.character(scaled_summary_t[1, ])
  288. scaled_summary_t
  289. scaled_summary_t <- scaled_summary_t[-1, ] # 删除第一行
  290. scaled_summary_t
  291. write.csv(scaled_summary_t,file="scaled_15disorder_structure summary.CSV",quote=F,row.names = T)
  292. scaled_summary<-read.csv("scaled_15disorder_structure summary.CSV",header=T,sep=",",row.names = 1)
  293. rownames(scaled_summary)
  294. colnames(scaled_summary)
  296. p2 <- pca_result$x[, c(1:4), drop = FALSE]
  297. df_pca_asso<-cbind(p2,scaled_summary)
  298. df_pca_asso
  299. help(rcorr)
  300. library(Hmisc)
  301. res<-rcorr(as.matrix(df_pca_asso),type="spearman")
  302. res
  303. R <- res$r
  304. P <- res$P
  305. vars <- colnames(df_pca_asso)
  306. vars
  307. library(tidyverse)
  309. n <- length(vars)
  310. upper_idx <- upper.tri(P, diag = FALSE)
  312. p_unique <- P[upper_idx]
  313. p_unique
  314. p_unique <- p_unique[!is.na(p_unique)]
  315. p_unique
  316. # FDR
  317. p_unique_fdr <- p.adjust(p_unique, method = "fdr")
  318. P_fdr <- matrix(NA, nrow = n, ncol = n)
  319. colnames(P_fdr) <- rownames(P_fdr) <- vars
  320. P_fdr[upper_idx] <- p_unique_fdr
  321. P_fdr[lower.tri(P_fdr)] <- t(P_fdr)[lower.tri(P_fdr)]
  323. diag(P_fdr) <- NA
  324. R
  325. P_fdr
  326. df_plot <- R %>%
  327. as.data.frame() %>%
  328. rownames_to_column("Var1") %>%
  329. pivot_longer(-Var1, names_to = "Var2", values_to = "cor") %>%
  330. left_join(
  331. P_fdr %>%
  332. as.data.frame() %>%
  333. rownames_to_column("Var1") %>%
  334. pivot_longer(-Var1, names_to = "Var2", values_to = "p"),
  335. by = c("Var1", "Var2")
  336. )
  338. df_plot$Var1 <- factor(df_plot$Var1, levels = vars)
  339. df_plot$Var2 <- factor(df_plot$Var2, levels = vars)
  340. df_plot
  341. ##dot plot
  342. df_plot <- df_plot %>%
  343. mutate(
  344. abs_cor = abs(cor),
  345. logp = -log10(p + 1e-10)
  346. )
  347. df_plot
  348. df_plot <- df_plot %>%
  349. mutate(
  350. p_fdr = as.vector(P_fdr[cbind(
  351. match(Var1, vars),
  352. match(Var2, vars)
  353. )]),
  354. sig = case_when(
  355. p_fdr < 0.001 ~ "***",
  356. p_fdr < 0.01 ~ "**",
  357. p_fdr < 0.05 ~ "*",
  358. TRUE ~ ""
  359. )
  360. )
  361. df_plot
  362. library(ggplot2)
  364. # bubble plot
  366. ggplot(df_plot, aes(Var1, Var2)) +
  367. geom_point(aes(size = logp,
  368. color = cor)) +
  369. geom_text(aes(label = sig), size = 4, vjust = 0.7) + # 星号标注
  370. scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
  371. limits = range(df_plot$cor, na.rm = TRUE)) +
  372. scale_size(range = c(0, 10)) +
  373. scale_alpha(range = c(0.2, 1)) +
  374. theme_minimal()+
  375. theme(
  376. axis.text.x = element_text(angle = 45, hjust = 1),
  377. panel.grid = element_blank()
  378. ) +
  379. labs(
  380. size = "-log10(p)",
  381. color = "r",
  382. x = NULL, y = NULL,
  383. title = "Cor: PC1-4/ structures"
  384. )
  387. ## PC1 and variables
  388. library(dplyr)
  389. library(ggplot2)
  391. # ## PC and developmental stages----------------------------------------------------
  393. scaled<-read.csv("combine_r_rank_biserial-10structure-scaled.CSV",header=TRUE,sep=",")
  394. str(scaled)
  395. head(scaled)
  396. scaled$dev_stage <- factor(scaled$dev_stage,
  397. levels = c("Early_prenatal" ,"Early_mid_prenatal","Late_mid_prenatal","Late_prenatal", "Early_infancy","Late_infancy", "Early_childhood", "Late_childhood", "Adolescence",
  398. "Young_adulthood","Middle_adulthood"))
  400. ##disorder
  401. disor<-rownames(pca_result$x)
  402. disor
  404. scaled_summary <- scaled %>%
  405. group_by(dev_stage) %>%
  406. summarise(across(disor, mean, na.rm = TRUE))
  407. scaled_summary_t<-t(scaled_summary)
  408. scaled_summary_t
  409. colnames(scaled_summary_t)
  410. colnames(scaled_summary_t) <- as.character(scaled_summary_t[1, ])
  411. scaled_summary_t
  412. scaled_summary_t <- scaled_summary_t[-1, ]
  413. scaled_summary_t
  414. write.csv(scaled_summary_t,file="scaled_15disorder_dev_stage summary.CSV",quote=F,row.names = T)
  415. scaled_summary<-read.csv("scaled_15disorder_dev_stage summary.CSV",header=T,sep=",",row.names = 1)
  416. rownames(scaled_summary)
  417. colnames(scaled_summary)
  419. p2 <- pca_result$x[, c(1:4), drop = FALSE]
  420. df_pca_asso<-cbind(p2,scaled_summary)
  421. df_pca_asso
  423. library(Hmisc)
  424. res<-rcorr(as.matrix(df_pca_asso),type = "spearman")
  425. res
  426. R <- res$r
  427. P <- res$P
  428. vars <- colnames(df_pca_asso)
  429. vars
  430. library(tidyverse)
  432. n <- length(vars)
  434. upper_idx <- upper.tri(P, diag = FALSE)
  436. p_unique <- P[upper_idx]
  437. p_unique
  438. p_unique <- p_unique[!is.na(p_unique)]
  439. p_unique
  440. # FDR
  441. p_unique_fdr <- p.adjust(p_unique, method = "fdr")
  442. P_fdr <- matrix(NA, nrow = n, ncol = n)
  443. colnames(P_fdr) <- rownames(P_fdr) <- vars
  445. P_fdr[upper_idx] <- p_unique_fdr
  447. P_fdr[lower.tri(P_fdr)] <- t(P_fdr)[lower.tri(P_fdr)]
  449. diag(P_fdr) <- NA
  450. R
  451. P_fdr
  452. df_plot <- R %>%
  453. as.data.frame() %>%
  454. rownames_to_column("Var1") %>%
  455. pivot_longer(-Var1, names_to = "Var2", values_to = "cor") %>%
  456. left_join(
  457. P_fdr %>%
  458. as.data.frame() %>%
  459. rownames_to_column("Var1") %>%
  460. pivot_longer(-Var1, names_to = "Var2", values_to = "p"),
  461. by = c("Var1", "Var2")
  462. )
  464. df_plot$Var1 <- factor(df_plot$Var1, levels = vars)
  465. df_plot$Var2 <- factor(df_plot$Var2, levels = vars)
  466. df_plot
  468. df_plot <- df_plot %>%
  469. mutate(
  470. abs_cor = abs(cor),
  471. logp = -log10(p + 1e-10) # 避免 log(0)
  472. )
  473. df_plot
  474. df_plot <- df_plot %>%
  475. mutate(
  476. p_fdr = as.vector(P_fdr[cbind(
  477. match(Var1, vars),
  478. match(Var2, vars)
  479. )]),
  480. sig = case_when(
  481. p_fdr < 0.001 ~ "***",
  482. p_fdr < 0.01 ~ "**",
  483. p_fdr < 0.05 ~ "*",
  484. TRUE ~ ""
  485. )
  486. )
  488. ggplot(df_plot, aes(Var1, Var2)) +
  489. geom_point(aes(size = logp,
  490. color = cor)) +
  491. geom_text(aes(label = sig), size = 4, vjust = 0.7) + # 星号标注
  492. scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
  493. limits = range(df_plot$cor, na.rm = TRUE)) +
  494. scale_size(range = c(0, 10)) +
  495. scale_alpha(range = c(0.2, 1)) +
  496. theme_minimal()+
  497. theme(
  498. axis.text.x = element_text(angle = 45, hjust = 1),
  499. panel.grid = element_blank()
  500. ) +
  501. labs(
  502. size = "-log10(p)",
  503. color = "r",
  504. x = NULL, y = NULL,
  505. title = "Cor: PC1-4/ dev_stage"
  506. )
  509. # ### PC and single variable ----------------------------------------------------
  510. ## PC2 and late infancy
  511. library(dplyr)
  512. library(ggplot2)
  513. scaled<-read.csv("combine_r_rank_biserial-10structure-scaled.CSV",header=TRUE,sep=",")
  514. str(scaled)
  515. head(scaled)
  516. scaled$dev_stage2 <- ifelse(scaled$dev_stage == "Late_infancy",
  517. "Late_infancy",
  518. "other")
  520. scaled$dev_stage2 <- factor(scaled$dev_stage2,
  521. levels = c("Late_infancy", "other"))
  522. summary(scaled$dev_stage2)
  524. disor<-rownames(pca_result$x)
  525. disor
  527. scaled_summary <- scaled %>%
  528. group_by(dev_stage2) %>%
  529. summarise(across(disor, mean, na.rm = TRUE))
  530. scaled_summary_t<-t(scaled_summary)
  531. scaled_summary_t
  532. colnames(scaled_summary_t)
  533. colnames(scaled_summary_t) <- as.character(scaled_summary_t[1, ])
  534. scaled_summary_t
  535. scaled_summary_t <- scaled_summary_t[-1, ]
  536. scaled_summary_t
  537. write.csv(scaled_summary_t,file="scaled_15disorder_late infancy summary.CSV",quote=F,row.names = T)
  538. scaled_summary<-read.csv("scaled_15disorder_late infancy summary.CSV",header=T,sep=",",row.names = 1)
  539. rownames(scaled_summary)
  540. colnames(scaled_summary)
  541. as.data.frame(pca_result$x[,1])
  543. df_pca_asso<-cbind(pca_result$x[,1:4],scaled_summary)
  544. df_pca_asso$PC1
  546. library(Hmisc)
  547. res<-rcorr(as.matrix(df_pca_asso))
  548. res
  549. disor
  550. ##"Cor: PC2/ Late_infancy enrichment"
  551. library(ggpubr)
  552. ggplot(data = df_pca_asso, aes(x = PC2, y =Late_infancy )) +
  553. geom_point() +
  554. geom_smooth(method = "lm", se = TRUE, color = "red") +
  555. stat_cor(method = "spearman",label.x.npc = 0.1, label.y.npc = 1.0,) +
  556. labs(title = "Cor: PC2/ Late_infancy",
  557. x = "PC2",
  558. y = "Late_infancy enrichment") +
  559. theme_minimal()+
  560. geom_text_repel(
  561. aes(label = disor),
  562. size = 3,
  563. max.overlaps = 20,
  564. box.padding = 0.5,
  565. point.padding = 0.3,
  566. min.segment.length = 0.2,
  567. segment.color = "grey50"
  568. )
  571. ## PC4 and cerebellum
  572. library(dplyr)
  573. library(ggplot2)
  574. scaled<-read.csv("combine_r_rank_biserial-11structure-scaled.CSV",header=TRUE,sep=",")
  575. str(scaled)
  576. head(scaled)
  577. scaled$structure2 <- ifelse(scaled$structure == "cerebellum",
  578. "cerebellum",
  579. "other")
  581. scaled$structure2 <- factor(scaled$structure2,
  582. levels = c("cerebellum", "other"))
  583. summary(scaled$structure2)
  584. disor<-rownames(pca_result$x)
  585. disor
  586. scaled_summary <- scaled %>%
  587. group_by(structure2) %>%
  588. summarise(across(disor, mean, na.rm = TRUE))
  589. scaled_summary_t<-t(scaled_summary)
  590. scaled_summary_t
  591. colnames(scaled_summary_t)
  592. colnames(scaled_summary_t) <- as.character(scaled_summary_t[1, ])
  593. scaled_summary_t
  594. scaled_summary_t <- scaled_summary_t[-1, ]
  595. scaled_summary_t
  596. write.csv(scaled_summary_t,file="scaled_15disorder_cerebellum summary.CSV",quote=F,row.names = T)
  597. scaled_summary<-read.csv("scaled_15disorder_cerebellum summary.CSV",header=T,sep=",",row.names = 1)
  598. rownames(scaled_summary)
  599. colnames(scaled_summary)
  600. as.data.frame(pca_result$x[,1])
  601. df_pca_asso<-cbind(pca_result$x[,1:4],scaled_summary)
  602. df_pca_asso$PC1
  604. library(Hmisc)
  605. res<-rcorr(as.matrix(df_pca_asso))
  606. res
  607. disor
  608. ##"Cor: PC4/ 和cerebellum enrichment"
  609. library(ggpubr)
  610. ggplot(data = df_pca_asso, aes(x = PC4, y =cerebellum )) +
  611. geom_point() +
  612. geom_smooth(method = "lm", se = TRUE, color = "red") +
  613. stat_cor(method = "spearman",label.x.npc = 0.2, label.y.npc = 1.0,) +
  614. labs(title = "Cor: PC4/ cerebellum enrichment",
  615. x = "PC4",
  616. y = "cerebellum enrichment") +
  617. theme_minimal()+
  618. geom_text_repel(
  619. aes(label = disor),
  620. size = 3,
  621. max.overlaps = 20,
  622. box.padding = 0.5,
  623. point.padding = 0.3,
  624. min.segment.length = 0.2,
  625. segment.color = "grey50"
  626. )
  628. # ##ENIGMA-AHBA enrichment correlations ------------------------
  629. library("data.table")
  630. library("dplyr")
  631. library(ggplot2)
  632. library(ggpubr)
  633. brainregion <-read.csv("brain_region_Annotation-6-fine.CSV",header=TRUE,sep=",")
  634. colnames(brainregion)
  635. head(brainregion)
  636. disord <-read.csv("combine_r_rank_biserial_rpkm-ave.CSV",header=TRUE,sep=",")
  637. head(disord)
  639. Disor_region_anno<-dplyr::left_join(disord,brainregion,by="X")
  640. head(Disor_region_anno)
  642. Disor_region_anno_enig <- Disor_region_anno %>%
  643. filter(!is.na(Structure_comn)) %>%
  644. group_by(Structure_comn) %>%
  645. summarise(across(all_of(colnames(Disor_region_anno)[-c(1,17:21)]), ave))
  646. Disor_region_anno_enig<- unique(Disor_region_anno_enig, fromLast = TRUE)
  647. Disor_region_anno_enig
  648. write.csv(Disor_region_anno_enig, file="Disor_region_anno.CSV", quote = F,row.names = T)
  649. # ##ASD ------------------------------------------------------------------
  650. enig_sub <-read.csv("asd_meta-analysis_case-controls_CortThick.CSV",header=TRUE,sep=",")
  651. enig_tck <-read.csv("asd_meta-analysis_case-controls_SubVol.CSV",header=TRUE,sep=",")
  653. enig<-rbind(enig_sub[,c(2:3,9)],enig_tck[,c(2:3,9)])
  654. enig$Structure<-sub('L_','',enig$Structure,ignore.case = F)
  655. enig$Structure<-sub('R_','',enig$Structure,ignore.case = F)
  656. enig$Structure<-sub('LL','Ll',enig$Structure,ignore.case = F)
  657. enig$Structure<-sub('RL','Rl',enig$Structure,ignore.case = F)
  658. enig$Structure<-sub('L','',enig$Structure,ignore.case = F)
  659. enig$Structure<-sub('R','',enig$Structure,ignore.case = F)
  660. enig
  661. enig2 <- enig %>%
  662. group_by(Structure) %>%
  663. summarise( d_icv_ave = mean(d_icv, na.rm = TRUE))
  664. print(enig2,n=40)
  666. region_enig <-read.csv("brain region-enigma.CSV",header=TRUE,sep=",")
  667. region_enig
  668. enig_comn<-dplyr::left_join(region_enig,enig2,by="Structure")
  669. enig_comn
  671. enig_merge<-dplyr::left_join(enig_comn,Disor_region_anno_enig,by="Structure_comn")
  672. str(enig_merge)
  673. enig_merge2 <- enig_merge %>%
  674. filter(!is.na(Structure_comn)) %>%
  675. group_by(Structure_comn) %>%
  676. summarise(across(all_of(colnames(enig_merge)[c(4:19)]), ave))
  678. enig_merge2<- unique(enig_merge2, fromLast = TRUE)
  680. enig_merge3<-enig_merge2 %>%
  681. filter(!is.na(d_icv_ave))
  683. ##dotplot
  684. ggplot(data = enig_merge3, aes(x = d_icv_ave, y =ASD )) +
  685. geom_point() +
  686. geom_smooth(method = "lm", se = TRUE, color = "red") +
  687. stat_cor(method = "spearman",label.x.npc = 'left', label.y.npc = "top",) +
  688. labs(title = "Correlation: brain MRI changes/ risk gene overexpression",
  689. x = "ASD_case-controls-d_icv",
  690. y = "ASD gene overexpression (rank-biserial)") +
  691. theme_minimal()+
  692. geom_text(
  693. data = subset(enig_merge3, Structure_comn %in% c("pal", "thal", "precuneus","lingual","caud", "hippo",
  694. "amyg","cuneus","put")),
  695. aes(label = Structure_comn),
  696. vjust = -0.5, # 标签垂直偏移
  697. hjust = 0.1 # 标签水平偏移
  698. )
  701. # WGCNA -----------------------------------------------------------------
  702. # BiocManager::install("WGCNA")
  703. # BiocManager::install("preprocessCore")
  704. library("WGCNA")
  705. library("glue")
  706. ## export dir
  707. outdir <- "C:/WGCNA"
  708. dir.create(file.path(outdir), showWarnings = FALSE)
  710. data1<-read.table("expression matrix_RPKM.CSV",header=TRUE,sep=",",row.names=1)
  711. dim(data1)
  712. row_sums <- rowSums(data1)
  713. sort(row_sums)[1:10]
  714. # delete very low expression genes
  715. row_sums
  716. data2 <- data1[row_sums >10, ]
  717. dim(data2)
  718. data2<-t(data2)
  719. head(data2)
  720. ## Plot a sample dendrogram
  721. setwd("C:/WGCNA")
  722. sampleTree <- hclust(dist(data2), method = "average")
  723. pdf(file = "Sample_dendrogram_rm_outlier.pdf", width = 12, height = 9)
  724. par(cex = 0.6)
  725. par(mar = c(2,5,5,1))
  726. plot(sampleTree, main = "Sample clustering to detect outliers",
  727. sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
  728. abline(h = 100000, col = "red")
  729. dev.off()
  730. ## clust 1 contains the samples we want to keep
  731. clust <- cutreeStatic(sampleTree, cutHeight = 100000, minSize = 5)
  732. clust
  733. keepSamples <- (clust != 0)
  734. notkeepSamples <- clust[clust == 0]
  735. length(notkeepSamples)
  736. datExpr <- data2[keepSamples, ]
  737. str(datExpr)
  738. keepSamples
  739. write.csv(datExpr, file = "datExpr.csv")
  741. ## Determination of the soft-thresholding power====================================================================================
  742. powers <- c(1:20) # often 20 in practice
  743. sft <- pickSoftThreshold(datExpr, powerVector = powers, RsquaredCut = 0.85)
  744. sft
  745. pdf("Soft_Threshold.pdf")
  746. par(mfrow = c(1, 2))
  747. plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  748. xlab = "Soft Threshold (power)", ylab = "SFT, signed R^2",
  749. type = "n", main = paste("Scale independence"))
  750. text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  751. labels = powers, col = "red")
  752. abline(h = 0.85, col = "red") ## 阈值线
  753. plot(sft$fitIndices[, 1], sft$fitIndices[, 5], type = "n",
  754. xlab = "Soft Threshold (power)", ylab = "Mean Connectivity",
  755. main = paste("Mean connectivity"))
  756. text(sft$fitIndices[, 1], sft$fitIndices[, 5],
  757. labels = powers, col = "red")
  758. dev.off()
  760. ## Obtain an appropriate soft-threshold
  761. power <- sft$powerEstimate
  762. if (is.na(power))
  763. {
  764. nSamples <- ncol(datExpr)
  765. power <- ifelse(nSamples < 20, 9, ifelse(nSamples < 30, 8, ifelse(nSamples < 40, 7, 6)))
  766. }
  767. power
  769. #Module identification====================================================================================
  770. module_detect_dir <- "C:/WGCNA"
  771. dir.create(file.path(module_detect_dir), showWarnings = FALSE)
  772. enableWGCNAThreads(nThreads = 2)
  773. str(datExpr)
  774. nGenes <- ncol(datExpr)
  775. block_thred <- 40000
  776. if(nGenes <= block_thred){
  777. net <- blockwiseModules(datExpr, corType = "pearson", maxBlockSize = 1.1 * nGenes,
  778. networkType = "unsigned", power = power, minModuleSize = 40,
  779. mergeCutHeight = 0.25, numericLabels = TRUE, saveTOMs = TRUE,
  780. pamRespectsDendro = FALSE, saveTOMFileBase = glue("{module_detect_dir}/TOM"))
  781. } else {
  782. net <- blockwiseModules(datExpr, corType = "pearson", maxBlockSize = block_thred,
  783. networkType = "unsigned", power = power, minModuleSize = 40,
  784. mergeCutHeight = 0.25, numericLabels = TRUE, saveTOMs = TRUE,
  785. pamRespectsDendro = FALSE, saveTOMFileBase = glue("{module_detect_dir}/TOM"))
  786. }
  788. moduleLables <- net$colors
  789. moduleColors <- labels2colors(moduleLables)
  790. moduleLables
  791. # #Extract the whole-genome module color table--------------------------------------
  792. str(moduleLables)
  793. head(moduleLables)
  794. genome_module<-cbind(moduleLables,moduleColors)
  795. genome_module<-as.data.frame(genome_module)
  796. factor_modules<-factor(genome_module$moduleColors)
  797. dummy_modules<-model.matrix(~factor(factor_modules)-1)
  798. str(dummy_modules)
  799. genome_module_dummy<-cbind(genome_module,dummy_modules)
  800. colnames(genome_module_dummy) <- gsub("factor\\(factor_modules\\)", "", colnames(genome_module_dummy))
  801. write.csv(genome_module_dummy, file = glue("{module_detect_dir}/genome_module_dummy.csv"), row.names = TRUE)
  804. ## save the results
  805. MEs <- net$MEs
  806. geneTree <- net$dendrograms[[1]]
  807. save(MEs, moduleLables, moduleColors, geneTree, net,
  808. file = glue("{module_detect_dir}/networkConstruction_and_modules.RData"))
  810. ## Plot the clustering tree and co-expression modules
  811. blocks <- length(unique(net$blocks))
  812. if(blocks == 1){
  813. pdf(glue("{module_detect_dir}/Gene_Cluster_Dendrogram.pdf"))
  814. plotDendroAndColors(net$dendrogram[[1]], moduleColors[net$blockGenes[[1]]],
  815. groupLabels = "Module colors", dendroLabels = FALSE,
  816. hang = 0.03, addGuide = TRUE, guideHang = 0.05)
  817. dev.off()
  818. } else {
  819. pdf(glue("{module_detect_dir}/Gene_Cluster_Dendrogram.pdf"))
  820. for (i in 1:blocks){
  821. plotDendroAndColors(net$dendrogram[[i]], moduleColors[net$blockGenes[[i]]],
  822. groupLabels = "Module colors",
  823. main = glue("Gene dendrogram and module colors in block {i}"),
  824. dendroLabels = FALSE,
  825. hang = 0.03, addGuide = TRUE, guideHang = 0.05)
  826. }
  827. dev.off()
  828. }
  831. # Co-expression network visualization
  832. dissTOM <- 1 - TOMsimilarityFromExpr(datExpr, power = power)
  833. set.seed(123)
  834. select <- sample(nGenes, size = as.numeric(400))
  835. selectTOM <- dissTOM[select, select]
  836. selectColors <- moduleColors[select]
  837. selectTree <- hclust(as.dist(selectTOM), method = "average")
  838. plotDiss = selectTOM^7
  839. diag(plotDiss) = NA
  840. pdf(glue("{outdir}/Network_heatmap_plot_selectGenes.pdf"))
  841. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes", col = heat.colors(12))
  842. dev.off()
  844. # compute module eigengenes
  845. MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
  846. MEs0
  847. MEsCol <- orderMEs(MEs0)
  848. str(MEsCol)
  849. write.csv(MEsCol,file=glue("{module_detect_dir}/mescol.csv"), row.names = T, quote = FALSE)
  850. MEsCol <- read.csv(glue("{module_detect_dir}/mescol.csv"), sep=",", row.names = 1)
  851. MEsCol
  853. #=## Phenotypic data association analysis - ====================================================================================
  854. ## Read phenotypic data
  855. data1 <- read.csv("brian regions-wgcna.csv", sep=",", row.names = 1)
  856. allTraits<-data1[,1:3]
  857. str(allTraits)
  858. Sample <- rownames(datExpr)
  859. str(Sample)
  860. traitRows <- match(Sample, rownames(allTraits))
  861. traitRows
  862. datTraits <- data.frame(allTraits[traitRows, ])
  863. table(rownames(datTraits) == rownames(datExpr))
  864. dev_stage<-factor(datTraits$dev_stage)
  865. coarse_structure<-factor(datTraits$coarse_structure)
  866. dummy_dev<-model.matrix(~factor(dev_stage)-1)
  867. dummy_structure<-model.matrix(~factor(coarse_structure)-1)
  868. str(dummy_dev)
  869. str(dummy_structure)
  871. # regression analysis for 41 MEs -----------------------------------------------------
  873. library(dplyr)
  874. str(MEsCol)
  875. head(dummy_structure)
  876. MEsColall<-cbind(MEsCol,dummy_structure,dummy_dev)
  877. write.csv(MEsColall,file=glue("{module_detect_dir}/MEsColall-21Var.csv"), row.names = T, quote = FALSE)
  878. PValues<-data.frame(matrix(nrow = 18, ncol = length(MEsCol)))
  879. BS<-data.frame(matrix(nrow = 18, ncol = length(MEsCol)))
  880. str(BS)
  881. colnames(MEsCol)
  882. for (i in 1:length(MEsCol))
  883. {
  884. formula <- as.formula(paste(colnames(MEsCol)[i], "~ dummy_structure+ dummy_dev"))
  885. model <- lm(formula, data = MEsColall)
  886. summary(model)
  887. # 进行其他处理，如汇总结果等
  888. PValue<-summary(model)$ coefficients[, "Pr(>|t|)"]
  889. B<-summary(model)$coefficients[,"Estimate"]
  890. PValues[,i]<-unname(PValue)
  891. BS[,i]<-unname(B)
  892. }
  893. str(PValues)
  894. head(PValues)
  895. colnames(PValues)<-colnames(MEsCol)
  896. colnames(BS)<-colnames(MEsCol)
  897. rownames(PValues)<-names(PValue)
  898. rownames(BS)<-names(PValue)
  899. names(PValue)
  900. head(BS)
  901. tBS<-t(BS)
  902. tPValues<-t(PValues)
  903. colnames(tPValues) <- gsub("dummy_devfactor\\(dev_stage\\)", "", colnames(tPValues))
  904. colnames(tPValues) <- gsub("dummy_structurefactor\\(coarse_structure\\)", "", colnames(tPValues))
  905. colnames(tBS) <- gsub("dummy_devfactor\\(dev_stage\\)", "", colnames(tBS))
  906. colnames(tBS) <- gsub("dummy_structurefactor\\(coarse_structure\\)", "", colnames(tBS))
  907. write.csv(tBS, file = glue("{module_detect_dir}/structure—dev-covariate-regression_B-18var.CSV"),row.names = T)
  908. write.csv(tPValues, file = glue("{module_detect_dir}/structure—dev-covariate-regression_PV-18var.csv"), row.names = T, quote = FALSE)
  911. ## Select the modules with most significant spatiotemporal enrichment (at least one regression coefficient greater than 0.05).
  912. library("WGCNA")
  913. library("glue")
  914. tBS<-read.csv(glue("{module_detect_dir}/structure—dev-covariate-regression_B-18var.CSV"), check.names = FALSE, sep=",", row.names = 1)
  915. tPValues<-read.csv(glue("{module_detect_dir}/structure—dev-covariate-regression_PV-18var.CSV"), check.names = FALSE,sep=",", row.names = 1)
  916. tBS<- as.matrix(tBS)
  917. tPValues<- as.matrix(tPValues)
  918. colnames(tBS)
  919. rownames(tBS)
  920. #(Intercept)
  921. tBS2<-tBS[,c("(Intercept)","Neocortex","Striatum","Ganglionic_eminence","Amygdala",
  922. "Hippocampus","Brain_stem","Cerebellum","Early prenatal", "Early mid-prenatal", "Late mid-prenatal","Late prenatal","Early infancy","Late infancy","Early childhood","Late childhood","Adolescence","Middle adulthood")]
  923. tPValues2<-tPValues[,c("(Intercept)","Neocortex","Striatum","Ganglionic_eminence","Amygdala",
  924. "Hippocampus","Brain_stem","Cerebellum","Early prenatal", "Early mid-prenatal", "Late mid-prenatal","Late prenatal","Early infancy","Late infancy","Early childhood","Late childhood","Adolescence","Middle adulthood")]
  925. colnames(tBS2)[1] <- "Tha | Young Adult"
  926. colnames(tPValues2)[1] <- "Tha | Young Adult"
  927. head(tPValues2)
  930. ##FDR ------------------------------------------------
  931. p_adj_by_disease <- apply(tPValues2, 2, p.adjust, method = "BH")
  933. p_adj_by_disease
  934. write.csv(p_adj_by_disease, file="structure—dev-covariate-regression_PV-18var-P-41M-FDR.CSV", quote = F,row.names = T)
  936. ##module with MAX B＞0.05 --------------------------------------------------
  938. threshold<-0.05
  939. threshold
  940. tBS2_sig <- tBS2[apply(tBS2, 1, max) > threshold, ]
  941. str(tBS2_sig)
  942. rownames(tBS2_sig)
  943. p_adj_fdr<-p_adj_by_disease[ rownames(tBS2_sig),]
  944. write.csv(tBS2_sig, file = glue("{module_detect_dir}/structure—dev-covariate-regression_B-18var_B大于0.05ME.CSV"),row.names = T)
  945. write.csv(p_adj_fdr, file = glue("{module_detect_dir}/structure—dev-covariate-regression_PV-18var_B大于0.05ME-fdr.csv"), row.names = T, quote = FALSE)
  947. ### textMatrix
  948. textMatrix <-signif(tBS2_sig, 1)
  950. textMatrix[abs(tBS2_sig) <= 0.05] <- ""
  951. textMatrix[p_adj_fdr >= 0.05] <- ""
  952. textMatrix
  953. dim(textMatrix) <- dim(tBS2_sig)
  955. # heatmap
  956. pdf(glue("{outdir}/Module-trait_relationship_dev_str_covari-18var-B0.05-fdr0.05.pdf"), width = 8, height =10)
  957. par(mar = c(9, 9, 3, 3),cex.main = 2)
  958. labeledHeatmap(Matrix = tBS2_sig, xLabels = colnames(tBS2_sig), yLabels = rownames(tBS2_sig),
  959. ySymbols = rownames(tBS2_sig), colorLabels = FALSE, colors = blueWhiteRed(50),
  960. cex.lab=1.0,
  961. textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.6, zlim = c(-0.25, 0.25),
  962. main = paste("Module-trait relationships"),
  963. clustRow = TRUE,
  964. rowDendro = TRUE,
  966. verticalSeparator.x=8)
  967. dev.off()
  970. # Chi-square test between module genes and the 15 merged risk gene sets -----------------------------------------
  971. #Create a loop to batch integrate risk gene data into genome module data
  972. library("data.table")
  973. riskgene<-read.table("15 disorder_merged gene.CSV",header=TRUE,sep=",",row.names=1)
  974. str(riskgene)
  975. braingene<-fread("genome_module_dummy.CSV",header=TRUE,sep=",")
  976. braingene<-braingene[,-c(2,3)]
  977. colnames(braingene)[1]<-"gene"
  978. colnames(braingene)[2:42]<-paste("ME",colnames(braingene)[2:42],sep="")
  979. colnames(braingene)
  980. str(braingene)
  981. colname=c(colnames(riskgene))
  982. str(riskgene)
  983. dim(riskgene)
  984. dim(braingene)
  985. braingeneXrisk<-braingene
  987. for(i in 1:15){
  988. v2<-colname[i]
  989. riskgenei<-data.table(c(riskgene[,i]),"1")
  990. names(riskgenei) <- c("gene", v2)
  991. str(riskgenei)
  992. braingeneXrisk<-dplyr::left_join(braingeneXrisk,riskgenei,by="gene")
  993. }
  994. braingeneXrisk[is.na(braingeneXrisk)]<-0
  995. str(braingeneXrisk)
  996. braingeneXrisk<-braingeneXrisk[!duplicated(braingeneXrisk,by="gene"),]
  997. str(braingeneXrisk)
  999. ##Chi-square test
  1000. library(gmodels)
  1001. library("data.table")
  1002. library("dplyr")
  1003. str(braingeneXrisk)
  1004. colnames(braingeneXrisk)
  1005. setDT(braingeneXrisk, keep.rownames=FALSE, key=NULL, check.names=FALSE)
  1006. variables<-c(colnames(braingeneXrisk))
  1007. variables
  1008. ncol(braingene)
  1009. ncol(braingeneXrisk)
  1010. p_values<-data.frame()
  1011. OR<-data.frame()
  1012. Risk_Positive_prop<-data.frame()
  1013. imax<-ncol(braingene)
  1014. imax
  1015. jmax<-ncol(braingeneXrisk)-ncol(braingene)
  1016. jmax
  1017. c(imax+1:jmax)
  1018. variables[2:ncol(braingene)]
  1019. variables[ncol(braingene)+1:jmax]
  1020. braingeneXrisk[[1]]
  1022. for (i in 2:imax)
  1023. {for (j in c(imax+1:jmax))
  1024. {result<-CrossTable(braingeneXrisk[[i]], braingeneXrisk[[j]],expected = F, format = "SAS",fisher = T, prop.c = F,prop.t = F,prop.chisq = F)
  1025. p_values[i-1,j-ncol(braingene)]<-result$fisher.gt$p.value
  1026. OR[i-1,j-ncol(braingene)]<-result$fisher.gt$estimate[[1]]
  1027. Risk_Positive_prop[i-1,j-ncol(braingene)]<-result$prop.col[2,2]
  1028. }
  1029. }
  1030. rownames(p_values)<-c(variables[2:imax])
  1031. colnames(p_values)<-c(variables[imax+1:jmax])
  1032. rownames(OR)<-c(variables[2:imax])
  1033. colnames(OR)<-c(variables[imax+1:jmax])
  1034. rownames(Risk_Positive_prop)<-c(variables[2:imax])
  1035. colnames(Risk_Positive_prop)<-c(variables[imax+1:jmax])
  1036. rownames(Risk_Positive_prop)
  1038. write.csv(p_values, file="P_Values-41ME.CSV", quote = F,row.names = T)
  1039. write.csv(OR,file="OR-41ME.CSV", quote=F,row.names = T)
  1044. ##Conduct enrichment analysis with the enrich function in the clusterProfiler package
  1045. library(dplyr)
  1046. library(stringr)
  1047. library(rlang)
  1048. library(ggplot2)
  1049. require(clusterProfiler)
  1050. library(clusterProfiler)
  1051. library(biomaRt)
  1052. library(AnnotationHub)
  1053. library(AnnotationDbi)
  1054. library(ggstar)
  1055. library("org.Hs.eg.db")
  1056. library("ggtext")
  1058. # enriChKEGG---------------
  1059. # KEGG analysis of significantly enriched modules---------------------------------------------------------------
  1060. sig_module<-read.table("sig modules-16.CSV",header=TRUE,sep=",",row.names=1)
  1061. colnames(sig_module)
  1062. posi_gene_Entrez<-vector("list", ncol(sig_module))
  1063. names(posi_gene_Entrez)<-colnames(sig_module)
  1064. colnames(sig_module)
  1066. for (j in 1:16)
  1067. {
  1068. posi_sig_module<-sig_module[,j]
  1069. length(posi_sig_module)
  1070. posi_sig_module <- posi_sig_module[!posi_sig_module==""]
  1071. length(posi_sig_module)
  1072. posi_sig_module
  1073. entrezid_all = mapIds(x = org.Hs.eg.db,
  1074. keys = posi_sig_module,
  1075. keytype = "SYMBOL",
  1076. column = "ENTREZID")
  1077. entrezid_all
  1078. entrezid_all = na.omit(entrezid_all)
  1079. entrezid_all = data.frame(entrezid_all)
  1080. posi_gene_Entrez[[j]]<- entrezid_all[,1]
  1081. }
  1082. summary(posi_gene_Entrez)
  1084. save(posi_gene_Entrez, file ='KEGG_sig module gene.RData')
  1086. load('KEGG_sig module gene.RData')
  1087. summary(posi_gene_Entrez)
  1088. ###qvalueCutoff=0.05###
  1089. ###Compare multiple enrichment results using the compareCluster function in clusterProfiler
  1090. KEGG_result<-compareCluster(posi_gene_Entrez, fun="enrichKEGG", organism="hsa",qvalueCutoff=0.05)
  1091. table(KEGG_result@compareClusterResult$Cluster)
  1092. str(KEGG_result)
  1093. KEGG_result@compareClusterResult$Description<-str_replace_all(KEGG_result@compareClusterResult$Description, " - Homo sapiens \\(.*?\\)", "")
  1094. str(KEGG_result@compareClusterResult$Description)
  1095. KEGG_result@compareClusterResult$Log_qvalue<- -log(KEGG_result@compareClusterResult$qvalue,10)
  1096. write.csv(KEGG_result,"KEGG_sig module gene-q0.05.csv")
  1097. dotpl<- dotplot(KEGG_result, x="Cluster",color= "Log_qvalue",
  1098. showCategory =3,
  1099. by="Log_qvalue",
  1100. size ="GeneRatio",
  1101. font.size = 10,
  1102. label_format = 30)+
  1103. theme(axis.text.x = element_text(angle = 45, hjust = 1,size=8),axis.text.y = element_text(hjust = 1,size=8),plot.title = element_text(size = 10, hjust = 0.5))+
  1104. labs(title="Enrichment analysis of module genes for KEGG pathways",x = '', y = '')
  1105. ggsave(dotpl,filename="KEGG_sig module gene-q0.05.png", width =6, height =7,units = "in", dpi = 300, bg = "transparent")
  1108. ##enrichGO--export go terms top 3significant
  1109. for (j in 1:16)
  1110. {
  1111. summary(posi_gene_Entrez)
  1112. enrichGO_result<-enrichGO(posi_gene_Entrez[[j]], OrgDb=org.Hs.eg.db, keyType = "ENTREZID", ont = "ALL",
  1113. pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.05,
  1114. minGSSize =3, maxGSSize = 500,
  1115. readable = FALSE,
  1116. pool = FALSE
  1117. )
  1118. str(enrichGO_result)
  1120. go <- enrichGO_result
  1121. str(go)
  1122. go <- go[order(go$ONTOLOGY, go$qvalue, decreasing = F),]
  1123. go_top5 <- go %>%
  1124. group_by(ONTOLOGY) %>%
  1125. slice_head(n = 3) %>%
  1126. ungroup()
  1127. go_top5 <- go_top5[order(go_top5$ONTOLOGY, go_top5$qvalue, decreasing = T),]
  1128. go_top5
  1130. go_top5$Description <- sapply(go_top5$Description, function(x) paste(strwrap(x, width = 30), collapse = "\n"))
  1131. str(go_top5)
  1132. go_top5$ShortDesc <- ifelse(nchar(go_top5$Description) > 40,
  1133. paste0(substr(go_top5$Description, 1, 37), "..."),
  1134. go_top5$Description)
  1135. go_top5$ShortDesc
  1136. go_top5$ShortDesc <- sub("acetylcholine", "Ach", go_top5$ShortDesc) # 提取首个词
  1137. barpl<-ggplot(go_top5, aes(x = reorder(ShortDesc, -qvalue),y= -log10(qvalue))) +
  1138. geom_col(aes(fill = ONTOLOGY), width = 0.5, show.legend =F) +
  1139. facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y') +
  1140. theme(axis.text = element_text(size = 12), panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
  1141. scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  1142. coord_flip() +
  1143. labs(title=colnames(sig_module)[j],x = '', y = '-Log10(q)' )
  1144. ggsave(barpl,filename=paste0('enrich.go_sig module gene',colnames(sig_module)[j],'-q0.05.png'), width = 4,height=4,dpi=300, units="in", bg = "transparent")
  1145. }

R code for data analysis.R, under CC-BY-4.0 · at the source

Overview

Authors: Weiqing Liu1, Tomomi Shimogori2
ORCID iDs: Weiqing Liu, Tomomi Shimogori
  1. Clinical Research Center for Mental Disorders, Shanghai Pudong New Area Mental Health Center, School of Medicine, Tongji University,Shanghai, China
  2. Laboratory for Molecular Mechanisms of Brain Development, RIKEN Center for Brain Science,Wako, Japan
Institutions: Tongji University (China); RIKEN Center for Brain Science (Japan)
Journal: Communications biology, volume 9, issue 1, article 634
Dates: received 23 September 2025; accepted 31 March 2026; published online 13 May 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s42003-026-10045-x · PMID 42129508 · PMCID PMC13172025 · OpenAlex W4414931447
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.21203/rs.3.rs-7667905/v1
Status: code verified
Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Connectivity, Spectral & time-frequency
Keywords: Development, Disease genetics
MeSH: Brain*, Genetic Predisposition to Disease*, Mental Disorders*, Transcriptome*, Gene Expression Profiling, Genome-Wide Association Study, Humans, Neurodevelopment (* major topic)
Topic: Health, Environment, Cognitive Aging (Health, Toxicology and Mutagenesis, Environmental Science), according to OpenAlex
Funding: Key Specialty Construction Project of Shanghai Pudong Health Commission (PWZzk2022-19), Shanghai Pudong New Area Health Committee Disciplinary Leader Program (PWRd2021-06), Shanghai Pudong Science & Technology Development Foundation (PKJ2025-Y24),Open Project of the Key Laboratory of Neuropsychiatric Disorders and Developmental Genetics, Ministry of Education (2025GDND02), and the New Quality Clinical Specialty Program of High-end Medical Disciplinary Construction in Shanghai Pudong New Area (2025-PWXZ-08),JSPS KAKENHI (18F18736); National Institutes of Health grants (U01MH124619), the Funding Program by MEXT Grants-in-Aid for Scientific Research on Innovative Areas “Cell type census of adaptive neuronal circuits: biological mechanisms of structural and functional organization (Adaptive Circuit Census, ACC)” (KAKENHI 21H05240), and the Brain Mapping by Integrated Neurotechnologies for Disease Studies (Brain/MINDS) project from the Japan Agency for Medical Research and Development (JP20dm0207001) and Brain/MINDS2.0 (JP23wm0625001), and JSPS KAKENHI (18F18736)
Citations: cited by 1 paper (Europe PMC); 70 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

liuwq2010/Neuropsychiatric-risk-gene-spatiotemporal-expression

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: aaec55250e394a7553d4aad991fa7f86cc0f0653, 21 September 2025
Languages: R (1)
Size: 2 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: license file
Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: clusterProfiler (1 file), ComplexHeatmap (1 file), data.table (1 file), easystats (1 file), ggplot2 (1 file), ggpubr (1 file), tidyverse (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
2 files
At the source: github.com/liuwq2010/Neuropsychiatric-risk-gene-spatiotemporal-expression

Zenodo 19097352

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Languages: R (1)
Size: 1 file, 1 script
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: clusterProfiler (1 file), ComplexHeatmap (1 file), data.table (1 file), easystats (1 file), ggplot2 (1 file), ggpubr (1 file), rstatix (1 file), tidyverse (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
1 file
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s42003-026-10045-x.

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;
  • 2 scripts, each with its path and the digest of its content;
  • 15 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

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s42003-026-10045-x.

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, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 2 keywords, 8 MeSH terms, 2 funders, 67 references.

Cite

This paper

Liu, W., & Shimogori, T. (2026). Spatiotemporal brain transcriptomics reveal risk gene hot-spots in major neuropsychiatric disorders. Communications biology, 9(1), 634. https://doi.org/10.1038/s42003-026-10045-x

BibTeX

@article{liu2026spatiotemporal,
author = {Liu, Weiqing and Shimogori, Tomomi},
title = {{Spatiotemporal brain transcriptomics reveal risk gene hot-spots in major neuropsychiatric disorders}},
journal = {Communications biology},
year = {2026},
month = may,
volume = {9},
number = {1},
pages = {634},
publisher = {Nature Publishing Group},
issn = {2399-3642},
doi = {10.1038/s42003-026-10045-x},
url = {https://doi.org/10.1038/s42003-026-10045-x},
pmid = {42129508},
pmcid = {PMC13172025}
}

RIS

TY - JOUR
AU - Liu, Weiqing
AU - Shimogori, Tomomi
TI - Spatiotemporal brain transcriptomics reveal risk gene hot-spots in major neuropsychiatric disorders
T2 - Communications biology
J2 - Commun Biol
PY - 2026
DA - 2026/05/13
VL - 9
IS - 1
SP - 634
SN - 2399-3642
PB - Nature Publishing Group
DO - 10.1038/s42003-026-10045-x
UR - https://doi.org/10.1038/s42003-026-10045-x
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42003-026-10045-x",
"type": "article-journal",
"title": "Spatiotemporal brain transcriptomics reveal risk gene hot-spots in major neuropsychiatric disorders",
"container-title": "Communications biology",
"author": [
{
"family": "Liu",
"given": "Weiqing"
},
{
"family": "Shimogori",
"given": "Tomomi"
}
],
"container-title-short": "Commun Biol",
"volume": "9",
"issue": "1",
"page": "634",
"DOI": "10.1038/s42003-026-10045-x",
"PMID": "42129508",
"PMCID": "PMC13172025",
"ISSN": "2399-3642",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s42003-026-10045-x",
"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.1371/journal.pcbi.1014422 [code]
Deciphering cell type-specific causal genetic effects on brain imaging-derived phenotypes and disorders with single-cell Mendelian randomization.
Journal: PLoS computational biology
In common: clusterProfiler, data.table, ggplot2, 1 other tool, genetics / omics, cellular / molecular, 10 references
[2] doi:10.1038/s41593-026-02247-7 [code]
Transcriptomic and phenotypic convergence of neurodevelopmental disorder risk genes in vitro and in vivo.
Journal: Nature neuroscience
In common: genetics / omics, cellular / molecular, 13 references
[3] doi:10.1038/s41380-026-03497-4 [code]
Transcriptome-informed brain cartography of polygenic risk and association with brain structure in major psychiatric disorders.
Journal: Molecular psychiatry
In common: clusterProfiler, ComplexHeatmap, ggplot2, 1 other tool, genetics / omics, cellular / molecular, 8 references
[4] doi:10.1016/j.celrep.2026.117500 [code]
Spatio-molecular gene expression reflects dorsal anterior cingulate cortex structure and function in the human brain.
Journal: Cell reports
In common: easystats, ComplexHeatmap, ggpubr, 2 other tools, genetics / omics, cellular / molecular, 7 references
[5] doi:10.1111/adb.70179 [code]
Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala.
Journal: Addiction biology
In common: rstatix, clusterProfiler, ComplexHeatmap, 3 other tools, genetics / omics, cellular / molecular, 3 references
[6] 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: WGCNA, data.table, ggplot2, 1 other tool, genetics / omics, cellular / molecular, 5 references
[7] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: WGCNA, clusterProfiler, ComplexHeatmap, 4 other tools, 2 references
[8] doi:10.1038/s41467-026-75193-4 [code]
Multi-ancestry gene expression models amplify transcriptome-wide association study discovery and validation.
Journal: Nature communications
In common: ComplexHeatmap, ggpubr, data.table, 2 other tools, genetics / omics, cellular / molecular, 4 references
[9] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: WGCNA, rstatix, clusterProfiler, 5 other tools, genetics / omics
[10] doi:10.1016/j.xhgg.2026.100652 [code]
CRISPR-engineered deletion of POGZ alters transcription factor binding at promoters of genes involved in synaptic signaling.
Journal: HGG advances
In common: WGCNA, ggpubr, data.table, 2 other tools, cellular / molecular, 4 references

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.