cluster_colors_AB <- c("#8B63A6","#8B63A6","#d7301f","#d7301f","#0570b0","#0570b0") # purple, red, blue
cluster_colors_glia <- c("#53B848","#53B848","#E8C51E","#E8C51E","#364eA1","#364eA1") # green, yellow, blue

cluster_colors <- c("#0A7B83","#7ABAF2","#FFD265","#595959","#2AA876","#F19C65","#8B63A6","#607D8B","#ff5252","#4C1B1B","#fb9a99")
new.cluster.ids <- c( "L5", "L3-5" ,"WM", "L1-5 Low UMI", "L2-3", "L6", "L1/Astrocytes",  "Meninges", "Blood Vessel", "Interneuron", "Mixed glia" )
names(cluster_colors) = new.cluster.ids

Load and combine stratification metadata

load(paste0(work_dir,"/../cell2location/UpdatedMetadataFeb23/ST_metadata_032723.withGenes.RData"))
# Denis classification
library(readxl)
data_prefix = "new_strats_062623"
denis_strat = fread(paste0(work_dir,"/../cell2location/stratifications/new_strats_062623.csv"), sep = ",") %>% 
  as.data.frame()
#names(denis_strat) # "no_abeta_NEW"   "low_abeta_NEW"  "high_abeta_NEW" "low_RGN_neg"    "low_RGN_pos"    "high_RGN_neg"   "high_RGN_pos" 

# Add stratifications to the metadata
denis_strat_m = denis_strat %>% 
  pivot_longer(cols = everything(), names_to = "strat", values_to = "barcode") %>% 
  distinct()
denis_strat_m$strat = gsub("_NEW","",denis_strat_m$strat)
denis_strat_m$strat = gsub("abeta","Abeta",denis_strat_m$strat)
denis_strat_m = bind_rows(denis_strat_m, 
                          data.frame(strat = "RGN_neg", 
                                     barcode = unique(c(denis_strat$low_RGN_neg, denis_strat$high_RGN_neg))))
denis_strat_m = bind_rows(denis_strat_m, 
                          data.frame(strat = "RGN_pos", 
                                     barcode = unique(c(denis_strat$low_RGN_pos, denis_strat$high_RGN_pos))))

metadata_with_strat = metadata %>% inner_join(denis_strat_m, by = c("V1"="barcode"))
### Filter by GM only
metadata_with_strat = metadata_with_strat[metadata_with_strat$brainRegion == "GM",]

#table(metadata_with_strat$strat)
message("Total number of spots: ", nrow(metadata_with_strat))
message("Stratification groups: ", 
        paste(names(table(metadata_with_strat$strat)), table(metadata_with_strat$strat), sep = " = ", collapse = "; "))

Prepare gene expression data

load(paste0(work_dir,"/230927/RGN_NICHES_PseudoBulk_GM.RData"))

Make gene background

niches_db = NICHES::LoadOmniPath("human")
ligand_background = unique(na.omit(unlist(as.list(niches_db$source.subunits))))
receptors_background = unique(na.omit(unlist(as.list(niches_db$target.subunits))))
lr_background = unique(c(ligand_background,receptors_background))

res_df = niches_df_norm$mat_counts %>% rownames_to_column("ligand_receptor")
# Get gene background
lr_sep = "—"
res_df$ligand_receptor = gsub("\"","",res_df$ligand_receptor)
res_df$ligand_receptor = gsub("HLA-","HLA/",res_df$ligand_receptor)
res_df$ligand_receptor = gsub("EXOC3-","EXOC3/",res_df$ligand_receptor)
res_df$ligand_receptor = gsub("PIK3CD-","PIK3CD/",res_df$ligand_receptor)
res_df$ligand_receptor = gsub("IGKV1-","IGKV1/",res_df$ligand_receptor)
res_df$ligand_receptor = gsub("KRTAP4-","KRTAP4/",res_df$ligand_receptor)

fix_gene_names <- function(genelist){
  genelist = gsub("HLA/","HLA-",genelist)
  genelist = gsub("EXOC3/","EXOC3-",genelist)
  genelist = gsub("PIK3CD/","PIK3CD-",genelist)
  genelist = gsub("IGKV1/","IGKV1-",genelist)
  genelist = gsub("KRTAP4/","KRTAP4-",genelist)
  return(genelist)
}

background_genes = unique(unlist(str_split(res_df$ligand_receptor,lr_sep)))
background_genes = unique(unlist(str_split(background_genes,"-")))
background_genes = fix_gene_names(background_genes)

library(HGNChelper)
background_genes = background_genes[which(utf8::utf8_valid(background_genes))]
#background_genes.check = checkGeneSymbols(background_genes)
#background_genes.check[!background_genes.check$Approved,] # 7 invalid
# Get expression matrices (pseudo-bulk)
metadata_pseudo = distinct(niches_df_norm$long[,c("sample_section","Bnum","SpotSet","age","ad_reagan")]) %>% 
  mutate(rowid = paste0(sample_section,"_",SpotSet)) %>% 
  column_to_rownames("rowid")

metadata_4deg_RGNpos_RGNneg = metadata_pseudo %>% filter(SpotSet %in% c("RGN_neg","RGN_pos"))
metadata_4deg_RGNpos_RGNneg$SpotSet = factor(metadata_4deg_RGNpos_RGNneg$SpotSet, levels = c("RGN_neg","RGN_pos"))

metadata_4deg_RGNpos_RGNneg_low = metadata_pseudo %>% filter(SpotSet %in% c("low_RGN_neg","low_RGN_pos"))
metadata_4deg_RGNpos_RGNneg_low$SpotSet = factor(metadata_4deg_RGNpos_RGNneg_low$SpotSet, levels = c("low_RGN_neg","low_RGN_pos"))

metadata_4deg_RGNpos_RGNneg_high = metadata_pseudo %>% filter(SpotSet %in% c("high_RGN_neg","high_RGN_pos"))
metadata_4deg_RGNpos_RGNneg_high$SpotSet = factor(metadata_4deg_RGNpos_RGNneg_high$SpotSet, levels = c("high_RGN_neg","high_RGN_pos"))
res_mean_RGNpos_RGNneg = analisys_wrapper(
  counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_RGNpos_RGNneg)] , 
  metadata_4deg = metadata_4deg_RGNpos_RGNneg, evcodes=T)

res_mean_RGNpos_RGNneg_low = analisys_wrapper(
  counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_RGNpos_RGNneg_low)] , 
  metadata_4deg = metadata_4deg_RGNpos_RGNneg_low, evcodes=T)

