OSCR

Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators.

Code ↔ Paper

7 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 7 matches
  1. [1] § Results › Transcriptomic profiling reveals diverse gene targeting effects on microglial state ↔ image_analysis/Crispr30_publication_final.Rmd, lines 393–469 · score 0.86 · HLA DMB, SLC11A1, FCGR3B, NCKAP1L, ARHGAP25, BIN2
  2. [2] § Results › Transcriptomic profiling reveals diverse gene targeting effects on microglial state ↔ drug_seq/generate_string_plot_final.Rmd, lines 32–48 · score 0.77 · HLA DMB, SLC11A1, NCKAP1L, ARHGAP25, BIN2, TLR6
  3. [3] § Results › Functional programs induced by targeted genes ↔ drug_seq/03_01_analyze_deg_go.R, lines 441–500 · score 0.66 · cell adhesion, complement activation, protein production, cytokine, GO, chemokine
  4. [4] § Results › Functional programs induced by targeted genes ↔ image_analysis/Crispr30_publication_final.Rmd, lines 393–469 · score 0.59 · CLEC7A, FCGR2A, CD14, IRF8, SYK, TREM2
  5. [5] § Results › Functional programs induced by targeted genes ↔ drug_seq/generate_string_plot_final.Rmd, lines 32–48 · score 0.58 · CLEC7A, FCGR2A, CD14, IRF8, SYK, TREM2
  6. [6] § Results › Transcriptomic profiling reveals diverse gene targeting effects on microglial state ↔ drug_seq/03_01_analyze_deg_go.R, lines 113–229 · score 0.52 · log2 fold change, log2FC, Asterisks, Pairwise, pvalue, DEGs
  7. [7] § Results › Functional programs induced by targeted genes ↔ drug_seq/03_01_analyze_deg_go.R, lines 113–229 · score 0.51 · log2 fold change, log2FC, Asterisks, GO, heatmaps, DEG

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,454 lines · 49 KB · no license · 3 matches

  1. LIB='/cluster/tufts/patralab/rbator01/R_libs/4.4.0'
  2. .libPaths(c("",LIB))
  4. library(pheatmap)
  5. library(tidyverse)
  6. library(ComplexHeatmap)
  7. library(openxlsx)
  8. library(ggforce)
  9. library(data.table)
  10. library(clusterProfiler)
  11. library('org.Hs.eg.db')
  12. library(data.table)
  13. library(circlize)
  14. library(grid)
  15. library(forcats)
  16. library(stringr)
  17. library(janitor)
  18. library(DESeq2)
  19. library(ggbreak)
  21. select = dplyr::select
  22. rename = dplyr::rename
  23. filter=dplyr::filter
  25. # paths and setup ----
  26. setwd('/cluster/tufts/patralab/rbator01/perlis_lab/cspr_apr25/')
  28. # DEG thresholds
  29. padj_thresh = 0.1
  30. praw_thresh = 0.1
  32. # load metadata and counts ----
  33. meta = read.xlsx("analysis/metadata_batch1/241210_Joy30_Metadata_format.xlsx")
  34. ctrl_counts = read.xlsx("analysis/nct1_expression/nct1_cpm_ko_genes.xlsx")
  36. rm_samples = fread("analysis/de/removed_samples_control_batch_q1q3_filtered_batch1_batch2.csv") %>%
  37. mutate(tmp = name) %>%
  38. separate(tmp, into=c("ko_gene","well","batch"), sep="\\|")
  40. n_rm = rm_samples %>%
  41. filter(batch == "b1") %>%
  42. group_by(ko_gene) %>%
  43. mutate(n_rm = n()) %>%
  44. ungroup() %>%
  45. as.data.frame()
  47. # deg results b1 -----
  48. b1 = fread("analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv")
  49. head(b1)
  51. res_sig_b1 = b1 %>%
  52. filter(pvalue < praw_thresh)
  54. # ko_gene in gene_name at praw < praw_thresh
  55. ko_in_praw_df_b1 <- res_sig_b1 %>%
  56. group_by(ko_gene) %>%
  57. summarise(ko_in_praw_deg_b1 = as.integer(ko_gene %in% gene_name), .groups = "drop")%>%
  58. distinct()
  60. # ko_gene in gene_name at padj < padj_thresh
  61. ko_in_padj_df_b1 <- res_sig_b1 %>%
  62. filter(padj < padj_thresh) %>%
  63. group_by(ko_gene) %>%
  64. summarise(ko_in_padj_deg_b1 = as.integer(ko_gene %in% gene_name), .groups = "drop")%>%
  65. distinct()
  67. # count DEGs
  68. t_b1 = as.data.frame(table((b1 %>% filter(padj < padj_thresh))$ko_gene))
  69. colnames(t_b1) = c("ko_gene","ndeg_padj_0.1_b1")
  71. top_b1 <- b1 %>%
  72. group_by(ko_gene) %>%
  73. arrange(padj) %>%
  74. slice_min(order_by = padj, n = 10, with_ties = FALSE) %>%
  75. select(ko_gene, gene_name) %>%
  76. rename(top_b1 = gene_name)
  78. top_b1_collapsed <- aggregate(top_b1 ~ ko_gene, data = top_b1, FUN = function(x) paste(x, collapse = ","))
  80. # assemble summary workbook ---------
  82. final = ctrl_counts %>%
  83. rename(ko_gene=gene_name) %>%
  84. mutate(coverage_pass = ifelse(mean_cpm > 5, 1,0)) %>%
  85. full_join(n_rm %>% dplyr::select(ko_gene, n_rm), by="ko_gene") %>%
  86. left_join(t_b1, by = "ko_gene") %>%
  87. left_join(ko_in_padj_df_b1, by = "ko_gene") %>%
  88. left_join(ko_in_praw_df_b1, by = "ko_gene") %>%
  89. left_join(top_b1_collapsed, by = "ko_gene")
  92. # Replace NA with 0
  93. final <- final %>%
  94. mutate(across(starts_with("ko_in_padj_deg_"), ~ if_else(is.na(.), 0L, .))) %>%
  95. mutate(across(starts_with("ko_in_praw_deg_"), ~ if_else(is.na(.), 0L, .))) %>%
  96. mutate(across(starts_with("ndeg_padj_0.1_"), ~ if_else(is.na(.), 0L, as.integer(.))))
  98. final = final %>%
  99. distinct()
  100. write.xlsx(final, "analysis/de/summary_batch1_q1q3_filter_rmlow_22aug25.xlsx")
  102. # Omics Fig 1A: DEG target overlap, batch 1 only =======
  103. dds = readRDS("analysis/processed/umi.trim.1mm_all_dedup.counts.dds.rds")
  104. genes_rm = c('RBFOX3','FCGR3B','PLD4','FCAR','NTC1')
  106. meta = colData(dds) %>%
  107. as.data.frame() %>%
  108. filter(!(KO.Target %in% genes_rm))
  110. all_targets = unique(meta$KO.Target)
  111. all_targets = factor(all_targets, levels = all_targets)
  113. ## Helper: build KO–KO log2FC heatmap with asterisks for significant entries
  114. make_heat <- function(deg_file,
  115. title_text,
  116. out_pdf,
  117. lfc_type = c("raw", "shrunken"),
  118. sig_col = c("padj", "pvalue"),
  119. sig_thresh = 0.1) {
  121. lfc_type <- match.arg(lfc_type) # "raw" or "shrunken"
  122. sig_col <- match.arg(sig_col) # "padj" or "pvalue"
  124. # Append info to title for clarity
  125. #title_text <- paste0(title_text, " ", lfc_type, " (", sig_col, " < ", sig_thresh, ")")
  127. # read & standardize column names
  128. deg <- data.table::fread(deg_file) %>%
  129. dplyr::rename(gene = gene_name)
  131. # ensure chosen significance column is present
  132. if (!sig_col %in% colnames(deg)) {
  133. stop(sprintf("Column '%s' not found in '%s'.", sig_col, deg_file))
  134. }
  136. # choose LFC column
  137. lfc_col <- if (lfc_type == "raw") "log2FoldChange_raw" else "log2FoldChange"
  139. # keep targets, aggregate duplicates:
  140. # - LFC columns: mean
  141. # - star_flag: TRUE if any chosen sig_col < sig_thresh
  142. agg <- deg %>%
  143. dplyr::select(gene, ko_gene, log2FoldChange_raw, log2FoldChange, pvalue, padj) %>%
  144. dplyr::filter(gene %in% all_targets, ko_gene %in% all_targets) %>%
  145. dplyr::group_by(gene, ko_gene) %>%
  146. dplyr::summarise(
  147. log2FoldChange_raw = mean(log2FoldChange_raw, na.rm = TRUE),
  148. log2FoldChange = mean(log2FoldChange, na.rm = TRUE),
  149. star_flag = any(.data[[sig_col]] < sig_thresh, na.rm = TRUE),
  150. .groups = "drop"
  151. ) %>%
  152. dplyr::mutate(dplyr::across(c(log2FoldChange_raw, log2FoldChange),
  153. ~ ifelse(is.nan(.x), NA_real_, .x)))
  155. # complete to full target grid; impute chosen LFC to 0
  156. heat_long <- agg %>%
  157. tidyr::complete(
  158. gene = all_targets,
  159. ko_gene = all_targets
  160. ) %>%
  161. dplyr::mutate(
  162. value = .data[[lfc_col]],
  163. value = tidyr::replace_na(value, 0),
  164. star_flag = tidyr::replace_na(star_flag, FALSE),
  165. gene = factor(gene, levels = all_targets),
  166. ko_gene = factor(ko_gene, levels = all_targets)
  167. )
  169. # value matrix
  170. heat_mat <- heat_long %>%
  171. dplyr::select(gene, ko_gene, value) %>%
  172. tidyr::pivot_wider(names_from = ko_gene, values_from = value) %>%
  173. tibble::column_to_rownames("gene") %>%
  174. as.matrix()
  176. # asterisk (logical) matrix
  177. star_mat <- heat_long %>%
  178. dplyr::select(gene, ko_gene, star_flag) %>%
  179. tidyr::pivot_wider(names_from = ko_gene, values_from = star_flag) %>%
  180. tibble::column_to_rownames("gene") %>%
  181. as.matrix()
  183. ## cluster columns (correlation distance)
  184. d_col <- as.dist(1 - stats::cor(heat_mat, use = "pairwise.complete.obs"))
  185. col_clust <- stats::hclust(d_col, method = "average")
  186. col_order_names <- colnames(heat_mat)[col_clust$order]
  188. # row order aligned to clustered columns (append non-overlapping rows)
  189. common <- intersect(col_order_names, rownames(heat_mat))
  190. first_part <- match(common, rownames(heat_mat))
  191. row_order <- c(first_part, setdiff(seq_len(nrow(heat_mat)), first_part))
  193. row_gp <- grid::gpar(fontfamily = "Arial", fontsize = 8)
  194. col_gp <- grid::gpar(fontfamily = "Arial", fontsize = 8)
  195. title_gp <- grid::gpar(fontfamily = "Arial", fontsize = 12, fontface = "bold")
  196. legend_gp <- grid::gpar(fontfamily = "Arial")
  198. ht <- ComplexHeatmap::Heatmap(
  199. heat_mat,
  200. name = "log2FC", #ifelse(lfc_type == "raw", "log2FC_raw", "log2FC_shrunken"),
  201. col = circlize::colorRamp2(
  202. c(min(heat_mat, na.rm = TRUE), 0, max(heat_mat, na.rm = TRUE)),
  203. c("blue", "white", "red")
  204. ),
  205. cluster_columns = as.dendrogram(col_clust),
  206. cluster_rows = FALSE,
  207. row_order = row_order,
  208. show_row_names = TRUE,
  209. show_column_names = TRUE,
  210. row_names_gp = row_gp,
  211. column_names_gp = col_gp,
  212. column_title_gp = title_gp,
  213. heatmap_legend_param = list(
  214. labels_gp = legend_gp,
  215. title_gp = grid::gpar(fontfamily = "Arial", fontface = "bold")
  216. ),
  217. cell_fun = function(j, i, x, y, width, height, fill) {
  218. if (!is.na(star_mat[i, j]) && star_mat[i, j]) {
  219. grid::grid.text("*", x = x, y = y,
  220. gp = grid::gpar(fontfamily = "Arial", fontsize = 10, fontface = "bold"))
  221. }
  222. }
  223. )
  225. # Use Cairo to embed Arial cleanly in the PDF
  226. Cairo::CairoPDF(out_pdf, height = 5, width = 6, family = "Arial")
  227. on.exit(dev.off(), add = TRUE)
  228. ComplexHeatmap::draw(ht)
  229. }
  231. # KO–KO overlap using raw LFC and p-value < 0.1
  233. make_heat(
  234. deg_file = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  235. title_text = "",
  236. out_pdf = "analysis/de/plots/deg_target_overlap_b1_lfcraw_pval_0.1.pdf",
  237. lfc_type = "raw",
  238. sig_col = "pvalue",
  239. sig_thresh = 0.1 )
  242. make_heat(
  243. deg_file = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  244. title_text = "",
  245. out_pdf = "analysis/de/plots/deg_target_overlap_b1_lfcraw_padj_0.1.pdf",
  246. lfc_type = "raw",
  247. sig_col = "padj",
  248. sig_thresh = 0.1 )
  251. # Barplot NUMBER DEG =======
  253. tab <- b1 %>%
  254. group_by(ko_gene) %>%
  255. summarise(n = sum(padj < 0.1, na.rm = TRUE), .groups = "drop")
  257. # Color/category per KO from the *diagonal* row (ko_gene == gene_name)
  258. diag_map <- b1 %>%
  259. filter(ko_gene == gene_name) %>%
  260. transmute(
  261. ko_gene,
  262. p_raw = pvalue,
  263. lfc = dplyr::coalesce(log2FoldChange, log2FoldChange_raw),
  264. direction = case_when(
  265. is.finite(p_raw) & p_raw < 0.1 & is.finite(lfc) & lfc > 0 ~ "sig_up",
  266. is.finite(p_raw) & p_raw < 0.1 & is.finite(lfc) & lfc < 0 ~ "sig_down",
  267. TRUE ~ "other"
  268. )
  269. )
  271. tab <- tab %>%
  272. left_join(diag_map[, c("ko_gene", "direction")], by = "ko_gene") %>%
  273. mutate(direction = tidyr::replace_na(direction, "other"),
  274. direction = factor(direction, levels = c("sig_up", "sig_down", "other")))
  277. # Find outlier and compute a sensible break
  278. top_row <- tab %>% slice_max(n, n = 1, with_ties = FALSE)
  279. top_val <- top_row$n
  280. max_others <- max(tab$n[tab$ko_gene != top_row$ko_gene], na.rm = TRUE)
  282. gap <- max(5, floor(0.05 * (top_val - max_others))) # ~5% gap or ≥5
  283. lower_brk <- max_others + gap
  284. upper_brk <- top_val - gap
  286. # Axis ticks: show normal ticks below break + a single tick at 1100 above break
  287. lower_ticks <- pretty(c(0, max_others))
  288. upper_tick <- 1100
  289. all_breaks <- unique(c(lower_ticks, upper_tick))
  291. p <- ggplot(tab, aes(x = fct_reorder(ko_gene, n), y = n, fill = direction)) +
  292. geom_col(color = "black") +
  293. coord_flip(clip = "off") +
  294. labs(
  295. x = "Target Gene",
  296. y = "n (DEGs with padj < 0.1)",
  297. title = "",
  298. fill = "DEG status of Target Gene"
  299. ) +
  300. scale_fill_manual(
  301. values = c(sig_up = "red", sig_down = "blue", other = "black"),
  302. labels = c(
  303. sig_up = "p < 0.1 & log2FC > 0",
  304. sig_down = "p < 0.1 & log2FC < 0",
  305. other = "otherwise"
  306. )
  307. ) +
  308. # Broken axis
  309. ggbreak::scale_y_break(c(lower_brk, upper_brk)) +
  310. # Only label lower ticks + one tick (1100) in the upper segment
  311. scale_y_continuous(
  312. breaks = all_breaks,
  313. labels = function(b) ifelse(b == upper_tick, "1100", b)
  314. ) +
  315. theme_minimal(base_size = 12, base_family = "Arial") +
  316. theme(
  317. text = element_text(family = "Arial"),
  318. plot.title = element_text(face = "bold"),
  319. legend.position = "right",
  320. plot.margin = margin(5.5, 20, 5.5, 5.5)
  321. ) +
  322. # ensure the 1100 tick is in range even if top_val < 1100
  323. expand_limits(y = max(top_val, upper_tick) * 1.02)
  325. systemfonts::match_font("Arial")
  327. show(p)
  328. ggsave("analysis/de/plots/deg_per_ko_gene.pdf", p,width = 7.5, height = 5, device = cairo_pdf, family = "Arial",onefile = FALSE)
  330. # Same barplot but all blue (no categories in legend)
  331. p <- ggplot(tab, aes(x = fct_reorder(ko_gene, n), y = n)) +
  332. geom_col(color = "black", fill = "blue") +
  333. coord_flip(clip = "off") +
  334. labs(
  335. x = "Target Gene",
  336. y = "n (DEGs with padj < 0.1)",
  337. title = ""
  338. ) +
  339. ggbreak::scale_y_break(c(lower_brk, upper_brk)) +
  340. scale_y_continuous(
  341. breaks = all_breaks,
  342. labels = function(b) ifelse(b == upper_tick, "1100", b)
  343. ) +
  344. theme_minimal(base_size = 12, base_family = "Arial") +
  345. theme(
  346. text = element_text(family = "Arial"),
  347. plot.title = element_text(face = "bold"),
  348. legend.position = "none", # hide legend
  349. plot.margin = margin(5.5, 20, 5.5, 5.5)
  350. ) +
  351. expand_limits(y = max(top_val, upper_tick) * 1.02)
  353. print(p)
  354. ggsave("analysis/de/plots/deg_per_ko_gene_blue.pdf", p,
  355. width = 5, height = 5, device = cairo_pdf, family = "Arial", onefile = FALSE)
  358. # GO ENRICHMENT for targets that have >10 DEG =================
  360. b1_filter = fread("analysis/de/all_res_vs_ntc1_padj_0.1_10deg_batch1_q1q3_filter_rmlow.tsv")
  363. ck<- compareCluster(geneCluster = ens ~ ko_gene,
  364. data = b1_filter,
  365. OrgDb = org.Hs.eg.db,
  366. keyType="ENSEMBL",
  367. fun = "enrichGO",
  368. ont="BP",
  369. readable=T)
  371. saveRDS(ck, "analysis/de/go_bp_batch1_padj_0.1_10deg.rds")
  372. write.xlsx(ck, "analysis/de/go_bp_batch1_padj_0.1_10deg.xlsx")
  374. ck<- compareCluster(geneCluster = ens ~ ko_gene,
  375. data = b1_filter,
  376. OrgDb = org.Hs.eg.db,
  377. keyType="ENSEMBL",
  378. fun = "enrichGO",
  379. ont="CC",
  380. readable=T)
  382. saveRDS(ck, "analysis/de/go_cc_batch1_padj_0.1_10deg.rds")
  383. write.xlsx(ck, "analysis/de/go_cc_batch1_padj_0.1_10deg.xlsx")
  385. ck<- compareCluster(geneCluster = ens ~ ko_gene,
  386. data = b1_filter,
  387. OrgDb = org.Hs.eg.db,
  388. keyType="ENSEMBL",
  389. fun = "enrichGO",
  390. ont="MF",
  391. readable=T)
  393. saveRDS(ck, "analysis/de/go_mf_batch1_padj_0.1_10deg.rds")
  394. write.xlsx(ck, "analysis/de/go_mf_batch1_padj_0.1_10deg.xlsx")
  396. ## Quick sanity plots for GO enrichments (BP / MF / CC)
  398. ck = readRDS("analysis/de/go_bp_batch1_padj_0.1.rds")
  399. dotplot(ck, show=3)+
  400. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  401. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  403. ck = readRDS("analysis/de/go_mf_batch1_padj_0.1.rds")
  404. dotplot(ck, show=3)+
  405. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  406. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  408. ck = readRDS("analysis/de/go_cc_batch1_padj_0.1.rds")
  409. dotplot(ck, show=3)+
  410. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  411. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  413. ## GO BP: filter by keyword families ----
  415. keywords <- c(
  416. "phago", "autophag",
  417. "neuro",
  418. "macrophage", "activation",
  419. "adhesion", "attachment", "integrin", "actin",
  420. "chemotaxis", "chemokine", "cytokine",
  421. "inflamm", "stress",
  422. "lyso",
  423. "metabolic", "oxydat",
  424. "mitochondria"
  425. )
  427. ck = readRDS("analysis/de/go_bp_batch1_padj_0.1_10deg.rds")
  429. # Collapse keywords into regex pattern
  430. pattern <- paste(keywords, collapse = "|")
  432. # Filter the dataframe based on keyword match in the specified column
  433. ck_filter = ck %>%
  434. filter(grepl(pattern, ck@compareClusterResult$Description, ignore.case = TRUE))
  436. dotplot(ck_filter, show=3)+
  437. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  438. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  441. # Select terms from collaborators for initial combination -----
  442. #BP
  443. ck = readRDS("analysis/de/go_bp_batch1_padj_0.1_10deg.rds")
  445. s = read.xlsx("analysis/de/go_bp_batch1_padj_0.1_10deg_sds.xlsx") %>%
  446. mutate(select = ifelse(select == 1, 's',select)) %>%
  447. filter(select == 's') %>%
  448. select(ID, select)
  450. l = read.xlsx("analysis/de/go_bp_batch1_padj_0.1_10deg_LM.xlsx") %>%
  451. filter(!is.na(LM)) %>%
  452. rename(select = LM) %>%
  453. mutate(select = 'l') %>%
  454. select(ID, select)
  456. c = read.xlsx("analysis/de/go_bp_batch1_padj_0.1_10deg.xlsx") %>%
  457. filter(!is.na(CAMILLA_list)) %>%
  458. rename(select = CAMILLA_list ) %>%
  459. select(ID, select)
  462. all_select = rbind(s,l,c) %>%
  463. distinct()
  465. dedup_select <- all_select %>%
  466. mutate(
  467. is_num = grepl("^\\d+$", select),
  468. # convert only the numeric entries; letters -> NA (no warning)
  469. num_val = as.integer(replace(select, !is_num, NA_character_)),
  470. letter_score = case_when(
  471. !is_num & select == "s" ~ 0L, # "s" beats other letters
  472. !is_num ~ 1L,
  473. TRUE ~ NA_integer_
  474. ),
  475. key1 = if_else(is_num, 0L, 1L), # numbers (0) beat letters (1)
  476. key2 = if_else(is_num, num_val, letter_score) # among numbers: smallest wins; among letters: "s" wins
  477. ) %>%
  478. arrange(ID, key1, key2) %>%
  479. distinct(ID, .keep_all = TRUE) %>% # one row per GO ID
  480. select(ID, select)
  482. ck_select = ck@compareClusterResult %>%
  483. left_join(dedup_select, by="ID") %>%
  484. filter(p.adjust < 0.05)
  487. lab_map <- c(
  488. `1` = "Phagocytosis / Complement activation",
  489. `2` = "Cell adhesion",
  490. `3` = "Chemotaxis",
  491. `4` = "PH regulation",
  492. `5` = "Chemokines/Cytokines",
  493. `6` = "Protein production",
  494. l = "liam",
  495. s = "steve"
  496. )
  498. ck_select$select_label <- lab_map[as.character(ck_select$select)]
  500. write.xlsx(ck_select, "analysis/de/go_bp_batch1_padj_0.1_10deg_allselect.xlsx")
  502. # MF ------
  503. ck = readRDS("analysis/de/go_mf_batch1_padj_0.1")
  505. l = read.xlsx("analysis/de/go_mf_batch1_padj_0.1_10deg_LM.xlsx") %>%
  506. filter(!is.na(LM)) %>%
  507. rename(select = LM) %>%
  508. mutate(select = 'l') %>%
  509. select(ID, select)
  511. c = read.xlsx("analysis/de/go_mf_batch1_padj_0.1_10deg.xlsx") %>%
  512. filter(!is.na(Camilla_list)) %>%
  513. rename(select = Camilla_list ) %>%
  514. select(ID,select)
  516. all_select = rbind(l,c) %>%
  517. distinct()
  519. dedup_select <- all_select %>%
  520. mutate(
  521. is_num = grepl("^\\d+$", select),
  522. # convert only the numeric entries; letters -> NA (no warning)
  523. num_val = as.integer(replace(select, !is_num, NA_character_)),
  524. letter_score = case_when(
  525. !is_num & select == "s" ~ 0L, # "s" beats other letters
  526. !is_num ~ 1L,
  527. TRUE ~ NA_integer_
  528. ),
  529. key1 = if_else(is_num, 0L, 1L), # numbers (0) beat letters (1)
  530. key2 = if_else(is_num, num_val, letter_score) # among numbers: smallest wins; among letters: "s" wins
  531. ) %>%
  532. arrange(ID, key1, key2) %>%
  533. distinct(ID, .keep_all = TRUE) %>% # one row per GO ID
  534. select(ID, select)
  536. ck_select = ck %>%
  537. filter(ID %in% dedup_select$ID)
  539. ck_select_sim <- simplify(ck_select, cutoff = 0.7, by = "p.adjust", select_fun = min)
  541. dotplot(ck_select_sim, show=200)+
  542. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  543. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  545. ck_select_sim_df = ck_select_sim@compareClusterResult %>%
  546. left_join(dedup_select, by="ID")
  548. lab_map <- c(
  549. `1` = "Phagocytosis / Complement activation",
  550. `2` = "Cell adhesion",
  551. `3` = "Chemotaxis",
  552. `4` = "PH regulation",
  553. `5` = "Chemokines/Cytokines",
  554. `6` = "Protein production",
  555. l = "liam",
  556. s = "steve"
  557. )
  559. ck_select_sim_df$select_label <- lab_map[as.character(ck_select_sim_df$select)]
  561. write.xlsx(ck_select_sim_df, "analysis/de/go_mf_batch1_padj_0.1_10deg_allselect.xlsx")
  563. # CC -----
  564. ck = readRDS("analysis/de/go_cc_batch1_padj_0.1_10deg.rds")
  566. l = read.xlsx("analysis/de/go_cc_batch1_padj_0.1_10deg_LM.xlsx") %>%
  567. filter(!is.na(LM)) %>%
  568. rename(select = LM) %>%
  569. mutate(select = 'l') %>%
  570. select(ID, select)
  572. c = read.xlsx("analysis/de/go_cc_batch1_padj_0.1_10deg.xlsx") %>%
  573. filter(!is.na(Camilla_list)) %>%
  574. rename(select = Camilla_list ) %>%
  575. select(ID,select)
  578. all_select = rbind(l,c) %>%
  579. distinct()
  581. dedup_select <- all_select %>%
  582. mutate(
  583. is_num = grepl("^\\d+$", select),
  584. # convert only the numeric entries; letters -> NA (no warning)
  585. num_val = as.integer(replace(select, !is_num, NA_character_)),
  586. letter_score = case_when(
  587. !is_num & select == "s" ~ 0L, # "s" beats other letters
  588. !is_num ~ 1L,
  589. TRUE ~ NA_integer_
  590. ),
  591. key1 = if_else(is_num, 0L, 1L), # numbers (0) beat letters (1)
  592. key2 = if_else(is_num, num_val, letter_score) # among numbers: smallest wins; among letters: "s" wins
  593. ) %>%
  594. arrange(ID, key1, key2) %>%
  595. distinct(ID, .keep_all = TRUE) %>% # one row per GO ID
  596. select(ID, select)
  598. ck_select = ck %>%
  599. filter(ID %in% dedup_select$ID)
  601. ck_select_sim <- simplify(ck_select, cutoff = 0.7, by = "p.adjust", select_fun = min)
  603. dotplot(ck_select_sim, show=200)+
  604. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  605. scale_y_discrete(labels = ~ str_wrap(as.character(.x), 100))
  607. ck_select_sim_df = ck_select_sim@compareClusterResult %>%
  608. left_join(dedup_select, by="ID")
  610. lab_map <- c(
  611. `1` = "Phagocytosis / Complement activation",
  612. `2` = "Cell adhesion",
  613. `3` = "Chemotaxis",
  614. `4` = "PH regulation",
  615. `5` = "Chemokines/Cytokines",
  616. `6` = "Protein production",
  617. l = "liam",
  618. s = "steve"
  619. )
  621. ck_select_sim_df$select_label <- lab_map[as.character(ck_select_sim_df$select)]
  623. write.xlsx(ck_select_sim_df, "analysis/de/go_cc_batch1_padj_0.1_10deg_allselect.xlsx")
  626. # make a plot of the CC individual results with collaborator selections
  628. # ---------- helpers ----------
  629. safe_min <- function(x) if (all(is.na(x))) NA_real_ else min(x, na.rm = TRUE)
  630. strip_go_brackets <- function(x) sub("\\s*\\[[^\\]]+\\]$", "", x)
  632. # choose one 'select' per y_id:
  633. # - if any numbers → pick the smallest number
  634. # - else letters only → prefer "s", otherwise first other letter
  635. choose_select <- function(x) {
  636. xo <- as.character(x)
  637. nums <- suppressWarnings(as.integer(xo))
  638. if (any(!is.na(nums))) return(as.character(min(nums, na.rm = TRUE)))
  639. if (any(xo == "s")) return("s")
  640. setdiff(unique(xo), "s")[1]
  641. }
  643. # Start from ck_select_sim_df; map columns ----------
  645. base_df <- ck_select_sim_df %>%
  646. transmute(
  647. comparison = Cluster,
  648. GO_ID = ID,
  649. GO_Term = Description,
  650. padj = p.adjust,
  651. Significant_genes_count = Count,
  652. select
  653. )
  655. # OPTIONAL: union of Top-10 terms per comparison by padj
  656. top_ids <- base_df %>%
  657. group_by(comparison) %>%
  658. slice_min(order_by = padj, n = 100, with_ties = FALSE) %>%
  659. pull(GO_ID) %>% unique()
  660. base_df <- base_df %>% filter(GO_ID %in% top_ids)
  662. # Build plot_df ----------
  663. eps <- 1e-300
  664. cap <- 10
  666. # y_id is now JUST the term
  667. plot_df <- base_df %>%
  668. mutate(
  669. mlog10 = pmin(-log10(pmax(padj, eps)), cap),
  670. y_id = GO_Term
  671. ) %>%
  672. group_by(y_id, comparison) %>%
  673. summarise(
  674. mlog10 = if (all(is.na(mlog10))) NA_real_ else max(mlog10, na.rm = TRUE),
  675. Significant_genes_count = if (all(is.na(Significant_genes_count)))
  676. NA_real_ else max(Significant_genes_count, na.rm = TRUE),
  677. .groups = "drop"
  678. ) %>%
  679. # annotation: collapse select per DESCRIPTION
  680. left_join(
  681. base_df %>%
  682. transmute(y_id = GO_Term, select) %>%
  683. group_by(y_id) %>%
  684. summarise(select = choose_select(select), .groups = "drop"),
  685. by = "y_id"
  686. )
  688. # order y by best padj across comparisons
  689. y_levels <- base_df %>%
  690. group_by(GO_Term) %>%
  691. summarise(best_padj = safe_min(padj), .groups = "drop") %>%
  692. arrange(best_padj) %>%
  693. pull(GO_Term) %>% unique()
  695. plot_df <- plot_df %>% mutate(y_id = factor(y_id, levels = y_levels))
  697. # DOTPLOT of −log10(padj) ----------
  698. gg_dot <- ggplot(plot_df, aes(x = comparison, y = y_id)) +
  699. geom_point(aes(size = Significant_genes_count, color = mlog10), na.rm = TRUE) +
  700. scale_size_continuous(name = "Significant genes", range = c(1.8, 8)) +
  701. scale_color_viridis_c(name = expression(-log[10](padj)), option = "C",
  702. direction = 1, na.value = "grey85") +
  703. labs(title = "Selected GO Terms",
  704. subtitle = sprintf("Colored by −log10(padj), capped at %s;\n size = Significant_genes_count", cap),
  705. x = "Comparison", y = NULL) +
  706. theme_bw() +
  707. theme(axis.text.y = element_text(size = 9),
  708. axis.text.x = element_text(angle = 45, hjust = 1))
  709. print(gg_dot)
  711. # FINAL plot with all annotations. =====
  713. # Target order from functional centered on NEG
  714. well_summary=read.xlsx('20250415_Joshua_Bowen_8511/2025822_functionaldata_celllevel_exp31_only_to_show_ipa_go.xlsx')
  716. exp31 = well_summary %>%
  717. filter(Experiment == "EXP31")
  719. control <- "NEG"
  721. average_by_experiment_treatment_well <- function(cell_data) {
  722. cell_data %>%
  723. group_by(Experiment, Treatment, Metadata_Well) %>%
  724. summarise(
  725. across(AreaShape_Area:Cell_Level_Phago_Index, mean, na.rm = TRUE),
  726. .groups = "drop"
  727. )
  728. }
  730. neg_mean <- well_summary %>%
  731. filter(Treatment == control) %>%
  732. summarise(m = mean(Cell_Level_Phago_Index, na.rm = TRUE), .groups = "drop") %>%
  733. pull(m)
  735. centered_order_tbl <- well_summary %>%
  736. group_by(Treatment) %>%
  737. summarise(mean_clpi = mean(Cell_Level_Phago_Index, na.rm = TRUE), .groups = "drop") %>%
  738. mutate(delta_clpi = mean_clpi - neg_mean) %>%
  739. filter(Treatment != control) %>%
  740. arrange(desc(delta_clpi))
  742. # read GO dataframe ------
  744. df_ann <- openxlsx::read.xlsx("analysis/de/go_bp_batch1_padj_0.1_10deg_allselect_rebeccaedit.xlsx")
  746. df <- df_ann %>%
  747. janitor::clean_names() %>%
  748. mutate(
  749. select_label = as.character(select_label),
  750. select_label = stringr::str_squish(select_label),
  751. select_label = dplyr::na_if(select_label, "") # treat "" as NA
  752. ) %>%
  753. filter(!is.na(select_label)) %>%
  754. mutate(
  755. padj_safe = pmax(p_adjust, 1e-300) # guard tiny zeros
  756. )
  758. df_plot0 <- df %>%
  759. transmute(
  760. Pathway = description,
  761. Target = ko_gene,
  762. Count = count,
  763. Group = select_label,
  764. colval = padj_safe
  765. )
  767. # Keep only targets present in well_summary order
  768. go_targets <- unique(df_plot0$Target)
  770. target_levels <- centered_order_tbl %>%
  771. filter(Treatment %in% go_targets) %>%
  772. pull(Treatment)
  774. df_plot0 <- df_plot0 %>%
  775. filter(Target %in% target_levels)
  777. # count filter
  778. df_plot0 <- df_plot0 %>% filter(Count > 2)
  780. # ----------------------------
  781. # Collapse "X" with response/regulation variants
  782. # Keep, per base X: most treatments, then strongest mean -log10(padj),
  783. # then preference: exact < response < cellular response < positive < negative < regulation
  784. # ----------------------------
  785. strip_variants_to_base <- function(s) {
  786. s %>%
  787. # response variants
  788. stringr::str_replace(regex("^cellular\\s+response\\s+to\\s+", ignore_case = TRUE), "") %>%
  789. stringr::str_replace(regex("^response\\s+to\\s+", ignore_case = TRUE), "") %>%
  790. # regulation variants
  791. stringr::str_replace(regex("^positive\\s+regulation\\s+of\\s+", ignore_case = TRUE), "") %>%
  792. stringr::str_replace(regex("^negative\\s+regulation\\s+of\\s+", ignore_case = TRUE), "") %>%
  793. stringr::str_replace(regex("^regulation\\s+of\\s+", ignore_case = TRUE), "") %>%
  794. stringr::str_squish()
  795. }
  797. variant_type <- function(s) {
  798. if (stringr::str_detect(s, regex("^cellular\\s+response\\s+to\\s+", TRUE))) return("resp_cellular")
  799. if (stringr::str_detect(s, regex("^response\\s+to\\s+", TRUE))) return("response")
  800. if (stringr::str_detect(s, regex("^positive\\s+regulation\\s+of\\s+", TRUE))) return("pos")
  801. if (stringr::str_detect(s, regex("^negative\\s+regulation\\s+of\\s+", TRUE))) return("neg")
  802. if (stringr::str_detect(s, regex("^regulation\\s+of\\s+", TRUE))) return("reg")
  803. return("exact")
  804. }
  806. df_plot0 <- df_plot0 %>%
  807. mutate(
  808. base_name = strip_variants_to_base(Pathway),
  809. var_type = vapply(Pathway, variant_type, character(1))
  810. )
  812. # metrics per GO term across all targets
  813. path_metrics <- df_plot0 %>%
  814. group_by(Pathway, base_name) %>%
  815. summarise(
  816. n_treat = n_distinct(Target),
  817. mean_sig = mean(colval, na.rm = TRUE), # smaller is more significant
  818. .groups = "drop"
  819. )
  821. # preference tier: exact (1) < response (2) < cellular response (3) < positive (4) < negative (5) < regulation (6)
  822. path_metrics <- path_metrics %>%
  823. mutate(tier = dplyr::case_when(
  824. grepl("^cellular\\s+response\\s+to\\s+", Pathway, ignore.case = TRUE) ~ 3L,
  825. grepl("^response\\s+to\\s+", Pathway, ignore.case = TRUE) ~ 2L,
  826. grepl("^positive\\s+regulation\\s+of\\s+", Pathway, ignore.case = TRUE) ~ 4L,
  827. grepl("^negative\\s+regulation\\s+of\\s+", Pathway, ignore.case = TRUE) ~ 5L,
  828. grepl("^regulation\\s+of\\s+", Pathway, ignore.case = TRUE) ~ 6L,
  829. TRUE ~ 1L
  830. ))
  832. winners <- path_metrics %>%
  833. arrange(base_name, desc(n_treat), mean_sig, tier, Pathway) %>%
  834. group_by(base_name) %>%
  835. dplyr::slice(1) %>%
  836. ungroup()
  838. # keep only the winning GO term for each base X
  839. df_plot0 <- df_plot0 %>%
  840. semi_join(winners %>% dplyr::select(Pathway), by = "Pathway") %>%
  841. mutate(Group = factor(Group, levels = unique(Group)))
  843. # Keep only certain "metabolic process" categories
  845. keep_processes <- c(
  846. "ATP metabolic process",
  847. "regulation of lipid metabolic process",
  848. "superoxide metabolic process"
  849. )
  851. df_plot0 <- df_plot0 %>%
  852. filter(
  853. # keep all terms that don't end in "metabolic process"
  854. !stringr::str_detect(Pathway, "metabolic process$") |
  855. # OR are one of the specific ones you care about
  856. Pathway %in% keep_processes
  857. )
  859. # Keep union of top 10 per target
  860. n_keep <- 12
  861. top10_per_treat <- df_plot0 %>%
  862. group_by(Target) %>%
  863. arrange(colval, Pathway, .by_group = TRUE) %>% # <- ASC padj
  864. slice_head(n = n_keep) %>%
  865. ungroup()
  867. # Collect the unique GO terms that are top-10 in at least one target
  868. keepers <- unique(top10_per_treat$Pathway)
  870. # Filter the *full* table to those GO terms for *all* Treatments where they exist
  871. df_plot0 <- df_plot0 %>%
  872. dplyr::filter(Pathway %in% keepers)
  874. # Y ordering by max -log10(padj) within Group + dividers
  876. ord <- df_plot0 %>%
  877. group_by(Group, Pathway) %>%
  878. summarise(ordkey = suppressWarnings(min(colval, na.rm = TRUE)), .groups = "drop") %>% # <- MIN padj
  879. mutate(ordkey = ifelse(is.infinite(ordkey), NA_real_, ordkey)) %>%
  880. group_by(Group) %>%
  881. arrange(ordkey, Pathway, .by_group = TRUE) %>% # <- ASC padj
  882. mutate(y_in_group = row_number()) %>%
  883. ungroup()
  885. if (nrow(ord) == 0) stop("No rows left to plot after filters and collapsing.")
  887. sizes <- ord %>%
  888. dplyr::count(Group, name = "n") %>%
  889. dplyr::mutate(offset = dplyr::lag(cumsum(n), default = 0))
  891. ord <- ord %>% left_join(sizes, by = "Group") %>%
  892. mutate(y = y_in_group + offset)
  894. y_top <- max(ord$y) + 0.6
  895. y_bottom <- min(ord$y) - 0.4
  897. # positions for right-side labels and horizontal dividers
  898. label_pos <- ord %>%
  899. dplyr::group_by(Group) %>%
  900. dplyr::summarise(y = (min(y) + max(y))/2, .groups = "drop")
  902. # compute one boundary after each group, then drop ONLY the very last boundary by y
  903. dividers <- ord %>%
  904. dplyr::group_by(Group) %>%
  905. dplyr::summarise(y = max(y) + 0.5, .groups = "drop") %>%
  906. dplyr::filter(y < max(y)) # <- robust: keep all but the largest boundary
  908. df_plot <- df_plot0 %>%
  909. left_join(ord %>% dplyr::select(Group, Pathway, y), by = c("Group","Pathway"))
  911. write.xlsx(df_plot, "analysis/de/go_bp_batch1_padj_0.1_10deg_allselect_rebeccaedit_finalcatstoplot.xlsx")
  912. # Plot (x order = ΔCLPI-descending; size=count; color=-log10(padj))
  913. left_tag <- ".LEFT"
  914. right_tag <- ".RIGHT"
  915. x_levels <- target_levels
  917. df_plot <- df_plot %>%
  918. mutate(Target2 = factor(Target, levels = c(left_tag, x_levels, right_tag)))
  920. label_pos <- label_pos %>% mutate(x_lab = right_tag)
  922. p <- ggplot(df_plot, aes(x = Target2, y = y, size = Count, colour = colval)) +
  923. geom_hline(data = dividers, aes(yintercept = y),
  924. inherit.aes = FALSE, linewidth = 0.4, colour = "grey70") +
  925. geom_point(alpha = 0.9) +
  926. geom_text(data = label_pos, aes(x = x_lab, y = y, label = Group),
  927. inherit.aes = FALSE, hjust = 0, vjust = 0.5, fontface = "bold") +
  928. scale_x_discrete(
  929. drop = FALSE,
  930. limits = c(left_tag, x_levels, right_tag),
  931. breaks = x_levels,
  932. labels = x_levels,
  933. expand = expansion(add = c(0, 0.65))
  934. ) +
  935. scale_y_reverse(
  936. limits = c(y_top, y_bottom),
  937. breaks = ord$y,
  938. labels = ord$Pathway,
  939. expand = c(0, 0)
  940. ) +
  941. scale_size_area(max_size = 11, name = "Gene Count") +
  942. # clusterProfiler-like red↔blue (red = more significant)
  943. scale_colour_gradient(
  944. name = "padj",
  945. low = "red", # smaller padj = red
  946. high = "blue",
  947. oob = scales::squish,
  948. guide = guide_colorbar(reverse = TRUE) # <- lower numbers appear at the top
  949. ) +
  950. labs(x = "Target Gene", y = NULL) +
  951. coord_cartesian(clip = "off") +
  952. theme_bw(base_size = 12) +
  953. theme(
  954. plot.margin = margin(5.5, 70, 5.5, 5.5),
  955. axis.text.x = element_text(angle = 45, hjust = 1),
  956. legend.position = "right"
  957. )
  959. p
  961. # make a dotplot annotated with a boxplot --------
  962. # Pick a global font family
  963. base_family <- "Arial"
  965. # Keep the SAME x limits as dotplot (including sentinel columns)
  966. x_limits_full <- c(left_tag, x_levels, right_tag)
  968. # Compute Treatment-level mean (across wells), subtract NEG mean
  969. .delta_bar_data <- function(data, feature, control = "NEG") {
  970. # Only keep NEG + treatments that appear in the dotplot order
  971. present <- intersect(c(control, x_levels), unique(data$Treatment))
  973. means <- data %>%
  974. filter(Treatment %in% present) %>%
  975. group_by(Treatment) %>%
  976. summarise(mean_val = mean(.data[[feature]], na.rm = TRUE), .groups = "drop")
  978. neg_mean <- means %>% filter(Treatment == control) %>% pull(mean_val)
  979. if (length(neg_mean) == 0 || is.na(neg_mean)) {
  980. stop("NEG control not found in the provided data for feature: ", feature)
  981. }
  983. # Delta = gene mean - NEG mean, for non-NEG treatments in the dotplot order
  984. bars <- tibble(Treatment = x_levels) %>%
  985. left_join(means, by = "Treatment") %>%
  986. mutate(delta = mean_val - neg_mean)
  988. # Build a full x vector including sentinels (no bars for them)
  989. df <- tibble(Treatment2 = factor(x_limits_full, levels = x_limits_full)) %>%
  990. left_join(
  991. bars %>%
  992. transmute(Treatment2 = factor(Treatment, levels = x_limits_full),
  993. delta),
  994. by = "Treatment2"
  995. ) %>%
  996. mutate(sign = case_when(
  997. is.na(delta) ~ "zero", # nothing drawn
  998. delta > 0 ~ "pos",
  999. delta < 0 ~ "neg",
  1000. TRUE ~ "zero"
  1001. ))
  1002. df
  1003. }
  1005. # Use Arial in the two delta bar panels
  1006. make_delta_bar_layer <- function(data, feature, panel_title, control = "NEG") {
  1007. df <- .delta_bar_data(data, feature, control = control)
  1009. ggplot(df, aes(x = Treatment2, y = delta, fill = sign)) +
  1010. geom_col(width = 0.7, na.rm = TRUE) +
  1011. geom_hline(yintercept = 0, linewidth = 0.4, color = "grey50") +
  1012. scale_fill_manual(values = c(pos = "red", neg = "blue", zero = "grey80"), guide = "none") +
  1013. scale_x_discrete(drop = FALSE, limits = x_limits_full, breaks = NULL, expand = expansion(add = c(0, 0.65))) +
  1014. labs(title = panel_title, x = NULL, y = NULL) +
  1015. theme_minimal(base_size = 12, base_family = base_family) + # <- Arial
  1016. theme(
  1017. plot.title = element_text(hjust = 0, face = "bold", size = 11, margin = margin(b = 2)),
  1018. axis.text.x = element_blank(),
  1019. axis.title.x = element_blank(),
  1020. axis.ticks.x = element_blank(),
  1021. legend.position = "none",
  1022. plot.margin = margin(5.5, 70, 0, 5.5)
  1023. )
  1024. }
  1026. p <- p +
  1027. theme_bw(base_size = 12, base_family = base_family) + # <- Arial
  1028. theme(plot.margin = margin(5.5, 70, 5.5, 5.5))
  1030. p <- p +
  1031. geom_text(
  1032. data = label_pos,
  1033. aes(x = x_lab, y = y, label = Group),
  1034. inherit.aes = FALSE,
  1035. hjust = 0, vjust = 0.5,
  1036. fontface = "bold",
  1037. family = base_family,
  1038. size = 3
  1039. )
  1041. # Rebuild your two annotation bars with plotmath titles using Δ (delta)
  1042. p_bar_pi <- make_delta_bar_layer(
  1043. exp31, "Cell_Level_Phago_Index",
  1044. panel_title = expression(Delta~" well averaged PI"),
  1045. control = "NEG"
  1046. )
  1047. p_bar_ecc <- make_delta_bar_layer(
  1048. exp31, "AreaShape_Eccentricity",
  1049. panel_title = expression(Delta~" well-averaged Eccentricity"),
  1050. control = "NEG"
  1051. )
  1053. # make sure label_pos has the x column for the right margin labels
  1054. label_pos <- label_pos %>% mutate(x_lab = right_tag)
  1056. # 1) Remove any existing text/label geoms from the dotplot
  1057. strip_text_layers <- function(g) {
  1058. g$layers <- Filter(function(lyr) !inherits(lyr$geom, c("GeomText", "GeomLabel")), g$layers)
  1059. g
  1060. }
  1061. p <- strip_text_layers(p)
  1063. # Add ONE right-side label layer (Arial + smaller size)
  1064. base_family <- "Arial"
  1065. p <- p +
  1066. geom_text(
  1067. data = label_pos,
  1068. aes(x = x_lab, y = y, label = Group),
  1069. inherit.aes = FALSE,
  1070. hjust = 0, vjust = 0.5,
  1071. fontface = "bold",
  1072. family = base_family,
  1073. size = 3.2
  1074. )
  1076. p <- p +
  1077. theme(
  1078. axis.text.x = element_text(
  1079. angle = 45, hjust = 1, vjust = 1,
  1080. family = "Arial", size = 12
  1081. ),
  1082. # give a bit more bottom margin so angled labels don't clip
  1083. plot.margin = margin(5.5, 70, 10, 5.5)
  1084. )
  1086. # combine and save
  1087. combined <- p_bar_pi / p_bar_ecc / p + patchwork::plot_layout(heights = c(0.7, 0.7, 5))
  1089. dir.create("analysis/de", recursive = TRUE, showWarnings = FALSE)
  1090. dev <- if (capabilities("cairo")) grDevices::cairo_pdf else grDevices::pdf
  1091. ggplot2::ggsave(
  1092. filename = "analysis/de/go_final_dotplot_barplot.pdf",
  1093. plot = combined,
  1094. width = 10, height = 14, units = "in",
  1095. device = dev
  1096. )
  1098. # SINGLE CATEGORY gene heatmap ======
  1099. plot_group_heatmaps <- function(df,
  1100. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1101. group = "Lysosome",
  1102. out_dir = "analysis/plots/single_cat/",
  1103. alpha = 0.1,
  1104. lfc_thresh = 0.2,
  1105. centered_order_tbl = NULL,
  1106. min_genes_after_deg = 1,
  1107. min_treats_after_filter = 1,
  1108. cluster_rows = TRUE,
  1109. cluster_cols = TRUE,
  1110. plot_height = 17,
  1111. plot_width = 5,
  1112. base_family = "Arial",
  1113. lfc_limits = c(-4, 2)
  1114. ) {
  1116. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv"
  1117. out_dir = "analysis/plots/single_cat/"
  1118. alpha = 0.1
  1119. lfc_thresh = 0.2
  1120. min_genes_after_deg = 1
  1121. min_treats_after_filter = 1
  1122. cluster_rows = TRUE
  1123. cluster_cols = TRUE
  1124. base_family = "Arial"
  1125. lfc_limits = c(-4, 2)
  1126. df = df_phago
  1127. group = "Phagocytosis"
  1128. centered_order_tbl = centered_order_tbl
  1129. plot_width = sz$width_in
  1130. plot_height = sz$height_in
  1132. # --- Safety & inputs ---
  1133. stopifnot(all(c("select_label","description","gene_id","ko_gene") %in% names(df)))
  1134. if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
  1136. # Read DEG (only once)
  1137. deg <- data.table::fread(deg_path) %>%
  1138. dplyr::mutate(gene = gene_name)
  1140. # Terms in this group
  1141. terms <- df %>%
  1142. dplyr::filter(select_label == group) %>%
  1143. dplyr::pull(description) %>%
  1144. unique()
  1145. if (length(terms) == 0) stop(sprintf("No terms found with select_label == '%s'.", group))
  1147. # Helper to sanitize filenames
  1148. sanitize <- function(x) gsub("[^A-Za-z0-9_\\-]+", "_", x)
  1150. # PDF device that embeds fonts
  1151. .pdf_device <- if (capabilities("cairo")) grDevices::cairo_pdf else grDevices::pdf
  1153. # Plot each term
  1154. out_files <- list()
  1155. for (term_desc in terms) {
  1156. message("Plotting [", group, "]: ", term_desc)
  1158. # genes for this term
  1159. genes <- df %>%
  1160. dplyr::filter(select_label == group, description == term_desc) %>%
  1161. dplyr::pull(gene_id) %>%
  1162. stringr::str_split("/", simplify = FALSE) %>%
  1163. unlist(use.names = FALSE) %>%
  1164. stringr::str_trim() %>%
  1165. unique()
  1166. genes <- genes[genes != ""]
  1167. if (!length(genes)) { warning("No genes parsed for: ", term_desc); next }
  1169. # treatments that actually list this term in df
  1170. treats_in_term_df <- df %>%
  1171. dplyr::filter(select_label == group, description == term_desc) %>%
  1172. dplyr::pull(ko_gene) %>%
  1173. unique()
  1174. if (!length(treats_in_term_df)) { warning("No treatments listed for this term in df: ", term_desc); next }
  1176. # subset DEG
  1177. deg_sel <- deg %>%
  1178. dplyr::filter(gene %in% genes) %>%
  1179. dplyr::select(gene, log2FoldChange, padj, ko_gene)
  1180. if (!nrow(deg_sel)) { warning("No DEG rows matched for: ", term_desc); next }
  1182. # wide -> matrices
  1183. lfc_wide <- deg_sel %>%
  1184. dplyr::select(gene, ko_gene, log2FoldChange) %>%
  1185. tidyr::pivot_wider(names_from = ko_gene, values_from = log2FoldChange)
  1186. p_wide <- deg_sel %>%
  1187. dplyr::select(gene, ko_gene, padj) %>%
  1188. tidyr::pivot_wider(names_from = ko_gene, values_from = padj)
  1190. lfc_raw <- lfc_wide %>% tibble::column_to_rownames("gene") %>% as.matrix()
  1191. p_raw <- p_wide %>% tibble::column_to_rownames("gene") %>% as.matrix()
  1193. # drop cols not present in df for this term
  1194. lfc_raw <- lfc_raw[, colnames(lfc_raw) %in% treats_in_term_df, drop = FALSE]
  1195. p_raw <- p_raw[, colnames(p_raw) %in% treats_in_term_df, drop = FALSE]
  1197. # drop all-NA columns, then order by ΔCLPI (only remaining)
  1198. nonempty_cols <- colSums(!is.na(lfc_raw)) > 0
  1199. lfc_raw <- lfc_raw[, nonempty_cols, drop = FALSE]
  1200. p_raw <- p_raw[, nonempty_cols, drop = FALSE]
  1202. desired_order <- intersect(centered_order_tbl$Treatment, colnames(lfc_raw))
  1203. lfc_raw <- lfc_raw[, desired_order, drop = FALSE]
  1204. p_raw <- p_raw[, desired_order, drop = FALSE]
  1205. if (ncol(lfc_raw) < min_treats_after_filter) { warning("Too few treatments after filtering for: ", term_desc); next }
  1207. # fill NA (after pruning)
  1208. lfc_mat <- lfc_raw; lfc_mat[is.na(lfc_mat)] <- 0
  1209. p_mat <- p_raw; p_mat[is.na(p_mat)] <- 1
  1211. # drop genes with no data across remaining treatments
  1212. nonempty_rows <- rowSums(!is.na(lfc_raw)) > 0
  1213. lfc_mat <- lfc_mat[nonempty_rows, , drop = FALSE]
  1214. p_mat <- p_mat[nonempty_rows, , drop = FALSE]
  1215. if (nrow(lfc_mat) < min_genes_after_deg) { warning("Too few genes with data left for: ", term_desc); next }
  1217. heatmap_max <- suppressWarnings(max(lfc_mat, na.rm = TRUE)); if (!is.finite(heatmap_max)) heatmap_max <- 1
  1218. heatmap_min <- suppressWarnings(min(lfc_mat, na.rm = TRUE)); if (!is.finite(heatmap_min)) heatmap_min <- -1
  1220. # star overlay (use Arial)
  1221. cell_fun_star <- function(j, i, x, y, w, h, fill) {
  1222. if (!is.na(p_mat[i, j]) && p_mat[i, j] < alpha &&
  1223. !is.na(lfc_mat[i, j]) && abs(lfc_mat[i, j]) > lfc_thresh) {
  1224. grid::grid.text("*", x, y, gp = grid::gpar(fontfamily = base_family))
  1225. }
  1226. }
  1227. heatmap_min <- lfc_limits[1]
  1228. heatmap_max <- lfc_limits[2]
  1230. # keep white at 0 (diverging)
  1231. col_fun <- circlize::colorRamp2(
  1232. c(heatmap_min, 0, heatmap_max),
  1233. c("darkblue", "white", "red")
  1234. )
  1236. # legend ticks (top→bottom shows 2, 0, -2, -4). If yours renders bottom→top,
  1237. # swap to at = c(-4, -2, 0, 2) to match your preference.
  1238. legend_at <- c(2, 0, -2, -4)
  1239. legend_labels <- c("2", "0", "-2", "-4")
  1241. pdf_file <- file.path(out_dir, paste0(sanitize(group), "__", sanitize(term_desc), ".pdf"))
  1242. .pdf_device(pdf_file, height = plot_height, width = plot_width, family = base_family)
  1244. view(lfc_mat)
  1245. ht <- ComplexHeatmap::Heatmap(
  1246. lfc_mat,
  1247. name = "Log2FC",
  1248. col = col_fun,
  1249. cluster_columns = cluster_cols,
  1250. cluster_rows = cluster_rows,
  1251. column_names_gp = grid::gpar(fontsize = 10, fontfamily = base_family),
  1252. row_names_gp = grid::gpar(fontsize = 10, fontfamily = base_family),
  1253. heatmap_legend_param = list(
  1254. title = "Log2FoldChange\nKO vs. Control",
  1255. title_gp = grid::gpar(fontfamily = base_family, fontface = "bold"),
  1256. labels_gp = grid::gpar(fontfamily = base_family),
  1257. at = legend_at,
  1258. labels = legend_labels
  1259. ),
  1260. column_title = term_desc,
  1261. column_title_gp = grid::gpar(fontfamily = base_family, fontface = "bold"),
  1262. cell_fun = cell_fun_star
  1263. )
  1265. ComplexHeatmap::draw(ht)
  1266. grDevices::dev.off()
  1267. message("Saved: ", pdf_file)
  1268. out_files[[term_desc]] <- pdf_file
  1269. }
  1270. invisible(out_files)
  1271. }
  1273. estimate_heatmap_dims <- function(df_plot, deg_path, group, centered_order_tbl) {
  1274. # genes for this group/term set (split "gene_id" on '/')
  1275. genes <- df_plot %>%
  1276. filter(select_label == group) %>%
  1277. pull(gene_id) %>%
  1278. str_split("/", simplify = FALSE) %>%
  1279. unlist(use.names = FALSE) %>%
  1280. str_trim() %>%
  1281. unique()
  1282. genes <- genes[genes != ""]
  1284. # treatments mentioned for these terms
  1285. treats_in_term_df <- df_plot %>%
  1286. filter(select_label == group) %>%
  1287. pull(ko_gene) %>%
  1288. unique()
  1290. # read DEG and subset to these genes/treatments
  1291. deg <- data.table::fread(deg_path) %>% mutate(gene = gene_name)
  1292. deg_sel <- deg %>%
  1293. filter(gene %in% genes, ko_gene %in% treats_in_term_df) %>%
  1294. select(gene, ko_gene, log2FoldChange)
  1296. # wide, then drop all-NA cols and keep centered order
  1297. lfc_wide <- deg_sel %>%
  1298. pivot_wider(names_from = ko_gene, values_from = log2FoldChange)
  1299. if (nrow(lfc_wide) == 0) return(list(nr = 0L, nc = 0L))
  1301. lfc_raw <- lfc_wide %>% column_to_rownames("gene") %>% as.matrix()
  1302. nonempty_cols <- colSums(!is.na(lfc_raw)) > 0
  1303. lfc_raw <- lfc_raw[, nonempty_cols, drop = FALSE]
  1305. desired_order <- intersect(centered_order_tbl$Treatment, colnames(lfc_raw))
  1306. lfc_raw <- lfc_raw[, desired_order, drop = FALSE]
  1308. # drop genes with all NA across remaining treatments
  1309. nonempty_rows <- rowSums(!is.na(lfc_raw)) > 0
  1310. lfc_raw <- lfc_raw[nonempty_rows, , drop = FALSE]
  1312. list(nr = nrow(lfc_raw), nc = ncol(lfc_raw))
  1313. }
  1315. heatmap_device_size <- function(nr, nc,
  1316. cell_w_mm = 6,
  1317. cell_h_mm = 6.5,
  1318. row_names_mm = 35,
  1319. col_names_mm = 12,
  1320. row_dend_mm = 4,
  1321. col_dend_mm = 4,
  1322. legend_mm = 18,
  1323. gap_mm = 6) {
  1324. mm_to_in <- function(mm) mm / 25.4
  1325. heatmap_w_in <- mm_to_in(cell_w_mm * nc)
  1326. heatmap_h_in <- mm_to_in(cell_h_mm * nr)
  1327. width_in <- heatmap_w_in + mm_to_in(row_names_mm + row_dend_mm + gap_mm + legend_mm)
  1328. height_in <- heatmap_h_in + mm_to_in(col_names_mm + col_dend_mm + gap_mm)
  1329. list(width_in = width_in, height_in = height_in)
  1330. }
  1333. # Example: Phagocytosis panel
  1334. df_phago <- df %>% filter(description == "phagocytosis")
  1335. head(df_phago)
  1337. dims <- estimate_heatmap_dims(df_phago, deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1338. group = "Phagocytosis", centered_order_tbl = centered_order_tbl)
  1340. sz <- heatmap_device_size(dims$nr, dims$nc)
  1342. files_phago <- plot_group_heatmaps(
  1343. df_phago,
  1344. group = "Phagocytosis",
  1345. centered_order_tbl = centered_order_tbl,
  1346. plot_width = sz$width_in,
  1347. plot_height = sz$height_in
  1348. )
  1350. # Lysosome terms
  1351. df_plot = df %>%
  1352. filter(description == "lysosome organization")
  1354. dims <- estimate_heatmap_dims(df_plot,
  1355. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1356. group = "Lysosome",
  1357. centered_order_tbl = centered_order_tbl)
  1359. sz <- heatmap_device_size(dims$nr, dims$nc)
  1360. files_lyso <- plot_group_heatmaps(
  1361. df_plot,
  1362. group = "Lysosome",
  1363. centered_order_tbl = centered_order_tbl,
  1364. plot_width = sz$width_in,
  1365. plot_height = sz$height_in
  1366. )
  1368. # Activation terms
  1369. df_plot = df %>%
  1370. filter(description == "macrophage activation")
  1373. head(df_plot)
  1374. dims <- estimate_heatmap_dims(df_plot,
  1375. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1376. group = "Inflammation/Activation",
  1377. centered_order_tbl = centered_order_tbl)
  1379. sz <- heatmap_device_size(dims$nr, dims$nc)
  1380. files_act <- plot_group_heatmaps(
  1381. df_plot,
  1382. group = "Inflammation/Activation",
  1383. centered_order_tbl = centered_order_tbl,
  1384. plot_width = sz$width_in,
  1385. plot_height = sz$height_in
  1386. )
  1388. # Actin terms
  1389. df_plot = df %>%
  1390. filter(description == "regulation of cell shape")
  1391. dims <- estimate_heatmap_dims(df_plot,
  1392. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1393. group = "Actin/Cell Shape",
  1394. centered_order_tbl = centered_order_tbl)
  1396. sz <- heatmap_device_size(dims$nr, dims$nc)
  1397. files_act <- plot_group_heatmaps(
  1398. df_plot,
  1399. group = "Actin/Cell Shape",
  1400. centered_order_tbl = centered_order_tbl,
  1401. plot_width = sz$width_in,
  1402. plot_height = sz$height_in
  1403. )
  1405. # Actin term
  1406. df_plot = df %>%
  1407. filter(description == "actin crosslink formation")
  1408. dims <- estimate_heatmap_dims(df_plot,
  1409. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1410. group = "Actin/Cell Shape",
  1411. centered_order_tbl = centered_order_tbl)
  1413. sz <- heatmap_device_size(dims$nr, dims$nc)
  1414. files_act <- plot_group_heatmaps(
  1415. df_plot,
  1416. group = "Actin/Cell Shape",
  1417. centered_order_tbl = centered_order_tbl,
  1418. plot_width = sz$width_in,
  1419. plot_height = sz$height_in
  1420. )
  1422. df_plot = df %>%
  1423. filter(description == "actin filament bundle assembly")
  1424. dims <- estimate_heatmap_dims(df_plot,
  1425. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1426. group = "Actin/Cell Shape",
  1427. centered_order_tbl = centered_order_tbl)
  1429. sz <- heatmap_device_size(dims$nr, dims$nc)
  1430. files_act <- plot_group_heatmaps(
  1431. df_plot,
  1432. group = "Actin/Cell Shape",
  1433. centered_order_tbl = centered_order_tbl,
  1434. plot_width = sz$width_in,
  1435. plot_height = sz$height_in
  1436. )
  1439. df_plot = df %>%
  1440. filter(description == "regulation of actin cytoskeleton organization")
  1441. dims <- estimate_heatmap_dims(df_plot,
  1442. deg_path = "analysis/de/all_res_vs_ntc1_batch1_q1q3_filter_rmlow.tsv",
  1443. group = "Actin/Cell Shape",
  1444. centered_order_tbl = centered_order_tbl)
  1446. sz <- heatmap_device_size(dims$nr, dims$nc)
  1447. files_act <- plot_group_heatmaps(
  1448. df_plot,
  1449. group = "Actin/Cell Shape",
  1450. centered_order_tbl = centered_order_tbl,
  1451. plot_width = sz$width_in,
  1452. plot_height = sz$height_in
  1453. )

