ST NICHES
A-beta LR comparison (GM only)
Rush Alzheimer’s Disease
Center
Chicago, IL 60612
Chicago, IL 60612
2023-12-12
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.idsLoad 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))))
high_Abeta_RGN_pos = unique(intersect(denis_strat$high_abeta_NEW, unique(denis_strat_m$barcode[denis_strat_m$strat=="RGN_pos"])))
denis_strat_m = bind_rows(denis_strat_m, data.frame(strat = "high_Abeta_RGN_pos",
barcode = high_Abeta_RGN_pos))
high_Abeta_RGN_neg = unique(intersect(denis_strat$high_abeta_NEW, unique(denis_strat_m$barcode[denis_strat_m$strat=="RGN_neg"])))
denis_strat_m = bind_rows(denis_strat_m, data.frame(strat = "high_Abeta_RGN_neg",
barcode = high_Abeta_RGN_neg))
low_Abeta_RGN_pos = unique(intersect(denis_strat$low_abeta_NEW, unique(denis_strat_m$barcode[denis_strat_m$strat=="RGN_pos"])))
denis_strat_m = bind_rows(denis_strat_m, data.frame(strat = "low_Abeta_RGN_pos",
barcode = low_Abeta_RGN_pos))
low_Abeta_RGN_neg = unique(intersect(denis_strat$low_abeta_NEW, unique(denis_strat_m$barcode[denis_strat_m$strat=="RGN_neg"])))
denis_strat_m = bind_rows(denis_strat_m, data.frame(strat = "low_Abeta_RGN_neg",
barcode = low_Abeta_RGN_neg))
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 = "; "))Load NICHES
niches_csv = list.files("/pastel/Github_scripts/SpeakEasy_ST/docs/cell_comms/niches_output/param_k4", pattern = ".csv", full.names = T)
niches_df = data.frame()
for (i in 1:length(niches_csv)){
filen = niches_csv[[i]]
samplen = gsub("(.*)/(.*)_k4\\.csv","\\2",filen)
samplen = gsub("_A1","_1",samplen)
samplen = gsub("_B1","_2",samplen)
samplen = gsub("_C1","_3",samplen)
samplen = gsub("_D1","_4",samplen)
df = fread(filen)
df_l = df %>%
pivot_longer(cols = -ligand_receptor, names_to = "barcode", values_to = "estimate") %>%
mutate(sample_section = samplen)
df_l = df_l[df_l$estimate>0,]
niches_df = bind_rows(niches_df, df_l)
}
save(niches_df, file = paste0(work_dir,"/NICHES/niches_df.RData"))Prepare gene expression data
load(paste0(work_dir,"/niches_df.RData"))
# > head(niches_df)
# ligand_receptor barcode estimate sample_section
# 1 CALM1—TRPC3 B16002_B16002_A1_GAAACATAGGAAACAG-1 1.223412 B16002_1
# 2 CALM1—TRPC3 B16002_B16002_A1_GTGAGGAGCGGTTGAG-1 1.717043 B16002_1
# 3 CALM1—TRPC3 B16002_B16002_A1_ATTAATTCGGTCACTC-1 1.641994 B16002_1
# 4 CALM1—TRPC3 B16002_B16002_A1_AACCTCGCTTTAGCCC-1 2.422880 B16002_1
# 5 CALM1—TRPC3 B16002_B16002_A1_CCCGCAAATAATCATC-1 3.321304 B16002_1
# 6 CALM1—TRPC3 B16002_B16002_A1_GTCCGAGAGCAATCAT-1 3.010470 B16002_1
niches_df_norm = suppressWarnings(get_niches_pseudo_bulk_norm2(niches_df = niches_df,
metadata.filt = metadata_with_strat,
pseudoBulkby = "strat",
selected_barcodes = NA,
pseudo_bulk_operation = "median", # Data used for the normalization
min_num_spots_to_include = 10,
normalization = "cpm",
na2zero = T))
save(niches_df_norm, file = paste0(work_dir,"/230927/RGN_NICHES_PseudoBulk_GM.RData"))Make gene background
niches_db = NICHES::LoadOmniPath("human")
ligand_complex_background = niches_db$source.subunits %>% as.data.frame() %>%
unite(source.complex, everything(), na.rm = T) %>% distinct() %>% remove_rownames()
receptor_complex_background = niches_db$target.subunits %>% as.data.frame() %>%
unite(target.subunits, everything(), na.rm = T) %>% distinct() %>% remove_rownames()
paste_sep <- function(x,y){return(paste(x,y, sep = "_"))}
lr_background_complex_background = outer(ligand_complex_background$source.complex,
receptor_complex_background$target.subunits,
FUN = "paste_sep")
dim(lr_background_complex_background) <- NULL
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 invalidanalisys_wrapper <- function(counts, metadata_4deg, p_val_cutoff = 0.05, evcodes=F){
########################################################################
# Filter genes and samples
#counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_RGNpos_RGNneg_high)]
metadata_4deg = as.data.frame(metadata_4deg, check.names = F)
counts_cpm <- as.data.frame(edgeR::cpm(counts))
keep.exp <- rowSums(counts_cpm > 1) >= (0.5 * ncol(counts_cpm) )
counts = counts[keep.exp,]
counts_cpm = counts_cpm[keep.exp,]
column2drop = colSums(counts)==0
message("Dropping ", sum(column2drop), " sections with 0 library size")
counts = counts[,which(!column2drop)]
metadata_4deg = metadata_4deg[colnames(counts),]
# Transform into log cpm
counts_log_cpm <- as.data.frame(edgeR::cpm(counts, log = TRUE))
########################################################################
# PCA
pca <- prcomp(t(counts_log_cpm), scale. = TRUE)
variances <- pca$sdev^2
explained <- variances / sum(variances)
pca_data <- cbind(metadata_4deg, pca$x[, 1:5])
# PC1 v. PC2
pc1vpc2 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = SpotSet)) +
geom_text(aes(label = SpotSet), size = 8, alpha = 0.5) +
labs(x = sprintf("PC%d (%.2f%%)", 1, round(explained[1] * 100, 2)),
y = sprintf("PC%d (%.2f%%)", 2, round(explained[2] * 100, 2))) +
theme(legend.position = "none")
pc1vpc2_points <- ggplot(pca_data, aes(x = PC1, y = PC2, color = SpotSet)) +
geom_point(size = rel(3)) +
scale_color_discrete(name = "Spot stratification") +
labs(x = sprintf("PC%d (%.2f%%)", 1, round(explained[1] * 100, 2)),
y = sprintf("PC%d (%.2f%%)", 2, round(explained[2] * 100, 2))) +
theme(legend.position = c(0.35, 0.8),
legend.background = element_rect(color = "black", size = 1,
linetype = 1))
pc1vpc2_points
########################################################################
# DE
design <- model.matrix(~ SpotSet, data = metadata_4deg)
colnames(design)[1] <- "Intercept"
y <- DGEList(counts)
y <- calcNormFactors(y)
v1 <- voom(y, design, normalize.method = "quantile")
corfit1 <- duplicateCorrelation(v1, design, block = metadata_4deg$Bnum)
# corfit1$consensus
v2 <- voom(y, design, normalize.method = "quantile",
block = metadata_4deg$Bnum, correlation = corfit1$consensus)
corfit2 <- duplicateCorrelation(v2, design, block = metadata_4deg$Bnum)
# corfit2$consensus
fit <- lmFit(v2, design, block = metadata_4deg$Bnum,
correlation = corfit2$consensus)
test <- eBayes(fit)
result <- topTable(test, coef = 2, number = nrow(counts))
########################################################################
# VOLCANO
res_df_All = result
res_df_All$ligand_receptor = rownames(res_df_All)
res_df_All$P.adj = result$adj.P.Val
res_df_All = res_df_All[order(res_df_All$P.adj),]
res_df_All = na.omit(res_df_All)
topList = res_df_All[p.adjust(res_df_All$P.adj, method = "fdr") <= 0.05,"ligand_receptor"]
matches <- unique(grep(paste(ADgene_list,collapse="|"), topList, value=TRUE))
genes2show = unique(c(res_df_All[res_df_All$logFC<0,]$ligand_receptor[1:5],
res_df_All[res_df_All$logFC>0,]$ligand_receptor[1:5],
matches))
p_volcano <- suppressWarnings(plot_volcano2(res = res_df_All, show_top_n_genes = 20, title = "RGN differential Ligand-Receptor",
adjP = T,
betaCol = "logFC",
pCol = "P.Value",
labelCol = "ligand_receptor",
genes2show = genes2show))
########################################################################
# ENRICHMENT
enrich_res_LR = run_enrichments(res_df_All,
plot_title_label = "RGN differential L-R",
beta_col = "logFC",
pval_col = "P.adj",
method = "ligand_receptor",
p_val_cutoff = p_val_cutoff,
ordered_query = F,
source2include = c("GO:CC","REAC","GO:MF","GO:BP"),
evcodes = evcodes,
top_n = 10,
background_genes = lr_background)
enrich_res_LR$gprofile_df$model = "LR"
enrich_res_L = run_enrichments(res_df_All, method = "ligand",
plot_title_label = "RGN differential L-R",
beta_col = "logFC",
pval_col = "P.adj",
p_val_cutoff = p_val_cutoff,
ordered_query = F,
source2include = c("GO:CC","REAC","GO:MF","GO:BP"),
evcodes = evcodes,
top_n = 10,
background_genes = ligand_background)
enrich_res_L$gprofile_df$model = "L"
enrich_res_R = run_enrichments(res_df_All, method = "receptor",
plot_title_label = "RGN differential L-R",
beta_col = "logFC",
pval_col = "P.adj",
p_val_cutoff = p_val_cutoff,
ordered_query = F,
source2include = c("GO:CC","REAC","GO:MF","GO:BP"),
evcodes = evcodes,
top_n = 10,
background_genes = receptors_background)
enrich_res_R$gprofile_df$model = "R"
enrich_res = bind_rows(enrich_res_LR$gprofile_df,enrich_res_L$gprofile_df,enrich_res_R$gprofile_df)
top_paths = enrich_res[enrich_res$source %in% c("GO:CC","REAC","GO:MF","GO:BP"),]
top_paths = top_paths[top_paths$isTop,
c("term_name","logP","model","Direction")]
top_paths$term_name = top_paths$term_name %>% str_trunc(50)
top_paths = top_paths[top_paths$logP > -log10(0.05),]
lr_test_plot = ggplot(top_paths, aes(x = model, y = term_name)) +
geom_tile(aes(fill = logP), colour="black") +
geom_text(aes(label = round(logP,digits = 1))) +
scale_fill_distiller(name = expression(-log[10](P-value)), palette = "RdPu", direction = 1) +
facet_grid(.~Direction, scales = "free_x", space = "free_x") +
theme_classic() + xlab('') + ylab('') +
theme(strip.text.x = element_text(size = 12, color = "black"), axis.text = element_text(size = 12, color = "black"))
########################################################################
# TOP GENES
top_3 = rownames(res_df_All)[1:3]
counts_long = counts[top_3,] %>%
rownames_to_column("ligand_receptor") %>%
pivot_longer(cols = -ligand_receptor, names_to = "sampleID", values_to = "counts")
counts_long$sample_section = gsub("(.*?)_(.*?)_(.*)","\\1_\\2",counts_long$sampleID)
counts_long$SpotSet = gsub("(.*?)_(.*?)_(.*)","\\3",counts_long$sampleID)
top1 = ggplot(counts_long, aes(x = SpotSet, y = counts)) +
geom_boxplot(outlier.shape = NA) + ggbeeswarm::geom_quasirandom() + facet_wrap(~ligand_receptor) + theme_bw() +
labs(title = "Counts")
counts_long = counts_cpm[top_3,] %>%
rownames_to_column("ligand_receptor") %>%
pivot_longer(cols = -ligand_receptor, names_to = "sampleID", values_to = "counts")
counts_long$sample_section = gsub("(.*?)_(.*?)_(.*)","\\1_\\2",counts_long$sampleID)
counts_long$SpotSet = gsub("(.*?)_(.*?)_(.*)","\\3",counts_long$sampleID)
top2 = ggplot(counts_long, aes(x = SpotSet, y = counts)) +
geom_boxplot(outlier.shape = NA) + ggbeeswarm::geom_quasirandom() + facet_wrap(~ligand_receptor) + theme_bw() +
labs(title = "CPM")
counts_long = as.data.frame(v2$E[top_3,]) %>%
rownames_to_column("ligand_receptor") %>%
pivot_longer(cols = -ligand_receptor, names_to = "sampleID", values_to = "counts")
counts_long$sample_section = gsub("(.*?)_(.*?)_(.*)","\\1_\\2",counts_long$sampleID)
counts_long$SpotSet = gsub("(.*?)_(.*?)_(.*)","\\3",counts_long$sampleID)
top3 = ggplot(counts_long, aes(x = SpotSet, y = counts)) +
geom_boxplot(outlier.shape = NA) + ggbeeswarm::geom_quasirandom() + facet_wrap(~ligand_receptor) + theme_bw() +
labs(title = "Quantile-voom")
top_genes_plot = suppressWarnings(ggarrange(top1,top2,top3, ncol = 1, nrow = 3))
########################################################################
res_list = list(
counts = counts,
counts_cpm = counts_cpm,
counts_log_cpm = counts_log_cpm,
counts_voom = as.data.frame(v2$E),
PCA_raw = pc1vpc2_points,
DE_result = res_df_All,
p_volcano = p_volcano,
enrich_res = enrich_res,
enrich_res_LR = enrich_res_LR,
enrich_res_L = enrich_res_L,
enrich_res_R = enrich_res_R,
lr_test_plot = lr_test_plot,
top_genes_plot = top_genes_plot
)
return(res_list)
}# 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_noAbeta_lowAbeta = metadata_pseudo %>% filter(SpotSet %in% c("no_Abeta","low_Abeta"))
metadata_4deg_noAbeta_lowAbeta$SpotSet = factor(metadata_4deg_noAbeta_lowAbeta$SpotSet, levels = c("no_Abeta","low_Abeta"))
metadata_4deg_noAbeta_highAbeta = metadata_pseudo %>% filter(SpotSet %in% c("no_Abeta","high_Abeta"))
metadata_4deg_noAbeta_highAbeta$SpotSet = factor(metadata_4deg_noAbeta_highAbeta$SpotSet, levels = c("no_Abeta","high_Abeta"))
metadata_4deg_lowAbeta_highAbeta = metadata_pseudo %>% filter(SpotSet %in% c("low_Abeta","high_Abeta"))
metadata_4deg_lowAbeta_highAbeta$SpotSet = factor(metadata_4deg_lowAbeta_highAbeta$SpotSet, levels = c("high_Abeta","low_Abeta"))
metadata_4deg_lowAbeta_highAbeta_RGN_pos = metadata_pseudo %>% filter(SpotSet %in% c("low_Abeta_RGN_pos","high_Abeta_RGN_pos"))
metadata_4deg_lowAbeta_highAbeta_RGN_pos$SpotSet = factor(metadata_4deg_lowAbeta_highAbeta_RGN_pos$SpotSet, levels = c("high_Abeta_RGN_pos","low_Abeta_RGN_pos"))
metadata_4deg_lowAbeta_highAbeta_RGN_neg = metadata_pseudo %>% filter(SpotSet %in% c("low_Abeta_RGN_neg","high_Abeta_RGN_neg"))
metadata_4deg_lowAbeta_highAbeta_RGN_neg$SpotSet = factor(metadata_4deg_lowAbeta_highAbeta_RGN_neg$SpotSet, levels = c("high_Abeta_RGN_neg","low_Abeta_RGN_neg"))res_mean_noAbeta_lowAbeta = analisys_wrapper(
counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_noAbeta_lowAbeta)] ,
metadata_4deg = metadata_4deg_noAbeta_lowAbeta, evcodes = T)
res_mean_noAbeta_highAbeta = analisys_wrapper(
counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_noAbeta_highAbeta)] ,
metadata_4deg = metadata_4deg_noAbeta_highAbeta, p_val_cutoff = 0.1, evcodes = T)
res_mean_lowAbeta_highAbeta = analisys_wrapper(
counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_lowAbeta_highAbeta)] ,
metadata_4deg = metadata_4deg_lowAbeta_highAbeta, evcodes = T)
res_mean_lowAbeta_highAbeta_RGN_pos = analisys_wrapper(
counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_lowAbeta_highAbeta_RGN_pos)] ,
metadata_4deg = metadata_4deg_lowAbeta_highAbeta_RGN_pos, evcodes = T)
res_mean_lowAbeta_highAbeta_RGN_neg = analisys_wrapper(
counts = niches_df_norm$mat_mean[,rownames(metadata_4deg_lowAbeta_highAbeta_RGN_neg)] ,
metadata_4deg = metadata_4deg_lowAbeta_highAbeta_RGN_neg, evcodes = T)
save(res_mean_noAbeta_lowAbeta, file = paste0(work_dir,"/",run_id,"/res_mean_noAbeta_lowAbeta_GM.RData"))
save(res_mean_noAbeta_highAbeta, file = paste0(work_dir,"/",run_id,"/res_mean_noAbeta_highAbeta_GM.RData"))
save(res_mean_lowAbeta_highAbeta, file = paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_GM.RData"))
save(res_mean_lowAbeta_highAbeta_RGN_pos, file = paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_RGN_pos_GM.RData"))
save(res_mean_lowAbeta_highAbeta_RGN_neg, file = paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_RGN_neg_GM.RData"))load(paste0(work_dir,"/",run_id,"/res_mean_noAbeta_lowAbeta_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_noAbeta_highAbeta_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_RGN_pos_GM.RData"))
load(paste0(work_dir,"/",run_id,"/res_mean_lowAbeta_highAbeta_RGN_neg_GM.RData"))NO/LOW/HIGH
res_df_model = res_mean_lowAbeta_highAbeta$DE_result %>% mutate(model = "LOWvsHIGH")
res_df_model = res_df_model %>% bind_rows(res_mean_noAbeta_lowAbeta$DE_result %>% mutate(model = "LOWvsNO"))
res_df_model = res_df_model %>% bind_rows(res_mean_noAbeta_highAbeta$DE_result %>% mutate(model = "HIGHvsNO"))
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)
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_noAbeta_lowAbeta$p_volcano + labs(title = "LOWvsNO"),
res_mean_noAbeta_highAbeta$p_volcano + labs(title = "HIGHvsNO"),
res_mean_lowAbeta_highAbeta$p_volcano + labs(title = "LOWvsHIGH"),
#res_mean_lowAbeta_highAbeta_RGN_pos$p_volcano + labs(title = "LOWvsHIGH (glia+)"),
#res_mean_lowAbeta_highAbeta_RGN_neg$p_volcano + labs(title = "LOWvsHIGH (glia-)"),
ncol = 3)
Top gene-pairs
ggarrange(
res_mean_noAbeta_lowAbeta$top_genes_plot,
res_mean_noAbeta_highAbeta$top_genes_plot,
res_mean_lowAbeta_highAbeta$top_genes_plot,
#res_mean_lowAbeta_highAbeta_RGN_pos$top_genes_plot,
#res_mean_lowAbeta_highAbeta_RGN_neg$top_genes_plot,
ncol = 3,
labels = c(
"LvsN"
,"HvsN"
,"LvsH"
#,"LvsH (glia+)"
#,"LvsH (glia-)"
),
hjust = -8)
Enrichments
# The adjusted p-values are reported in the gprofiler results.
enrichments_combined_df = bind_rows(
res_mean_noAbeta_lowAbeta$enrich_res %>% mutate(contrast = "LOW-NO"),
res_mean_noAbeta_highAbeta$enrich_res %>% mutate(contrast = "HIGH-NO"),
res_mean_lowAbeta_highAbeta$enrich_res %>% mutate(contrast = "LOW-HIGH")
#,res_mean_lowAbeta_highAbeta_RGN_pos$enrich_res %>% mutate(contrast = "LOW-HIGH (glia+)")
#,res_mean_lowAbeta_highAbeta_RGN_neg$enrich_res %>% mutate(contrast = "LOW-HIGH (glia-)")
) %>% 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_Abeta_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_noAbeta_lowAbeta$enrich_res_LR$enrich_p + labs(title = "LOWvsNO"),
res_mean_noAbeta_highAbeta$enrich_res_LR$enrich_p + labs(title = "HIGHvsNO"),
res_mean_lowAbeta_highAbeta$enrich_res_LR$enrich_p + labs(title = "LOWvsHIGH")
#,res_mean_lowAbeta_highAbeta_RGN_pos$enrich_res_LR$enrich_p + labs(title = "LOWvsHIGH (glia+)")
#,res_mean_lowAbeta_highAbeta_RGN_neg$enrich_res_LR$enrich_p + labs(title = "LOWvsHIGH (glia-)")
,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(
"LOW-NO"
,"LOW-HIGH"
#,"LOW-HIGH (glia+)"
#,"LOW-HIGH (glia-)"
#"HIGH-NO"
)]
enrich_selected_highlight = enrich_selected_highlight[,c(
"LOW-NO"
,"LOW-HIGH"
#,"LOW-HIGH (glia+)"
#,"LOW-HIGH (glia-)"
#"HIGH-NO"
)]
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-NO"] = "Low vs No Aβ"
#colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="HIGH-NO"] = "High vs No Aβ"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH"] = "Low vs High Aβ"
#colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH (glia+)"] = "Low vs High Aβ (glia+)"
#colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH (glia-)"] = "Low vs High Aβ (glia-)"
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"))
}
},
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),
#column_names_max_height = max_text_width(colnames(enrich_selected_highlight), gp = gpar(fontsize = 12)),
row_names_max_width = max_text_width(rownames(enrich_selected_highlight), gp = gpar(fontsize = 9))))
Cairo::CairoPDF(file = paste0(work_dir,"/Figures/",run_id,"_Abeta_LR_enrichment_heatmap.pdf"), family = "Arial", width = 4.5, height = 5)
enrich_heatmap
dev.off()## png
## 2Comparison between No vs Low vs High
res_combined = res_mean_noAbeta_lowAbeta$DE_result[,c("ligand_receptor","t","P.adj")] %>%
inner_join(res_mean_lowAbeta_highAbeta$DE_result[,c("ligand_receptor","t","P.adj")],
suffix = c(".no_low",".low_high"),
by = "ligand_receptor")
res_combined = na.omit(res_combined)
res_combined$color = "NS"
res_combined$color[res_combined$P.adj.no_low<=0.05 & res_combined$P.adj.low_high<=0.05] = "Both Signf"
res_combined$color[res_combined$P.adj.no_low<=0.05 & res_combined$P.adj.low_high>0.05] = "Low vs No AB differential only"
res_combined$color[res_combined$P.adj.no_low>0.05 & res_combined$P.adj.low_high<=0.05] = "Low vs High AB differential only"
res_combined$color = factor(res_combined$color,
levels = c("Both Signf",
"Low vs No AB differential only",
"Low vs High AB differential only",
"NS"))
res_combined = res_combined[order(rowMeans(res_combined[,c("P.adj.low_high","P.adj.no_low")])), ]
topList = res_combined[res_combined$color == "Both Signf","ligand_receptor"]
topList2 = res_combined[res_combined$color == "Low vs No AB differential only","ligand_receptor"]
topList3 = res_combined[res_combined$color == "Low vs High AB differential only","ligand_receptor"]
completeList = res_combined[res_combined$color != "NS","ligand_receptor"]
matches <- unique(grep(paste(ADgene_list,collapse="|"), unique(c(topList,topList2,topList3)), value=TRUE))
genes2show = unique(c(matches))
enrichments_combined_df_signif = enrichments_combined_df[enrichments_combined_df$p_value<=0.05,]
genes_l1 = paste(unlist(enrichments_combined_df_signif[grep("apoptotic",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l1 = unique(str_split(genes_l1,",")[[1]])
print(paste0("apoptotic: ",paste0(genes_l1,collapse = ";")))## [1] "apoptotic: SLC7A5;THBS1;TNFRSF11B;PRNP;VEGFB;KDR;ST6GAL1;SLIT2;SEMA3E;NRP1;FAP;TGFB1;WNT7A;HRAS;GRM4;ARF6;APP;ITGB1;FAS;CRKL;PLAUR;BMP2;UBB;TLR2;HSP90AA1;ERBB2;HSP90AB1;TGFB2;NGF;BEX3;ATP2A2;DCC;CLU;TREM2;IL6R;HSPA1A;GNAI2;FLNA;IGF1;LRP2;TGFBR2;SEMA4D;CAV1;APOE;CXCR4;SMAD3;RIPK1;AXL;KREMEN1;CD74;CD44;TJP1;PIK3R1;CD47;MIF;CALR;ANGPTL4;ADORA1;RTN4;APLP1;HMGB1;ACKR3;HNRNPK;IL6ST;GSTP1;TRAF2;TLR4;MEGF10;TGFBR1;AGT;TNFRSF1B;AR;EGFR;EPHA7;LRP6;HSP90B1;ITGAV;RAMP2;ITGA1;SERPINE1;DDR2;AQP1;CX3CL1;CX3CR1;GRN;RYR2;PDE1A;BMP7;ACVR1;KCNQ3;LTF;LGALS9;PPP2R2B;ROBO2;TNFRSF19;SORT1;DNAJA3;MAGED1"genes_l2 = paste(unlist(enrichments_combined_df_signif[grep("mapk",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l2 = unique(str_split(genes_l2,",")[[1]])
print(paste0("MAPK: ",paste0(genes_l2,collapse = ";")))## [1] "MAPK: F2R;GNAI2;C5AR1;THBS1;NTRK3;FN1;BMP2;LPAR3;VEGFB;KDR;ERBB2;F2RL1;BMP4;DRD2;APP;GRM4;CD81;IL6R;PDGFD;HRAS;HMGB1;TNFRSF19;IGF1;MIF;TIMP3;CD74;TGFB2;FGF2;CAV1;PTPRC;HAVCR2;FAS;IGFBP3;APOE;TEK;AR;NRP1;DAG1;IGFBP4;FZD8;HTR2C;AGT;ADRA2A;EPHA7;NRXN1;ADRB2;CRKL;IGFBP6;CSPG4;ADCYAP1;RAMP3;EPHA4;TGFB3;TLR4;CD4;PDGFC;RIPK1;PDGFRB;TREM2;FLT1;ITGA1;WNT7A;CD44;LILRB4;IGF1R;GHR;LGALS9;SEMA7A;CDON;FGF1;ACKR3;CSF1R;GSTP1;TRAF2;FBLN1;TF;MBP;ITGAV"genes_l3 = paste(unlist(enrichments_combined_df_signif[grep("extracellular matrix ",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l3 = unique(str_split(genes_l3,",")[[1]])
print(paste0("extracellular matrix: ",paste0(genes_l3,collapse = ";")))## [1] "extracellular matrix: LAMB3;ITGA3;ITGB1;LAMC1;ITGA7;LAMA4;LAMC3;LAMB1;ITGAV;ITGB8;LAMA3;DAG1;BMP4;BMP1;NCAM1;ITGA6;LRP4;LAMA2;BMP2;LAMA5;COL9A3;ITGA2B;PTPRS;SDC4;COL9A1;TNXB;ITGA5;COL7A1;ADAM10;BMP7;SDC2;F11R;JAM2;JAM3;PDGFB;COL4A2;TGFB3;LAMB2;BCAN;ITGB4;TGFB2;COL9A2;ADAM17;NID1;MMP2;EFEMP2;COL4A1;CD47;LTBP1;ITGB5;ITGB2;A2M;COL18A1;NRXN1;DDR1;COL11A1;COL4A5;ITGA9;ADAM15;CD151;HSPG2;COL1A2;VWF;DCN;KDR;TIMP1;AGRN;PECAM1;ICAM1;VCAN;COL6A2;COL5A1;SDC3;MMP15;FN1;FBN1;FGF2;DDR2;ADAM9;NCSTN;FBLN1;VTN;TNR;ICAM2;SPARC;PDGFA;NCAN;PSEN1;APP;P4HB;LAMA1;COL6A1;BGN;COL5A3;COL5A2;NTN4;TIMP2;MMP16;COL16A1;CD44;LRP1;TGFB1;NOTCH1;SMAD3;ANTXR1;TIE1;FGFR4;TNFRSF1A;MELTF"genes_l4 = paste(unlist(enrichments_combined_df_signif[grep("cell motility",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l4 = unique(str_split(genes_l4,",")[[1]])
print(paste0("cell motility: ",paste0(genes_l4,collapse = ";")))## [1] "cell motility: ITGA3;ITGB1;THBS4;LAMA4;FGF7;FGFR1;MET;PLXNB1;SEMA4D;ACVR1B;BMPR1A;BMPR2;ACVR1C;LAMB1;ITGAV;LAMA3;DAG1;PDGFRB;HBEGF;BMP4;TGFBR1;ITGA6;SEMA3F;PLXNA3;LAMA2;BMP2;LAMA5;EGFR;ACVRL1;CALR;LRP1;SEMA4A;NRP1;PLXNA4;ITGA2B;ENG;SDC4;SERPINE2;TNXB;WNT5B;CD81;ITGA5;TGFB1;FGF22;CSF1;ADAM10;MDK;PDGFC;WNT5A;ACVR1;SEMA4G;PLXNB2;SEMA4B;MCAM;PLXNA2;BMP7;JAM2;SEMA3A;EPHB2;JAG1;JAM3;SEMA5B;PLXNA1;PDGFB;S1PR1;IL6R;VEGFD;TNFRSF14;NRP2;NOTCH1;FGF18;EPHA4;SEMA5A;TGFB2;LGALS9;FGF17;THY1;FGF9;IGF1;IGF1R;APOE;SEMA3C;ADAM17;MMP2;EFNA1;CD9;INSR;TGFBR2;CD47;FLT1;FZD4;UNC5D;EPHA2;GNAI2;SORL1;RET;SMAD3;PTPRG;TLR4;ERBB4;HGF;CCN1;CAV1;INSL3;ADORA1;FLNA;TJP1;CLDN5;CSF1R;SEMA6D;ADAM15;CD151;CADM4;PTPRM;GRN;RTN4;TREM2;CXCL12;CXCL16;SHH;SEMA4F;NTN1;UNC5C;WNT7A;CD74;STX4;JAK2;DCN;KDR;SEMA3D;PLXND1;VEGFB;TIMP1;CD200;PTPRJ;PECAM1;SEMA3B;LRIG2;ICAM1;PGF;ACKR3;FN1;FGF2;SEMA6A;DDR2;ADAM9;TIE1;ADGRG1;MSN;PTN;ADRA2A;PTPRK;PLXNC1;SEMA7A;FBLN1;FADD;ROBO1;NTRK3;CD99L2;RACK1;RELN;SCARB1;VSIR;AQP1;PDGFRA;SLIT2;EDN1;C1QBP;FLRT2;VTN;PTPRT;TNR;SPARC;PDGFA;PIK3R1;VCL;APP;ARF6;TRADD;LAMA1;LGALS3;PODXL;SEMA6C;SFRP1;VEGFA;LAMC1;MERTK;EPHA3;DLL4;FZD3;ITGA7;DCC;EFNB1;CD63;SDC2;ITGB4;BAMBI;PTK7;EDNRB;AXL;IL16;EPHA8;C5;DDR1;EFNB2;SDC3;EGF;KIAA0319L;ITGB2;F11R;SEMA4C;CDH2;L1CAM;FGF1;EPHB4;TYRO3;ITGB5;IL34;EPHB3;GPC1;NEO1;SEMA6B;PLAT;PTPRF;AFDN;SELL;PIK3CA;CX3CL1;FLRT3;VCAN;NTNG2;ITGA9;APOA1;FAP;PRKCZ;NRCAM;CORT;EPHB1;PSEN1;CHL1;ADGRB1;GPC4;HTR6;GREM1;TMIGD3;CCN2;FGF13;RECK;SIRPA;CCN3;NTNG1;TNFSF12;GFRA1;EDN3;IRAK4;GAS6;NTN4;CNTN2;SELPLG;LRP8;PODXL2;EDNRA;LGALS8;MAPK1;SLIT1;CD44;LYN;TBC1D24;ECM1;HRAS;FLT4;SFTPD;IL33;COL5A1;ARF4;S100A9;FYN;FGFR4;GPC6;APCDD1;THBS1;DRD2;F2RL1;ROBO4;SEMA3E;ADORA3;SERPINF1;CRKL;NBL1;ITGAX;CXCL14;CXCR4;PLAU;MIF;CHRD;TF;HMGB1;BSG;AGT;IGFBP5;SERPINE1;APOD;CX3CR1;MYLK;C5AR1;MDGA1;GPC3;ARTN;HSP90AA1;IGFBP6;ARPC5;VCAM1;CLEC14A;MEGF10;ITGA1;RPS19;F2R;DPP4;PDGFD;IFITM1;PTPRC;IGFBP3;TEK;ITGA4;C3AR1;PTAFR;TNC;PRSS3;VSTM4;CSPG4;LRP5;PTPRO;PLTP"genes_l5 = paste(unlist(enrichments_combined_df_signif[grep("neurogenesis",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l5 = unique(str_split(genes_l5,",")[[1]])
print(paste0("neurogenesis: ",paste0(genes_l5,collapse = ";")))## [1] "neurogenesis: ITGA3;ITGB1;MET;PLXNB1;SEMA4D;TGFB1;EPHA3;EFNA5;DLL4;ITGA6;WNT7A;FZD3;LRP6;FGFR1;ALCAM;FLNA;GRM5;BMP4;BMPR1A;BMPR2;DCC;DSCAM;NRP1;EPHB2;CSF1;GDF7;EFNB1;ERBB4;SDC2;SEMA4A;PLXNA1;LAMA2;ITGB4;BMP2;LAMB2;DAG1;NLGN3;NRXN3;PTK7;EDNRB;AXL;FGFR2;FN1;WNT10B;EPHA8;EFNA3;BMP7;LAMB1;LRP1;EPHA4;NLGN1;NOTCH3;FGF2;EPHA7;DLL1;EGFR;GRM7;CXCL12;SDC4;DDR1;BMPR1B;ROBO3;NCAM1;EFNB2;EFNB3;NCAN;TREM2;SEMA6D;SLITRK2;PTPRD;SEMA4G;SEMA5A;OPRM1;PTPRS;KIAA0319L;UNC5A;THY1;WNT2B;NOTCH1;SLIT3;ROBO2;CNTN5;SEMA3A;SEMA4C;PLXNB2;LDLR;FLRT2;LRRC4C;PLXNA4;CDH2;L1CAM;PLXNA2;S1PR1;CNTN1;VEGFD;MDK;SEMA3F;ARF6;EPHA2;PLXND1;TYRO3;NRP2;SEMA4F;MAG;PTPRZ1;WNT5A;LAMC3;GDF11;IL34;FARP2;EPHB3;NLGN4X;NRXN1;SLITRK3;RIMS1;ROBO1;ERBB3;MMP2;SLIT2;GPC1;PTN;NTN1;NEO1;SHH;HHIP;APP;PARD3;NECTIN1;SEMA6B;PTPRG;JAG2;NOTCH2;WNT5B;PTPRF;NTM;NEGR1;ADAM17;SEMA6C;PTPRM;NFASC;UCN;JAG1;EFNA1;EPHA10;UNC5D;RET;STX3;SHANK1;CX3CL1;NLGN2;C3;FZD4;TGFBR1;AGRN;FLRT3;VCAN;TNR;LRP4;NTNG2;LGI1;UNC5C;ADCY1;PLXNA3;NPR2;ERBB2;RGMA;LGI4;VLDLR;LRTM2;RTN4R;RTN4;LEPR;JAK2;SEMA4B;PRKCZ;CALR;NRCAM;SCN1B;WNT3;ADGRV1;SORL1;TGM2;PTCH1;EPHB1;DLL3;RYK;SEMA3C;SEMA5B;TNXB;IL15RA;CSF1R;VEGFA;PSEN1;NCSTN;FZD1;EDN1;ADAM22;KCNJ10;ADM;CHL1;SERPINE2;FZD8;BCAN;ADGRB1;RTN4RL2;SLITRK4;EPHA6;SEMA3D;EPHA5;TLR4;PAPLN;IGF1R;SFRP1;WNT7B;DKK1;LRIG2;ATXN10;CD9;ADORA2A;GPR37;CTF1;IL6ST;FGF13;NTNG1;MXRA8;VAPA;GFRA1;DSCAML1;CHRM1;EDN3;GAS6;NTN4;APOE;CNTN2;CTHRC1;SOCS2;SEMA6A;NTRK3;TGFB2;CNTN4;UNC5B;TNFRSF21;CNTNAP2;CNTN6;LAMA1;DDR2;TCTN1;LRP8;PLXNC1;SEMA7A;SEMA3B;CNR1;EDNRA;SLITRK1;OPCML;APLP2;NDP;TLR2;MAPK1;SLIT1;WNT2;LYN;TBC1D24;BMP6;CNTNAP1;IL33;ARF4;GRN;S100A9;FYN;VTN;ANOS1;CEP290;KREMEN1;BDNF;ADGRG1;FGFR4;RELN;ANXA2;ADCYAP1;LRP2;VCL;FLOT1;RSPO2;APCDD1;IL1RAPL1;SLC7A5;DRD2;ROBO4;MDGA1;SEMA3E;SERPINF1;CRKL;NBL1;UBB;ARTN;HSP90AA1;HSP90AB1;NGF;CLU;RIMS2;LPAR3;VCAM1;SPP1;NPTN;CXCR4;S1PR5;SDK2;RTN4RL1;APLP1;HMGB1;CDON;OMG;TNFRSF1B;ITGA1;APOD;ISLR2;CX3CR1;SLC45A3;B2M;KCNQ3;SCARF1;C5AR1;C1QA;ITGA4;CSPG4;NPY;PTPRO;TNC;ETV5"genes_l6 = paste(unlist(enrichments_combined_df_signif[grep("synapse",tolower(enrichments_combined_df_signif$term_name)),]$intersection),collapse = ",")
genes_l6 = unique(str_split(genes_l6,",")[[1]])
print(paste0("synapse: ",paste0(genes_l6,collapse = ";")))## [1] "synapse: ITGA3;ITGB1;PLXNB1;SEMA4D;LAMA5;EFNA5;WNT7A;LRFN5;FLNA;GRM5;DSCAM;NRP1;EPHB2;PTPRT;ERBB4;ADAM10;SEMA4A;LAMB2;DAG1;COL4A5;COL4A1;NLGN3;NRXN3;PTK7;SHANK2;ADGRL1;THBS2;EPHA4;NLGN1;NRXN2;EPHA7;EFNB2;EFNB3;TREM2;SLITRK2;PTPRD;PTPRS;LRRC4B;ROBO2;CNTN5;PLXNB2;LRRTM2;FLRT2;LRRC4C;PLXNA4;CDH2;L1CAM;SEMA3F;ARF6;PLXND1;NRP2;LRRN1;WNT5A;EPHB3;NLGN4X;NRXN1;SLITRK3;LRFN4;NTN1;APP;NECTIN1;PTPRF;NEGR1;EFNA1;SHANK1;ADGRL2;CX3CL1;NLGN2;C3;AGRN;FLRT3;CDH6;TNR;LRP4;NTNG2;ERBB2;LRTM2;LRRC4;NRCAM;EPHB1;RYK;SPTBN2;C1QB;PSEN1;FZD1;LGI2;BCAN;ADGRB1;SLITRK4;GPC4;IGF1R;WNT7B;DKK1;FGF13;NTNG1;APOE;CNTN2;INSR;LRP8;LRRN3;PLXNC1;SLITRK1;GRIN2B;TLR2;IGSF21;SLIT1;CNTNAP1;IL1RAP;ARF4;SORT1;FYN;ARF1;GRIN2C;ADRA2A;SV2A;ITGA5;SEMA4F;NRG1;GRM2;PLAT;STX3;FZD4;LGI1;ADCY1;LIN7C;RTN4R;RTN4;JAK2;SEMA4B;PRKCZ;KCND2;GRM1;ADAM22;NETO2;ADORA2A;ADAM23;CHRM1;STX4;NPTX2;NPTXR;CNR1;HTR2A;LYN;HRAS;PDYN;BMPR2;DCC;CADM1;OPRD1;NECTIN3;ADORA1;GRM7;SEMA4C;CACNA1C;IGSF11;LRFN1;MPDZ;ATP1A3;FABP5;EZR;BDNF;RELN;IL1RAPL1;DRD2;PRNP;MDGA1;SEMA3E;CRKL;HSPA8;NPTN;SDK2;AMIGO2;ICAM5;ELFN1;CX3CR1;C5AR1;DNAJA3;GPC6;FLOT1;SCN8A;FZD3;LAMA4;EFNB1;LAMA2;GNA11;FGFR2;SV2B;NCAN;CNTN1;PTPRZ1;RIMS1;GNA15;CHRM3;CADM2;GPC1;PTN;UCN;KCNJ4;GABBR2;VCAN;ADAM11;UNC5C;GNAI2;NPR2;CORT;PTPRA;NCSTN;KCNJ10;SERPINE2;SV2C;HTR6;RPH3A;GPR37;DSCAML1;CNTNAP2;FGF12;OLFM2;CNTN6;LGI3;STX1A;TTYH1;MAPK1;TBC1D24;SLC3A2;HTR2C;GRIN2D;ADORA3;GAD1;GRM4;PLD1;FXYD6;ADCY8;NGF;ATP2A2;PORCN;CLU;RPS27A;DYSF;RIMS2;LPAR3;UBA52;PLD2;GRM3;SLC16A7;GNB1;MBP;KCNQ3;KCNQ5;RPS19;GLRA2;F2R;ITGAM;C1QA;PTPRO;TNC;SLC17A7;CRH;LAMP1;S1PR2;SNX14;NPY;TFR2;ITGA8;ADORA2B;YWHAQ;NTRK3"genes_from_pathways = unique(c(genes_l1,genes_l2,genes_l3,genes_l4,genes_l5,genes_l6))
matches_from_pathways <- unique(grep(paste(genes_from_pathways,collapse="|"), completeList, value=TRUE))
genes2show = unique(c(matches_from_pathways))res_combined$gene2label = ""
res_combined$gene2label[res_combined$ligand_receptor %in% genes2show] =
res_combined$ligand_receptor[res_combined$ligand_receptor %in% genes2show]
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)],";apoptosis")
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)],";MAPK")
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)],";extracellular matrix")
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 motility")
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)],";neurogenesis")
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)],";synapse")
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")
res_combined$gene2label = gsub("—","-",res_combined$gene2label)gene_list = c(
"APP—DCC", # ;apoptosis;MAPK;extracellular matrix;cell motility;neurogenesis;synapse;AD genes
"LRPAP1—SORL1", # ;cell motility;neurogenesis;AD genes
"APOE—SORL1", # ;apoptosis;MAPK;cell motility;neurogenesis;synapse;AD genes
"CD99—CD99L2", # ;cell motility;neurogenesis;AD genes
"C1QA—CSPG4", # ;MAPK;cell motility;neurogenesis;synapse;AD genes
"LRRC4B—PTPRF", # ;cell motility;neurogenesis;synapse;AD genes
"LRRC4B—PTPRS", # ;extracellular matrix;neurogenesis;synapse;AD genes
"TGFA—ADAM17", # ;extracellular matrix;cell motility;neurogenesis;AD genes
"HLA-DRA—CD63", # ;cell motility;AD genes
"HLA-C—CD81",
"HLA-A—APLP2",
"SPP1—ITGA9", # ;extracellular matrix;cell motility;neurogenesis
"ADAM10—IL6R", # ;apoptosis;MAPK;extracellular matrix;cell motility;synapse
"TYROBP—SEMA6D",
"APOE—LDLR",
"SPP1—ITGAV-ITGB5"
)
#res_combined_sel[res_combined_sel$gene2label %in% gene_list,c("ligand_receptor","source")]Selecting genes to show in the scatter plot…
topList1 = res_combined[res_combined$color == "Both Signf","ligand_receptor"]
topList2 = res_combined[res_combined$color == "Low vs No AB differential only","ligand_receptor"]
topList3 = res_combined[res_combined$color == "Low vs High AB differential only","ligand_receptor"]
topList4 = res_combined[res_combined$color == "Low vs High AB differential only" & res_combined$t.low_high>0,"ligand_receptor"]
topList5 = res_combined[res_combined$color != "NS" & res_combined$t.no_low>0 & res_combined$t.low_high>0,"ligand_receptor"]
topList6 = res_combined[res_combined$color == "Low vs No AB differential only" & res_combined$t.no_low>0,"ligand_receptor"]
genes2show = unique(c(gene_list,topList1[1:2],topList2[1:5],topList3[1:5],topList4[1:6],topList5[1:2],topList6[1:2]))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)
signif_colors = setNames(rev(c("grey","#3B4992FF","#EE0000FF","#631879FF")), levels(res_combined$color))
signif_colors = setNames(rev(c("grey","#d7301f","#0570b0","#8B63A6")), levels(res_combined$color))
(scatter_p = ggplot(res_combined, aes(x = t.no_low, y = t.low_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 = 3,
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.84, 0.15),
legend.title.align = 1,
legend.background=element_rect(fill = alpha("white", 0)),
legend.key.size = unit(0.45, "cm")) +
labs(x = "Low vs No AB differential (t-statistics)",
y = "Low vs High AB Differential (t-statistics)",
title = "",
color = "Significance"))
save(res_combined, file = paste0(work_dir,"/Figures/",run_id,"_Abeta_LR_no_low_high_scatter.RData"))
pdf(file = paste0(work_dir,"/Figures/",run_id,"_Abeta_LR_no_low_high_scatter.pdf"), width = 8, height = 6)
scatter_p
dev.off()## png
## 2LOW/HIGH - Glia +/-
res_df_model = res_mean_lowAbeta_highAbeta$DE_result %>% mutate(model = "LOWvsHIGH")
# res_df_model = res_df_model %>% bind_rows(res_mean_noAbeta_lowAbeta$DE_result %>% mutate(model = "LOWvsNO"))
# res_df_model = res_df_model %>% bind_rows(res_mean_noAbeta_highAbeta$DE_result %>% mutate(model = "HIGHvsNO"))
res_df_model = res_df_model %>% bind_rows(res_mean_lowAbeta_highAbeta_RGN_pos$DE_result %>% mutate(model = "LOWvsHIGH (glia+)"))
res_df_model = res_df_model %>% bind_rows(res_mean_lowAbeta_highAbeta_RGN_neg$DE_result %>% mutate(model = "LOWvsHIGH (glia-)"))
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)
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_noAbeta_lowAbeta$p_volcano + labs(title = "LOWvsNO"),
# res_mean_noAbeta_highAbeta$p_volcano + labs(title = "HIGHvsNO"),
res_mean_lowAbeta_highAbeta$p_volcano + labs(title = "LOWvsHIGH"),
res_mean_lowAbeta_highAbeta_RGN_pos$p_volcano + labs(title = "LOWvsHIGH (glia+)"),
res_mean_lowAbeta_highAbeta_RGN_neg$p_volcano + labs(title = "LOWvsHIGH (glia-)"),
ncol = 3)
Top gene-pairs
ggarrange(
# res_mean_noAbeta_lowAbeta$top_genes_plot,
# res_mean_noAbeta_highAbeta$top_genes_plot,
res_mean_lowAbeta_highAbeta$top_genes_plot,
res_mean_lowAbeta_highAbeta_RGN_pos$top_genes_plot,
res_mean_lowAbeta_highAbeta_RGN_neg$top_genes_plot,
ncol = 3,
labels = c(
# "LvsN",
# "HvsN",
"LvsH"
,"LvsH (glia+)"
,"LvsH (glia-)"
),
hjust = -8)
Enrichments
enrichments_combined_df = bind_rows(
# res_mean_noAbeta_lowAbeta$enrich_res %>% mutate(contrast = "LOW-NO"),
# res_mean_noAbeta_highAbeta$enrich_res %>% mutate(contrast = "HIGH-NO"),
res_mean_lowAbeta_highAbeta$enrich_res %>% mutate(contrast = "LOW-HIGH")
,res_mean_lowAbeta_highAbeta_RGN_pos$enrich_res %>% mutate(contrast = "LOW-HIGH (glia+)")
,res_mean_lowAbeta_highAbeta_RGN_neg$enrich_res %>% mutate(contrast = "LOW-HIGH (glia-)")
) %>% 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_Abeta_selected.xlsx"))
enrichments_combined_df$highlight = enrichments_combined_df$term_name%in%enrich_selected$term_name[enrich_selected$highlight=="T"]
createDT(enrichments_combined_df[enrichments_combined_df$p_value<0.05,])absmax <- function(x) { x[which.max( abs(x) )]}
absmin <- function(x) { x[which.min( abs(x) )]}
enrich_selected_highlight_pval = enrichments_combined_df[enrichments_combined_df$highlight,] %>%
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,] %>%
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(
#"LOW-NO",
"LOW-HIGH",
"LOW-HIGH (glia+)"
,"LOW-HIGH (glia-)"
#"HIGH-NO"
)]
enrich_selected_highlight = enrich_selected_highlight[,c(
#"LOW-NO",
"LOW-HIGH",
"LOW-HIGH (glia+)"
,"LOW-HIGH (glia-)"
#"HIGH-NO"
)]
#colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-NO"] = "Low vs No Aβ"
#colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="HIGH-NO"] = "High vs No Aβ"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH"] = "Low vs High Aβ"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH (glia+)"] = "Low vs High Aβ (glia+)"
colnames(enrich_selected_highlight)[colnames(enrich_selected_highlight)=="LOW-HIGH (glia-)"] = "Low vs High Aβ (glia-)"
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"))
}
},
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),
#column_names_max_height = max_text_width(colnames(enrich_selected_highlight), gp = gpar(fontsize = 12)),
row_names_max_width = max_text_width(rownames(enrich_selected_highlight), gp = gpar(fontsize = 9))))
Cairo::CairoPDF(file = paste0(work_dir,"/Figures/",run_id,"_Abeta_glia_LR_enrichment_heatmap.pdf"), family = "Arial", width = 4.7, height = 5)
enrich_heatmap
dev.off()## png
## 2Comparison between Low vs High (glia+/glia-)
res_combined =
res_mean_lowAbeta_highAbeta_RGN_pos$DE_result[,c("ligand_receptor","t","P.adj")] %>%
inner_join(res_mean_lowAbeta_highAbeta_RGN_neg$DE_result[,c("ligand_receptor","t","P.adj")],
suffix = c(".glia_pos",".glia_neg"),
by = "ligand_receptor")
res_combined = na.omit(res_combined)
res_combined$color = "NS"
res_combined$color[res_combined$P.adj.glia_pos<=0.05 & res_combined$P.adj.glia_neg<=0.05] = "Both Signf"
res_combined$color[res_combined$P.adj.glia_pos<=0.05 & res_combined$P.adj.glia_neg>0.05] = "Low vs High AB (glia-high) only"
res_combined$color[res_combined$P.adj.glia_pos>0.05 & res_combined$P.adj.glia_neg<=0.05] = "Low vs High AB (glia-low) only"
res_combined$color = factor(res_combined$color,
levels = c("Both Signf",
"Low vs High AB (glia-high) only",
"Low vs High AB (glia-low) only",
"NS"))
res_combined = res_combined[order(rowMeans(res_combined[,c("P.adj.glia_pos","P.adj.glia_neg")])), ]topList1 = res_combined[res_combined$color == "Both Signf","ligand_receptor"]
topList2 = res_combined[res_combined$color == "Low vs High AB (glia-high) only","ligand_receptor"]
topList3 = res_combined[res_combined$color == "Low vs High AB (glia-low) only","ligand_receptor"]
topList4 = res_combined[res_combined$color == "Low vs High AB (glia-low) only" & res_combined$t.glia_neg>0,"ligand_receptor"]
topList5 = res_combined[res_combined$color == "Both Signf" & res_combined$t.glia_neg>0,"ligand_receptor"]
topList6 = res_combined[res_combined$color != "NS" & res_combined$t.glia_neg>0 & res_combined$t.glia_pos<0,"ligand_receptor"]
topList7 = res_combined[res_combined$color == "Low vs High AB (glia-low) only" & res_combined$t.glia_neg>0 & res_combined$t.glia_pos>0,"ligand_receptor"]
topList8 = res_combined[res_combined$color == "Low vs High AB (glia-high) only" & res_combined$t.glia_pos>0,"ligand_receptor"]
# matches <- unique(grep(paste(ADgene_list,collapse="|"), unique(c(topList1,topList2,topList3)), value=TRUE))
# genes2show = unique(c(matches))
gene_list = gene_list[gene_list %in% topList1]
genes2show = unique(c(gene_list,topList1[1:3],topList2[1:3],topList3[1:3],topList4[1:3],topList5[1:3],topList6[1:3],topList7[1:5],topList8[1:3]))
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)signif_colors = setNames(rev(c("grey","#3B4992FF","#EE0000FF","#631879FF")), levels(res_combined$color))
signif_colors = setNames(rev(c("grey","#d7301f","#0570b0","#8B63A6")), levels(res_combined$color))
(scatter_p = ggplot(res_combined, aes(x = t.glia_neg, y = t.glia_pos)) +
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 = 5,
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.15),
legend.title.align = 1,
legend.background=element_rect(fill = alpha("white", 0)),
legend.key.size = unit(0.45, "cm")) +
labs(x = "Low vs High AB (glia-low) (t-statistics)",
y = "Low vs High AB (glia-high) (t-statistics)",
title = "",
color = "Significance"))
save(res_combined, file = paste0(work_dir,"/Figures/",run_id,"_Abeta_LR_low_high_scatter.RData"))
pdf(file = paste0(work_dir,"/Figures/",run_id,"_Abeta_LR_low_high_scatter.pdf"), width = 8, height = 6)
scatter_p
dev.off()## png
## 2Session info
## 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.4 factoextra_1.0.7
## [4] FactoMineR_2.9 bestNormalize_1.9.1 sva_3.42.0
## [7] genefilter_1.76.0 mgcv_1.9-0 nlme_3.1-153
## [10] edgeR_3.36.0 ggfortify_0.4.16 gridExtra_2.3
## [13] paletteer_1.5.0 gtools_3.9.5 RColorBrewer_1.1-3
## [16] circlize_0.4.15 ComplexHeatmap_2.10.0 ggstatsplot_0.12.1
## [19] furrr_0.3.1 future_1.33.0 variancePartition_1.24.1
## [22] BiocParallel_1.28.3 limma_3.50.3 lme4_1.1-35.1
## [25] Matrix_1.6-3 yardstick_1.2.0 workflowsets_1.0.1
## [28] workflows_1.1.3 tune_1.1.2 rsample_1.2.0
## [31] recipes_1.0.8 parsnip_1.1.1 modeldata_1.2.0
## [34] infer_1.0.5 dials_1.2.0 scales_1.3.0
## [37] broom_1.0.5 tidymodels_1.1.1 reshape2_1.4.4
## [40] kableExtra_1.3.4 rstatix_0.7.2 NICHES_1.0.0
## [43] downloadthis_0.3.3 HGNChelper_0.8.1 data.table_1.14.8
## [46] ggeasy_0.1.4 ggthemes_4.2.4 ggsci_3.0.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.4 tidyr_1.3.0
## [58] tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0
## [61] Seurat_5.0.0 SeuratObject_5.0.0 sp_2.1-1
##
## loaded via a namespace (and not attached):
## [1] ica_1.0-3 svglite_2.1.2 class_7.3-19
## [4] foreach_1.5.2 lmtest_0.9-40 crayon_1.5.2
## [7] rbibutils_2.2.16 MASS_7.3-60 backports_1.4.1
## [10] rlang_1.1.2 XVector_0.34.0 ROCR_1.0-11
## [13] irlba_2.3.5.1 nloptr_2.0.3 rjson_0.2.21
## [16] bit64_4.0.5 glue_1.6.2 pheatmap_1.0.12
## [19] rngtools_1.5.2 sctransform_0.4.1 vipor_0.4.5
## [22] pbkrtest_0.5.2 parallel_4.1.2 spatstat.sparse_3.0-3
## [25] AnnotationDbi_1.56.2 BiocGenerics_0.40.0 dotCall64_1.1-0
## [28] spatstat.geom_3.2-7 tidyselect_1.2.0 fitdistrplus_1.1-11
## [31] XML_3.99-0.14 zoo_1.8-12 butcher_0.3.3
## [34] xtable_1.8-4 RcppHNSW_0.5.0 magrittr_2.0.3
## [37] evaluate_0.23 zlibbioc_1.40.0 Rdpack_2.6
## [40] cli_3.6.2 doRNG_1.8.6 rstudioapi_0.15.0
## [43] miniUI_0.1.1.1 DiceDesign_1.9 logger_0.2.2
## [46] bslib_0.5.1 rpart_4.1-15 aod_1.3.2
## [49] polynom_1.4-1 fastDummies_1.7.3 shiny_1.7.5.1
## [52] xfun_0.41 clue_0.3-65 parameters_0.21.3
## [55] cluster_2.1.2 caTools_1.18.2 KEGGREST_1.34.0
## [58] venn_1.11 listenv_0.9.0 Biostrings_2.62.0
## [61] png_0.1-8 ipred_0.9-14 zeallot_0.1.0
## [64] withr_2.5.2 bitops_1.0-7 cellranger_1.1.0
## [67] plyr_1.8.9 hardhat_1.3.0 coda_0.19-4
## [70] pillar_1.9.0 gplots_3.1.3 GlobalOptions_0.1.2
## [73] cachem_1.0.8 multcomp_1.4-25 scatterplot3d_0.3-44
## [76] GetoptLong_1.0.5 vctrs_0.6.5 ellipsis_0.3.2
## [79] generics_0.1.3 lava_1.7.3 tools_4.1.2
## [82] beeswarm_0.4.0 munsell_0.5.0 emmeans_1.8.9
## [85] fastmap_1.1.1 compiler_4.1.2 abind_1.4-5
## [88] httpuv_1.6.12 plotly_4.10.3 GenomeInfoDbData_1.2.7
## [91] prodlim_2023.08.28 ggpp_0.5.5 lattice_0.20-45
## [94] deldir_1.0-9 utf8_1.2.4 later_1.3.1
## [97] jsonlite_1.8.7 pbapply_1.7-2 carData_3.0-5
## [100] estimability_1.4.1 lazyeval_0.2.2 promises_1.2.1
## [103] car_3.1-2 doParallel_1.0.17 goftest_1.2-3
## [106] checkmate_2.3.0 spatstat.utils_3.0-4 reticulate_1.34.0
## [109] rmarkdown_2.25 sandwich_3.0-2 cowplot_1.1.1
## [112] webshot_0.5.5 Rtsne_0.16 Biobase_2.54.0
## [115] uwot_0.1.16 igraph_1.5.1 survival_3.5-7
## [118] yaml_2.3.7 systemfonts_1.0.5 htmltools_0.5.7
## [121] memoise_2.0.1 locfit_1.5-9.8 IRanges_2.28.0
## [124] viridisLite_0.4.2 digest_0.6.33 RhpcBLASctl_0.23-42
## [127] rappdirs_0.3.3 leaps_3.1 mime_0.12
## [130] bayestestR_0.13.1 spam_2.10-0 RSQLite_2.3.3
## [133] future.apply_1.11.0 blob_1.2.4 S4Vectors_0.32.4
## [136] lhs_1.1.6 labeling_0.4.3 splines_4.1.2
## [139] rematch2_2.1.2 Cairo_1.6-1 RCurl_1.98-1.12
## [142] hms_1.1.3 colorspace_2.1-0 ggbeeswarm_0.7.2
## [145] shape_1.4.6 nnet_7.3-16 sass_0.4.7
## [148] Rcpp_1.0.11 RANN_2.6.1 mvtnorm_1.2-3
## [151] GPfit_1.0-8 multcompView_0.1-9 fansi_1.0.6
## [154] tzdb_0.4.0 parallelly_1.36.0 R6_2.5.1
## [157] ggridges_0.5.4 lifecycle_1.0.4 statsExpressions_1.5.2
## [160] datawizard_0.9.0 curl_5.1.0 ggsignif_0.6.4
## [163] minqa_1.2.6 leiden_0.4.3.1 jquerylib_0.1.4
## [166] snakecase_0.11.1 RcppAnnoy_0.0.21 TH.data_1.1-2
## [169] iterators_1.0.14 spatstat.explore_3.2-5 gower_1.0.1
## [172] htmlwidgets_1.6.2 polyclip_1.10-6 crosstalk_1.2.0
## [175] timechange_0.2.0 rvest_1.0.3 globals_0.16.2
## [178] flashClust_1.01-2 insight_0.19.6 patchwork_1.1.3
## [181] spatstat.random_3.2-1 progressr_0.14.0 codetools_0.2-18
## [184] matrixStats_1.1.0 prettyunits_1.2.0 RSpectra_0.16-1
## [187] GenomeInfoDb_1.30.1 correlation_0.8.4 gtable_0.3.4
## [190] DBI_1.1.3 stats4_4.1.2 highr_0.10
## [193] tensor_1.5 httr_1.4.7 KernSmooth_2.23-20
## [196] vroom_1.6.4 stringi_1.8.2 progress_1.2.2
## [199] farver_2.1.1 annotate_1.72.0 magick_2.8.1
## [202] timeDate_4022.108 DT_0.30 xml2_1.3.5
## [205] admisc_0.33 boot_1.3-28 OmnipathR_3.9.8
## [208] scattermore_1.2 bit_4.0.5 spatstat.data_3.0-3
## [211] pkgconfig_2.0.3 knitr_1.45