res_mean_RGNpos_RGNneg_high = analisys_wrapper(
  counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_RGNpos_RGNneg_high)], 
  metadata_4deg = metadata_4deg_RGNpos_RGNneg_high, evcodes=T)

save(res_mean_RGNpos_RGNneg, file = paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_GM.RData"))
save(res_mean_RGNpos_RGNneg_low, file = paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_low_GM.RData"))
save(res_mean_RGNpos_RGNneg_high, file = paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_high_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_low_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_RGNpos_RGNneg_high_GM.RData"))

Correlation between DEs

res_df_model = res_mean_RGNpos_RGNneg$DE_result %>% mutate(model = "glia-high vs. glia-low")
res_df_model = res_df_model %>% bind_rows(res_mean_RGNpos_RGNneg_low$DE_result %>% mutate(model = "glia-high vs. glia-low (low Abeta)"))
res_df_model = res_df_model %>% bind_rows(res_mean_RGNpos_RGNneg_high$DE_result %>% mutate(model = "glia-high vs. glia-low (high Abeta)"))

dups = res_df_model %>%
  dplyr::group_by(ligand_receptor, model) %>%
  dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
  dplyr::filter(n > 1L)
res_df_model = res_df_model[!res_df_model$ligand_receptor %in% unique(dups$ligand_receptor), ]
tmp_df = res_df_model %>% pivot_wider(id_cols = ligand_receptor, names_from = model, values_from = t)

cor_df = as.data.frame(cor(tmp_df[,-1],use = "pairwise.complete.obs", method = "pearson"))
breaksList = seq(-1, 1, by = 0.01)

pheatmap::pheatmap(cor_df, color=colorRampPalette(c(pal_aaas()(2)[2],"white",pal_aaas()(2)[1]))(length(breaksList)), 
                   display_numbers = T, number_color = "black",
                   main = "Correlation (t-statistics)", breaks = breaksList)

Complete analysis

PCA

res_mean_RGNpos_RGNneg$PCA_raw + res_mean_RGNpos_RGNneg_low$PCA_raw + res_mean_RGNpos_RGNneg_high$PCA_raw

DE

rownames(res_df_model) = NULL
res_df_model = res_df_model[order(res_df_model$adj.P.Val),]
createDT(res_df_model)
ggarrange(res_mean_RGNpos_RGNneg$p_volcano + labs(title = "glia-high vs. glia-low"),
  res_mean_RGNpos_RGNneg_low$p_volcano + labs(title = "glia-high vs. glia-low (low)"),
  res_mean_RGNpos_RGNneg_high$p_volcano + labs(title = "glia-high vs. glia-low (high)"), ncol = 3)

Top gene-pairs

ggarrange(res_mean_RGNpos_RGNneg$top_genes_plot,
  res_mean_RGNpos_RGNneg_low$top_genes_plot,
  res_mean_RGNpos_RGNneg_high$top_genes_plot, ncol = 3, labels = c("All","Low","High"), hjust = -8)

Enrichments

enrichments_combined_df = bind_rows(res_mean_RGNpos_RGNneg$enrich_res %>% mutate(contrast = "RGNpos-RNGneg"),
res_mean_RGNpos_RGNneg_low$enrich_res %>% mutate(contrast = "lowRGNpos-lowRNGneg"),
res_mean_RGNpos_RGNneg_high$enrich_res %>% mutate(contrast = "highRGNpos-highRNGneg")) %>% 
  arrange(-logP) %>%
  select(contrast, Direction, source, term_name, p_value, precision, recall, model, intersection) %>% filter(model == "LR") 

enrich_selected = read_xlsx(paste0(work_dir,"/../ST_NICHES_RGN_selected.xlsx"))
enrich_selected = enrich_selected$term_name[enrich_selected$highlight=="T"]

createDT(enrichments_combined_df[enrichments_combined_df$p_value<0.05,])
suppressWarnings(ggarrange(res_mean_RGNpos_RGNneg$enrich_res_LR$enrich_p + labs(title = "glia-high vs. glia-low"),
  res_mean_RGNpos_RGNneg_low$enrich_res_LR$enrich_p + labs(title = "glia-high vs. glia-low (low Abeta)"),
  res_mean_RGNpos_RGNneg_high$enrich_res_LR$enrich_p + labs(title = "glia-high vs. glia-low (high Abeta)"), nrow = 3))
#top_10_UP = enrichments_combined_df[enrichments_combined_df$Direction == "UP",]$term_name[1:5]
#top_10_DOWN = enrichments_combined_df[enrichments_combined_df$Direction == "DOWN",]$term_name[1:5]
#enrich_selected = unique(c(enrich_selected,top_10_UP,top_10_DOWN))

enrichments_combined_df$highlight = enrichments_combined_df$term_name %in% enrich_selected

absmax <- function(x) { x[which.max( abs(x) )]}
absmin <- function(x) { x[which.min( abs(x) )]}

# up_down_mtx = enrichments_combined_df[enrichments_combined_df$highlight,] %>% 
#   filter(p_value <= 0.05) %>%
#   mutate(signed_logP = ifelse(Direction=="DOWN",log10(p_value),-log10(p_value))) %>%
#   dplyr::select(term_name,contrast,signed_logP,Direction) %>% distinct() %>%
#   pivot_wider(id_cols = term_name, names_from = c(contrast,Direction), values_from = signed_logP) %>%
#   column_to_rownames("term_name")
# rg = max(abs(up_down_mtx),na.rm = T)
# pheatmap(up_down_mtx, cluster_rows = F, cluster_cols = F, 
#          color=colorRampPalette(c("navy", "white", "darkred"))(100),
#          breaks = seq(-rg, rg, length.out = 100))

enrich_selected_highlight_pval = enrichments_combined_df[enrichments_combined_df$highlight,] %>%
  filter(p_value <= 0.05) %>%
  dplyr::select(term_name,contrast,p_value) %>% distinct() %>%
  pivot_wider(id_cols = term_name, names_from = contrast, values_from = p_value, values_fn = absmin) %>% 
  column_to_rownames("term_name")
enrich_selected_highlight_pval[is.na(enrich_selected_highlight_pval)] <- 2

enrich_selected_highlight = enrichments_combined_df[enrichments_combined_df$highlight,] %>% 
  filter(p_value <= 0.05) %>%
  mutate(signed_logP = ifelse(Direction=="DOWN",log10(p_value),-log10(p_value))) %>%
  dplyr::select(term_name,contrast,signed_logP) %>% distinct() %>%
  pivot_wider(id_cols = term_name, names_from = contrast, values_from = signed_logP, values_fn = absmax) %>% 
  column_to_rownames("term_name")

enrich_selected_highlight_pval = enrich_selected_highlight_pval[,c("lowRGNpos-lowRNGneg","highRGNpos-highRNGneg","RGNpos-RNGneg")]
enrich_selected_highlight = enrich_selected_highlight[,c("lowRGNpos-lowRNGneg","highRGNpos-highRNGneg","RGNpos-RNGneg")]
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="RGNpos-RNGneg"] = "glia-high vs. glia-low (All)"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="lowRGNpos-lowRNGneg"] = "glia-high vs. glia-low (low AB)"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="highRGNpos-highRNGneg"] = "glia-high vs. glia-low (high AB)"
enrich_selected_highlight = as.matrix(enrich_selected_highlight)