03_01_analyze_deg_go.R at commit 7ea1019, no license · at the source

Overview

Authors: Joy E. Horng1,2, Liam T. McCrea1, Rebecca E. Batorsky3, Joshua J. Bowen1, Camilla Boschian1,4, Yoonjae Song1, Roy H. Perlis1,4,5, Steven D. Sheridan1,4,5
ORCID iDs: Joshua J. Bowen, Roy H. Perlis, Steven D. Sheridan
  1. Center for Genomic Medicine and Department of Psychiatry, Massachusetts General Hospital,Boston, MA USA
  2. Department of Molecular Biology, Massachusetts General Hospital,Boston, MA USA
  3. Tufts Institute for Artificial Intelligence, Tufts University,Medford, MA USA
  4. Department of Psychiatry, Harvard Medical School,Boston, MA USA
  5. Harvard Stem Cell Institute,Cambridge, MA USA
Institutions: Massachusetts General Hospital (United States); Tufts University (United States); Harvard University (United States); Harvard Stem Cell Institute (United States)
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology, volume 51, issue 9, pages 1718-1727
Dates: received 15 November 2025; accepted 31 March 2026; published online 16 April 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41386-026-02406-1 · PMID 41992032 · PMCID PMC13389400 · OpenAlex W4416337181
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.11.14.688490
Status: code verified
Categories: genetics / omics (modality), human (organism), schizophrenia / psychosis (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions
Keywords: Cellular neuroscience, Genetics of the nervous system
MeSH: Microglia*, Schizophrenia*, Gene Expression Profiling, Genomics, Humans, Phagocytosis (* major topic)
Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex
Funding: HHS | NIH | National Institute of Mental Health (NIMH) (R01MH120227, R01MH131687)
Citations: not cited yet (Europe PMC); 75 references in the paper

Abstract

Microglia are increasingly recognized as key regulators of neural circuit development and putative contributors to the pathophysiology of neuropsychiatric disorders such as schizophrenia (SCZ). However, the functional impact of SCZ-associated genes in microglia remains largely unexplored. Here, we performed an arrayed CRISPR targeting screen of 30 SCZ-associated genes predicted to be differentially expressed in human microglia-like cells. Target genes were prioritized based on post-mortem transcriptomic relevance and predicted ontology-based roles in phagocytosis pathways. We quantified phagocytic activity and morphological changes following gene targeting using high-content confocal imaging. Key targets, including CYFIP1, MSR1, TREM2, SYK, ITGB2, ITGAM, and IRF8, modulated phagocytosis and altered morphological properties consistent with activation states, validating their functional roles in microglia. To elucidate transcriptional impact, we further applied a multiplexed RNA sequencing platform across gene targets. These analyses revealed gene-specific transcriptional signatures, implicating divergent pathways related to phagocytic, activation, cytoskeletal, and lysosomal function. Together, these findings demonstrate the utility of CRISPR-based functional genomics in characterizing microglia function and identifying new target genes and mechanisms that may underlie their contributions to SCZ pathophysiology.

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 7 matches between paragraphs and lines of code.

mgb-cqh/crispr30

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 7ea101924cd85986ebddc043eb9d661e32c241b3, 21 November 2025
Languages: R (6), Shell (5)
Size: 14 files, 11 scripts
Software Heritage: not archived
Found in: the text, “RNA-seq data analysis”
Holds: README, 2 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (6 files), data.table (5 files), clusterProfiler (4 files), DESeq2 (4 files), pheatmap (4 files), ComplexHeatmap (3 files), ggpubr (3 files), SAMtools (2 files), broom (1 file), circlize (1 file), emmeans (1 file), FastQC (1 file), ggplot2 (1 file), igraph (1 file), limma (1 file), lme4 (1 file), multcomp (1 file), patchwork (1 file), STAR (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
12 files
At the source: github.com/mgb-cqh/crispr30

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

Gene expression data will be made available for download from the NCBI Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE319209) upon publication. Additional data that support the findings of this study will be available from the corresponding authors upon reasonable request.

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

Versions

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

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 2 keywords, 6 MeSH terms, 1 funder, 74 references.

Cite

This paper

Horng, J. E., McCrea, L. T., Batorsky, R. E., Bowen, J. J., Boschian, C., Song, Y., Perlis, R. H., & Sheridan, S. D. (2026). Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators. Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology, 51(9), 1718-1727. https://doi.org/10.1038/s41386-026-02406-1

BibTeX

@article{horng2026functional,
author = {Horng, Joy E. and McCrea, Liam T. and Batorsky, Rebecca E. and Bowen, Joshua J. and Boschian, Camilla and Song, Yoonjae and Perlis, Roy H. and Sheridan, Steven D.},
title = {{Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators}},
journal = {Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology},
year = {2026},
month = apr,
volume = {51},
number = {9},
pages = {1718--1727},
publisher = {Nature Publishing Group},
issn = {0893-133X},
doi = {10.1038/s41386-026-02406-1},
url = {https://doi.org/10.1038/s41386-026-02406-1},
pmid = {41992032},
pmcid = {PMC13389400}
}

RIS

TY - JOUR
AU - Horng, Joy E.
AU - McCrea, Liam T.
AU - Batorsky, Rebecca E.
AU - Bowen, Joshua J.
AU - Boschian, Camilla
AU - Song, Yoonjae
AU - Perlis, Roy H.
AU - Sheridan, Steven D.
TI - Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators
T2 - Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
J2 - Neuropsychopharmacology
PY - 2026
DA - 2026/04/16
VL - 51
IS - 9
SP - 1718
EP - 1727
SN - 0893-133X
PB - Nature Publishing Group
DO - 10.1038/s41386-026-02406-1
UR - https://doi.org/10.1038/s41386-026-02406-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41386-026-02406-1",
"type": "article-journal",
"title": "Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators",
"container-title": "Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology",
"author": [
{
"family": "Horng",
"given": "Joy E."
},
{
"family": "McCrea",
"given": "Liam T."
},
{
"family": "Batorsky",
"given": "Rebecca E."
},
{
"family": "Bowen",
"given": "Joshua J."
},
{
"family": "Boschian",
"given": "Camilla"
},
{
"family": "Song",
"given": "Yoonjae"
},
{
"family": "Perlis",
"given": "Roy H."
},
{
"family": "Sheridan",
"given": "Steven D."
}
],
"container-title-short": "Neuropsychopharmacology",
"volume": "51",
"issue": "9",
"page": "1718-1727",
"DOI": "10.1038/s41386-026-02406-1",
"PMID": "41992032",
"PMCID": "PMC13389400",
"ISSN": "0893-133X",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41386-026-02406-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
16
]
]
}
}

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.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: FastQC, multcomp, STAR, 16 other tools, genetics / omics, 4 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: multcomp, limma, DESeq2, 10 other tools, cellular / molecular, 4 references
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: SAMtools, limma, broom, 12 other tools, cellular / molecular
[4] doi:10.1038/s41467-026-74753-y [code]
A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: STAR, SAMtools, limma, 11 other tools, cellular / molecular, 1 reference
[5] doi:10.1016/j.isci.2026.116412 [code]
KOLF2.1J iTF-Microglia: A standardized platform to study microglial transcriptional regulatory networks in CNS disease.
Journal: iScience
In common: SAMtools, DESeq2, circlize, 7 other tools, genetics / omics, 6 references
[6] 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: STAR, SAMtools, limma, 9 other tools, genetics / omics, 3 references
[7] 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: limma, DESeq2, igraph, 10 other tools, genetics / omics, cellular / molecular, 2 references
[8] doi:10.1038/s41593-026-02300-5 [code]
Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: FastQC, SAMtools, DESeq2, 9 other tools, genetics / omics, cellular / molecular, 1 reference
[9] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: FastQC, STAR, SAMtools, 8 other tools, genetics / omics, 3 references
[10] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: STAR, SAMtools, limma, 8 other tools, cellular / molecular, 3 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.