(enrich_heatmap = ComplexHeatmap::Heatmap(enrich_selected_highlight,
                        cell_fun = function(j, i, x, y, width, height, fill) {
                               if((enrich_selected_highlight_pval[i, j] <= 0.001)) {
                                 grid.text("***", x, y, gp = gpar(fontsize = 8, fontfamily = "sans"))
                               } else if((enrich_selected_highlight_pval[i, j] <= 0.01)) {
                                 grid.text("**", x, y, gp = gpar(fontsize = 8, fontfamily = "sans"))
                               } else if((enrich_selected_highlight_pval[i, j] <= 0.05)) {
                                 grid.text("*", x, y, gp = gpar(fontsize = 8, fontfamily = "sans"))
                               } else if((enrich_selected_highlight_pval[i, j] <= 0.1)) {
                                 grid.text(".", x, y, gp = gpar(fontsize = 8, fontfamily = "sans"))
                               }
                          },
                        na_col = "lightgrey", name = "Signed log10(P)",
                        col = colorRamp2(c(-max(abs(enrich_selected_highlight),na.rm = T), -2,0,2,
                                           max(abs(enrich_selected_highlight),na.rm = T)), 
                                         rev(brewer.pal(n=7, name="RdBu"))[c(1,3,4,5,7)]),
                        cluster_columns = F, cluster_rows = F, row_names_gp = gpar(fontsize = 9),
                        column_names_gp = gpar(rotate = 45),
                        row_names_side = "left", show_row_names = T,
                        column_names_rot = 45,    rect_gp = gpar(col = "white", lwd = 0.1),
                        row_names_max_width = max_text_width(rownames(enrich_selected_highlight), gp = gpar(fontsize = 9))))

pdf(file = paste0(work_dir,"/Figures/",figDir_id,"_GliaNet_LR_enrichment_heatmap.pdf"), width = 5, height = 7)
enrich_heatmap
dev.off()
## png 
##   2

Comparison between contrasts

## Get only LR pairs with expression in example feature plot (B16002_4)
LR_with_expression_in_example = niches_df_norm$long %>% filter(sample_section == "B16002_4") %>%
  select(ligand_receptor, estimate_fraction) %>%
  distinct() %>% filter(estimate_fraction > 100)
LR_with_expression_in_example = unique(LR_with_expression_in_example$ligand_receptor)
res_combined = res_mean_RGNpos_RGNneg_low$DE_result[,c("ligand_receptor","t","P.adj")] %>%
  inner_join(res_mean_RGNpos_RGNneg_high$DE_result[,c("ligand_receptor","t","P.adj")], 
             suffix = c(".low",".high"), 
             by = "ligand_receptor") 

res_combined = na.omit(res_combined)
res_combined$color = "NS"
res_combined$color[res_combined$P.adj.low<=0.05 & res_combined$P.adj.high<=0.05] = "Both Signf"
res_combined$color[res_combined$P.adj.low<=0.05 & res_combined$P.adj.high>0.05] = "glia-high vs. glia-low (low AB) only"
res_combined$color[res_combined$P.adj.low>0.05 & res_combined$P.adj.high<=0.05] = "glia-high vs. glia-low (high AB) only"
res_combined$color = factor(res_combined$color, 
                            levels = c("Both Signf",
                                       "glia-high vs. glia-low (low AB) only",
                                       "glia-high vs. glia-low (high AB) only",
                                       "NS"))


res_combined = res_combined[order(rowMeans(res_combined[,c("P.adj.high","P.adj.low")])), ]

topList = res_combined[res_combined$color == "Both Signf","ligand_receptor"]
completeList = res_combined[res_combined$color != "NS","ligand_receptor"]
completeList = unique(completeList[completeList %in% LR_with_expression_in_example])

matches <- unique(grep(paste(ADgene_list,collapse="|"), completeList, value=TRUE))

enrichments_combined_df_signif = enrichments_combined_df[enrichments_combined_df$p_value<=0.05,]
genes_l1 = paste(unlist(enrichments_combined_df_signif[grep("autophagy",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l1 = unique(str_split(genes_l1,",")[[1]])
print(paste0("autophagy: ",paste0(genes_l1,collapse = ";")))
## [1] "autophagy: UBA52;RPS27A;UBB;UBC;HSPA8;VIM;HSP90AA1;GFAP"
genes_l2 = paste(unlist(enrichments_combined_df_signif[grep("ion binding",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l2 = unique(str_split(genes_l2,",")[[1]])
print(paste0("ion binding: ",paste0(genes_l2,collapse = ";")))
## [1] "ion binding: SPARC;SMAP1;LDLR;ACTR2;APOE;ZYX;RPS27A;APP;VLDLR;CALM3;ADCY8;CALM2;ACVR1;NOTCH1;GPC1;COL18A1;GNAS;PRNP;ARF1;NOTCH4;TGFBR1;BMPR1B;PKM;EGFR;SPP1;GNAI2;TYROBP;CALR;ITGA3;HSPA8;TF;ABCA1;S100A1;RYR1;SMAD3;HSP90AA1;GAS6;MERTK;TGFBR2;KCNAB1;ATP1A3;AGRN;EFEMP1;SLIT1;GRM7;RIPK1;TIMP3;CD320;ADAM17;GSTP1;NOTCH2;TIMP2;PSAP;KIT;CACNA1C;RTN4;NELL2;CD4;TRAF2;ADCY9;RTN4R;RYR2;NRXN3"
genes_l3 = paste(unlist(enrichments_combined_df_signif[grep("vesicle",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l3 = unique(str_split(genes_l3,",")[[1]])
print(paste0("vesicle: ",paste0(genes_l3,collapse = ";")))
## [1] "vesicle: SPARC;CTTN;LRPAP1;LDLR;ACTR2;APOE;UBA52;ZYX;GFRA1;MIF;CD44;RPS27A;VWF;APP;TSPAN15;UBB;CALM3;ADCY8;CALM2;UBC;AQP1;PDGFB;ARRB2;NOTCH1;GPC1;COL18A1;GNAS;CNTN1;PRNP;ARF1;TGFBR1;ACKR3;PKM;EGFR;SPP1;GNAI2;PLXNA1;TYROBP;CALR;ITGA3;ARPC5;HSPA8;TF;ABCA1;PTPN13;VIM;RYR1;HSP90AA1;GAS6;B2M;SERPINF1;WNT3;ATP1A3;AGRN;CTNNB1;EFEMP1;CD74;RIPK1;CD53;FCGRT;PDGFC;TIMP3;ADCYAP1R1;GSTP1;TIMP2;CLU;GLG1;PSAP;GSTO1;PDGFA;CD63;KIT;RTN4RL1;CD4;TRAF2;OPRM1;CD82;RTN4R;ACKR1;FLRT2;CD81;FN1;FGFRL1;CD9"
genes_l4 = paste(unlist(enrichments_combined_df_signif[grep("cell adhesion",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l4 = unique(str_split(genes_l4,",")[[1]])
print(paste0("cell adhesion: ",paste0(genes_l4,collapse = ";")))
## [1] "cell adhesion: PLXNB1;SEMA4D;ITGB1;EPHA7;NOTCH1;JAG2;DSCAM;PTPRT;RET;TNFSF13B;TFRC;LRFN5;LAMA2;ITGAV;ITGB8;FBN1;EFNA5;PTPRM;DCC;PLXNA1;LAMB3;DSCAML1;SHH;EFNB1;TJP1;FLRT3;ITGA6;NRG1;EGFR;CNTN5;COL7A1;TGFB1;ERBB3;CD9;EPHB2;NECTIN3;PVR;PLXNB2;ITGB4;LAMA5;NRXN1;FYN;THY1;LAMB1;CD22;IGSF21;WNT5A;PTPRK;PDGFB;S1PR1;NRP2;NFASC;CD47;SIRPA;TNFRSF14;LRRC4;CDH6;SDC4;CADM4;CADM2;CADM3;LAMA4;LGALS9;FARP2;PLXNA4;CD58;CSF1;CELSR3;CNTN1;NLGN3;SEMA6A;THBS2;RTN4;LRRC4B;ENTPD1;LAMC1;PTPRS;NLGN1;DAG1;PTPRD;NRCAM;NRXN3;EPHB1;EFNB2;EPHB6;CCN3;BCAN;ERBB2;PTPRA;CNTN2;LAMA3;MELTF;NRXN2;CD151;NTN4;DDR1;ADGRL1;TNFSF9;GP1BB;L1CAM;CLDN11;BMP7;NLGN4X;COL18A1;ITGA3;CNTN4;APP;ADORA2A;AZGP1;NECTIN1;UNC5D;EPHA4;TYK2;JAK1;IGSF11;NECTIN2;EPHB3;EFNB3;PRKCZ;PLXNA2;AFDN;NEGR1;ITGB5;ROBO2;NRP1;ADAM10;UNC5C;FLRT2;ICOSLG;NTM;BMP4;NLGN2;NCAM1;NTNG2;STX4;ITGA7;CNTNAP1;VCAN;TNR;SCN1B;LAMB2;GPC4;FGFRL1;ADAM15;F11R;JAM3;CADM1;MAG;SIGLEC8;LSAMP;EFNA1;TNFRSF21;KITLG;SLITRK3;PTPRF;THBS3;CHL1;APLP1;MOG;EPHA3;CD81;ARF6;CNTNAP2;ADAM22;CDH13;NCAM2;ALCAM;OPCML;PAPLN;CD200;ADAM23;PDGFRA;VEGFA;LRRC4C;SEMA5A;PARD3;SPON1;SLC7A1;CD4;FLOT1;AXL;THBS4;OMG;ANXA1;LRRN2;CXCL12;RGMB;JAM2;CD63;PIK3R1;ENG;NTNG1;COL16A1;ANOS1;MDK;ANXA2;ROBO1;MCAM;SLITRK1;COL5A3;NID1;ROBO3;COL6A1;EFEMP2;STX3;LRP6;LRFN4;CCL28;CD83;LAMC3;WNT10B;PLXNC1;NOTCH4;PECAM1;ADAM9;PLXNA3;ATP5F1B;PIK3CA;EPHA8;IL1RAPL1;EPHA2;DLL1;SPON2;ADGRB1;NTN1;SLITRK2;NEO1;MFGE8;PLXND1;PSEN1;SH2B3;VCAM1;EZR;ADAM17;ITGB2;ICAM1;IGF1;THBS1;ITGAM;KDR;FN1;SFRP2;ITGAX;FGF2;HMGB1;CALR;ITGA2B;FAP;PTPRZ1;ITGA1;ITGA8;ITGA5;TJP2;VWF;SPP1;ITGA9;BCAM;ICAM5;CD44;COL6A2;CD93;B2M;CXCR4;SERPINE1;CD74;TNC;BMP2;PRNP;DDR2;SELPLG;BMP6;MDGA1;SIGLEC10;SELL;TGFB2;CLDN5;RELN;CHRD;SDK2;PCDHB7;TEK;BMPR2;SPTBN2;ADAM11;SPTAN1;FGF1;NOTCH3;PTN;GFRA1;SEMA7A"
genes_l5 = paste(unlist(enrichments_combined_df_signif[grep("cell migration",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l5 = unique(str_split(genes_l5,",")[[1]])
print(paste0("cell migration: ",paste0(genes_l5,collapse = ";")))
## [1] "cell migration: ADAM17;GNAI2;F2R;ICAM1;IGFBP5;IGF1;THBS1;SEMA6D;KDR;LGALS9;SPARC;ENG;FN1;SFRP2;ITGAX;FLT1;PDGFD;CXCL14;CXCR4;FGF2;HMGB1;CALR;ITGA2B;SERPINE1;ITGAV;DCN;CD74;CD81;TNC;NBL1;BMP2;F2RL1;DDR2;LRP1;JAM3;IGFBP3;LAMA4;MIF;SDC4;DRD2;HGF;ITGA5;IL6R;PDGFB;APP;LAMA2;RTN4;TGFB2;PDGFA;CLDN5;ITGB1;CD47;LAMA5;IFITM1;FGFR1;C5AR1;TGFB1;RELN;CHRD;ITGA6;TLR4;TEK;VCAM1;CD44;ITGB2;VCAN;LAMC1;EDNRA;FAP;F11R;SELPLG;ITGA1;MDGA1;GPC3;SELL;ITGA9;MET;PLXNB1;SEMA4D;FGF7;NOTCH1;PTPRT;RET;TIE1;PTPRM;SEMA3F;PLXNA1;ACVR1B;SHH;TJP1;SEMA4B;EGFR;PTPRG;HBEGF;CD9;EPHB2;SEMA4G;PLXNB2;THY1;LAMB1;WNT5A;PTPRK;S1PR1;NRP2;SEMA4F;TNFRSF14;PLXNA4;CSF1;ERBB4;SEMA6A;CXCL16;ADORA1;FGF17;DAG1;PDGFC;PDGFRB;CCN3;BMPR2;BMPR1A;LAMA3;CD151;BMP7;ITGA3;GRN;UNC5D;EPHA4;SLIT2;PLXNA2;AFDN;NRP1;ADAM10;UNC5C;FLRT2;BMP4;NTNG2;CORT;STX4;TNR;SEMA5B;ADAM15;NRG3;EFNA1;LRIG2;TREM2;IL34;CSF1R;KITLG;NTRK3;ARF6;SEMA6C;CDH13;SEMA3A;SEMA4A;CD200;PDGFRA;VEGFA;SEMA5A;FGF1;THBS4;ANXA1;CXCL12;JAM2;PIK3R1;TGFBR1;NTNG1;IGF1R;MDK;FGF9;ROBO1;MCAM;PTN;SCARB1;ECM1;CCL28;PLXNC1;PECAM1;ADAM9;PLXNA3;GPI;SEMA3B;SORL1;ATP5F1B;FGF18;INSR;EPHA2;WNT5B;SEMA7A;ACVR1C;ADGRB1;NTN1;TBC1D24;MYLK;PLXND1;SEMA6B;EFNB1;DCC;GPC1;FLRT3;ITGB4;FYN;SIRPA;FZD3;NRCAM;EPHB1;EFNB2;CNTN2;NTN4;DDR1;L1CAM;EPHB3;PRKCZ;FGF13;ITGB5;SDC2;ITGA7;GPC4;PTPRF;CHL1;EPHA3;SDC3;AXL;HTR6;CD63;PIK3CA;GFRA1;BAMBI;EPHA8;NEO1;PSEN1;TGFBR3"
genes_l6 = paste(unlist(enrichments_combined_df_signif[grep("myelination",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l6 = unique(str_split(genes_l6,",")[[1]])
print(paste0("myelination: ",paste0(genes_l6,collapse = ";")))
## [1] "myelination: GPC1;TGFB1;CD9;ITGB4;NFASC;CNTN1;DAG1;ERBB2;CNTN2;SCN8A;PTPRZ1;CNTNAP1;JAM3;MAG;TNFRSF21;ADAM22;PARD3;LGI4;JAM2;PTN;NCSTN"
genes_fig5 = c("AQP4","SPARC","C3AR1","SYT2","GABRD","PVALB","SCN1B","EEF1A1","RORB","ALDH1A3","TRIM22","SEMA3C","CD68","SERPINE2","ZNF561")

genes_from_pathways = unique(c(genes_l1,genes_l2,genes_l3,genes_l4,genes_l5,genes_l6,genes_fig5))
# top_n = 2
# genes_from_pathways = unique(c(genes_l1[1:top_n],genes_l2[1:top_n],genes_l3[1:top_n],genes_l4[1:top_n],genes_l5[1:top_n],genes_l6[1:top_n],genes_fig5[1:top_n]))
matches_from_pathways <- unique(grep(paste(genes_from_pathways,collapse="|"), completeList, value=TRUE))
 
genes2show = unique(c(matches_from_pathways))

#gene_list = c("SPARC—ENG","GFAP—APLNR","RPH3A—NRXN1","APP—TSPAN15","SEMA6D—PLXNC1","CLU—TREM2","APOE—LDLR","DSCAM—DCC")
gene_list = c(
  "SPARC—ENG", # ;vesicle;cell adhesion;cell migration;Fig5
  "GFAP—APLNR", # ;vesicle;cell adhesion;cell migration;AD genes
  "RPH3A—NRXN1", # ;cell adhesion
  "APP—TSPAN15", # ;vesicle;AD genes
  "SEMA6D—PLXNC1", # ;vesicle;cell adhesion;cell migration;AD genes
  "CLU—TREM2", # ;ion binding;cell adhesion;cell migration;AD genes
  "APOE—LDLR", # ;vesicle;cell adhesion;cell migration;AD genes
  "DSCAM—DCC", # ;cell adhesion;cell migration
  "PLXNA1—TREM2", # ;vesicle;cell adhesion;cell migration;myelination;AD genes
  "TGFB1—ITGAV", # ;cell adhesion;myelination
  "TNFRSF14—LRFN5", # ;autophagy;cell adhesion;cell migration
  "ADAM15—ITGAV", # ;autophagy;ion binding;cell adhesion;cell migration
  "GNB3—GABBR2",
  "RIMS1—SLC17A7",
  "CXCL16—GRM7",
  "APLN—GRM7",
  "CCL28—GRM7",
  "CORT—GRM7"
  )

topList1 = res_combined[res_combined$color == "Both Signf","ligand_receptor"]
topList2 = res_combined[res_combined$color == "glia-high vs. glia-low (low AB) only","ligand_receptor"]
topList3 = res_combined[res_combined$color == "glia-high vs. glia-low (high AB) only","ligand_receptor"]
topList4 = res_combined[res_combined$color != "NS" & res_combined$t.low<0 & res_combined$t.high>0,"ligand_receptor"]
topList5 = res_combined[res_combined$color == "glia-high vs. glia-low (low AB) only" & res_combined$t.low>0,"ligand_receptor"]
genes2show = unique(c(gene_list,topList1[1:2],topList2[1:3],topList3[1:3],topList4[1:6],topList5[1]))
res_combined$gene2label = ""
res_combined$gene2label[res_combined$ligand_receptor %in% genes2show] = res_combined$ligand_receptor[res_combined$ligand_receptor %in% genes2show]
res_combined$gene2label = gsub("—","-",res_combined$gene2label)

cluster_colors_AB <- c("#8B63A6","#8B63A6","#d7301f","#d7301f","#0570b0","#0570b0") # purple, red, blue
cluster_colors_glia <- c("#53B848","#53B848","#E8C51E","#E8C51E","#364eA1","#364eA1") # green, yellow, blue

signif_colors = setNames(rev(c("grey","#364eA1","#E8C51E","#53B848")), levels(res_combined$color))

(scatter_p = ggplot(res_combined, aes(x = t.low, y = t.high)) +
  geom_abline(color = "darkgray", linetype = "dashed") +
  geom_point(alpha = 0.25, aes(color = color), 
             show.legend = T) +
  geom_text_repel(aes(label = gene2label, color = color), 
                  show.legend = F, segment.alpha = .3,
                  max.overlaps = Inf, size = 3.5, force = 10, 
                  fontface = "bold.italic") +
  scale_color_manual(values = signif_colors) +
  stat_cor(method = "pearson", size = 4.5) +
  theme_test(base_size = 13) +
  guides(color = 
           guide_legend(label.position = "left", 
                        label.hjust = 1, 
                        override.aes = list(size = 4, alpha = 1, fill = NA))) +
  theme(legend.position = c(0.8, 0.12),
        legend.title.align = 1, 
        legend.background=element_rect(fill = alpha("white", 0)),
        legend.key.size = unit(0.45, "cm")) +
labs(x = "glia-high vs. glia-low (low AB) (t-statistics)", 
     y = "glia-high vs. glia-low (high AB) (t-statistics)", 
     title = "",
     color = "Significance"))

save(res_combined, file = paste0(work_dir,"/Figures/",figDir_id,"_GliaNet_LR_combined_results.RData"))

pdf(file = paste0(work_dir,"/Figures/",figDir_id,"_GliaNet_LR_low_high_scatter.pdf"), width = 8, height = 6)
scatter_p
dev.off()
## png 
##   2

Examples of LR pairs

res_combined_sel = res_combined %>% filter(gene2label != "")
res_combined_sel$source = ""
res_combined_sel$source[grep(paste(genes_l1,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l1,collapse="|"),res_combined_sel$gene2label)],";autophagy")
res_combined_sel$source[grep(paste(genes_l2,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l2,collapse="|"),res_combined_sel$gene2label)],";ion binding")
res_combined_sel$source[grep(paste(genes_l3,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l3,collapse="|"),res_combined_sel$gene2label)],";vesicle")
res_combined_sel$source[grep(paste(genes_l4,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l4,collapse="|"),res_combined_sel$gene2label)],";cell adhesion")
res_combined_sel$source[grep(paste(genes_l5,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l5,collapse="|"),res_combined_sel$gene2label)],";cell migration")
res_combined_sel$source[grep(paste(genes_l6,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_l6,collapse="|"),res_combined_sel$gene2label)],";myelination")
res_combined_sel$source[grep(paste(genes_fig5,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(genes_fig5,collapse="|"),res_combined_sel$gene2label)],";Fig5")
res_combined_sel$source[grep(paste(ADgene_list,collapse="|"),res_combined_sel$gene2label)] = paste0(res_combined_sel$source[grep(paste(ADgene_list,collapse="|"),res_combined_sel$gene2label)],";AD genes")
createDT(res_combined_sel)
gene_list2 = c("GFAP-APLNR",
                    "SPARC-ENG",
                    "CLU-TREM2",
                    "RPH3A-NRXN1",
                    "APP-TSPAN15",
                    "PLXNA1-TREM2",
                    "DSCAM-DCC",
                    "TGFB1-ITGAV",
                    "APOE-LDLR",
               "GNB3—GABBR2",
               "RIMS1—SLC17A7",
               "CXCL16—GRM7")
gene_list2 = gsub("-","—",gene_list2)

library(tidyverse)
library(data.table)
library(reshape2)
prep4boxplot = function(counts_voom,contrast = ""){
  pseudo_bulk_mtx = counts_voom %>% as.data.frame()
  mtx_1_long = setNames(melt(as.matrix(pseudo_bulk_mtx)),c("gene","sample","expression")) %>% 
    mutate(strat = gsub("(.*?)_(.*?)_(.*)","\\3",sample),
           sample_section = gsub("(.*?)_(.*?)_(.*)_(.*)","\\1_\\2",sample))
  mtx_1_long$contrast = contrast
  return(mtx_1_long)
}
mtx_voom_long = bind_rows(prep4boxplot(res_mean_RGNpos_RGNneg$counts_voom, contrast = "glia-high vs. glia-low (combined)"),
  prep4boxplot(res_mean_RGNpos_RGNneg_low$counts_voom, contrast = "glia-high vs. glia-low (low AB)"),
  prep4boxplot(res_mean_RGNpos_RGNneg_high$counts_voom, contrast = "glia-high vs. glia-low (high AB)"))

mtx_voom_long2plot_Pvalues = res_df_model %>%
  dplyr::filter(ligand_receptor %in% gene_list2) %>%
  select(ligand_receptor,P.adj,model) %>%
  pivot_wider(names_from = model,values_from = P.adj)
colnames(mtx_voom_long2plot_Pvalues)[colnames(mtx_voom_long2plot_Pvalues) == "glia-high vs. glia-low"] = "combined"
colnames(mtx_voom_long2plot_Pvalues)[colnames(mtx_voom_long2plot_Pvalues) == "glia-high vs. glia-low (high Abeta)"] = "high AB"
colnames(mtx_voom_long2plot_Pvalues)[colnames(mtx_voom_long2plot_Pvalues) == "glia-high vs. glia-low (low Abeta)"] = "low AB"
mtx_voom_long2plot_Pvalues_long = mtx_voom_long2plot_Pvalues %>%
  pivot_longer(cols = -ligand_receptor, names_to = "contrast", values_to = "P.adj")

mtx_voom_long2plot = mtx_voom_long %>% 
  dplyr::filter(gene %in% gene_list2) 

mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="RGN_neg"] <- "glia-low (combined)"
mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="RGN_pos"] <- "glia-high (combined)"
mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="low_RGN_neg"] <- "glia-low (low AB)"
mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="low_RGN_pos"] <- "glia-high (low AB)"
mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="high_RGN_neg"] <- "glia-low (high AB)"
mtx_voom_long2plot$strat[mtx_voom_long2plot$strat=="high_RGN_pos"] <- "glia-high (high AB)"

mtx_voom_long2plot$contrast = gsub("(.*) \\((.*)\\)","\\2",mtx_voom_long2plot$contrast)
mtx_voom_long2plot$contrast_strat = paste0(mtx_voom_long2plot$contrast,"_",mtx_voom_long2plot$strat)

mtx_voom_long2plot$contrast_strat_labels = gsub("(.*)_(.*) \\((.*)","\\2",mtx_voom_long2plot$contrast_strat)
contrast_strat_labels = mtx_voom_long2plot$contrast_strat_labels
names(contrast_strat_labels) = mtx_voom_long2plot$contrast_strat

cluster_colors_AB <- c("#8B63A6","#8B63A6","#d7301f","#d7301f","#0570b0","#0570b0") # purple, red, blue
cluster_colors_glia <- c("#53B848","#53B848","#E8C51E","#E8C51E","#364eA1","#364eA1") # green, yellow, blue

signif_colors = setNames(rev(c("#364eA1","#E8C51E","#53B848")), unique(mtx_voom_long2plot$contrast))

mtx_voom_long2plot = mtx_voom_long2plot %>% 
  left_join(mtx_voom_long2plot_Pvalues_long, by = c("contrast","gene"="ligand_receptor")) 
# convert pvalues to asterisk code
mtx_voom_long2plot$signif = ifelse(mtx_voom_long2plot$P.adj <= 0.001,"***",
                                   ifelse(mtx_voom_long2plot$P.adj <= 0.01,"**",
                                          ifelse(mtx_voom_long2plot$P.adj <= 0.05,"*","")))

mtx_voom_long2plot_labels = mtx_voom_long2plot %>% 
  dplyr::group_by(gene,contrast,contrast_strat) %>% 
  dplyr::summarize(p.adj = min(P.adj)) %>%
  mutate(signif = ifelse(p.adj <= 0.001,"***",
                         ifelse(p.adj <= 0.01,"**",
                                ifelse(p.adj <= 0.05,"*",""))))
mtx_voom_long2plot$gene = gsub("—","-",mtx_voom_long2plot$gene)
mtx_voom_long2plot_labels$gene = gsub("—","-",mtx_voom_long2plot_labels$gene)
mtx_voom_long2plot_labels$gene = factor(mtx_voom_long2plot_labels$gene, 
          levels = gsub("—","-",gene_list2))
mtx_voom_long2plot$gene = factor(mtx_voom_long2plot$gene, 
          levels = gsub("—","-",gene_list2))

# mtx_voom_long2plot_labels$`.y.` = "expression"
# mtx_voom_long2plot_labels$group1 = "glia-high"
# mtx_voom_long2plot_labels$group1 = "glia-low"
# mtx_voom_long2plot_labels$p = mtx_voom_long2plot_labels$p.adj
# mtx_voom_long2plot_labels$`p.signif` = mtx_voom_long2plot_labels$signif
# mtx_voom_long2plot_labels$contrast_strat = gsub("(.*)_(.*) \\((.*)","\\2",mtx_voom_long2plot_labels$contrast_strat)
# mtx_voom_long2plot_labels = unique(mtx_voom_long2plot_labels)

mtx_voom_long2plot$contrast = factor(mtx_voom_long2plot$contrast, 
                                     levels = c("combined","low AB","high AB"))
mtx_voom_long2plot$contrast_strat = factor(mtx_voom_long2plot$contrast_strat, 
                                           levels = c("combined_glia-high (combined)",
                                                      "combined_glia-low (combined)",
                                                      "low AB_glia-high (low AB)",
                                                      "low AB_glia-low (low AB)",
                                                      "high AB_glia-high (high AB)",
                                                      "high AB_glia-low (high AB)"))

(p = ggplot(mtx_voom_long2plot, 
            aes(x = contrast_strat, y = expression, fill = contrast)) +
  geom_boxplot(outlier.shape = NA) + 
  ggbeeswarm::geom_quasirandom(alpha = .2) + 
  # stat_pvalue_manual(mtx_voom_long2plot_labels, label = "p.signif") +
  geom_text(data = mtx_voom_long2plot_labels %>% filter(grepl("glia-high",contrast_strat)),
            aes(label = signif,
                x = stage(contrast_strat, after_stat = x + 0.45),
                y = Inf),
            color = "black",
            hjust = 0.5, vjust = 1, size = 20/.pt) +
  scale_x_discrete(labels = contrast_strat_labels) + 
  scale_y_continuous(breaks = scales::breaks_extended(3)) +
  scale_fill_manual(values = signif_colors) +
  facet_wrap(~gene, ncol = 3, scales = "free") + 
  theme_classic(base_size = 20) +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=20, colour="black",face = "italic"),
        strip.background = element_rect(colour="white", fill="white")) +
  easy_rotate_x_labels(angle = 20, side = "right") +
  labs(y = "Normalized expression", x = "", fill = "") +
  easy_move_legend(to = "bottom") + 
  easy_change_legend(what = "direction", to = "horizontal"))

pdf(file = paste0(work_dir,"/Figures/",figDir_id,"_GliaNet_LR_boxplots.pdf"), width = 14, height = 12)
p
dev.off()
## png 
##   2

Session info

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] readxl_1.4.3             ggrepel_0.9.5            factoextra_1.0.7        
##  [4] FactoMineR_2.11          bestNormalize_1.9.1      sva_3.42.0              
##  [7] genefilter_1.76.0        mgcv_1.9-1               nlme_3.1-153            
## [10] edgeR_3.36.0             ggfortify_0.4.17         gridExtra_2.3           
## [13] paletteer_1.6.0          gtools_3.9.5             RColorBrewer_1.1-3      
## [16] circlize_0.4.16          ComplexHeatmap_2.15.4    ggstatsplot_0.12.4      
## [19] furrr_0.3.1              future_1.34.0            variancePartition_1.24.1
## [22] BiocParallel_1.28.3      limma_3.50.3             lme4_1.1-35.1           
## [25] Matrix_1.6-5             yardstick_1.3.1          workflowsets_1.1.0      
## [28] workflows_1.1.4          tune_1.2.1               rsample_1.2.1           
## [31] recipes_1.1.0            parsnip_1.2.1            modeldata_1.4.0         
## [34] infer_1.0.7              dials_1.2.1              scales_1.3.0            
## [37] broom_1.0.6              tidymodels_1.2.0         reshape2_1.4.4          
## [40] kableExtra_1.4.0         rstatix_0.7.2            NICHES_1.0.0            
## [43] downloadthis_0.4.0       HGNChelper_0.8.14        data.table_1.15.4       
## [46] ggeasy_0.1.4             ggthemes_5.1.0           ggsci_3.2.0             
## [49] ggpubr_0.6.0             janitor_2.2.0            lubridate_1.9.3         
## [52] forcats_1.0.0            stringr_1.5.1            dplyr_1.1.4             
## [55] purrr_1.0.2              readr_2.1.5              tidyr_1.3.1             
## [58] tibble_3.2.1             ggplot2_3.5.1            tidyverse_2.0.0         
## [61] Seurat_5.1.0             SeuratObject_5.0.2       sp_2.1-4                
## 
## loaded via a namespace (and not attached):
##   [1] pbapply_1.7-2          lattice_0.20-45        vctrs_0.6.5           
##   [4] blob_1.2.4             survival_3.7-0         prodlim_2024.06.25    
##   [7] spatstat.data_3.1-2    later_1.3.2            nloptr_2.1.1          
##  [10] DBI_1.2.3              rappdirs_0.3.3         uwot_0.2.2            
##  [13] zlibbioc_1.40.0        htmlwidgets_1.6.4      mvtnorm_1.2-5         
##  [16] GlobalOptions_0.1.2    leaps_3.2              leiden_0.4.3.1        
##  [19] parallel_4.1.2         pbkrtest_0.5.2         irlba_2.3.5.1         
##  [22] Rcpp_1.0.13            KernSmooth_2.23-20     DT_0.33               
##  [25] promises_1.3.0         magick_2.8.4           statsExpressions_1.5.5
##  [28] RSpectra_0.16-2        RhpcBLASctl_0.23-42    digest_0.6.36         
##  [31] png_0.1-8              sctransform_0.4.1      cowplot_1.1.3         
##  [34] pkgconfig_2.0.3        venn_1.12              spatstat.random_3.3-1 
##  [37] gower_1.0.1            ggbeeswarm_0.7.2       estimability_1.5.1    
##  [40] iterators_1.0.14       minqa_1.2.7            GPfit_1.0-8           
##  [43] reticulate_1.38.0      spam_2.10-0            beeswarm_0.4.0        
##  [46] GetoptLong_1.0.5       xfun_0.46              bslib_0.8.0           
##  [49] zoo_1.8-12             tidyselect_1.2.1       ica_1.0-3             
##  [52] viridisLite_0.4.2      rlang_1.1.4            jquerylib_0.1.4       
##  [55] glue_1.7.0             lhs_1.2.0              matrixStats_1.3.0     
##  [58] emmeans_1.10.3         multcompView_0.1-10    lava_1.8.0            
##  [61] ggsignif_0.6.4         bayestestR_0.14.0      labeling_0.4.3        
##  [64] httpuv_1.6.15          class_7.3-19           TH.data_1.1-2         
##  [67] annotate_1.72.0        jsonlite_1.8.8         XVector_0.34.0        
##  [70] bit_4.0.5              mime_0.12              systemfonts_1.1.0     
##  [73] gplots_3.1.3.1         stringi_1.8.4          insight_0.20.2        
##  [76] processx_3.8.4         spatstat.sparse_3.1-0  scattermore_1.2       
##  [79] spatstat.explore_3.3-1 rbibutils_2.2.16       hardhat_1.4.0         
##  [82] bitops_1.0-8           cli_3.6.3              Rdpack_2.6            
##  [85] RSQLite_2.3.7          pheatmap_1.0.12        correlation_0.8.5     
##  [88] timechange_0.3.0       rstudioapi_0.16.0      butcher_0.3.4         
##  [91] locfit_1.5-9.10        listenv_0.9.1          miniUI_0.1.1.1        
##  [94] BiocGenerics_0.40.0    lifecycle_1.0.4        timeDate_4032.109     
##  [97] munsell_0.5.1          cellranger_1.1.0       caTools_1.18.2        
## [100] codetools_0.2-18       websocket_1.4.2        coda_0.19-4.1         
## [103] Biobase_2.54.0         GenomeInfoDb_1.30.1    vipor_0.4.7           
## [106] lmtest_0.9-40          admisc_0.35            ggpp_0.5.8-1          
## [109] xtable_1.8-4           ROCR_1.0-11            flashClust_1.01-2     
## [112] scatterplot3d_0.3-44   abind_1.4-5            farver_2.1.2          
## [115] parallelly_1.38.0      RANN_2.6.1             RcppAnnoy_0.0.22      
## [118] goftest_1.2-3          logger_0.3.0           patchwork_1.2.0       
## [121] cluster_2.1.2          future.apply_1.11.2    zeallot_0.1.0         
## [124] prettyunits_1.2.0      ggridges_0.5.6         igraph_2.0.3          
## [127] parameters_0.22.1      spatstat.utils_3.0-5   htmltools_0.5.8.1     
## [130] yaml_2.3.10            utf8_1.2.4             plotly_4.10.4         
## [133] XML_3.99-0.17          withr_3.0.0            fitdistrplus_1.2-1    
## [136] bit64_4.0.5            rngtools_1.5.2         effectsize_0.8.9      
## [139] splitstackshape_1.4.8  doRNG_1.8.6            multcomp_1.4-26       
## [142] foreach_1.5.2          Biostrings_2.62.0      progressr_0.14.0      
## [145] chromote_0.2.0         memoise_2.0.1          evaluate_0.24.0       
## [148] tzdb_0.4.0             ps_1.7.7               curl_5.2.1            
## [151] fansi_1.0.6            fastDummies_1.7.3      highr_0.11            
## [154] tensor_1.5             polynom_1.4-1          checkmate_2.3.2       
## [157] aod_1.3.3              cachem_1.1.0           deldir_2.0-4          
## [160] rjson_0.2.21           clue_0.3-65            tools_4.1.2           
## [163] sass_0.4.9             sandwich_3.1-0         magrittr_2.0.3        
## [166] RCurl_1.98-1.16        car_3.1-2              xml2_1.3.6            
## [169] httr_1.4.7             rmarkdown_2.27         boot_1.3-28           
## [172] globals_0.16.3         R6_2.5.1               nnet_7.3-16           
## [175] RcppHNSW_0.6.0         progress_1.2.3         KEGGREST_1.34.0       
## [178] shape_1.4.6.1          rematch2_2.1.2         splines_4.1.2         
## [181] snakecase_0.11.1       carData_3.0-5          colorspace_2.1-1      
## [184] generics_0.1.3         stats4_4.1.2           pillar_1.9.0          
## [187] GenomeInfoDbData_1.2.7 plyr_1.8.9             dotCall64_1.1-1       
## [190] gtable_0.3.5           rvest_1.0.4            knitr_1.48            
## [193] IRanges_2.28.0         fastmap_1.2.0          crosstalk_1.2.1       
## [196] Cairo_1.6-2            doParallel_1.0.17      datawizard_0.12.2     
## [199] AnnotationDbi_1.56.2   backports_1.5.0        S4Vectors_0.32.4      
## [202] ipred_0.9-15           spatstat.univar_3.0-0  vroom_1.6.5           
## [205] hms_1.1.3              Rtsne_0.17             shiny_1.8.1.1         
## [208] OmnipathR_3.9.8        polyclip_1.10-7        DiceDesign_1.10       
## [211] lazyeval_0.2.2         crayon_1.5.3           MASS_7.3-60.0.1       
## [214] svglite_2.1.3          rpart_4.1-15           compiler_4.1.2        
## [217] spatstat.geom_3.3-2