Long-term effects of adolescent risperidone treatment on the mouse cortex.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Materials and methods › SnRNA-seq statistics ↔ Antipsychotic_github.Rmd, lines 931–1055 · score 0.85 · FDR threshold, Cellular Component, EnrichGO, Biological Process, Molecular Function, CC
[2] § Results › SnRNA-seq analysis shows transcriptional changes associated with synaptic function in risperidone-treated mouse cortex ↔ Antipsychotic_github.Rmd, lines 2263–2404 · score 0.72 · log2FC, Upregulated genes, downregulated genes, GO enrichment, log10 adjusted, MF
[3] § Results › High-dimensional WGCNA identifies key gene modules associated with risperidone treatment ↔ Antipsychotic_github.Rmd, lines 3495–3592 · score 0.64 · co expression network, module assignment, hub gene, glutamatergic neurons, Dendogram, WGCNA
[4] § Results › High-dimensional WGCNA identifies key gene modules associated with risperidone treatment ↔ Antipsychotic_github.Rmd, lines 3384–3493 · score 0.61 · M10, GO terms, hdWGCNA, M8, M9, M6

Paper

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

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 5,888 lines · 195 KB · no license · 4 matches

---
title: "Antipsychotic_manuscript"
output: html_notebook
---
Libraries
```{r}
library(Seurat)
library(DESeq2)
library(clusterProfiler)
library(enrichR)
library(GeneOverlap)
library(WGCNA)
library(hdWGCNA)
library(ensembldb)
library(org.Mm.eg.db)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(tidyr)
library(Matrix)
library(cowplot)
library(patchwork)
library(tidyverse)
library(magrittr)
library(igraph)
library(pheatmap)
library(RColorBrewer)
library(openxlsx)
library(writexl)
```
Load h5 files
```{r}
counts_r1 <- Read_CellBender_h5_Mat(file_name = "~/Risperidone_Proj/Hahn_AA01_1/Hahn_AA01/count/Hahn_AA01_Risperdal_1/outs/cellbender_output/cellbender_filtered_feature_bc_matrix_filtered.h5")
counts_r2 <- Read_CellBender_h5_Mat(file_name = "~/Risperidone_Proj/Hahn_AA01_1/Hahn_AA01/count/Hahn_AA01_Risperdal_2/outs/cellbender_output/cellbender_filtered_feature_bc_matrix_filtered.h5")
counts_s1 <- Read_CellBender_h5_Mat(file_name = "~/Risperidone_Proj/Hahn_AA01_1/Hahn_AA01/count/Hahn_AA01_Saline_1/outs/cellbender_output/cellbender_filtered_feature_bc_matrix_filtered.h5")
counts_s2 <- Read_CellBender_h5_Mat(file_name = "~/Risperidone_Proj/Hahn_AA01_1/Hahn_AA01/count/Hahn_AA01_Saline_2/outs/cellbender_output/cellbender_filtered_feature_bc_matrix_filtered.h5")
```
Create Seurat object
```{r}
SO_r1 <- CreateSeuratObject(counts = counts_r1, project = "risperdal", min.cells = 3, min.features = 200)
SO_r2 <- CreateSeuratObject(counts = counts_r2, project = "risperdal", min.cells = 3, min.features = 200)
SO_s1 <- CreateSeuratObject(counts = counts_s1, project = "saline", min.cells = 3, min.features = 200)
SO_s2 <- CreateSeuratObject(counts = counts_s2, project = "saline", min.cells = 3, min.features = 200)
```
Merge Seurat objects
```{r}
Mcortex <- merge(SO_r1, y = c(SO_r2, SO_s1, SO_s2), add.cell.ids = c("R1", "R2", "S1", "S2"), project = "Mcortex")
# Create joint count matrix
Mcortex <- JoinLayers(Mcortex)
Layers(Mcortex[["RNA"]])
```
Quality Control
```{r}
# Add number of genes per UMI for each cell to metadata
Mcortex$log10GenesPerUMI <- log10(Mcortex$nFeature_RNA) / log10(Mcortex$nCount_RNA)
# Compute percent mito ratio
Mcortex$mitoRatio <- PercentageFeatureSet(object = Mcortex, pattern = "^mt-")
Mcortex$mitoRatio <- [email hidden]$mitoRatio / 100
# Visualize mitoRatio as a violin plot
VlnPlot(Mcortex, features = c("mitoRatio"), ncol = 3, raster=FALSE)
# Visualize the number of cell counts per sample
[email hidden] %>%
mutate(orig.ident = factor(orig.ident, levels = c("saline", "risperdal"))) %>%
ggplot(aes(x = orig.ident, fill = orig.ident)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.3, size = 3.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
ggtitle("NCells")
# Visualize the number UMIs/transcripts per cell
[email hidden] %>%
ggplot(aes(color=orig.ident, x=nCount_RNA, fill= orig.ident)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
ylab("Cell density") +
geom_vline(xintercept = 500)
# Visualize the distribution of genes detected per cell via histogram
[email hidden] %>%
ggplot(aes(color=orig.ident, x=nFeature_RNA, fill= orig.ident)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
geom_vline(xintercept = 300)
# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
[email hidden] %>%
ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=mitoRatio)) +
geom_point() +
scale_colour_gradient(low = "gray90", high = "black") +
stat_smooth(method=lm, aes(color = mitoRatio)) +
scale_x_log10() +
scale_y_log10() +
theme_classic() +
geom_vline(xintercept = 500) +
geom_hline(yintercept = 250) +
facet_wrap(~orig.ident)
# Visualize the distribution of mitochondrial gene expression detected per cell
[email hidden] %>%
ggplot(aes(color=orig.ident, x=mitoRatio, fill=orig.ident)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
geom_vline(xintercept = 0.2)
# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
[email hidden] %>%
ggplot(aes(x=log10GenesPerUMI, color = orig.ident, fill=orig.ident)) +
geom_density(alpha = 0.2) +
theme_classic() +
geom_vline(xintercept = 0.8)
#Cell-Level Filtering
# Filter out low quality reads using selected thresholds - these will change with experiment
Mcortex_QC <- subset(x = Mcortex,
subset= (nCount_RNA >= 500) &
(nFeature_RNA >= 250) &
(log10GenesPerUMI > 0.80) &
(mitoRatio < 0.20))
# Gene-Level Filtering
# Extract counts
counts_cell <- GetAssayData(object = Mcortex_QC, layer = "counts")
# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts_cell > 0
# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10
# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts_cell[keep_genes, ]
# Reassign to filtered Seurat object
Mcortex_QC <- CreateSeuratObject(filtered_counts, meta.data = [email hidden])
# Visualize Cell counts after filtering
[email hidden] %>%
mutate(orig.ident = factor(orig.ident, levels = c("saline", "risperdal"))) %>%
ggplot(aes(x = orig.ident, fill = orig.ident)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.3, size = 3.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
ggtitle("NCells post-filter")
```
Dimensionality reduction and normalization
```{r}
# Normalize and scale the data
Mcortex_QC <- SCTransform(Mcortex_QC, verbose = FALSE)
# Perform PCA
Mcortex_QC <- RunPCA(Mcortex_QC, features = VariableFeatures(object = Mcortex_QC))
```
Re-organize metadata
```{r}
# Create new metadata called Tx, and change saline to Control
Mcortex_QC$Tx <- ifelse(Mcortex_QC$orig.ident == "saline", "Control", "Risperidone")
# Make it into a factor for consistent ordering in plot
Mcortex_QC$Tx <- factor(Mcortex_QC$Tx, levels = c("Control", "Risperidone"))
```
Add metadata for predicted sex
```{r}
# read CSV file with predicted sex information
sex_info <- read.csv("~/Risperidone_Proj/Hahn_AA01_1/Hahn_AA01/aggr/Hahn_AA01/outs/predicted_sex.csv")
# Add prefix based on suffix
sex_info <- sex_info %>%
mutate(
new_barcode = paste0(
case_when(
str_detect(X, "-1$") ~ "R1_",
str_detect(X, "-2$") ~ "R2_",
str_detect(X, "-3$") ~ "S1_",
str_detect(X, "-4$") ~ "S2_",
TRUE ~ NA_character_
),
str_remove(X, "-[1-4]$"),
"-1"
)
)
# Check the new barcode
View(sex_info)
# Create a named vector for easy lookup
barcode_to_sex <- setNames(sex_info$Predicted_Sex, sex_info$new_barcode)
# Perform matching between Mcortex_QC and sex_info
Predicted_Sex_vec <- barcode_to_sex[match(colnames(Mcortex_QC), names(barcode_to_sex))]
# Find the barcodes from sex_info that are not found in Mcortex_QC
unmatched_barcodes <- sex_info[is.na(match(sex_info$new_barcode, colnames(Mcortex_QC))), ]
# Save the unmatched barcodes to a separate CSV file for record-keeping
write.csv(unmatched_barcodes, "unmatched_barcodes_from_RisperdalProj4.csv", row.names = FALSE)
# Find the shared barcodes between Mcortex_QC and sex_info
shared_barcodes <- intersect(colnames(Mcortex_QC), sex_info$new_barcode)
# Extract the Predicted_Sex values for the shared barcodes from sex_info
Predicted_Sex_shared <- sex_info$Predicted_Sex[match(shared_barcodes, sex_info$new_barcode)]
# Create a vector for Predicted_Sex with NA for unmatched barcodes
Predicted_Sex_vec <- rep(NA, ncol(Mcortex_QC)) # Default to NA for all cells
names(Predicted_Sex_vec) <- colnames(Mcortex_QC) # Ensure the names are the same as the column names of Mcortex_QC
# Assign the Predicted_Sex values to the shared barcodes
Predicted_Sex_vec[shared_barcodes] <- Predicted_Sex_shared
# Add the Predicted_Sex metadata column to Mcortex_QC
Mcortex_QC <- AddMetaData(Mcortex_QC, metadata = Predicted_Sex_vec, col.name = "Predicted_Sex")
# Check new metadata
View([email hidden])
```
Clustering and UMAP visualization
```{r}
# Determine the dimensionality of the dataset
ElbowPlot(Mcortex_QC)
Mcortex_QC <- FindNeighbors(Mcortex_QC, dims = 1:20)
Mcortex_QC <- FindClusters(Mcortex_QC, resolution = 0.5)
# Run non-linear dimensional reduction (UMAP)
Mcortex_QC <- RunUMAP(Mcortex_QC, dims = 1:20)
# visualize UMAP
DimPlot(Mcortex_QC, reduction = "umap", label = TRUE, raster = FALSE, pt.size = 0.1, alpha = 0.5) + ggtitle("")
```
Find cluster markers
```{r}
# find markers for every cluster compared to all remaining cells, report only the positive ones
cluster.markers <- FindAllMarkers(Mcortex_QC, only.pos = TRUE)
cluster.markers <- cluster.markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
arrange(desc(avg_log2FC))
write.xlsx(cluster.markers, file = "cluster_markers_Proj4_all.xlsx")
```
Assign cell type identity to clusters
```{r}
# Define cluster ID's
cluster.id2 <- c(
"0" = "Glutamatergic L2/3 IT",
"1" = "Astrocytes",
"2" = "Glutamatergic L6 CT",
"3" = "Glutamatergic L5 IT",
"4" = "Glutamatergic (unspecified)",
"5" = "GABAergic Vip",
"6" = "Glutamatergic L5 IT",
"7" = "Glutamatergic L2/3 IT",
"8" = "Glutamatergic L6 IT",
"9" = "GABAergic Pvalb",
"10" = "Oligodendrocytes",
"11" = "Astrocytes",
"12" = "Glutamatergic L5/6 NP",
"13" = "GABAergic Sst",
"14" = "Glutamatergic L5 ET",
"15" = "Undetermined",
"16" = "Oligodendrocytes",
"17" = "GABAergic Lamp5",
"18" = "Glutamatergic (unspecified)",
"19" = "Glutamatergic L6 CT",
"20" = "Glutamatergic L5/6 NP",
"21" = "Endothelial cells",
"22" = "Oligodendrocytes",
"23" = "Glutamatergic L6 IT Car3",
"24" = "Glutamatergic (unspecified)",
"25" = "Undetermined",
"26" = "Oligodendrocytes"
)
Mcortex_QC$celltype2 <- unname(cluster.id2[clusters_char])
# Visualize on UMAP
DimPlot(Mcortex_QC, reduction = "umap", group.by = 'celltype2', label = TRUE, raster = FALSE, repel = TRUE, pt.size = 0.1, alpha = 0.5) + ggtitle("")
```
Count the number of cells per cell type
```{r}
# NCell by cell type
[email hidden] %>%
ggplot(aes(x = celltype1, fill = celltype1)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.3, size = 3) +
theme_classic() +
theme(
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold", margin = margin(b = 20)),
legend.position = "none"
) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
ggtitle("NCells by cell type")
```
DE analysis (Pseudobulk) - Glut (all glutamatergic neurons)
```{r}
# Get the expression matrix
expr_mat <- GetAssayData(Mcortex_QC, assay = "SCT", layer = "data")
Glut_cells <- colnames(Mcortex_QC)[Mcortex_QC$celltype1 == "Glutamatergic neurons"]
# Subset Seurat object for all glutamatergic neurons
Mcortex_QC_Glut <- subset(Mcortex_QC, cells = Glut_cells)
# Filter out all genes with count sum <10
Mcortex_QC_Glut <- subset(Mcortex_QC_Glut, features = rownames(Mcortex_QC_Glut)[Matrix::rowSums(GetAssayData(Mcortex_QC_Glut, assay = "RNA", layer = "counts")) >= 10])
# Extract sample_id from barcode names
Mcortex_QC_Glut$sample_id <- sapply(strsplit(Cells(Mcortex_QC_Glut), "_"), `[`, 1)
# Aggregate counts
expr_mat_Glut <- AggregateExpression(
Mcortex_QC_Glut,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Glut <- data.frame(
sample_id = colnames(expr_mat_Glut),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Glut)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Glut) <- sample_info_Glut$sample_id
# Build DESeq2 object
dds_Glut <- DESeqDataSetFromMatrix(
countData = expr_mat_Glut,
colData = sample_info_Glut,
design = ~ condition
)
# Run DE analysis
dds_Glut <- DESeq(dds_Glut)
# Get results
res_Glut <- results(dds_Glut, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Glut) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Glut.csv"),
row.names = TRUE
)
```
Glut DEG - heatmap
```{r}
# Count the number of DEGs with padj < 0.05
sum(res_Glut$padj < 0.05, na.rm = TRUE) #9092
# Select all DEGs with padj < 0.05 and remove NA
res_Glut_filtered <- res_Glut[!is.na(res_Glut$padj) & res_Glut$padj < 0.05, ]
# Then order and extract gene names
top_genes_Glut <- rownames(res_Glut_filtered[order(res_Glut_filtered$padj), ])
# Perform Variance-Stabilizing Transformation (VST)
vst_data_Glut <- vst(dds_Glut, blind = TRUE)
# Ensure genes exist in the dataset
common_genes_Glut <- intersect(top_genes_Glut, rownames(assay(vst_data_Glut)))
# Extract the VST-transformed expression values for these genes
heatmap_matrix_Glut <- assay(vst_data_Glut)[common_genes_Glut, ]
# Extract metadata for annotations
annotation_Glut <- sample_info_Glut %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Glut$condition <- factor(annotation_Glut$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Glut <- rownames(annotation_Glut)[order(annotation_Glut$condition)]
heatmap_matrix_Glut <- heatmap_matrix_Glut[, ordered_samples_Glut]
annotation_colors_Glut <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
#Generate Heatmap
pheatmap(heatmap_matrix_Glut,
scale = "row", # Normalize each gene across samples
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Glut,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "",
)
```
Glut DEG - volcano
```{r}
# Convert DESeq2 results to a data frame
res_df_Glut <- as.data.frame(res_Glut)
res_df_Glut$gene <- rownames(res_df_Glut)
# Remove any NA values
res_df_Glut <- na.omit(res_df_Glut)
# Define significance threshold
padj_threshold <- 0.05
log2FoldChange_threshold <- 0.5
# Add a new column to categorize genes for coloring
res_df_Glut$Significance <- "Not Significant"
res_df_Glut$Significance[res_df_Glut$padj < padj_threshold & abs(res_df_Glut$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_Glut <- c("Kcnh1", "Kcnq2", "Kcnq3", "Kcnq5", "Kcnd3", "Kcnj6", "Kcnj9", "Grin2a", "Grin2b", "Gabra5", "Gabrb1", "Gabrb3", "Gabrg3", "Gabrd", "Shank1", "Shank2", "Dlg1", "Dlg2", "Dlg4")
# Generate Volcano plot
volcano_plot_Glut <- ggplot(res_df_Glut, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_df_Glut %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_df_Glut %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_df_Glut %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Glut, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "Glut",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 20)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
# Print the plot
print(volcano_plot_Glut)
```
DE analysis (Pseudobulk) - PV-IN
```{r}
# Subset cell type
Mcortex_QC_Pvalb <- subset(Mcortex_QC, subset = celltype1 == "GABAergic Pvalb")
# Filter out all genes with count sum <10
Mcortex_QC_Pvalb <- subset(Mcortex_QC_Pvalb, features = rownames(Mcortex_QC_Pvalb)[Matrix::rowSums(GetAssayData(Mcortex_QC_Pvalb, assay = "RNA", layer = "counts")) >= 10])
# Extract sample_id from barcode names
Mcortex_QC_Pvalb$sample_id <- sapply(strsplit(Cells(Mcortex_QC_Pvalb), "_"), `[`, 1)
# Aggregate counts
expr_mat_Pvalb <- AggregateExpression(
Mcortex_QC_Pvalb,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Pvalb <- data.frame(
sample_id = colnames(expr_mat_Pvalb),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Pvalb)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Pvalb) <- sample_info_Pvalb$sample_id
# Build DESeq2 object
dds_Pvalb <- DESeqDataSetFromMatrix(
countData = expr_mat_Pvalb,
colData = sample_info_Pvalb,
design = ~ condition
)
# Run DE analysis
dds_Pvalb <- DESeq(dds_Pvalb)
# Get results
res_Pvalb <- results(dds_Pvalb, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Pvalb) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_PV-IN.csv"),
row.names = TRUE
)
```
PV-IN DEG - heatmap
```{r}
# Count the number of DEGs with padj < 0.05
sum(res_Pvalb$padj < 0.05, na.rm = TRUE) #2298
# Select all DEGs with padj < 0.05 and remove NA
res_Pvalb_filtered <- res_Pvalb[!is.na(res_Pvalb$padj) & res_Pvalb$padj < 0.05, ]
# Then order and extract gene names
top_genes_Pvalb <- rownames(res_Pvalb_filtered[order(res_Pvalb_filtered$padj), ])
# Perform Variance-Stabilizing Transformation (VST)
vst_data_Pvalb <- vst(dds_Pvalb, blind = TRUE)
# Ensure genes exist in the dataset
common_genes_Pvalb <- intersect(top_genes_Pvalb, rownames(assay(vst_data_Pvalb)))
# Extract the VST-transformed expression values for these genes
heatmap_matrix_Pvalb <- assay(vst_data_Pvalb)[common_genes_Pvalb, ]
# Extract metadata for annotations
annotation_Pvalb <- sample_info_Pvalb %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Pvalb$condition <- factor(annotation_Pvalb$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Pvalb <- rownames(annotation_Pvalb)[order(annotation_Pvalb$condition)]
heatmap_matrix_Pvalb <- heatmap_matrix_Pvalb[, ordered_samples_Pvalb]
# Change the colors manually
annotation_colors <- list(
condition = c(
"Control" = "#7ef29d",
"Risperidone" = "#0f68a9"))
annotation_colors_Pvalb <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
#Generate Heatmap
pheatmap(heatmap_matrix_Pvalb,
scale = "row", # Normalize each gene across samples
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Pvalb,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
```
PV-IN DEG - volcano
```{r}
# Convert DESeq2 results to a data frame
res_df_Pvalb <- as.data.frame(res_Pvalb)
res_df_Pvalb$gene <- rownames(res_df_Pvalb)
# Remove any NA values
res_df_Pvalb <- na.omit(res_df_Pvalb)
# Add a new column to categorize genes for coloring
res_df_Pvalb$Significance <- "Not Significant"
res_df_Pvalb$Significance[res_df_Pvalb$padj < padj_threshold & abs(res_df_Pvalb$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_Pvalb <- c("Kcna2", "Kcnc1", "Kcnc3", "Kcnd3", "Kcnq2", "Kcnq3", "Grin2a", "Grin2b", "Dlg1", "Shank1", "Shank2", "Gabrg3", "Gabra1", "Gabrb1")
# Generate Volcano plot
volcano_plot_Pvalb <- ggplot(res_df_Pvalb, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_df_Pvalb %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_df_Pvalb %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_df_Pvalb %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Pvalb, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 15,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "PV-IN",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
# Print the plot
print(volcano_plot_Pvalb)
```
DE analysis (Pseudobulk) - Sst-IN
```{r}
# Subset cell type
Mcortex_QC_Sst <- subset(Mcortex_QC, subset = celltype1 == "GABAergic Sst")
# Filter out all genes with count sum <10
Mcortex_QC_Sst <- subset(Mcortex_QC_Sst, features = rownames(Mcortex_QC_Sst)[Matrix::rowSums(GetAssayData(Mcortex_QC_Sst, assay = "RNA", layer = "counts")) >= 10])
# Extract sample_id from barcode names
Mcortex_QC_Sst$sample_id <- sapply(strsplit(Cells(Mcortex_QC_Sst), "_"), `[`, 1)
# Aggregate counts
expr_mat_Sst <- AggregateExpression(
Mcortex_QC_Sst,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Sst <- data.frame(
sample_id = colnames(expr_mat_Sst),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Sst)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Sst) <- sample_info_Sst$sample_id
# Build DESeq2 object
dds_Sst <- DESeqDataSetFromMatrix(
countData = expr_mat_Sst,
colData = sample_info_Sst,
design = ~ condition
)
# Run DE analysis
dds_Sst <- DESeq(dds_Sst)
# Get results
res_Sst <- results(dds_Sst, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Sst) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Sst-IN.csv"),
row.names = TRUE
)
```
Sst-IN DEG - heatmap
```{r}
# Count the number of DEGs with padj < 0.05
sum(res_Sst$padj < 0.05, na.rm = TRUE) #1305
# Select all DEGs with padj < 0.05 and remove NA
res_Sst_filtered <- res_Sst[!is.na(res_Sst$padj) & res_Sst$padj < 0.05, ]
# Then order and extract gene names
top_genes_Sst <- rownames(res_Sst_filtered[order(res_Sst_filtered$padj), ])
# Perform Variance-Stabilizing Transformation (VST)
vst_data_Sst <- vst(dds_Sst, blind = TRUE)
# Ensure genes exist in the dataset
common_genes_Sst <- intersect(top_genes_Sst, rownames(assay(vst_data_Sst)))
# Extract the VST-transformed expression values for these genes
heatmap_matrix_Sst <- assay(vst_data_Sst)[common_genes_Sst, ]
# Extract metadata for annotations
annotation_Sst <- sample_info_Sst %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Sst$condition <- factor(annotation_Sst$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Sst <- rownames(annotation_Sst)[order(annotation_Sst$condition)]
heatmap_matrix_Sst <- heatmap_matrix_Sst[, ordered_samples_Sst]
annotation_colors_Sst <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
#Generate Heatmap
pheatmap(heatmap_matrix_Sst,
scale = "row", # Normalize each gene across samples
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_PN,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "",
)
```
Sst-IN DEG - volcano
```{r}
# Convert DESeq2 results to a data frame
res_df_Sst <- as.data.frame(res_Sst)
res_df_Sst$gene <- rownames(res_df_Sst)
# Remove any NA values
res_df_Sst <- na.omit(res_df_Sst)
# Add a new column to categorize genes for coloring
res_df_Sst$Significance <- "Not Significant"
res_df_Sst$Significance[res_df_Sst$padj < padj_threshold & abs(res_df_Sst$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_Sst <- c("Kcnb1", "Kcnc1", "Kcnd3", "Kcnq3", "Kcnq2", "Kcnh7", "Grin2a", "Grin2b", "Dlg1", "Shank1", "Shank2", "Gabrg3")
# Generate Zoomed-in Volcano plot
volcano_plot_Sst_zoom <- ggplot(res_df_Sst, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_df_Sst %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 2) +
geom_point(data = res_df_Sst %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 2) +
geom_point(data = res_df_Sst %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 2) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Sst, gene, "")),
size = 5,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
coord_cartesian(xlim = c(-4, 4), ylim = c(0, 100)) + # 👈 Zoomed view
theme_minimal() +
labs(title = "",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(legend.position = "right",
axis.text = element_text(size = 14),
axis.title = element_text(size = 16))
# Print the plot
print(volcano_plot_Sst_zoom)
```
Graph specific DEG expression change by treatment
```{r}
## Glut
# To plot log2FC in DEG results (not the raw count)
genes_to_plot <- c("Drd2", "Drd3", "Ddc", "Domt", "Htr2a", "Htr1f", "Htr4", "Htr2c", "Htf5a", "Htr7", "Htr6", "Htr1d", "Htr2b", "Maoa")
# Subset the gene of interest
res_df_Glut_plot <- res_df_Glut[res_df_Glut$gene %in% genes_to_plot, c("gene", "log2FoldChange", "padj")]
# Create 2-line p-value label column
res_df_Glut_plot$label <- sprintf("p =\n%.2g", res_df_Glut_plot$padj)
# Create the bar plot for log2FC
ggplot(res_df_Glut_plot, aes(x = gene, y = log2FoldChange, fill = log2FoldChange > 0)) +
geom_col() +
scale_fill_manual(values = c("steelblue", "firebrick")) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_text_repel(
aes(label = label),
size = 4,
nudge_y = 0.5, # pushes labels slightly above/below bars
direction = "y", # keeps labels aligned vertically
segment.color = "grey50",
box.padding = 0.3,
point.padding = 0.2
) +
labs(
title = "Glutamatergic neurons saline vs. risperidone",
x = NULL,
y = "Log2 Fold Change"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14))
## PV-IN
# To plot log2FC in DEG results (not the raw count)
genes_to_plot <- c("Drd2", "Htr2a")
# Subset the gene of interest
res_df_Pvalb_plot <- res_df_Pvalb[res_df_Pvalb$gene %in% genes_to_plot, c("gene", "log2FoldChange", "padj")]
# Create 2-line p-value label column
res_df_Pvalb_plot$label <- sprintf("p =\n%.2g", res_df_Pvalb_plot$padj)
# Create the bar plot for log2FC
ggplot(res_df_Pvalb_plot, aes(x = gene, y = log2FoldChange, fill = log2FoldChange > 0)) +
geom_col() +
scale_fill_manual(values = c("steelblue", "firebrick")) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_text_repel(
aes(label = label),
size = 4,
nudge_y = 0.5, # pushes labels slightly above/below bars
direction = "y", # keeps labels aligned vertically
segment.color = "grey50",
box.padding = 0.3,
point.padding = 0.2
) +
labs(
title = "PV-IN saline vs. risperidone",
x = NULL,
y = "Log2 Fold Change"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14))
## Sst-IN
# To plot log2FC in DEG results (not the raw count)
genes_to_plot <- c("Drd2")
# Subset the gene of interest
res_df_Sst_plot <- res_df_Sst[res_df_Sst$gene %in% genes_to_plot, c("gene", "log2FoldChange", "padj")]
# Create 2-line p-value label column
res_df_Sst_plot$label <- sprintf("p =\n%.2g", res_df_Sst_plot$padj)
# Create the bar plot for log2FC
ggplot(res_df_Sst_plot, aes(x = gene, y = log2FoldChange, fill = log2FoldChange > 0)) +
geom_col() +
scale_fill_manual(values = c("steelblue", "firebrick")) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_text_repel(
aes(label = label),
size = 4,
nudge_y = 0.5, # pushes labels slightly above/below bars
direction = "y", # keeps labels aligned vertically
segment.color = "grey50",
box.padding = 0.3,
point.padding = 0.2
) +
labs(
title = "Sst-IN saline vs. risperidone",
x = NULL,
y = "Log2 Fold Change"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14))
```
Venn Diagram for DEGs
```{r}
library(VennDiagram)
library(grid)
# Set adjusted FDR threshold
padj_cutoff <- 0.05
# Get significant gene names from the DEG lists from 3 cell types
sig_genes_Glut <- rownames(res_Glut[which(res_Glut$padj < padj_cutoff), ])
sig_genes_Pvalb <- rownames(res_Pvalb[which(res_Pvalb$padj < padj_cutoff), ])
sig_genes_Sst <- rownames(res_Sst[which(res_Sst$padj < padj_cutoff), ])
# Create and draw the Venn diagram
venn_plot <- venn.diagram(
x = list(
Glut = sig_genes_Glut,
PV = sig_genes_Pvalb,
Sst = sig_genes_Sst
),
filename = NULL, # Don't save to file
fill = c("skyblue", "pink", "yellow"),
alpha = 0.5,
cex = 1.5,
cat.cex = 0,
cat.pos = 0,
cat.dist = 0.05,
scaled = FALSE,
main = "Overlap of DEGs in 3 cell types"
)
grid.draw(venn_plot)
# Prepare a named list of vectors (each becomes a sheet)
DEG_threecell <- list(
DEGs_Glut = data.frame(Gene = sig_genes_Glut),
DEGs_PV = data.frame(Gene = sig_genes_Pvalb),
DEGs_Sst = data.frame(Gene = sig_genes_Sst),
Shared_Glut_PV = data.frame(Gene = intersect(sig_genes_Glut, sig_genes_Pvalb)),
Shared_Glut_Sst = data.frame(Gene = intersect(sig_genes_Glut, sig_genes_Sst)),
Shared_Sst_PV = data.frame(Gene = intersect(sig_genes_Sst, sig_genes_Pvalb)),
Shared_all = data.frame(
Gene = Reduce(intersect, list(
sig_genes_Glut,
sig_genes_Sst,
sig_genes_Pvalb
))
)
)
# Save to Excel file
write_xlsx(DEG_threecell, path = "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/DE_modified/Proj4_threecell_overlap.xlsx")
## EnrichGO for DEGs that are shared by all 3 cell types
# Run EnrichGO - BP
BP_GO_shared <- enrichGO(
gene = DEG_threecell$Shared_all$Gene,
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run EnrichGO - MF
MF_GO_shared <- enrichGO(
gene = DEG_threecell$Shared_all$Gene,
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL",
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run EnrichGO - CC
CC_GO_shared <- enrichGO(
gene = DEG_threecell$Shared_all$Gene,
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Save as Excel file
write.xlsx(
list(
BP = as.data.frame(BP_GO_shared),
MF = as.data.frame(MF_GO_shared),
CC = as.data.frame(CC_GO_shared)
),
file = "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/EnrichGO/Proj4/Proj4_threecelltype_overlap.xlsx"
)
# Plot dotplot
dotplot(BP_GO_shared, showCategory = 10, title = "GO Enrichment (Biological Process)") +
theme(
text = element_text(size = 14), # Adjust global text size
axis.text.x = element_text(size = 12), # Adjust x-axis font size
axis.text.y = element_text(size = 14), # Adjust y-axis font size
plot.title = element_text(size = 12) # Adjust title font size
)
dotplot(MF_GO_shared, showCategory = 10, title = "GO Enrichment (Molecular Function)") +
theme(
text = element_text(size = 14), # Adjust global text size
axis.text.x = element_text(size = 12), # Adjust x-axis font size
axis.text.y = element_text(size = 14), # Adjust y-axis font size
plot.title = element_text(size = 12) # Adjust title font size
)
dotplot(CC_GO_shared, showCategory = 10, title = "GO Enrichment (Cellular Component)") +
theme(
text = element_text(size = 14), # Adjust global text size
axis.text.x = element_text(size = 12), # Adjust x-axis font size
axis.text.y = element_text(size = 14), # Adjust y-axis font size
plot.title = element_text(size = 12) # Adjust title font size
)
```
DE analysis (Pseudobulk) - PV-IN - F vs. M
```{r}
# Subset male and female separately
Mcortex_QC_Pvalb_F <- subset(Mcortex_QC_Pvalb, Predicted_Sex %in% "Female")
Mcortex_QC_Pvalb_M <- subset(Mcortex_QC_Pvalb, Predicted_Sex %in% "Male")
## Female
# Filter out all genes with count sum less than 10
Mcortex_QC_Pvalb_F <- subset(Mcortex_QC_Pvalb_F, features = rownames(Mcortex_QC_Pvalb_F)[Matrix::rowSums(GetAssayData(Mcortex_QC_Pvalb_F, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Pvalb_F <- AggregateExpression(
Mcortex_QC_Pvalb_F,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Pvalb_F <- data.frame(
sample_id = colnames(expr_mat_Pvalb_F),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Pvalb_F)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Pvalb_F) <- sample_info_Pvalb_F$sample_id
# Build DESeq2 object
dds_Pvalb_F <- DESeqDataSetFromMatrix(
countData = expr_mat_Pvalb_F,
colData = sample_info_Pvalb_F,
design = ~ condition
)
# Run DE analysis
dds_Pvalb_F <- DESeq(dds_Pvalb_F)
# Get results
res_Pvalb_F <- results(dds_Pvalb_F, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Pvalb_F) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_PV-IN_F.csv"),
row.names = TRUE
)
## Male
# Filter out all genes with count sum less than 10
Mcortex_QC_Pvalb_M <- subset(Mcortex_QC_Pvalb_M, features = rownames(Mcortex_QC_Pvalb_M)[Matrix::rowSums(GetAssayData(Mcortex_QC_Pvalb_M, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Pvalb_M <- AggregateExpression(
Mcortex_QC_Pvalb_M,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Pvalb_M <- data.frame(
sample_id = colnames(expr_mat_Pvalb_M),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Pvalb_M)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Pvalb_M) <- sample_info_Pvalb_M$sample_id
# Build DESeq2 object
dds_Pvalb_M <- DESeqDataSetFromMatrix(
countData = expr_mat_Pvalb_M,
colData = sample_info_Pvalb_M,
design = ~ condition
)
# Run DE analysis
dds_Pvalb_M <- DESeq(dds_Pvalb_M)
# Get results
res_Pvalb_M <- results(dds_Pvalb_M, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Pvalb_M) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_PV-IN_M.csv"),
row.names = TRUE
)
```
PV-IN - F vs. M DEG Heatmap
```{r}
## Female
# Count the number of DEGs with padj < 0.05
sum(res_Pvalb_F$padj < 0.05, na.rm = TRUE) #701
# Select all DEGs with adjusted p-value < 0.05
sig_Pvalb_F <- res_Pvalb_F[!is.na(res_Pvalb_F$padj) & res_Pvalb_F$padj < 0.05, ]
# Then order and extract gene names
top_genes_Pvalb_F <- rownames(sig_Pvalb_F[order(sig_Pvalb_F$padj), ])
# Perform VST
vst_data_Pvalb_F <- vst(dds_Pvalb_F, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Pvalb_F <- intersect(top_genes_Pvalb_F, rownames(assay(vst_data_Pvalb_F)))
# Extract VST expression matrix
heatmap_matrix_Pvalb_F <- assay(vst_data_Pvalb_F)[common_genes_Pvalb_F, ]
# Extract metadata for annotations
annotation_Pvalb_F <- sample_info_Pvalb_F %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Pvalb_F$condition <- factor(annotation_Pvalb_F$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Pvalb_F <- rownames(annotation_Pvalb_F)[order(annotation_Pvalb_F$condition)]
heatmap_matrix_Pvalb_F <- heatmap_matrix_Pvalb_F[, ordered_samples_Pvalb_F]
annotation_colors_Pvalb_F <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Pvalb_F,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Pvalb_F,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
## Male
# Count the number of DEGs with padj < 0.05
sum(res_Pvalb_M$padj < 0.05, na.rm = TRUE) #528
# Select all DEGs with adjusted p-value < 0.05
sig_Pvalb_M <- res_Pvalb_M[!is.na(res_Pvalb_M$padj) & res_Pvalb_M$padj < 0.05, ]
# Then order and extract gene names
top_genes_Pvalb_M <- rownames(sig_Pvalb_M[order(sig_Pvalb_M$padj), ])
# Perform VST
vst_data_Pvalb_M <- vst(dds_Pvalb_M, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Pvalb_M <- intersect(top_genes_Pvalb_M, rownames(assay(vst_data_Pvalb_M)))
# Extract VST expression matrix
heatmap_matrix_Pvalb_M <- assay(vst_data_Pvalb_M)[common_genes_Pvalb_M, ]
# Extract metadata for annotations
annotation_Pvalb_M <- sample_info_Pvalb_M %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Pvalb_M$condition <- factor(annotation_Pvalb_M$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Pvalb_M <- rownames(annotation_Pvalb_M)[order(annotation_Pvalb_M$condition)]
heatmap_matrix_Pvalb_M <- heatmap_matrix_Pvalb_M[, ordered_samples_Pvalb_M]
annotation_colors_Pvalb_M <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Pvalb_M,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Pvalb_M,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
```
PV-IN - F vs. M DEG volcano
```{r}
# Convert DESeq2 results to a data frame
res_Pvalb_F_df <- as.data.frame(res_Pvalb_F)
res_Pvalb_F_df$gene <- rownames(res_Pvalb_F_df)
res_Pvalb_M_df <- as.data.frame(res_Pvalb_M)
res_Pvalb_M_df$gene <- rownames(res_Pvalb_M_df)
# Remove any NA values
res_Pvalb_F_df <- na.omit(res_Pvalb_F_df)
res_Pvalb_M_df <- na.omit(res_Pvalb_M_df)
# Add a new column to categorize genes for coloring
res_Pvalb_F_df$Significance <- "Not Significant"
res_Pvalb_F_df$Significance[res_Pvalb_F_df$padj < padj_threshold & abs(res_Pvalb_F_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
res_Pvalb_M_df$Significance <- "Not Significant"
res_Pvalb_M_df$Significance[res_Pvalb_M_df$padj < padj_threshold & abs(res_Pvalb_M_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_PV_F <- c("Kcnd3", "Kcnt2", "Shank1", "Shank2", "Dlg1", "Gria1", "Grid1", "Grik3", "Gabrg3")
genes_of_interest_PV_M <- c("Kcnd3", "Kcnh1", "Nrg1", "Shank1", "Shank2", "Grid1", "Grik3", "Grm8")
# Generate Volcano plot
volcano_plot_PV_F <- ggplot(res_Pvalb_F_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Pvalb_F_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Pvalb_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Pvalb_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_PV_F, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "PV-IN F",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
volcano_plot_PV_M <- ggplot(res_Pvalb_M_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Pvalb_M_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Pvalb_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Pvalb_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_PV_M, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "PV-IN M",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
# Print the plot
print(volcano_plot_PV_F)
print(volcano_plot_PV_M)
```
PV-IN DEG - F vs. M list comparison - Venn Diagram
```{r}
# Set adjusted p-value threshold
padj_cutoff <- 0.05
# Get significant gene names
sig_genes_Pvalb_F <- rownames(res_Pvalb_F[which(res_Pvalb_F$padj < padj_cutoff), ])
sig_genes_Pvalb_M <- rownames(res_Pvalb_M[which(res_Pvalb_M$padj < padj_cutoff), ])
# Prepare a named list of vectors (each becomes a sheet)
deg_lists <- list(
DEGs_Male = data.frame(Gene = sig_genes_Pvalb_M),
DEGs_Female = data.frame(Gene = sig_genes_Pvalb_F),
Shared_DEGs = data.frame(Gene = intersect(sig_genes_Pvalb_M, sig_genes_Pvalb_F)),
Male_only = data.frame(Gene = setdiff(sig_genes_Pvalb_M, sig_genes_Pvalb_F)),
Female_only = data.frame(Gene = setdiff(sig_genes_Pvalb_F, sig_genes_Pvalb_M))
)
# Save to Excel file
write_xlsx(deg_lists, path = "~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_PV-IN_FvM.xlsx")
```
DE analysis (Pseudobulk) - Glut - F vs. M
```{r}
# Subset male and female separately
Mcortex_QC_Glut_F <- subset(Mcortex_QC_Glut, Predicted_Sex %in% "Female")
Mcortex_QC_Glut_M <- subset(Mcortex_QC_Glut, Predicted_Sex %in% "Male")
## Female
# Filter out all genes with count sum less than 10
Mcortex_QC_Glut_F <- subset(Mcortex_QC_Glut_F, features = rownames(Mcortex_QC_Glut_F)[Matrix::rowSums(GetAssayData(Mcortex_QC_Glut_F, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Glut_F <- AggregateExpression(
Mcortex_QC_Glut_F,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Glut_F <- data.frame(
sample_id = colnames(expr_mat_Glut_F),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Glut_F)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Glut_F) <- sample_info_Glut_F$sample_id
# Build DESeq2 object
dds_Glut_F <- DESeqDataSetFromMatrix(
countData = expr_mat_Glut_F,
colData = sample_info_Glut_F,
design = ~ condition
)
# Run DE analysis
dds_Glut_F <- DESeq(dds_Glut_F)
# Get results
res_Glut_F <- results(dds_Glut_F, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Glut_F) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Glut_F.csv"),
row.names = TRUE
)
## Male
# Filter out all genes with count sum less than 10
Mcortex_QC_Glut_M <- subset(Mcortex_QC_Glut_M, features = rownames(Mcortex_QC_Glut_M)[Matrix::rowSums(GetAssayData(Mcortex_QC_Glut_M, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Glut_M <- AggregateExpression(
Mcortex_QC_Glut_M,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Glut_M <- data.frame(
sample_id = colnames(expr_mat_Glut_M),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Glut_M)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Glut_M) <- sample_info_Glut_M$sample_id
# Build DESeq2 object
dds_Glut_M <- DESeqDataSetFromMatrix(
countData = expr_mat_Glut_M,
colData = sample_info_Glut_M,
design = ~ condition
)
# Run DE analysis
dds_Glut_M <- DESeq(dds_Glut_M)
# Get results
res_Glut_M <- results(dds_Glut_M, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Glut_M) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Glut_M.csv"),
row.names = TRUE
)
```
Glut - F vs. M DEG Heatmap
```{r}
## Female
# Count the number of DEGs with padj < 0.05
sum(res_Glut_F$padj < 0.05, na.rm = TRUE) #5534
# Select all DEGs with adjusted p-value < 0.05
sig_Glut_F <- res_Glut_F[!is.na(res_Glut_F$padj) & res_Glut_F$padj < 0.05, ]
# Then order and extract gene names
top_genes_Glut_F <- rownames(sig_Glut_F[order(sig_Glut_F$padj), ])
# Perform VST
vst_data_Glut_F <- vst(dds_Glut_F, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Glut_F <- intersect(top_genes_Glut_F, rownames(assay(vst_data_Glut_F)))
# Extract VST expression matrix
heatmap_matrix_Glut_F <- assay(vst_data_Glut_F)[common_genes_Glut_F, ]
# Extract metadata for annotations
annotation_Glut_F <- sample_info_Glut_F %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Glut_F$condition <- factor(annotation_Glut_F$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Glut_F <- rownames(annotation_Glut_F)[order(annotation_Glut_F$condition)]
heatmap_matrix_Glut_F <- heatmap_matrix_Glut_F[, ordered_samples_Glut_F]
annotation_colors_Glut_F <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Glut_F,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Glut_F,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
## Male
# Count the number of DEGs with padj < 0.05
sum(res_Glut_M$padj < 0.05, na.rm = TRUE) #5643
# Select all DEGs with adjusted p-value < 0.05
sig_Glut_M <- res_Glut_M[!is.na(res_Glut_M$padj) & res_Glut_M$padj < 0.05, ]
# Then order and extract gene names
top_genes_Glut_M <- rownames(sig_Glut_M[order(sig_Glut_M$padj), ])
# Perform VST
vst_data_Glut_M <- vst(dds_Glut_M, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Glut_M <- intersect(top_genes_Glut_M, rownames(assay(vst_data_Glut_M)))
# Extract VST expression matrix
heatmap_matrix_Glut_M <- assay(vst_data_Glut_M)[common_genes_Glut_M, ]
# Extract metadata for annotations
annotation_Glut_M <- sample_info_Glut_M %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Glut_M$condition <- factor(annotation_Glut_M$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Glut_M <- rownames(annotation_Glut_M)[order(annotation_Glut_M$condition)]
heatmap_matrix_Glut_M <- heatmap_matrix_Glut_M[, ordered_samples_Glut_M]
annotation_colors_Glut_M <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Glut_M,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Glut_M,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
```
Glut - F vs. M DEG volcano
```{r}
# Convert DESeq2 results to a data frame
res_Glut_F_df <- as.data.frame(res_Glut_F)
res_Glut_F_df$gene <- rownames(res_Glut_F_df)
res_Glut_M_df <- as.data.frame(res_Glut_M)
res_Glut_M_df$gene <- rownames(res_Glut_M_df)
# Remove any NA values
res_Glut_F_df <- na.omit(res_Glut_F_df)
res_Glut_M_df <- na.omit(res_Glut_M_df)
# Add a new column to categorize genes for coloring
res_Glut_F_df$Significance <- "Not Significant"
res_Glut_F_df$Significance[res_Glut_F_df$padj < padj_threshold & abs(res_Glut_F_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
res_Glut_M_df$Significance <- "Not Significant"
res_Glut_M_df$Significance[res_Glut_M_df$padj < padj_threshold & abs(res_Glut_M_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_Glut_F <- c("Kcnd3", "Kcnq3", "Kcnt2", "Kcnb1", "Shank1", "Shank2", "Grin2a", "Grin2b", "Dlg1", "Gria1", "Grid1", "Grik3", "Nrg1", "Erbb4", "Gabrb1", "Gabra5")
genes_of_interest_Glut_M <- c("Kcnd3", "Kcnq3", "Kcnt2", "Nrg1", "Erbb4", "Shank1", "Shank2", "Gria1", "Grid1", "Grik3", "Grin2a")
# Generate Volcano plot
volcano_plot_Glut_F <- ggplot(res_Glut_F_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Glut_F_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Glut_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Glut_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Glut_F, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 25,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "Glut F",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
volcano_plot_Glut_M <- ggplot(res_Glut_M_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Glut_M_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Glut_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Glut_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Glut_M, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "Glut M",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
# Print the plot
print(volcano_plot_Glut_F)
print(volcano_plot_Glut_M)
```
Glut DEG - F vs. M list comparison
```{r}
# Get significant gene names
sig_genes_Glut_F <- rownames(res_Glut_F[which(res_Glut_F$padj < padj_cutoff), ])
sig_genes_Glut_M <- rownames(res_Glut_M[which(res_Glut_M$padj < padj_cutoff), ])
# Prepare a named list of vectors (each becomes a sheet)
deg_lists <- list(
DEGs_Male = data.frame(Gene = sig_genes_Glut_M),
DEGs_Female = data.frame(Gene = sig_genes_Glut_F),
Shared_DEGs = data.frame(Gene = intersect(sig_genes_Glut_M, sig_genes_Glut_F)),
Male_only = data.frame(Gene = setdiff(sig_genes_Glut_M, sig_genes_Glut_F)),
Female_only = data.frame(Gene = setdiff(sig_genes_Glut_F, sig_genes_Glut_M))
)
# Save to Excel file
write_xlsx(deg_lists, path = "~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Glut_FvM.xlsx")
```
DE analysis (Pseudobulk) - Sst-IN - F vs. M
```{r}
# Subset male and female separately
Mcortex_QC_Sst_F <- subset(Mcortex_QC_Sst, Predicted_Sex %in% "Female")
Mcortex_QC_Sst_M <- subset(Mcortex_QC_Sst, Predicted_Sex %in% "Male")
## Female
# Filter out all genes with count sum less than 10
Mcortex_QC_Sst_F <- subset(Mcortex_QC_Sst_F, features = rownames(Mcortex_QC_Sst_F)[Matrix::rowSums(GetAssayData(Mcortex_QC_Sst_F, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Sst_F <- AggregateExpression(
Mcortex_QC_Sst_F,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Sst_F <- data.frame(
sample_id = colnames(expr_mat_Sst_F),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Sst_F)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Sst_F) <- sample_info_Sst_F$sample_id
# Build DESeq2 object
dds_Sst_F <- DESeqDataSetFromMatrix(
countData = expr_mat_Sst_F,
colData = sample_info_Sst_F,
design = ~ condition
)
# Run DE analysis
dds_Sst_F <- DESeq(dds_Sst_F)
# Get results
res_Sst_F <- results(dds_Sst_F, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Sst_F) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Sst_F.csv"),
row.names = TRUE
)
## Male
# Filter out all genes with count sum less than 10
Mcortex_QC_Sst_M <- subset(Mcortex_QC_Sst_M, features = rownames(Mcortex_QC_Sst_M)[Matrix::rowSums(GetAssayData(Mcortex_QC_Sst_M, assay = "RNA", layer = "counts")) >= 10])
# Aggregate counts
expr_mat_Sst_M <- AggregateExpression(
Mcortex_QC_Sst_M,
group.by = c("Tx", "sample_id"),
assays = "RNA"
)$RNA
# Create colData (metadata about your samples)
sample_info_Sst_M <- data.frame(
sample_id = colnames(expr_mat_Sst_M),
condition = ifelse(grepl("Risperidone", colnames(expr_mat_Sst_M)), "Risperidone", "Control"),
stringsAsFactors = FALSE
)
rownames(sample_info_Sst_M) <- sample_info_Sst_M$sample_id
# Build DESeq2 object
dds_Sst_M <- DESeqDataSetFromMatrix(
countData = expr_mat_Sst_M,
colData = sample_info_Sst_M,
design = ~ condition
)
# Run DE analysis
dds_Sst_M <- DESeq(dds_Sst_M)
# Get results
res_Sst_M <- results(dds_Sst_M, contrast = c("condition", "Risperidone", "Control"))
# Save results
write.csv(
as.data.frame(res_Sst_M) %>%
subset(!is.na(padj)) %>%
arrange(padj, desc(abs(log2FoldChange))),
file = paste0("~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Sst_M.csv"),
row.names = TRUE
)
```
Sst-IN - F vs. M DEG Heatmap
```{r}
## Female
# Count the number of DEGs with padj < 0.05
sum(res_Sst_F$padj < 0.05, na.rm = TRUE) #339
# Select all DEGs with adjusted p-value < 0.05
sig_Sst_F <- res_Sst_F[!is.na(res_Sst_F$padj) & res_Sst_F$padj < 0.05, ]
# Then order and extract gene names
top_genes_Sst_F <- rownames(sig_Sst_F[order(sig_Sst_F$padj), ])
# Perform VST
vst_data_Sst_F <- vst(dds_Sst_F, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Sst_F <- intersect(top_genes_Sst_F, rownames(assay(vst_data_Sst_F)))
# Extract VST expression matrix
heatmap_matrix_Sst_F <- assay(vst_data_Sst_F)[common_genes_Sst_F, ]
# Extract metadata for annotations
annotation_Sst_F <- sample_info_Sst_F %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Sst_F$condition <- factor(annotation_Sst_F$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Sst_F <- rownames(annotation_Sst_F)[order(annotation_Sst_F$condition)]
heatmap_matrix_Sst_F <- heatmap_matrix_Sst_F[, ordered_samples_Sst_F]
annotation_colors_Sst_F <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Sst_F,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Sst_F,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
## Male
# Count the number of DEGs with padj < 0.05
sum(res_Sst_M$padj < 0.05, na.rm = TRUE) #269
# Select all DEGs with adjusted p-value < 0.05
sig_Sst_M <- res_Sst_M[!is.na(res_Sst_M$padj) & res_Sst_M$padj < 0.05, ]
# Then order and extract gene names
top_genes_Sst_M <- rownames(sig_Sst_M[order(sig_Sst_M$padj), ])
# Perform VST
vst_data_Sst_M <- vst(dds_Sst_M, blind = TRUE)
# Keep only genes that are present in the VST-transformed data
common_genes_Sst_M <- intersect(top_genes_Sst_M, rownames(assay(vst_data_Sst_M)))
# Extract VST expression matrix
heatmap_matrix_Sst_M <- assay(vst_data_Sst_M)[common_genes_Sst_M, ]
# Extract metadata for annotations
annotation_Sst_M <- sample_info_Sst_M %>% select(condition)
# Make sure your 'condition' column is a factor in desired order
annotation_Sst_M$condition <- factor(annotation_Sst_M$condition, levels = c("Control", "Risperidone"))
# Now reorder the columns of your heatmap matrix
ordered_samples_Sst_M <- rownames(annotation_Sst_M)[order(annotation_Sst_M$condition)]
heatmap_matrix_Sst_M <- heatmap_matrix_Sst_M[, ordered_samples_Sst_M]
annotation_colors_Sst_M <- list(
condition = c(
"Control" = brewer.pal(8, "Set2")[1],
"Risperidone" = brewer.pal(8, "Set2")[2]
)
)
# Generate heatmap
pheatmap(heatmap_matrix_Sst_M,
scale = "row",
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
annotation_col = annotation_Sst_M,
annotation_colors = annotation_colors,
border_color = NA,
fontsize = 8,
main = "")
```
Sst-IN - F vs. M DEG volcano
```{r}
# Convert DESeq2 results to a data frame
res_Sst_F_df <- as.data.frame(res_Sst_F)
res_Sst_F_df$gene <- rownames(res_Sst_F_df)
res_Sst_M_df <- as.data.frame(res_Sst_M)
res_Sst_M_df$gene <- rownames(res_Sst_M_df)
# Remove any NA values
res_Sst_F_df <- na.omit(res_Sst_F_df)
res_Sst_M_df <- na.omit(res_Sst_M_df)
# Add a new column to categorize genes for coloring
res_Sst_F_df$Significance <- "Not Significant"
res_Sst_F_df$Significance[res_Sst_F_df$padj < padj_threshold & abs(res_Sst_F_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
res_Sst_M_df$Significance <- "Not Significant"
res_Sst_M_df$Significance[res_Sst_M_df$padj < padj_threshold & abs(res_Sst_M_df$log2FoldChange) > log2FoldChange_threshold] <- "Significant"
# Generate volcano plot with specific genes of interest labeled
# Define genes of interest
genes_of_interest_Sst_F <- c("Kcnd3", "Kcnq3", "Kcnt2", "Shank1", "Shank2", "Gria1", "Grid1", "Grik3", "Gabra5")
genes_of_interest_Sst_M <- c("Kcnd3", "Kcnt2", "Kcnh1", "Shank1", "Shank2", "Grik3")
# Generate Volcano plot
volcano_plot_Sst_F <- ggplot(res_Sst_F_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Sst_F_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Sst_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Sst_F_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Sst_F, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "Sst-IN F",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
volcano_plot_Sst_M <- ggplot(res_Sst_M_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(data = res_Sst_M_df %>% dplyr::filter(Significance == "Not Significant"),
aes(x = log2FoldChange, y = -log10(padj)),
color = "gray", alpha = 0.5, size = 1) +
geom_point(data = res_Sst_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange > 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "red", alpha = 0.5, size = 1) +
geom_point(data = res_Sst_M_df %>% dplyr::filter(Significance == "Significant" & log2FoldChange < 0),
aes(x = log2FoldChange, y = -log10(padj)),
color = "blue", alpha = 0.5, size = 1) +
geom_label_repel(
aes(label = ifelse(gene %in% genes_of_interest_Sst_M, gene, "")),
size = 7,
box.padding = 0.75,
label.padding = 0.35,
fill = alpha("white", 0),
color = "black",
segment.color = "black",
segment.size = 0.8,
segment.alpha = 0.8,
min.segment.length = 0,
force = 5,
max.overlaps = Inf
) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "black") +
theme_minimal() +
labs(title = "Sst-IN M",
x = "Log2 Fold Change",
y = "-Log10 Adjusted P-Value") +
theme(
legend.position = "right",
axis.title.x = element_text(size = 16), # X-axis label size
axis.title.y = element_text(size = 16), # Y-axis label size
axis.text.x = element_text(size = 14), # X-axis tick labels
axis.text.y = element_text(size = 14), # Y-axis tick labels
plot.title = element_text(size = 20, face = "bold", hjust = 0.5, margin = margin(b = 10)),
plot.margin = margin(t = 10, r = 10, b = 10, l = 10) # adds space above the whole plot
)
# Print the plot
print(volcano_plot_Sst_F)
print(volcano_plot_Sst_M)
```
Sst-IN DEG - F vs. M list comparison
```{r}
# Get significant gene names
sig_genes_Sst_F <- rownames(res_Sst_F[which(res_Sst_F$padj < padj_cutoff), ])
sig_genes_Sst_M <- rownames(res_Sst_M[which(res_Sst_M$padj < padj_cutoff), ])
# Prepare a named list of vectors (each becomes a sheet)
deg_lists <- list(
DEGs_Male = data.frame(Gene = sig_genes_Sst_M),
DEGs_Female = data.frame(Gene = sig_genes_Sst_F),
Shared_DEGs = data.frame(Gene = intersect(sig_genes_Sst_M, sig_genes_Sst_F)),
Male_only = data.frame(Gene = setdiff(sig_genes_Sst_M, sig_genes_Sst_F)),
Female_only = data.frame(Gene = setdiff(sig_genes_Sst_F, sig_genes_Sst_M))
)
# Save to Excel file
write_xlsx(deg_lists, path = "~/Risperidone_Proj/R_analysis/DE_modified/Proj4_C_v_R_Sst_FvM.xlsx")
```
EnrichGO - PV-IN
```{r}
# Filter DEG list
res_Pvalb_enrichGO <- rownames(res_Pvalb[!is.na(res_Pvalb$padj) &
res_Pvalb$padj < 0.05 &
abs(res_Pvalb$log2FoldChange) > 0.5, ])
# Run EnrichGO - BP
BP_GO_Pvalb <- enrichGO(
gene = res_Pvalb_enrichGO,
OrgDb = org.Mm.eg.db, # Use org.Mm.eg.db for mouse
keyType = "SYMBOL",
ont = "BP", # "BP" (Biological Process), "MF" (Molecular Function), "CC" (Cellular Component)
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# View the first few enriched GO terms
head(BP_GO_Pvalb)
# Run EnrichGO - MF
MF_GO_Pvalb <- enrichGO(
gene = res_Pvalb_enrichGO,
OrgDb = org.Mm.eg.db, # Use org.Mm.eg.db for mouse
keyType = "SYMBOL",
ont = "MF", # "BP" (Biological Process), "MF" (Molecular Function), "CC" (Cellular Component)
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# View the first few enriched GO terms
head(MF_GO_Pvalb)
# Run EnrichGO - CC
CC_GO_Pvalb <- enrichGO(
gene = res_Pvalb_enrichGO,
OrgDb = org.Mm.eg.db, # Use org.Mm.eg.db for mouse
keyType = "SYMBOL",
ont = "CC", # "BP" (Biological Process), "MF" (Molecular Function), "CC" (Cellular Component)
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# View the first few enriched GO terms
head(CC_GO_Pvalb)
# Save as Excel file
write.xlsx(
list(
BP = as.data.frame(BP_GO_Pvalb),
MF = as.data.frame(MF_GO_Pvalb),
CC = as.data.frame(CC_GO_Pvalb)
),
file = "~/Risperidone_Proj/R_analysis/EnrichGO/Proj4_Pvalb_GO.xlsx"
)
```
EnrichGO - PV-IN (up/down separate)
```{r}
# Convert the rownames into a new column of gene symbols
res_df_Pvalb <- as.data.frame(res_Pvalb)
res_df_Pvalb$gene <- rownames(res_df_Pvalb)
# Upregulated genes (padj < 0.05 & log2FC > 0.5)
up_Pvalb <- res_df_Pvalb %>%
dplyr::filter(padj < 0.05 & log2FoldChange > 0.5) %>%
dplyr::pull(gene)
# Downregulated genes (padj < 0.05 & log2FC < -0.5)
down_Pvalb <- res_df_Pvalb %>%
dplyr::filter(padj < 0.05, log2FoldChange < -0.5) %>%
dplyr::pull(gene)
## BP
# Run GO enrichment
BP_up_Pvalb <- enrichGO(gene = up_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
BP_down_Pvalb <- enrichGO(gene = down_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_BP_up_Pvalb <- as.data.frame(BP_up_Pvalb)
df_BP_down_Pvalb <- as.data.frame(BP_down_Pvalb)
# Prepare subset for plotting (top 20 each)
df_BP_up_Pvalb_plot <- df_BP_up_Pvalb[1:min(20, nrow(df_BP_up_Pvalb)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_BP_down_Pvalb_plot <- df_BP_down_Pvalb[1:min(20, nrow(df_BP_down_Pvalb)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_BP_Pvalb <- bind_rows(df_BP_up_Pvalb_plot, df_BP_down_Pvalb_plot)
# Plot
ggplot(combined_BP_Pvalb, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- PV-IN (BP)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
# Plot - just downregulated
ggplot(df_BP_down_Pvalb_plot, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "PV-IN (BP)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
legend.position = "none",
axis.text.y = element_text(size = 14), # bigger GO term labels
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, face = "bold", hjust = 0.5)
)
## MF
# Run GO enrichment
MF_up_Pvalb <- enrichGO(gene = up_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
MF_down_Pvalb <- enrichGO(gene = down_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_MF_up_Pvalb <- as.data.frame(MF_up_Pvalb)
df_MF_down_Pvalb <- as.data.frame(MF_down_Pvalb)
# Prepare subset for plotting (top 20 each)
df_MF_up_Pvalb_plot <- df_MF_up_Pvalb[1:min(20, nrow(df_MF_up_Pvalb)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_MF_down_Pvalb_plot <- df_MF_down_Pvalb[1:min(20, nrow(df_MF_down_Pvalb)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_MF_Pvalb <- bind_rows(df_MF_up_Pvalb_plot, df_MF_down_Pvalb_plot)
# Plot
ggplot(combined_MF_Pvalb, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- PV-IN (MF)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
# Plot - just downregulated
ggplot(df_MF_down_Pvalb_plot, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "PV-IN (MF)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
legend.position = "none",
axis.text.y = element_text(size = 14), # bigger GO term labels
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, face = "bold", hjust = 0.5)
)
## CC
# Run GO enrichment
CC_up_Pvalb <- enrichGO(gene = up_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
CC_down_Pvalb <- enrichGO(gene = down_Pvalb, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_CC_up_Pvalb <- as.data.frame(CC_up_Pvalb)
df_CC_down_Pvalb <- as.data.frame(CC_down_Pvalb)
# Prepare subset for plotting (top 20 each)
df_CC_up_Pvalb_plot <- df_CC_up_Pvalb[1:min(20, nrow(df_CC_up_Pvalb)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_CC_down_Pvalb_plot <- df_CC_down_Pvalb[1:min(20, nrow(df_CC_down_Pvalb)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_CC_Pvalb <- bind_rows(df_CC_up_Pvalb_plot, df_CC_down_Pvalb_plot)
# Plot
ggplot(combined_CC_Pvalb, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- PV-IN (CC)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
# Save all enrichment results to Excel
# Create workbook
wb <- createWorkbook()
# Add BP data
addWorksheet(wb, "Up_BP")
addWorksheet(wb, "Down_BP")
writeData(wb, "Up_BP", df_BP_up_Pvalb)
writeData(wb, "Down_BP", df_BP_down_Pvalb)
addWorksheet(wb, "Up_MF")
addWorksheet(wb, "Down_MF")
writeData(wb, "Up_MF", df_MF_up_Pvalb)
writeData(wb, "Down_MF", df_MF_down_Pvalb)
addWorksheet(wb, "Up_CC")
addWorksheet(wb, "Down_CC")
writeData(wb, "Up_CC", df_CC_up_Pvalb)
writeData(wb, "Down_CC", df_CC_down_Pvalb)
# Save excel file
saveWorkbook(wb, "~/Risperidone_Proj/R_analysis/EnrichGO/Proj4_Pvalb_updown.xlsx", overwrite = TRUE)
```
EnrichGO - Glut (up/down separate)
```{r}
# Convert the rownames into a new column of gene symbols
res_df_Glut <- as.data.frame(res_Glut)
res_df_Glut$gene <- rownames(res_df_Glut)
# Upregulated genes (padj < 0.05 & log2FC > 0.5)
up_Glut <- res_df_Glut %>%
dplyr::filter(padj < 0.05, log2FoldChange > 0.5) %>%
dplyr::pull(gene)
# Downregulated genes (padj < 0.05 & log2FC < -0.5)
down_Glut <- res_df_Glut %>%
dplyr::filter(padj < 0.05, log2FoldChange < -0.5) %>%
dplyr::pull(gene)
## BP
# Run GO enrichment
BP_up_Glut <- enrichGO(gene = up_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
BP_down_Glut <- enrichGO(gene = down_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_BP_up_Glut <- as.data.frame(BP_up_Glut)
df_BP_down_Glut <- as.data.frame(BP_down_Glut)
# Prepare subset for plotting (top 20 each)
df_BP_up_Glut_plot <- df_BP_up_Glut[1:min(20, nrow(df_BP_up_Glut)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_BP_down_Glut_plot <- df_BP_down_Glut[1:min(20, nrow(df_BP_down_Glut)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_BP_Glut <- bind_rows(df_BP_up_Glut_plot, df_BP_down_Glut_plot)
# Plot
ggplot(combined_BP_Glut, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- Glut (BP)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
## MF
# Run GO enrichment
MF_up_Glut <- enrichGO(gene = up_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
MF_down_Glut <- enrichGO(gene = down_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_MF_up_Glut <- as.data.frame(MF_up_Glut)
df_MF_down_Glut <- as.data.frame(MF_down_Glut)
# Prepare subset for plotting (top 20 each)
df_MF_up_Glut_plot <- df_MF_up_Glut[1:min(20, nrow(df_MF_up_Glut)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_MF_down_Glut_plot <- df_MF_down_Glut[1:min(20, nrow(df_MF_down_Glut)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_MF_Glut <- bind_rows(df_MF_up_Glut_plot, df_MF_down_Glut_plot)
# Plot
ggplot(combined_MF_Glut, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- Glut (MF)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
# Plot - just downregulated
ggplot(df_MF_down_Glut_plot, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "Glut (MF)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
legend.position = "none",
axis.text.y = element_text(size = 14), # bigger GO term labels
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, face = "bold", hjust = 0.5)
)
## CC
# Run GO enrichment
CC_up_Glut <- enrichGO(gene = up_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
CC_down_Glut <- enrichGO(gene = down_Glut, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_CC_up_Glut <- as.data.frame(CC_up_Glut)
df_CC_down_Glut <- as.data.frame(CC_down_Glut)
# Prepare subset for plotting (top 20 each)
df_CC_up_Glut_plot <- df_CC_up_Glut[1:min(20, nrow(df_CC_up_Glut)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_CC_down_Glut_plot <- df_CC_down_Glut[1:min(20, nrow(df_CC_down_Glut)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_CC_Glut <- bind_rows(df_CC_up_Glut_plot, df_CC_down_Glut_plot)
# Plot
ggplot(combined_CC_Glut, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Term", y = "-log10 Adjusted p-value", title = "GO Enrichment- Glut (CC)") +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal()
# Save all enrichment results to Excel
# Create workbook
wb <- createWorkbook()
# Add BP data
addWorksheet(wb, "Up_BP")
addWorksheet(wb, "Down_BP")
writeData(wb, "Up_BP", df_BP_up_Glut)
writeData(wb, "Down_BP", df_BP_down_Glut)
addWorksheet(wb, "Up_MF")
addWorksheet(wb, "Down_MF")
writeData(wb, "Up_MF", df_MF_up_Glut)
writeData(wb, "Down_MF", df_MF_down_Glut)
addWorksheet(wb, "Up_CC")
addWorksheet(wb, "Down_CC")
writeData(wb, "Up_CC", df_CC_up_Glut)
writeData(wb, "Down_CC", df_CC_down_Glut)
# Save excel file
saveWorkbook(wb, "~/Risperidone_Proj/R_analysis/EnrichGO/Proj4_Glut_updown.xlsx", overwrite = TRUE)
```
EnrichGO - Sst-IN (up/down separate)
```{r}
# Convert the rownames into a new column of gene symbols
res_df_Sst <- as.data.frame(res_Sst)
res_df_Sst$gene <- rownames(res_df_Sst)
# Upregulated genes (padj < 0.05 & log2FC > 0.5)
up_Sst <- res_df_Sst %>%
dplyr::filter(padj < 0.05, log2FoldChange > 0.5) %>%
dplyr::pull(gene)
# Downregulated genes (padj < 0.05 & log2FC < -0.5)
down_Sst <- res_df_Sst %>%
dplyr::filter(padj < 0.05, log2FoldChange < -0.5) %>%
dplyr::pull(gene)
## BP
# Run GO enrichment
BP_up_Sst <- enrichGO(gene = up_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
BP_down_Sst <- enrichGO(gene = down_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_BP_up_Sst <- as.data.frame(BP_up_Sst)
df_BP_down_Sst <- as.data.frame(BP_down_Sst)
# Prepare subset for plotting (top 20 each)
df_BP_up_Sst_plot <- df_BP_up_Sst[1:min(20, nrow(df_BP_up_Sst)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_BP_down_Sst_plot <- df_BP_down_Sst[1:min(20, nrow(df_BP_down_Sst)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_BP_Sst <- bind_rows(df_BP_up_Sst_plot, df_BP_down_Sst_plot)
# Plot
ggplot(combined_BP_Sst, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(
x = "GO Term",
y = "-log10 Adjusted p-value",
title = "Sst-IN (BP)"
) +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
axis.text.y = element_text(size = 12), # increase GO term text size
axis.text.x = element_text(size = 10), # adjust x-axis text size
axis.title = element_text(size = 12) # increase axis title size
)
## MF
# Run GO enrichment
MF_up_Sst <- enrichGO(gene = up_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
MF_down_Sst <- enrichGO(gene = down_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "MF", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_MF_up_Sst <- as.data.frame(MF_up_Sst)
df_MF_down_Sst <- as.data.frame(MF_down_Sst)
# Prepare subset for plotting (top 20 each)
df_MF_up_Sst_plot <- df_MF_up_Sst[1:min(20, nrow(df_MF_up_Sst)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_MF_down_Sst_plot <- df_MF_down_Sst[1:min(20, nrow(df_MF_down_Sst)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_MF_Sst <- bind_rows(df_MF_up_Sst_plot, df_MF_down_Sst_plot)
# Plot
ggplot(combined_MF_Sst, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(
x = "GO Term",
y = "-log10 Adjusted p-value",
title = "Sst-IN (MF)"
) +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
axis.text.y = element_text(size = 12), # increase GO term text size
axis.text.x = element_text(size = 10), # adjust x-axis text size
axis.title = element_text(size = 12) # increase axis title size
)
## CC
# Run GO enrichment
CC_up_Sst <- enrichGO(gene = up_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
CC_down_Sst <- enrichGO(gene = down_Sst, OrgDb = org.Mm.eg.db, keyType = "SYMBOL",
ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
# Convert results to data frame
df_CC_up_Sst <- as.data.frame(CC_up_Sst)
df_CC_down_Sst <- as.data.frame(CC_down_Sst)
# Prepare subset for plotting (top 20 each)
df_CC_up_Sst_plot <- df_CC_up_Sst[1:min(20, nrow(df_CC_up_Sst)), ] %>%
mutate(Direction = "Up", log_padj = -log10(p.adjust))
df_CC_down_Sst_plot <- df_CC_down_Sst[1:min(20, nrow(df_CC_down_Sst)), ] %>%
mutate(Direction = "Down", log_padj = log10(p.adjust)) # keep positive for plotting
combined_CC_Sst <- bind_rows(df_CC_up_Sst_plot, df_CC_down_Sst_plot)
# Plot
ggplot(combined_CC_Sst, aes(x = reorder(Description, log_padj), y = log_padj, fill = Direction)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(
x = "GO Term",
y = "-log10 Adjusted p-value",
title = "Sst-IN (CC)"
) +
scale_fill_manual(values = c("Up" = "tomato", "Down" = "skyblue")) +
theme_minimal() +
theme(
axis.text.y = element_text(size = 12), # increase GO term text size
axis.text.x = element_text(size = 10), # adjust x-axis text size
axis.title = element_text(size = 12) # increase axis title size
)
# Save all enrichment results to Excel
# Create workbook
wb <- createWorkbook()
# Add BP data
addWorksheet(wb, "Up_BP")
addWorksheet(wb, "Down_BP")
writeData(wb, "Up_BP", df_BP_up_Sst)
writeData(wb, "Down_BP", df_BP_down_Sst)
addWorksheet(wb, "Up_MF")
addWorksheet(wb, "Down_MF")
writeData(wb, "Up_MF", df_MF_up_Sst)
writeData(wb, "Down_MF", df_MF_down_Sst)
addWorksheet(wb, "Up_CC")
addWorksheet(wb, "Down_CC")
writeData(wb, "Up_CC", df_CC_up_Sst)
writeData(wb, "Down_CC", df_CC_down_Sst)
# Save excel file
saveWorkbook(wb, "~/Risperidone_Proj/R_analysis/EnrichGO/Proj4_Sst_updown.xlsx", overwrite = TRUE)
```
hdWGCNA
```{r}
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345)
# optionally enable multithreading - allowing parallel execution with up to 8 working processes.
enableWGCNAThreads(nThreads = 8)
# Set up Seurat object for WGCNA - do not subset this seurat object after this has been run.
Mcortex_QC_WGCNA <- SetupForWGCNA(
Mcortex_QC,
gene_select = "fraction", # the gene selection approach, can also select variable or custom
fraction = 0.05, # genes expressed in 5% of cells. Fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "hdWGCNA" # the name of the hdWGCNA experiment
)
# Extract sample_id from barcode names for the Seurat object
Mcortex_QC_WGCNA$sample_id <- sapply(strsplit(Cells(Mcortex_QC_WGCNA), "_"), `[`, 1)
# construct metacells in each group
Mcortex_QC_WGCNA <- MetacellsByGroups(
seurat_obj = Mcortex_QC_WGCNA,
group.by = c("celltype1", "sample_id"), # specify the columns in [email hidden] to group by
reduction = 'pca', # select the dimensionality reduction to perform KNN on
k = 65, # nearest-neighbors parameter
max_shared = 10, # maximum number of shared cells between two metacells
ident.group = 'celltype1' # set the Idents of the metacell seurat object
)
# normalize metacell expression matrix:
Mcortex_QC_WGCNA <- NormalizeMetacells(Mcortex_QC_WGCNA)
# View what layers exist in the WGCNA Seurat object, and check which default assay.
DefaultAssay(Mcortex_QC_WGCNA)
DefaultAssay(Mcortex_QC_WGCNA) <- "SCT"
Layers(Mcortex_QC_WGCNA[["SCT"]])
saveRDS(Mcortex_QC_WGCNA, file='~/Risperidone_Proj/R_analysis/RDS_files/hdWGCNA_object.rds')
Mcortex_QC_WGCNA <- readRDS("~/Risperidone_Proj/R_analysis/RDS_files/hdWGCNA_object.rds")
```
hdWGCNA continued - analysis of PV-IN
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_WGCNA <- SetDatExpr(
Mcortex_QC_WGCNA,
group_name = "GABAergic Pvalb",
group.by = "celltype1",
assay = "SCT",
layer = "data"
)
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_WGCNA <- TestSoftPowers(
Mcortex_QC_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list <- PlotSoftPowers(Mcortex_QC_WGCNA)
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
# Construct co-expression network
Mcortex_QC_WGCNA <- ConstructNetwork(
Mcortex_QC_WGCNA,
tom_name = 'PV-IN') # name of the topoligical overlap matrix written to disk
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_WGCNA <- ScaleData(Mcortex_QC_WGCNA, features=VariableFeatures(Mcortex_QC_WGCNA))
Mcortex_QC_WGCNA <- ModuleEigengenes(
Mcortex_QC_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_WGCNA)
# To get non-harmonized module eigengenes (not run)
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=FALSE)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_WGCNA <- ModuleConnectivity(
Mcortex_QC_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Pvalb')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_WGCNA <- ResetModuleNames(
Mcortex_QC_WGCNA,
new_name = "PV-IN-M")
# Change color of modules
# get the module table
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- unique(modules$module)
# make a table of the module-color pairings
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
# print the dataframe
mod_colors_df
# load MetBrewer color scheme package
library(MetBrewer)
# get a table of just the module and it's unique color
mod_color_df <- GetModules(Mcortex_QC_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
# the number of unique modules (subtract 1 because the grey module stays grey):
n_mods <- nrow(mod_color_df) - 1
# using the "Signac" palette from metbrewer, selecting for the number of modules
new_colors <- paste0(met.brewer("Signac", n=n_mods))
# reset the module colors
Mcortex_QC_WGCNA <- ResetModuleColors(Mcortex_QC_WGCNA, new_colors)
# Plot dendogram
PlotDendrogram(Mcortex_QC_WGCNA, main='hdWGCNA Dendrogram for PV-IN')
# plot genes ranked by kME for each module
PlotKMEs(Mcortex_QC_WGCNA, ncol=3)
# get the module assignment table:
modules <- GetModules(Mcortex_QC_WGCNA) %>% subset(module != 'grey')
# show the first 6 columns:
head(modules[,1:6])
# save to Excel
write_xlsx(modules, "~/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_module_assignments.xlsx")
# get hub genes
hub_df <- GetHubGenes(Mcortex_QC_WGCNA, n_hubs = 100)
head(hub_df)
# save to Excel
write_xlsx(hub_df, "~/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_top100_hubgenes.xlsx")
```
Visualization of hdWGCNA results - PV-IN
```{r}
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
Mcortex_QC_WGCNA,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=3)
# plot module correlagram
ModuleCorrelogram(Mcortex_QC_WGCNA)
# Dot plot of hME's
# Get hME's and module tables again
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=TRUE)
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
[email hidden] <- cbind([email hidden], MEs)
# plot with Seurat's DotPlot function
p <- DotPlot(Mcortex_QC_WGCNA, features=mods, group.by = 'celltype1')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
# Individual module network plot
ModuleNetworkPlot(
Mcortex_QC_WGCNA,
outdir='ModuleNetworks_PV-IN', # new folder name
n_inner = 10, # number of genes in inner ring
n_outer = 15, # number of genes in outer ring
n_conns = Inf, # show all of the connections
plot_size=c(5,5), # larger plotting area
vertex.label.cex=1 # font size
)
# UMAP with network
# Increased allowed globals size
options(future.globals.maxSize = 20 * 1024^3) # allow up to 20 GB
# Run Module UMAP
Mcortex_QC_WGCNA <- RunModuleUMAP(
Mcortex_QC_WGCNA,
n_hubs = 5, # number of hub genes to include for the UMAP embedding
n_neighbors=10, # neighbors parameter for UMAP
min_dist=0.3 # min distance between points in UMAP space
)
# Plot UMAP of genes and their co-expression relationship
ModuleUMAPPlot(
Mcortex_QC_WGCNA,
edge.alpha=0.15,
sample_edges=TRUE, # If true, downsampling edges
edge_prop=0.1, # proportion of edges to sample (10% here)
label_hubs=1, # how many hub genes to plot per module?
vertex.label.cex = 0.5,
)
# hubgene network
HubGeneNetworkPlot(
Mcortex_QC_WGCNA,
n_hubs = 3, n_other=5,
edge_prop = 1,
mods = 'all',
vertex.label.cex = 0.75,
hub.vertex.size = 4,
other.vertex.size = 1
)
```
Correlation of hdGCNA results with treatment variables - PV-IN
```{r}
# --- Create a modified version of hdWGCNA::ModuleTraitCorrelation ---
ModuleTraitCorrelation_mod <- function(seurat_obj, traits, group.by = NULL, features = "hMEs",
cor_method = "pearson", subset_by = NULL, subset_groups = NULL,
wgcna_name = NULL, ...) {
if (is.null(wgcna_name)) {
wgcna_name <- seurat_obj@misc$active_wgcna
}
hdWGCNA:::CheckWGCNAName(seurat_obj, wgcna_name)
if (features == "hMEs") {
MEs <- hdWGCNA:::GetMEs(seurat_obj, TRUE, wgcna_name)
} else if (features == "MEs") {
MEs <- hdWGCNA:::GetMEs(seurat_obj, FALSE, wgcna_name)
} else if (features == "scores") {
MEs <- hdWGCNA:::GetModuleScores(seurat_obj, wgcna_name)
} else {
stop("Invalid feature selection. Valid choices: hMEs, MEs, scores, average")
}
if (!is.null(subset_by)) {
print("subsetting")
seurat_full <- seurat_obj
MEs <- MEs[[email hidden][[subset_by]] %in% subset_groups, ]
seurat_obj <- seurat_obj[, [email hidden][[subset_by]] %in% subset_groups]
}
if (sum(traits %in% colnames([email hidden])) != length(traits)) {
stop(paste("Some of the provided traits were not found in the Seurat obj:",
paste(traits[!(traits %in% colnames([email hidden]))], collapse = ", ")))
}
if (is.null(group.by)) {
group.by <- "temp_ident"
seurat_obj$temp_ident <- Idents(seurat_obj)
}
valid_types <- c("numeric", "factor", "integer")
data_types <- sapply(traits, function(x) class([email hidden][, x]))
if (!all(data_types %in% valid_types)) {
incorrect <- traits[!(data_types %in% valid_types)]
stop(paste0("Invalid data types for ", paste(incorrect, collapse = ", "),
". Accepted data types are numeric, factor, integer."))
}
if (any(data_types == "factor")) {
factor_traits <- traits[data_types == "factor"]
for (tr in factor_traits) {
warning(paste0("Trait ", tr, " is a factor with levels ",
paste0(levels([email hidden][, tr]), collapse = ", "),
". Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?"))
}
}
modules <- hdWGCNA:::GetModules(seurat_obj, wgcna_name)
mods <- levels(modules$module)
mods <- mods[mods != "grey"]
trait_df <- [email hidden][, traits, drop = FALSE]
# Modified: Correctly detect single-trait case
if (length(traits) == 1) {
trait_df <- data.frame(x = trait_df)
colnames(trait_df) <- traits
}
if (any(data_types == "factor")) {
factor_traits <- traits[data_types == "factor"]
for (tr in factor_traits) {
trait_df[, tr] <- as.numeric(trait_df[, tr])
}
}
cor_list <- list()
pval_list <- list()
fdr_list <- list()
temp <- Hmisc::rcorr(as.matrix(trait_df), as.matrix(MEs), type = cor_method)
cur_cor <- temp$r[traits, mods, drop = FALSE]
cur_p <- temp$P[traits, mods, drop = FALSE]
p_df <- reshape2::melt(cur_p)
if (length(traits) == 1) {
tmp <- rep(mods, length(traits))
tmp <- factor(tmp, levels = mods)
tmp <- tmp[order(tmp)]
p_df$Var1 <- traits
p_df$Var2 <- tmp
rownames(p_df) <- 1:nrow(p_df)
p_df <- dplyr::select(p_df, c(Var1, Var2, value))
}
p_df <- p_df %>%
dplyr::mutate(fdr = p.adjust(value, method = "fdr")) %>%
dplyr::select(c(Var1, Var2, fdr))
cur_fdr <- reshape2::dcast(p_df, Var1 ~ Var2, value.var = "fdr")
rownames(cur_fdr) <- cur_fdr$Var1
cur_fdr <- cur_fdr[, -1, drop = FALSE]
cor_list[["all_cells"]] <- cur_cor
pval_list[["all_cells"]] <- cur_p
fdr_list[["all_cells"]] <- cur_fdr
trait_df <- cbind(trait_df, [email hidden][, group.by])
colnames(trait_df)[ncol(trait_df)] <- "group"
MEs <- cbind(as.data.frame(MEs), [email hidden][, group.by])
colnames(MEs)[ncol(MEs)] <- "group"
if (class([email hidden][, group.by]) == "factor") {
group_names <- levels([email hidden][, group.by])
} else {
group_names <- levels(as.factor([email hidden][, group.by]))
}
trait_list <- dplyr::group_split(trait_df, group, .keep = FALSE)
ME_list <- dplyr::group_split(MEs, group, .keep = FALSE)
names(trait_list) <- group_names
names(ME_list) <- group_names
for (i in names(trait_list)) {
temp <- Hmisc::rcorr(as.matrix(trait_list[[i]]), as.matrix(ME_list[[i]]))
cur_cor <- temp$r[traits, mods, drop = FALSE]
cur_p <- temp$P[traits, mods, drop = FALSE]
p_df <- reshape2::melt(cur_p)
if (length(traits) == 1) {
tmp <- rep(mods, length(traits))
tmp <- factor(tmp, levels = mods)
tmp <- tmp[order(tmp)]
p_df$Var1 <- traits
p_df$Var2 <- tmp
rownames(p_df) <- 1:nrow(p_df)
p_df <- dplyr::select(p_df, c(Var1, Var2, value))
}
p_df <- p_df %>%
dplyr::mutate(fdr = p.adjust(value, method = "fdr")) %>%
dplyr::select(c(Var1, Var2, fdr))
cur_fdr <- reshape2::dcast(p_df, Var1 ~ Var2, value.var = "fdr")
rownames(cur_fdr) <- cur_fdr$Var1
cur_fdr <- cur_fdr[, -1, drop = FALSE]
cor_list[[i]] <- cur_cor
pval_list[[i]] <- cur_p
fdr_list[[i]] <- as.matrix(cur_fdr)
}
mt_cor <- list(cor = cor_list, pval = pval_list, fdr = fdr_list)
if (!is.null(subset_by)) {
seurat_full <- hdWGCNA:::SetModuleTraitCorrelation(seurat_full, mt_cor, wgcna_name)
seurat_obj <- seurat_full
} else {
seurat_obj <- hdWGCNA:::SetModuleTraitCorrelation(seurat_obj, mt_cor, wgcna_name)
}
seurat_obj
}
# convert treatment to factor
Mcortex_QC_WGCNA$Tx <- as.factor(Mcortex_QC_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c('Tx')
Mcortex_QC_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_WGCNA,
traits = cur_traits,
group.by='celltype1'
)
# get the mt-correlation results
mt_cor <- GetModuleTraitCorrelation(Mcortex_QC_WGCNA)
names(mt_cor)
names(mt_cor$cor)
head(mt_cor$cor$`GABAergic Pvalb`[,1:5])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor$cor[[celltype]]
pval_mat <- mt_cor$pval[[celltype]]
fdr_mat <- mt_cor$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, ""))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, ""), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.5,
combine=TRUE
)
```
hdWGCNA continued - enrichment analysis- PV-IN
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_WGCNA <- RunEnrichr(
Mcortex_QC_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df <- GetEnrichrTable(Mcortex_QC_WGCNA)
# look at the results
head(enrich_df)
# save to Excel
write_xlsx(enrich_df, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_enrichR_100.xlsx")
# make GO term bar plots - didn't re-run
EnrichrBarPlot(
Mcortex_QC_WGCNA,
outdir = "enrichr_plots", # name of output directory
n_terms = 10, # number of enriched terms to show (sometimes more are shown if there are ties)
plot_size = c(5,7), # width, height of the output .pdfs
logscale=TRUE # do you want to show the enrichment as a log scale?
)
# enrichr dotplot - did not re-run
EnrichrDotPlot(
Mcortex_QC_WGCNA,
mods = "all",
database = "GO_Molecular_Function_2023",
n_terms = 2,
term_size = 13,
p_adj = TRUE
) +
scale_color_stepsn(colors = rev(viridis::magma(256))) +
theme(
axis.text.x = element_text(size = 13), # X-axis tick label size
)
# Define the order of modules you want (e.g., numeric order)
module_order <- c("PV-IN-M1", "PV-IN-M2", "PV-IN-M3", "PV-IN-M4", "PV-IN-M5", "PV-IN-M6", "PV-IN-M7", "PV-IN-M8", "PV-IN-M9")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_clean <- enrich_df %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_clean <- enrich_df_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("PV-IN - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
```
hdWGCNA continued - analysis of Sst-IN
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_WGCNA <- SetDatExpr(
Mcortex_QC_WGCNA,
group_name = "GABAergic Sst",
group.by = "celltype1",
assay = "SCT",
layer = "data"
)
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_WGCNA <- TestSoftPowers(
Mcortex_QC_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list <- PlotSoftPowers(Mcortex_QC_WGCNA)
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
# Construct co-expression network
Mcortex_QC_WGCNA <- ConstructNetwork(
Mcortex_QC_WGCNA,
soft_power = 18,
tom_name = 'Sst-IN') # name of the topological overlap matrix written to disk
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_WGCNA <- ScaleData(Mcortex_QC_WGCNA, features=VariableFeatures(Mcortex_QC_WGCNA))
Mcortex_QC_WGCNA <- ModuleEigengenes(
Mcortex_QC_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_WGCNA)
# To get non-harmonized module eigengenes (not run)
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=FALSE)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_WGCNA <- ModuleConnectivity(
Mcortex_QC_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Sst')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_WGCNA <- ResetModuleNames(
Mcortex_QC_WGCNA,
new_name = "Sst-IN-M")
# Change color of modules
# get the module table
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- unique(modules$module)
# make a table of the module-color pairings
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
# print the dataframe
mod_colors_df
# load MetBrewer color scheme package
library(MetBrewer)
# get a table of just the module and it's unique color
mod_color_df <- GetModules(Mcortex_QC_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
# the number of unique modules (subtract 1 because the grey module stays grey):
n_mods <- nrow(mod_color_df) - 1
# using the "Signac" palette from metbrewer, selecting for the number of modules
new_colors <- paste0(met.brewer("Signac", n=n_mods))
# reset the module colors
Mcortex_QC_WGCNA <- ResetModuleColors(Mcortex_QC_WGCNA, new_colors)
# Plot dendogram
PlotDendrogram(Mcortex_QC_WGCNA, main='hdWGCNA Dendrogram for Sst-IN')
# plot genes ranked by kME for each module
PlotKMEs(Mcortex_QC_WGCNA, ncol=3)
# get the module assignment table:
modules <- GetModules(Mcortex_QC_WGCNA) %>% subset(module != 'grey')
# show the first 6 columns:
head(modules[,1:6])
# save to Excel
write_xlsx(modules, "~/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_module_assignments.xlsx")
# get hub genes
hub_df <- GetHubGenes(Mcortex_QC_WGCNA, n_hubs = 100)
head(hub_df)
# save to Excel
write_xlsx(hub_df, "~/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_top100_hubgenes.xlsx")
```
Visualization of hdWGCNA results - Sst-IN
```{r}
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
Mcortex_QC_WGCNA,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch the plot together with patchwork
wrap_plots(plot_list, ncol=4)
# plot module correlagram
ModuleCorrelogram(Mcortex_QC_WGCNA)
# Dot plot of hME's
# Get hME's and module tables again
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=TRUE)
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
[email hidden] <- cbind([email hidden], MEs)
# plot with Seurat's DotPlot function
p <- DotPlot(Mcortex_QC_WGCNA, features=mods, group.by = 'celltype1')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
# Individual module network plot
ModuleNetworkPlot(
Mcortex_QC_WGCNA,
outdir='ModuleNetworks_Sst-IN', # new folder name
n_inner = 10, # number of genes in inner ring
n_outer = 15, # number of genes in outer ring
n_conns = Inf, # show all of the connections
plot_size=c(5,5), # larger plotting area
vertex.label.cex=1 # font size
)
# UMAP with network
# Increased allowed globals size
options(future.globals.maxSize = 20 * 1024^3) # allow up to 20 GB
# Run Module UMAP
Mcortex_QC_WGCNA <- RunModuleUMAP(
Mcortex_QC_WGCNA,
n_hubs = 5, # number of hub genes to include for the UMAP embedding
n_neighbors=10, # neighbors parameter for UMAP
min_dist=0.3 # min distance between points in UMAP space
)
# Plot UMAP of genes and their co-expression relationship
ModuleUMAPPlot(
Mcortex_QC_WGCNA,
edge.alpha=0.15,
sample_edges=TRUE, # If true, downsampling edges
edge_prop=0.1, # proportion of edges to sample (10% here)
label_hubs=1, # how many hub genes to plot per module?
vertex.label.cex = 0.5,
)
# hubgene network
HubGeneNetworkPlot(
Mcortex_QC_WGCNA,
n_hubs = 3, n_other=5,
edge_prop = 1,
mods = 'all',
vertex.label.cex = 0.75,
hub.vertex.size = 4,
other.vertex.size = 1
)
```
Correlation of hdGCNA results with treatment variables - Sst-IN
```{r}
# convert treatment to factor
Mcortex_QC_WGCNA$Tx <- as.factor(Mcortex_QC_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c('Tx')
Mcortex_QC_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_WGCNA,
traits = cur_traits,
group.by='celltype1'
)
# get the mt-correlation results
mt_cor <- GetModuleTraitCorrelation(Mcortex_QC_WGCNA)
names(mt_cor)
names(mt_cor$cor)
head(mt_cor$cor$`GABAergic Sst`[,1:5])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor$cor[[celltype]]
pval_mat <- mt_cor$pval[[celltype]]
fdr_mat <- mt_cor$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, ""))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, ""), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.5,
combine=TRUE
)
```
hdWGCNA continued - enrichment analysis- Sst-IN
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_WGCNA <- RunEnrichr(
Mcortex_QC_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df <- GetEnrichrTable(Mcortex_QC_WGCNA)
# look at the results
head(enrich_df)
# save to Excel
write_xlsx(enrich_df, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_enrichR_100.xlsx")
# make GO term barplot - did not rerun
EnrichrBarPlot(
Mcortex_QC_WGCNA,
outdir = "enrichr_plots", # name of output directory
n_terms = 10, # number of enriched terms to show (sometimes more are shown if there are ties)
plot_size = c(5,7), # width, height of the output .pdfs
logscale=TRUE # do you want to show the enrichment as a log scale?
)
# enrichr dotplot - did not rerun
EnrichrDotPlot(
Mcortex_QC_WGCNA,
mods = "all",
database = "GO_Molecular_Function_2023",
n_terms = 2,
term_size = 13,
p_adj = TRUE
) +
scale_color_stepsn(colors = rev(viridis::magma(256))) +
theme(
axis.text.x = element_text(size = 12), # X-axis tick label size
)
# Define the order of modules you want (e.g., numeric order)
module_order <- c("Sst-IN-M1", "Sst-IN-M2", "Sst-IN-M3", "Sst-IN-M4", "Sst-IN-M5",
"Sst-IN-M6", "Sst-IN-M7", "Sst-IN-M8", "Sst-IN-M9", "Sst-IN-M10", "Sst-IN-M11")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_clean <- enrich_df %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 55),
# Keep numeric ordering of modules
module = factor(module, levels = module_order),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_clean <- enrich_df_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Sst-IN - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
```
hdWGCNA continued - analysis of Glut
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_WGCNA <- SetDatExpr(
Mcortex_QC_WGCNA,
group_name = "Glutamatergic neurons",
group.by = "celltype1",
assay = "SCT",
layer = "data"
)
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_WGCNA <- TestSoftPowers(
Mcortex_QC_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list <- PlotSoftPowers(Mcortex_QC_WGCNA)
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
# Construct co-expression network
Mcortex_QC_WGCNA <- ConstructNetwork(
Mcortex_QC_WGCNA,
tom_name = 'Glut') # name of the topological overlap matrix written to disk
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_WGCNA <- ScaleData(Mcortex_QC_WGCNA, features=VariableFeatures(Mcortex_QC_WGCNA))
Mcortex_QC_WGCNA <- ModuleEigengenes(
Mcortex_QC_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_WGCNA)
# To get non-harmonized module eigengenes (not run)
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=FALSE)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_WGCNA <- ModuleConnectivity(
Mcortex_QC_WGCNA,
group.by = 'celltype1', group_name = 'Glutamatergic neurons')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_WGCNA <- ResetModuleNames(
Mcortex_QC_WGCNA,
new_name = "Glut-M")
# Change color of modules
# get the module table
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- unique(modules$module)
# make a table of the module-color pairings
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
# print the dataframe
mod_colors_df
# load MetBrewer color scheme package
library(MetBrewer)
# get a table of just the module and it's unique color
mod_color_df <- GetModules(Mcortex_QC_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
# the number of unique modules (subtract 1 because the grey module stays grey):
n_mods <- nrow(mod_color_df) - 1
# using the "Signac" palette from metbrewer, selecting for the number of modules
new_colors <- paste0(met.brewer("Signac", n=n_mods))
# reset the module colors
Mcortex_QC_WGCNA <- ResetModuleColors(Mcortex_QC_WGCNA, new_colors)
# Plot dendogram
PlotDendrogram(Mcortex_QC_WGCNA, main='hdWGCNA Dendrogram for Glutamatergic neurons')
# plot genes ranked by kME for each module
PlotKMEs(Mcortex_QC_WGCNA, ncol=3)
# get the module assignment table:
modules <- GetModules(Mcortex_QC_WGCNA) %>% subset(module != 'grey')
# show the first 6 columns:
head(modules[,1:6])
# save to Excel
write_xlsx(modules, "~/Risperidone_Proj/R_analysis/hdWGCNA/Glut_module_assignments.xlsx")
# get hub genes
hub_df <- GetHubGenes(Mcortex_QC_WGCNA, n_hubs = 100)
head(hub_df)
# save to Excel
write_xlsx(hub_df, "~/Risperidone_Proj/R_analysis/hdWGCNA/Glut_top100_hubgenes.xlsx")
```
Visualization of hdWGCNA results - Glut
```{r}
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
Mcortex_QC_WGCNA,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch the plot together with patchwork
wrap_plots(plot_list, ncol=3)
# plot module correlagram
ModuleCorrelogram(Mcortex_QC_WGCNA)
# Dot plot of hME's
# Get hME's and module tables again
MEs <- GetMEs(Mcortex_QC_WGCNA, harmonized=TRUE)
modules <- GetModules(Mcortex_QC_WGCNA)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
[email hidden] <- cbind([email hidden], MEs)
# plot with Seurat's DotPlot function
p <- DotPlot(Mcortex_QC_WGCNA, features=mods, group.by = 'celltype1')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
# Individual module network plot
ModuleNetworkPlot(
Mcortex_QC_WGCNA,
outdir='ModuleNetworks_Glut', # new folder name
n_inner = 10, # number of genes in inner ring
n_outer = 15, # number of genes in outer ring
n_conns = Inf, # show all of the connections
plot_size=c(5,5), # larger plotting area
vertex.label.cex=1 # font size
)
# UMAP with network
# Increased allowed globals size
options(future.globals.maxSize = 20 * 1024^3) # allow up to 20 GB
# Run Module UMAP
Mcortex_QC_WGCNA <- RunModuleUMAP(
Mcortex_QC_WGCNA,
n_hubs = 5, # number of hub genes to include for the UMAP embedding
n_neighbors=10, # neighbors parameter for UMAP
min_dist=0.3 # min distance between points in UMAP space
)
# Plot UMAP of genes and their co-expression relationship
ModuleUMAPPlot(
Mcortex_QC_WGCNA,
edge.alpha=0.15,
sample_edges=TRUE, # If true, downsampling edges
edge_prop=0.1, # proportion of edges to sample (10% here)
label_hubs=1, # how many hub genes to plot per module?
vertex.label.cex = 0.5,
)
# hubgene network
HubGeneNetworkPlot(
Mcortex_QC_WGCNA,
n_hubs = 3, n_other=5,
edge_prop = 1,
mods = 'all',
vertex.label.cex = 0.75,
hub.vertex.size = 4,
other.vertex.size = 1
)
```
Correlation of hdGCNA results with treatment variables - Glut
```{r}
# convert treatment to factor
Mcortex_QC_WGCNA$Tx <- as.factor(Mcortex_QC_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c('Tx')
Mcortex_QC_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_WGCNA,
traits = cur_traits,
group.by='celltype1'
)
# get the mt-correlation results
mt_cor <- GetModuleTraitCorrelation(Mcortex_QC_WGCNA)
names(mt_cor)
names(mt_cor$cor)
head(mt_cor$cor$`Glutamatergic neurons`[,1:5])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor$cor[[celltype]]
pval_mat <- mt_cor$pval[[celltype]]
fdr_mat <- mt_cor$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, ""))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, ""), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.5,
combine=TRUE
)
```
hdWGCNA continued - enrichment analysis- Glut
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_WGCNA <- RunEnrichr(
Mcortex_QC_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df <- GetEnrichrTable(Mcortex_QC_WGCNA)
# look at the results
head(enrich_df)
# save to Excel
write_xlsx(enrich_df, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_enrichR_100.xlsx")
# make GO term bar plots - didn't re-run
EnrichrBarPlot(
Mcortex_QC_WGCNA,
outdir = "enrichr_plots", # name of output directory
n_terms = 10, # number of enriched terms to show (sometimes more are shown if there are ties)
plot_size = c(5,7), # width, height of the output .pdfs
logscale=TRUE # do you want to show the enrichment as a log scale?
)
# enrichr dotplot - did not re-run
EnrichrDotPlot(
Mcortex_QC_WGCNA,
mods = "all",
database = "GO_Molecular_Function_2023",
n_terms = 2,
term_size = 13,
p_adj = TRUE
) +
scale_color_stepsn(colors = rev(viridis::magma(256))) +
theme(
axis.text.x = element_text(size = 13), # X-axis tick label size
)
# Define the order of modules you want (e.g., numeric order)
module_order <- c("Glut-M1", "Glut-M2", "Glut-M3", "Glut-M4", "Glut-M5")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_clean <- enrich_df %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_clean <- enrich_df_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Glut - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
```
hdWGCNA - sex specific
Set-up
```{r}
# Subset Seurat object into male and female separately
Mcortex_QC_F <- subset(Mcortex_QC, Predicted_Sex %in% "Female")
Mcortex_QC_M <- subset(Mcortex_QC, Predicted_Sex %in% "Male")
# Filter out all genes with count sum less than 10
Mcortex_QC_F <- subset(Mcortex_QC_F, features = rownames(Mcortex_QC_F)[Matrix::rowSums(GetAssayData(Mcortex_QC_F, assay = "RNA", layer = "counts")) >= 10])
Mcortex_QC_M <- subset(Mcortex_QC_M, features = rownames(Mcortex_QC_M)[Matrix::rowSums(GetAssayData(Mcortex_QC_M, assay = "RNA", layer = "counts")) >= 10])
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345)
# optionally enable multithreading - allowing parallel execution with up to 8 working processes.
enableWGCNAThreads(nThreads = 32)
# Set up Seurat object for WGCNA - do not subset this seurat object after this has been run.
Mcortex_QC_F_WGCNA <- SetupForWGCNA(
Mcortex_QC_F,
gene_select = "fraction", # the gene selection approach, can also select variable or custom
fraction = 0.05, # genes expressed in 5% of cells. Fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "hdWGCNA_F" # the name of the hdWGCNA experiment
)
Mcortex_QC_M_WGCNA <- SetupForWGCNA(
Mcortex_QC_M,
gene_select = "fraction", # the gene selection approach, can also select variable or custom
fraction = 0.05, # genes expressed in 5% of cells. Fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "hdWGCNA_M" # the name of the hdWGCNA experiment
)
# Extract sample_id from barcode names for the Seurat object
Mcortex_QC_F_WGCNA$sample_id <- sapply(strsplit(Cells(Mcortex_QC_F_WGCNA), "_"), `[`, 1)
Mcortex_QC_M_WGCNA$sample_id <- sapply(strsplit(Cells(Mcortex_QC_M_WGCNA), "_"), `[`, 1)
# construct metacells in each group
Mcortex_QC_F_WGCNA <- MetacellsByGroups(
seurat_obj = Mcortex_QC_F_WGCNA,
group.by = c("celltype1", "sample_id"), # specify the columns in [email hidden] to group by
reduction = 'pca', # select the dimensionality reduction to perform KNN on
k = 20, # nearest-neighbors parameter (adjust based on cell number)
max_shared = 10, # maximum number of shared cells between two metacells
ident.group = 'celltype1') # set the Idents of the metacell seurat object
Mcortex_QC_M_WGCNA <- MetacellsByGroups(
seurat_obj = Mcortex_QC_M_WGCNA,
group.by = c("celltype1", "sample_id"), # specify the columns in [email hidden] to group by
reduction = 'pca', # select the dimensionality reduction to perform KNN on
k = 20, # nearest-neighbors parameter (adjust based on cell number)
max_shared = 10, # maximum number of shared cells between two metacells
ident.group = 'celltype1') # set the Idents of the metacell seurat object
# normalize metacell expression matrix:
Mcortex_QC_F_WGCNA <- NormalizeMetacells(Mcortex_QC_F_WGCNA)
Mcortex_QC_M_WGCNA <- NormalizeMetacells(Mcortex_QC_M_WGCNA)
# View what layers exist in the WGCNA Seurat object, and check which default assay.
DefaultAssay(Mcortex_QC_M_WGCNA)
DefaultAssay(Mcortex_QC_F_WGCNA) <- "SCT"
Layers(Mcortex_QC_M_WGCNA[["SCT"]])
saveRDS(Mcortex_QC_F_WGCNA, file='~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/RDS_files/Mcortex_QC_F_WGCNA.rds')
saveRDS(Mcortex_QC_M_WGCNA, file='~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/RDS_files/Mcortex_QC_M_WGCNA.rds')
Mcortex_QC_F_WGCNA <- readRDS("~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/RDS_files/Mcortex_QC_F_WGCNA.rds")
Mcortex_QC_M_WGCNA <- readRDS("~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/RDS_files/Mcortex_QC_M_WGCNA.rds")
```
hdWGCNA analysis of Glut - sex specific
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_F_WGCNA <- SetDatExpr(
Mcortex_QC_F_WGCNA,
group_name = "Glutamatergic neurons",
group.by = "celltype1",
assay = "SCT",
layer = "data")
Mcortex_QC_M_WGCNA <- SetDatExpr(
Mcortex_QC_M_WGCNA,
group_name = "Glutamatergic neurons",
group.by = "celltype1",
assay = "SCT",
layer = "data")
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_F_WGCNA <- TestSoftPowers(
Mcortex_QC_F_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
Mcortex_QC_M_WGCNA <- TestSoftPowers(
Mcortex_QC_M_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list_F <- PlotSoftPowers(Mcortex_QC_F_WGCNA)
plot_list_M <- PlotSoftPowers(Mcortex_QC_M_WGCNA)
# assemble with patchwork
wrap_plots(plot_list_F, ncol=2)
wrap_plots(plot_list_M, ncol=2)
# Construct co-expression network
Mcortex_QC_F_WGCNA <- ConstructNetwork(
Mcortex_QC_F_WGCNA,
tom_name = 'Glut_F') # name of the topological overlap matrix written to disk
Mcortex_QC_M_WGCNA <- ConstructNetwork(
Mcortex_QC_M_WGCNA,
tom_name = 'Glut_M')
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_F_WGCNA <- ScaleData(Mcortex_QC_F_WGCNA, features=VariableFeatures(Mcortex_QC_F_WGCNA))
Mcortex_QC_M_WGCNA <- ScaleData(Mcortex_QC_M_WGCNA, features=VariableFeatures(Mcortex_QC_M_WGCNA))
Mcortex_QC_F_WGCNA <- ModuleEigengenes(
Mcortex_QC_F_WGCNA,
group.by.vars="sample_id")
Mcortex_QC_M_WGCNA <- ModuleEigengenes(
Mcortex_QC_M_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_F_WGCNA)
hMEs <- GetMEs(Mcortex_QC_M_WGCNA)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_F_WGCNA <- ModuleConnectivity(
Mcortex_QC_F_WGCNA,
group.by = 'celltype1', group_name = 'Glutamatergic neurons')
Mcortex_QC_M_WGCNA <- ModuleConnectivity(
Mcortex_QC_M_WGCNA,
group.by = 'celltype1', group_name = 'Glutamatergic neurons')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_F_WGCNA@misc$hdWGCNA$wgcna_params
Mcortex_QC_M_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_F_WGCNA <- ResetModuleNames(
Mcortex_QC_F_WGCNA,
new_name = "F_Glut-M")
Mcortex_QC_M_WGCNA <- ResetModuleNames(
Mcortex_QC_M_WGCNA,
new_name = "M_Glut-M")
# Change color of modules - Female
library(MetBrewer)
modules <- GetModules(Mcortex_QC_F_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_F_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_F <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_F_WGCNA <- ResetModuleColors(Mcortex_QC_F_WGCNA, new_colors_F)
# Change color of modules - Male
modules <- GetModules(Mcortex_QC_M_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_M_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_M <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_M_WGCNA <- ResetModuleColors(Mcortex_QC_M_WGCNA, new_colors_M)
# Plot dendogram - did not run
PlotDendrogram(Mcortex_QC_F_WGCNA, main='Glut_F')
PlotDendrogram(Mcortex_QC_M_WGCNA, main='Glut_M')
# plot genes ranked by kME for each module - did not run
PlotKMEs(Mcortex_QC_F_WGCNA, ncol=3)
PlotKMEs(Mcortex_QC_M_WGCNA, ncol=3)
# get the module assignment table:
modules_F <- GetModules(Mcortex_QC_F_WGCNA) %>% subset(module != 'grey')
modules_M <- GetModules(Mcortex_QC_M_WGCNA) %>% subset(module != 'grey')
# save to Excel
write_xlsx(modules_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_module_assignments.xlsx")
write_xlsx(modules_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_module_assignments.xlsx")
# get hub genes
hub_df_F <- GetHubGenes(Mcortex_QC_F_WGCNA, n_hubs = 100)
hub_df_M <- GetHubGenes(Mcortex_QC_M_WGCNA, n_hubs = 100)
# save to Excel
write_xlsx(hub_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_top100_hubgenes.xlsx")
write_xlsx(hub_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_top100_hubgenes.xlsx")
```
hdWGCN module trait correlation - Glut - sex specific
```{r}
# convert treatment to factor
Mcortex_QC_F_WGCNA$Tx <- as.factor(Mcortex_QC_F_WGCNA$Tx)
Mcortex_QC_M_WGCNA$Tx <- as.factor(Mcortex_QC_M_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c("Tx") # single trait is fine now!
Mcortex_QC_F_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_F_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
Mcortex_QC_M_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_M_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
# get the mt-correlation results
mt_cor_F <- GetModuleTraitCorrelation(Mcortex_QC_F_WGCNA)
names(mt_cor_F)
names(mt_cor_F$cor)
head(mt_cor_F$cor$`Glutamatergic neurons`[,1:10])
mt_cor_M <- GetModuleTraitCorrelation(Mcortex_QC_M_WGCNA)
names(mt_cor_M)
names(mt_cor_M$cor)
head(mt_cor_M$cor$`Glutamatergic neurons`[,1:10])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor_F$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor_F$cor[[celltype]]
pval_mat <- mt_cor_F$pval[[celltype]]
fdr_mat <- mt_cor_F$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# For male
wb <- createWorkbook()
for (celltype in names(mt_cor_M$cor)) {
cor_mat <- mt_cor_M$cor[[celltype]]
pval_mat <- mt_cor_M$pval[[celltype]]
fdr_mat <- mt_cor_M$fdr[[celltype]]
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_F_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
PlotModuleTraitCorrelation(
Mcortex_QC_M_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
```
hdWGCNA enrichment analysis- Glut - sex specific
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_F_WGCNA <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F <- GetEnrichrTable(Mcortex_QC_F_WGCNA)
enrich_df_M <- GetEnrichrTable(Mcortex_QC_M_WGCNA)
# save to Excel
write_xlsx(enrich_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_enrichR.xlsx")
write_xlsx(enrich_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_enrichR.xlsx")
# Define the order of modules you want (e.g., numeric order)
module_order_F <- c("F_Glut-M1", "F_Glut-M2", "F_Glut-M3", "F_Glut-M4",
"F_Glut-M5", "F_Glut-M6", "F_Glut-M7", "F_Glut-M8",
"F_Glut-M9", "F_Glut-M10")
module_order_M <- c("M_Glut-M1", "M_Glut-M2", "M_Glut-M3", "M_Glut-M4",
"M_Glut-M5", "M_Glut-M6", "M_Glut-M7", "M_Glut-M8",
"M_Glut-M9", "M_Glut-M10")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_F_clean <- enrich_df_F %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean <- enrich_df_M %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean <- enrich_df_F_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean <- enrich_df_M_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Glut_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
q <- ggplot(enrich_df_M_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Glut_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(q)
## Re-run enrichR for max 100 genes saved separately, graph not made)
# perform enrichment tests
Mcortex_QC_F_WGCNA_100 <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA_100 <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F_100 <- GetEnrichrTable(Mcortex_QC_F_WGCNA_100)
enrich_df_M_100 <- GetEnrichrTable(Mcortex_QC_M_WGCNA_100)
# save to Excel
write_xlsx(enrich_df_F_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_enrichR_100.xlsx")
write_xlsx(enrich_df_M_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_enrichR_100.xlsx")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_F_clean_100 <- enrich_df_F_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean_100 <- enrich_df_M_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean_100 <- enrich_df_F_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean_100 <- enrich_df_M_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Glut_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
q <- ggplot(enrich_df_M_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Glut_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
print(q)
```
hdWGCNA analysis of PV-IN - sex specific
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_F_WGCNA <- SetDatExpr(
Mcortex_QC_F_WGCNA,
group_name = "GABAergic Pvalb",
group.by = "celltype1",
assay = "SCT",
layer = "data")
Mcortex_QC_M_WGCNA <- SetDatExpr(
Mcortex_QC_M_WGCNA,
group_name = "GABAergic Pvalb",
group.by = "celltype1",
assay = "SCT",
layer = "data")
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_F_WGCNA <- TestSoftPowers(
Mcortex_QC_F_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
Mcortex_QC_M_WGCNA <- TestSoftPowers(
Mcortex_QC_M_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list_F <- PlotSoftPowers(Mcortex_QC_F_WGCNA)
plot_list_M <- PlotSoftPowers(Mcortex_QC_M_WGCNA)
# assemble with patchwork
wrap_plots(plot_list_F, ncol=2)
wrap_plots(plot_list_M, ncol=2)
# Construct co-expression network
Mcortex_QC_F_WGCNA <- ConstructNetwork(
Mcortex_QC_F_WGCNA,
tom_name = 'PV-IN_F') # name of the topological overlap matrix written to disk
Mcortex_QC_M_WGCNA <- ConstructNetwork(
Mcortex_QC_M_WGCNA,
tom_name = 'PV-IN_M')
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_F_WGCNA <- ScaleData(Mcortex_QC_F_WGCNA, features=VariableFeatures(Mcortex_QC_F_WGCNA))
Mcortex_QC_M_WGCNA <- ScaleData(Mcortex_QC_M_WGCNA, features=VariableFeatures(Mcortex_QC_M_WGCNA))
Mcortex_QC_F_WGCNA <- ModuleEigengenes(
Mcortex_QC_F_WGCNA,
group.by.vars="sample_id")
Mcortex_QC_M_WGCNA <- ModuleEigengenes(
Mcortex_QC_M_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_F_WGCNA)
hMEs <- GetMEs(Mcortex_QC_M_WGCNA)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_F_WGCNA <- ModuleConnectivity(
Mcortex_QC_F_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Pvalb')
Mcortex_QC_M_WGCNA <- ModuleConnectivity(
Mcortex_QC_M_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Pvalb')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_F_WGCNA@misc$hdWGCNA$wgcna_params
Mcortex_QC_M_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_F_WGCNA <- ResetModuleNames(
Mcortex_QC_F_WGCNA,
new_name = "F_PV-IN-M")
Mcortex_QC_M_WGCNA <- ResetModuleNames(
Mcortex_QC_M_WGCNA,
new_name = "M_PV-IN-M")
# Change color of modules - Female
library(MetBrewer)
modules <- GetModules(Mcortex_QC_F_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_F_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_F <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_F_WGCNA <- ResetModuleColors(Mcortex_QC_F_WGCNA, new_colors_F)
# Change color of modules - Male
modules <- GetModules(Mcortex_QC_M_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_M_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_M <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_M_WGCNA <- ResetModuleColors(Mcortex_QC_M_WGCNA, new_colors_M)
# Plot dendogram
PlotDendrogram(Mcortex_QC_F_WGCNA, main='PV-IN_F')
PlotDendrogram(Mcortex_QC_M_WGCNA, main='PV-IN_M')
# plot genes ranked by kME for each module - did not run
PlotKMEs(Mcortex_QC_F_WGCNA, ncol=3)
PlotKMEs(Mcortex_QC_M_WGCNA, ncol=3)
# get the module assignment table:
modules_F <- GetModules(Mcortex_QC_F_WGCNA) %>% subset(module != 'grey')
modules_M <- GetModules(Mcortex_QC_M_WGCNA) %>% subset(module != 'grey')
# save to Excel
write_xlsx(modules_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_module_assignments.xlsx")
write_xlsx(modules_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_module_assignments.xlsx")
# get hub genes
hub_df_F <- GetHubGenes(Mcortex_QC_F_WGCNA, n_hubs = 100)
hub_df_M <- GetHubGenes(Mcortex_QC_M_WGCNA, n_hubs = 100)
# save to Excel
write_xlsx(hub_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_top100_hubgenes.xlsx")
write_xlsx(hub_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_top100_hubgenes.xlsx")
```
hdWGCN module trait correlation - PV-IN - sex specific
```{r}
# convert treatment to factor
Mcortex_QC_F_WGCNA$Tx <- as.factor(Mcortex_QC_F_WGCNA$Tx)
Mcortex_QC_M_WGCNA$Tx <- as.factor(Mcortex_QC_M_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c("Tx") # single trait is fine now!
Mcortex_QC_F_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_F_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
Mcortex_QC_M_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_M_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
# get the mt-correlation results
mt_cor_F <- GetModuleTraitCorrelation(Mcortex_QC_F_WGCNA)
names(mt_cor_F)
names(mt_cor_F$cor)
head(mt_cor_F$cor$`GABAergic Pvalb`[,1:5])
mt_cor_M <- GetModuleTraitCorrelation(Mcortex_QC_M_WGCNA)
names(mt_cor_M)
names(mt_cor_M$cor)
head(mt_cor_M$cor$`GABAergic Pvalb`[,1:5])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor_F$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor_F$cor[[celltype]]
pval_mat <- mt_cor_F$pval[[celltype]]
fdr_mat <- mt_cor_F$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# For male
wb <- createWorkbook()
for (celltype in names(mt_cor_M$cor)) {
cor_mat <- mt_cor_M$cor[[celltype]]
pval_mat <- mt_cor_M$pval[[celltype]]
fdr_mat <- mt_cor_M$fdr[[celltype]]
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_F_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
PlotModuleTraitCorrelation(
Mcortex_QC_M_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
```
hdWGCNA enrichment analysis- PV-IN - sex specific
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_F_WGCNA <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F <- GetEnrichrTable(Mcortex_QC_F_WGCNA)
enrich_df_M <- GetEnrichrTable(Mcortex_QC_M_WGCNA)
# save to Excel
write_xlsx(enrich_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_enrichR.xlsx")
write_xlsx(enrich_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_enrichR.xlsx")
# Define the order of modules you want (e.g., numeric order)
module_order_F <- c("F_PV-IN-M1", "F_PV-IN-M2", "F_PV-IN-M3", "F_PV-IN-M4",
"F_PV-IN-M5", "F_PV-IN-M6")
module_order_M <- c("M_PV-IN-M1", "M_PV-IN-M2", "M_PV-IN-M3", "M_PV-IN-M4",
"M_PV-IN-M5")
# Filter enrichment table for your database
db_to_plot <- "GO_Biological_Process_2023"
# Clean and organize enrichment table
enrich_df_F_clean <- enrich_df_F %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean <- enrich_df_M %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean <- enrich_df_F_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean <- enrich_df_M_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("PV-IN_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
q <- ggplot(enrich_df_M_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("PV-IN_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(q)
## Re-run enrichR for max 100 genes saved separately, graph not made)
# perform enrichment tests
Mcortex_QC_F_WGCNA_100 <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA_100 <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F_100 <- GetEnrichrTable(Mcortex_QC_F_WGCNA_100)
enrich_df_M_100 <- GetEnrichrTable(Mcortex_QC_M_WGCNA_100)
# save to Excel
write_xlsx(enrich_df_F_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_enrichR_100.xlsx")
write_xlsx(enrich_df_M_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_enrichR_100.xlsx")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_F_clean_100 <- enrich_df_F_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean_100 <- enrich_df_M_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean_100 <- enrich_df_F_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean_100 <- enrich_df_M_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("PV-IN_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
q <- ggplot(enrich_df_M_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("PV-IN_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
print(q)
```
hdWGCNA analysis of Sst-IN - sex specific
```{r}
# Set up the expression matrix for 3 cell types of interest
Mcortex_QC_F_WGCNA <- SetDatExpr(
Mcortex_QC_F_WGCNA,
group_name = "GABAergic Sst",
group.by = "celltype1",
assay = "SCT",
layer = "data")
Mcortex_QC_M_WGCNA <- SetDatExpr(
Mcortex_QC_M_WGCNA,
group_name = "GABAergic Sst",
group.by = "celltype1",
assay = "SCT",
layer = "data")
# Select soft-power threshold
# Test different soft powers:
Mcortex_QC_F_WGCNA <- TestSoftPowers(
Mcortex_QC_F_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
Mcortex_QC_M_WGCNA <- TestSoftPowers(
Mcortex_QC_M_WGCNA,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list_F <- PlotSoftPowers(Mcortex_QC_F_WGCNA)
plot_list_M <- PlotSoftPowers(Mcortex_QC_M_WGCNA)
# assemble with patchwork
wrap_plots(plot_list_F, ncol=2)
wrap_plots(plot_list_M, ncol=2)
# Construct co-expression network
Mcortex_QC_F_WGCNA <- ConstructNetwork(
Mcortex_QC_F_WGCNA,
tom_name = 'Sst-IN_F') # name of the topological overlap matrix written to disk
Mcortex_QC_M_WGCNA <- ConstructNetwork(
Mcortex_QC_M_WGCNA,
tom_name = 'Sst-IN_M')
# Compute harmonized module eingengenes
# need to re-run ScaleData because sample_id metadata added after SCTransform:
Mcortex_QC_F_WGCNA <- ScaleData(Mcortex_QC_F_WGCNA, features=VariableFeatures(Mcortex_QC_F_WGCNA))
Mcortex_QC_M_WGCNA <- ScaleData(Mcortex_QC_M_WGCNA, features=VariableFeatures(Mcortex_QC_M_WGCNA))
Mcortex_QC_F_WGCNA <- ModuleEigengenes(
Mcortex_QC_F_WGCNA,
group.by.vars="sample_id")
Mcortex_QC_M_WGCNA <- ModuleEigengenes(
Mcortex_QC_M_WGCNA,
group.by.vars="sample_id")
# Get module eigengenes (Harmonized by default)
hMEs <- GetMEs(Mcortex_QC_F_WGCNA)
hMEs <- GetMEs(Mcortex_QC_M_WGCNA)
# Compute eigengene-based module connectivity (kME)
Mcortex_QC_F_WGCNA <- ModuleConnectivity(
Mcortex_QC_F_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Sst')
Mcortex_QC_M_WGCNA <- ModuleConnectivity(
Mcortex_QC_M_WGCNA,
group.by = 'celltype1', group_name = 'GABAergic Sst')
# Look inside the hdWGCNA network parameters to check softpower beta
Mcortex_QC_F_WGCNA@misc$hdWGCNA$wgcna_params
Mcortex_QC_M_WGCNA@misc$hdWGCNA$wgcna_params
# rename the modules for convenience
Mcortex_QC_F_WGCNA <- ResetModuleNames(
Mcortex_QC_F_WGCNA,
new_name = "F_Sst-IN-M")
Mcortex_QC_M_WGCNA <- ResetModuleNames(
Mcortex_QC_M_WGCNA,
new_name = "M_Sst-IN-M")
# Change color of modules - Female
library(MetBrewer)
modules <- GetModules(Mcortex_QC_F_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_F_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_F <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_F_WGCNA <- ResetModuleColors(Mcortex_QC_F_WGCNA, new_colors_F)
# Change color of modules - Male
modules <- GetModules(Mcortex_QC_M_WGCNA)
mods <- unique(modules$module)
mod_colors_df <- dplyr::select(modules, c(module, color)) %>%
distinct %>% arrange(module)
rownames(mod_colors_df) <- mod_colors_df$module
mod_colors_df
mod_color_df <- GetModules(Mcortex_QC_M_WGCNA) %>%
dplyr::select(c(module, color)) %>%
distinct %>% arrange(module)
n_mods <- nrow(mod_color_df) - 1
new_colors_M <- paste0(met.brewer("Signac", n=n_mods))
Mcortex_QC_M_WGCNA <- ResetModuleColors(Mcortex_QC_M_WGCNA, new_colors_M)
# Plot dendogram
PlotDendrogram(Mcortex_QC_F_WGCNA, main='Sst-IN_F')
PlotDendrogram(Mcortex_QC_M_WGCNA, main='Sst-IN_M')
# plot genes ranked by kME for each module - did not run
PlotKMEs(Mcortex_QC_F_WGCNA, ncol=3)
PlotKMEs(Mcortex_QC_M_WGCNA, ncol=3)
# get the module assignment table:
modules_F <- GetModules(Mcortex_QC_F_WGCNA) %>% subset(module != 'grey')
modules_M <- GetModules(Mcortex_QC_M_WGCNA) %>% subset(module != 'grey')
# save to Excel
write_xlsx(modules_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_module_assignments.xlsx")
write_xlsx(modules_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_module_assignments.xlsx")
# get hub genes
hub_df_F <- GetHubGenes(Mcortex_QC_F_WGCNA, n_hubs = 100)
hub_df_M <- GetHubGenes(Mcortex_QC_M_WGCNA, n_hubs = 100)
# save to Excel
write_xlsx(hub_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_top100_hubgenes.xlsx")
write_xlsx(hub_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_top100_hubgenes.xlsx")
```
hdWGCN module trait correlation - Sst-IN - sex specific
```{r}
# convert treatment to factor
Mcortex_QC_F_WGCNA$Tx <- as.factor(Mcortex_QC_F_WGCNA$Tx)
Mcortex_QC_M_WGCNA$Tx <- as.factor(Mcortex_QC_M_WGCNA$Tx)
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Predicted_Sex is a factor with levels Female, Male, Unknown. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order? --> can't use it actually because only allows 2 variables. Need to subset the seurat object earlier on for just female and male, and then redo this.
## Warning in ModuleTraitCorrelation(Mcortex_QC_WGCNA, traits = cur_traits, :Trait Tx is a factor with levels Control, Risperidone. Levels will be converted to numeric IN THIS ORDER for the correlation, is this the expected order?
# list of traits to correlate
cur_traits <- c("Tx") # single trait is fine now!
Mcortex_QC_F_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_F_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
Mcortex_QC_M_WGCNA <- ModuleTraitCorrelation_mod(
Mcortex_QC_M_WGCNA,
traits = cur_traits,
group.by = "celltype1"
)
# get the mt-correlation results
mt_cor_F <- GetModuleTraitCorrelation(Mcortex_QC_F_WGCNA)
names(mt_cor_F)
names(mt_cor_F$cor)
head(mt_cor_F$cor$`GABAergic Sst`[,1:5])
mt_cor_M <- GetModuleTraitCorrelation(Mcortex_QC_M_WGCNA)
names(mt_cor_M)
names(mt_cor_M$cor)
head(mt_cor_M$cor$`GABAergic Sst`[,1:5])
# Create a new workbook
wb <- createWorkbook()
# Loop through each cell type / group
for (celltype in names(mt_cor_F$cor)) {
# Extract correlation, p-value, and FDR matrices
cor_mat <- mt_cor_F$cor[[celltype]]
pval_mat <- mt_cor_F$pval[[celltype]]
fdr_mat <- mt_cor_F$fdr[[celltype]]
# Convert to data.frames for writing
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
# Add worksheets
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
# Write data to the workbook
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
# Save the workbook
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# For male
wb <- createWorkbook()
for (celltype in names(mt_cor_M$cor)) {
cor_mat <- mt_cor_M$cor[[celltype]]
pval_mat <- mt_cor_M$pval[[celltype]]
fdr_mat <- mt_cor_M$fdr[[celltype]]
cor_df <- as.data.frame(cor_mat)
pval_df <- as.data.frame(pval_mat)
fdr_df <- as.data.frame(fdr_mat)
addWorksheet(wb, paste0(celltype, "_cor"))
addWorksheet(wb, paste0(celltype, "_pval"))
addWorksheet(wb, paste0(celltype, "_fdr"))
writeData(wb, paste0(celltype, "_cor"), cor_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_pval"), pval_df, rowNames = TRUE)
writeData(wb, paste0(celltype, "_fdr"), fdr_df, rowNames = TRUE)
}
output_path <- "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_moduletraitcor.xlsx"
saveWorkbook(wb, output_path, overwrite = TRUE)
# Plot heatmap for the correlation result
PlotModuleTraitCorrelation(
Mcortex_QC_F_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
PlotModuleTraitCorrelation(
Mcortex_QC_M_WGCNA,
label = 'fdr',
label_symbol = 'stars',
text_size = 3,
text_digits = 2,
text_color = 'black',
high_color = '#db2763',
mid_color = 'white',
low_color = '#3772ff',
plot_max = 0.6,
combine=TRUE
)
```
hdWGCNA enrichment analysis- Sst-IN - sex specific
```{r}
# define the enrichr databases to test
dbs <- c('GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023')
# perform enrichment tests
Mcortex_QC_F_WGCNA <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = Inf # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F <- GetEnrichrTable(Mcortex_QC_F_WGCNA)
enrich_df_M <- GetEnrichrTable(Mcortex_QC_M_WGCNA)
# save to Excel
write_xlsx(enrich_df_F, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_enrichR.xlsx")
write_xlsx(enrich_df_M, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_enrichR.xlsx")
# Define the order of modules you want (e.g., numeric order)
module_order_F <- c("F_Sst-IN-M1", "F_Sst-IN-M2", "F_Sst-IN-M3", "F_Sst-IN-M4",
"F_Sst-IN-M5", "F_Sst-IN-M6", "F_Sst-IN-M7")
module_order_M <- c("M_Sst-IN-M1", "M_Sst-IN-M2", "M_Sst-IN-M3", "M_Sst-IN-M4",
"M_Sst-IN-M5", "M_Sst-IN-M6", "M_Sst-IN-M7")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_F_clean <- enrich_df_F %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean <- enrich_df_M %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean <- enrich_df_F_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean <- enrich_df_M_clean %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Sst-IN_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
q <- ggplot(enrich_df_M_clean,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Sst-IN_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(q)
## Re-run enrichR for max 100 genes saved separately, graph not made)
# perform enrichment tests
Mcortex_QC_F_WGCNA_100 <- RunEnrichr(
Mcortex_QC_F_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
Mcortex_QC_M_WGCNA_100 <- RunEnrichr(
Mcortex_QC_M_WGCNA,
dbs=dbs,
max_genes = 100 # use max_genes = Inf to choose all genes
)
# retrieve the output table
enrich_df_F_100 <- GetEnrichrTable(Mcortex_QC_F_WGCNA_100)
enrich_df_M_100 <- GetEnrichrTable(Mcortex_QC_M_WGCNA_100)
# save to Excel
write_xlsx(enrich_df_F_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_enrichR_100.xlsx")
write_xlsx(enrich_df_M_100, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_enrichR_100.xlsx")
# Filter enrichment table for your database
db_to_plot <- "GO_Molecular_Function_2023"
# Clean and organize enrichment table
enrich_df_F_clean_100 <- enrich_df_F_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_F),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
enrich_df_M_clean_100 <- enrich_df_M_100 %>%
dplyr::filter(db == db_to_plot) %>%
group_by(module) %>%
arrange(Adjusted.P.value, desc(Combined.Score), .by_group = TRUE) %>%
slice_head(n = 2) %>%
ungroup() %>%
mutate(
# Remove "(GO:xxxx)" suffix
Term_clean = str_remove(Term, "\\s*\$GO:\\d+\$"),
# Wrap long GO terms to two lines if needed
Term_clean = str_wrap(Term_clean, width = 50),
# Keep numeric ordering of modules
module = factor(module, levels = module_order_M),
# log10 transform Combined.Score for better scaling
log10_combined = log10(Combined.Score)
)
# Reorder terms so they're grouped by module, with best p-values at top
enrich_df_F_clean_100 <- enrich_df_F_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
enrich_df_M_clean_100 <- enrich_df_M_clean_100 %>%
mutate(Term_clean = fct_reorder(Term_clean, as.numeric(module), .desc = TRUE))
# Plot
p <- ggplot(enrich_df_F_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Sst-IN_F - ", db_to_plot),
x = "Module",
y = "GO Term"
)
q <- ggplot(enrich_df_M_clean_100,
aes(x = module, y = Term_clean,
size = log10_combined)) +
# Main points: filled by -log10(FDR)
geom_point(
aes(fill = -log10(Adjusted.P.value),
shape = Adjusted.P.value <= 0.05),
color = "black", alpha = 0.9, stroke = 0.7
) +
# Set shapes manually: filled circle for sig, hollow for ns
scale_shape_manual(
values = c("TRUE" = 21, "FALSE" = 1),
labels = c("TRUE" = "FDR ≤ 0.05", "FALSE" = "FDR > 0.05"),
name = "Significance"
) +
scale_fill_stepsn(
colors = rev(viridis::magma(256)),
name = expression(-log[10]("FDR"))
) +
scale_size_continuous(
name = expression(log[10]("Enrichment"))
) +
theme_minimal(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
axis.text.y = element_text(size = 11),
legend.position = "right"
) +
labs(
title = paste0("Sst-IN_M - ", db_to_plot),
x = "Module",
y = "GO Term"
)
print(p)
print(q)
```
Finding corresponding modules between male and female
```{r}
library(readxl)
library(stringr)
# Function to read one Excel file and organize modules
read_module_file <- function(file_path) {
sheets <- excel_sheets(file_path)
# Extract cell type from file name (assumes format like "Glut_M_top100_hubgenes.xlsx")
cell_type <- str_extract(basename(file_path), "Glut|PV-IN|Sst-IN")
module_list <- list()
for (sh in sheets) {
df <- read_excel(file_path, sheet = sh)
# Auto-detect gene name column
gene_col <- grep("gene|symbol|Gene|Symbol", names(df), value = TRUE)[1]
genes <- unique(df[[gene_col]])
# Extract module id (e.g., M1, M2, ...)
module_id <- str_extract(sh, "M\\d+")
module_list[[module_id]] <- genes
}
return(module_list)
}
# ---- Read all files ----
# Modify paths if needed
male_files <- c("~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_M_top100_hubgenes_copy.xlsx",
"~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_M_top100_hubgenes_copy.xlsx",
"~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_M_top100_hubgenes_copy.xlsx")
female_files <- c("~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Glut_F_top100_hubgenes_copy.xlsx",
"~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/PV-IN_F_top100_hubgenes_copy.xlsx",
"~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Sst-IN_F_top100_hubgenes_copy.xlsx")
male_modules <- list()
female_modules <- list()
for (file in male_files) {
cell_type <- str_extract(file, "Glut|PV-IN|Sst-IN")
male_modules[[cell_type]] <- read_module_file(file)
}
for (file in female_files) {
cell_type <- str_extract(file, "Glut|PV-IN|Sst-IN")
female_modules[[cell_type]] <- read_module_file(file)
}
# Check the vectors
names(male_modules) # Should show: "Glut" "PV" "Sst"
names(male_modules$Glut) # Should show: "M1" "M2" "M3" ...
male_modules$Glut$M1[1:10] # First 10 hub genes
## Extract all genes from Seurat object, restrict to genes expressed in at least 5% of cells (to be set as background genes for Fisher's test)
Mcortex_QC <- readRDS ("~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/RDS_files/Mcortex_QC.rds")
counts <- GetAssayData(Mcortex_QC, slot = "counts") # raw counts
pct_cells <- Matrix::rowSums(counts > 0) / ncol(counts)
background_genes <- names(pct_cells[pct_cells > 0.05])
length(background_genes) #7404 (without filtering at least 5%, it was 24253 genes total)
## Compare Male vs. Female modules by hubgenes
library(pheatmap)
library(dplyr)
# Function to compute Jaccard similarity
jaccard <- function(x, y) {
length(intersect(x, y)) / length(union(x, y))
}
# Function to compute module comparison for a single cell type
compare_modules <- function(cell_type, male_modules, female_modules, background_genes) {
male_mods <- male_modules[[cell_type]]
female_mods <- female_modules[[cell_type]]
male_names <- names(male_mods)
female_names <- names(female_mods)
# Initialize matrices
overlap_mat <- matrix(0, nrow = length(male_names), ncol = length(female_names),
dimnames = list(male_names, female_names))
jaccard_mat <- overlap_mat
fisher_mat <- overlap_mat
# Compute metrics
for (m in male_names) {
for (f in female_names) {
genes_m <- male_mods[[m]]
genes_f <- female_mods[[f]]
# Overlap
overlap_mat[m, f] <- length(intersect(genes_m, genes_f))
# Jaccard
jaccard_mat[m, f] <- jaccard(genes_m, genes_f)
# Fisher exact test
a <- length(intersect(genes_m, genes_f))
b <- length(setdiff(genes_m, genes_f))
c <- length(setdiff(genes_f, genes_m))
d <- length(background_genes) - (a + b + c)
fisher_mat[m, f] <- fisher.test(matrix(c(a,b,c,d), nrow = 2))$p.value
}
}
# Identify best female match for each male module (by Jaccard)
best_matches <- apply(jaccard_mat, 1, function(x) {
female_names[which.max(x)]
})
# Return as list
return(list(
overlap = overlap_mat,
jaccard = jaccard_mat,
fisher_p = fisher_mat,
best_match = best_matches
))
}
# ---- Run comparison for all cell types ----
cell_types <- names(male_modules)
results <- list()
for (ct in cell_types) {
results[[ct]] <- compare_modules(ct, male_modules, female_modules, background_genes)
}
# View Glut results
results$Glut$overlap # Overlap counts
results$Glut$jaccard # Jaccard similarity
results$Glut$fisher_p # Fisher p-values
results$Glut$best_match # Best female module per male module
## Save the best match table
library(openxlsx)
wb <- createWorkbook()
for (ct in cell_types) {
df <- data.frame(
Male_Module = names(results[[ct]]$best_match),
Best_Female_Module = results[[ct]]$best_match,
stringsAsFactors = FALSE
)
addWorksheet(wb, ct)
writeData(wb, ct, df)
}
saveWorkbook(wb, "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Male_vs_Female_Module_BestMatches.xlsx", overwrite = TRUE)
## Plot Heatmap
library(pheatmap)
# Function to plot enhanced heatmap for one cell type
plot_module_heatmap <- function(cell_type, results, fisher_alpha = 0.05) {
# Extract matrices
jaccard_mat <- results[[cell_type]]$jaccard
overlap_mat <- results[[cell_type]]$overlap
fisher_mat <- results[[cell_type]]$fisher_p
# Create matrix of labels: overlap counts + optional significance
label_mat <- matrix("", nrow = nrow(overlap_mat), ncol = ncol(overlap_mat),
dimnames = dimnames(overlap_mat))
for (i in 1:nrow(overlap_mat)) {
for (j in 1:ncol(overlap_mat)) {
sig_star <- ifelse(fisher_mat[i,j] < fisher_alpha, "*", "")
label_mat[i,j] <- paste0(overlap_mat[i,j], sig_star)
}
}
# Plot heatmap
pheatmap(jaccard_mat,
display_numbers = label_mat, # overlay counts + significance
number_color = "black",
main = paste0(cell_type, " Male vs Female Module Similarity"),
cluster_rows = TRUE,
cluster_cols = TRUE,
fontsize_number = 10,
color = colorRampPalette(c("white", "steelblue"))(50),
border_color = "grey60",
fontsize = 12)
}
# ---- Plot for all cell types ---- (Heatmap color is Jaccard similarity in 0-1 scddale, numbers inside cells means overlap count, * indicates Fisher p-value <0.05, meaning significant overlap; Rows = male modules, Columns = female modules)
for (ct in names(results)) {
plot_module_heatmap(ct, results)
}
## Generate list of overlapping genes between male vs. female modules
get_common_genes <- function(cell_type, male_modules, female_modules) {
male_mods <- male_modules[[cell_type]]
female_mods <- female_modules[[cell_type]]
common_genes <- list()
for (m in names(male_mods)) {
for (f in names(female_mods)) {
common_genes[[paste0(m, "_vs_", f)]] <- intersect(male_mods[[m]], female_mods[[f]])
}
}
return(common_genes)
}
# List the common genes
common_genes_glut <- get_common_genes("Glut", male_modules, female_modules)
common_genes_PV <- get_common_genes("PV-IN", male_modules, female_modules)
common_genes_Sst <- get_common_genes("Sst-IN", male_modules, female_modules)
# Generate excel file with common gene list
library(openxlsx)
write_common_genes_excel <- function(male_modules, female_modules, file = "~/Hahn_Winter/HyowonChoi/RStudioServer/Risperidone_Proj/R_analysis/hdWGCNA/Male_vs_Female_Overlap_Genes.xlsx") {
wb <- createWorkbook()
for (ct in names(male_modules)) {
male_mods <- male_modules[[ct]]
female_mods <- female_modules[[ct]]
df_list <- list()
for (m in names(male_mods)) {
for (f in names(female_mods)) {
overlap <- intersect(male_mods[[m]], female_mods[[f]])
df_list[[paste0(m, "_vs_", f)]] <- paste(overlap, collapse = ", ")
}
}
df <- data.frame(
Module_Comparison = names(df_list),
Overlapping_Genes = unlist(df_list),
stringsAsFactors = FALSE
)
addWorksheet(wb, ct)
writeData(wb, ct, df)
}
saveWorkbook(wb, file, overwrite = TRUE)
}
write_common_genes_excel(male_modules, female_modules)
```
Checking for module redundancy/module similarity within sex and cell type
```{r}
# Compare modules within a single module list (e.g., male Glut)
compare_within_modules <- function(module_list, background_genes) {
mod_names <- names(module_list)
# Initialize matrices
overlap_mat <- matrix(0, nrow = length(mod_names), ncol = length(mod_names),
dimnames = list(mod_names, mod_names))
jaccard_mat <- overlap_mat
fisher_mat <- overlap_mat
for (m1 in mod_names) {
for (m2 in mod_names) {
genes1 <- module_list[[m1]]
genes2 <- module_list[[m2]]
# Overlap
a <- length(intersect(genes1, genes2))
b <- length(setdiff(genes1, genes2))
c <- length(setdiff(genes2, genes1))
d <- length(background_genes) - (a + b + c)
overlap_mat[m1, m2] <- a
jaccard_mat[m1, m2] <- a / length(union(genes1, genes2))
# Fisher exact test
fisher_mat[m1, m2] <- fisher.test(matrix(c(a, b, c, d), nrow = 2))$p.value
}
}
return(list(
overlap = overlap_mat,
jaccard = jaccard_mat,
fisher_p = fisher_mat
))
}
# Male Glut
male_glut_results <- compare_within_modules(male_modules$Glut, background_genes)
male_glut_results$overlap
male_glut_results$jaccard
male_glut_results$fisher_p
# Female Glut
female_glut_results <- compare_within_modules(female_modules$Glut, background_genes)
female_glut_results$overlap
female_glut_results$jaccard
female_glut_results$fisher_p
# Male PV-IN
male_PV_results <- compare_within_modules(male_modules[["PV-IN"]], background_genes)
male_PV_results$overlap
male_PV_results$jaccard
male_PV_results$fisher_p
# Female PV-IN
female_PV_results <- compare_within_modules(female_modules[["PV-IN"]], background_genes)
female_PV_results$overlap
female_PV_results$jaccard
female_PV_results$fisher_p
# Male Sst-IN
male_Sst_results <- compare_within_modules(male_modules[["Sst-IN"]], background_genes)
male_Sst_results$overlap
male_Sst_results$jaccard
male_Sst_results$fisher_p
# Female Sst-IN
female_Sst_results <- compare_within_modules(female_modules[["Sst-IN"]], background_genes)
female_Sst_results$overlap
female_Sst_results$jaccard
female_Sst_results$fisher_p
# Enhanced heatmap with overlap text
plot_within_heatmap <- function(results, cell_type, sex, fisher_alpha = 0.05) {
jaccard_mat <- results$jaccard
overlap_mat <- results$overlap
fisher_mat <- results$fisher_p
# Add overlap + significance
label_mat <- matrix("", nrow = nrow(overlap_mat), ncol = ncol(overlap_mat),
dimnames = dimnames(overlap_mat)

Antipsychotic_github.Rmd at commit b768f3a, no license · at the source

Overview

Authors: Abneil D Alicea-Pauneto^1,2, Hyowon Choi^1,2,3, Wenyu Zhang^1,2, Xinjian Li⁴, Laurence A Busque^1,2, Mingxuan Li⁵, Andrew Wu^1,2, Evan Zhang⁶, Adam D Marc^1,2, Jonathan Preall⁶, Kuan Hong Wang⁷, Robert E Featherstone⁸, Steven J Siegel⁸, Chang-Gyu Hahn^1,2,3, Karin E Borgmann-Winter^1,2,3

ORCID iDs: Kuan Hong Wang

Department of Neuroscience, Thomas Jefferson University, Philadelphia, PA USA
Vickie & Jack Farber Institute for Neuroscience, Thomas Jefferson University, Philadelphia, PA USA
Department of Psychiatry and Human Behavior, Sidney Kimmel Medical College, Thomas Jefferson University, Philadelphia, PA USA
Department of Neurology of the Second Affiliated Hospital, Interdisciplinary Institute of Neuroscience and Technology, Zhejiang University School of Medicine, Hangzhou, China
Department of Computer Science, Hertford College, University of Oxford, Oxford, UK
Cold Spring Harbor Laboratory Cancer Center, Cold Spring Harbor, NY USA
Department of Neuroscience, University of Rochester Medical Center, Rochester, NY USA
Department of Psychiatry and the Behavioral Sciences, University of Southern California, Los Angeles, CA USA

Institutions: Thomas Jefferson University (United States); Second Affiliated Hospital of Zhejiang University (China); Zhejiang University (China); University of Oxford (United Kingdom); Cancer Center at Cold Spring Harbor Laboratory; University of Rochester Medicine (United States); University of Rochester (United States); University of Southern California (United States)

Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology, volume 51, issue 10, pages 1844-1854

Dates: received 23 January 2026; accepted 23 June 2026; published online 8 July 2026; in print September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41386-026-02483-2 · PMID 42420423 · PMCID PMC13486667 · OpenAlex W7167711737

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), mouse (organism), schizophrenia / psychosis (population)

Methods: Statistics, Smoothing, state filtering, decompositions, fMRI & imaging, Single-unit activity, calcium imaging, Spectral & time-frequency

Keywords: Synaptic development, Cellular neuroscience, Genetics of the nervous system

MeSH: Antipsychotic Agents*, Cerebral Cortex*, Risperidone*, Animals, Anxiety, Female, Male, Mice, Mice, Inbred C57BL, Neurons (* major topic)

Topic: Neuroscience and Neuropharmacology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex

Funding: U.S. Department of Health & Human Services | NIH | National Institute of Mental Health (NIMH) (RO1MH138995, RO1MH116463); NIGMS NIH HHS (T32 GM008562); NIMH NIH HHS (R01 MH138995, R01 MH116463); U.S. Department of Health & Human Services | NIH | Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD) (T32GM008562)

Citations: not cited yet (Europe PMC); 56 references in the paper

Abstract

Since the introduction of second-generation antipsychotics, antipsychotics have been increasingly prescribed for children and adolescents, raising concerns about their long-term impact on neurodevelopment. Antipsychotics block dopaminergic and serotonergic receptors, potentially disrupting the maturation of neurocognitive processes, which is a public health concern. Previous studies have reported that adolescent antipsychotic treatment can cause persistent neurocognitive dysfunction in rodents, yet the neurobiological underpinnings remain unknown. To address this, we administered risperidone, a commonly used antipsychotic, to C57BL/6 mice during adolescence (3 to 6 weeks of age) and examined behavioral and neurobiological outcomes nine weeks post-treatment. Risperidone-treated mice exhibited subtle deficits in behavioral correlates of anxiety-like behavior. In vivo, two-photon calcium imaging of cortical neurons revealed a remarkable increase in the amplitude of calcium events with subtle sex-specific changes in the frequency, consistent with increased neuronal excitability. Single-nucleus RNA-sequencing (snRNA-seq) analyses showed widespread reductions in transcripts for voltage-sensitive and inwardly rectifying potassium channels in both pyramidal neurons and interneurons. Additionally, both cell types exhibited reduced Grin2a and Grin2b, as well as scaffolding proteins, indicative of weakened synaptic connectivity between excitatory and inhibitory neurons. Interestingly, we observed sex-dependent differences in the directionality of correlation between certain gene co-expression modules and risperidone treatment. Our results suggest that adolescent risperidone treatment induces lasting transcriptomic and functional changes associated with altered excitatory-inhibitory neuronal interactions that may underline cognitive and behavioral dysregulations.

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

Repository

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

Hahn-BorgmannWinter/Adolescent_antipsychotics

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: b768f3a632e45fccec74c82fd5084bdb451ccd50, 26 January 2026

Languages: R (1)

Size: 2 files, 1 script

Software Heritage: not archived

Found in: “Data availability”

Holds: README, 1 notebook

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: clusterProfiler (1 file), cowplot (1 file), DESeq2 (1 file), ggplot2 (1 file), igraph (1 file), patchwork (1 file), pheatmap (1 file), reshape2 (1 file), Seurat (1 file), tidyverse (1 file), WGCNA (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

2 files

Antipsychotic_github.Rmd, R, 5,888 lines, 4 matches
README.md, Text, 2 lines

At the source: github.com/Hahn-BorgmannWinter/Adolescent_antipsychotics

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

geo:GSE316091, at NCBI GEO; found in “Data availability”

Data availability

The sequence and raw count matrices generated and analyzed for snRNA-seq experiments are deposited to the Gene Expression Omnibus, GSE316091 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE316091) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE316091). All code used to generate results for snRNA-seq analysis is available at: https://github.com/Hahn-BorgmannWinter/Adolescent_antipsychotics.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 15 authors, 3 keywords, 10 MeSH terms, 4 funders, 56 references.

Cite

This paper

Alicea-Pauneto, A. D., Choi, H., Zhang, W., Li, X., Busque, L. A., Li, M., Wu, A., Zhang, E., Marc, A. D., Preall, J., Wang, K. H., Featherstone, R. E., Siegel, S. J., Hahn, C.-G., & Borgmann-Winter, K. E. (2026). Long-term effects of adolescent risperidone treatment on the mouse cortex. Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology, 51(10), 1844-1854. https://doi.org/10.1038/s41386-026-02483-2

BibTeX

@article{aliceapauneto2026long, author = {Alicea-Pauneto, Abneil D and Choi, Hyowon and Zhang, Wenyu and Li, Xinjian and Busque, Laurence A and Li, Mingxuan and Wu, Andrew and Zhang, Evan and Marc, Adam D and Preall, Jonathan and Wang, Kuan Hong and Featherstone, Robert E and Siegel, Steven J and Hahn, Chang-Gyu and Borgmann-Winter, Karin E}, title = {{Long-term effects of adolescent risperidone treatment on the mouse cortex}}, journal = {Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology}, year = {2026}, month = jul, volume = {51}, number = {10}, pages = {1844--1854}, publisher = {Nature Publishing Group}, issn = {0893-133X}, doi = {10.1038/s41386-026-02483-2}, url = {https://doi.org/10.1038/s41386-026-02483-2}, pmid = {42420423}, pmcid = {PMC13486667} }

RIS

TY - JOUR AU - Alicea-Pauneto, Abneil D AU - Choi, Hyowon AU - Zhang, Wenyu AU - Li, Xinjian AU - Busque, Laurence A AU - Li, Mingxuan AU - Wu, Andrew AU - Zhang, Evan AU - Marc, Adam D AU - Preall, Jonathan AU - Wang, Kuan Hong AU - Featherstone, Robert E AU - Siegel, Steven J AU - Hahn, Chang-Gyu AU - Borgmann-Winter, Karin E TI - Long-term effects of adolescent risperidone treatment on the mouse cortex T2 - Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology J2 - Neuropsychopharmacology PY - 2026 DA - 2026/07/08 VL - 51 IS - 10 SP - 1844 EP - 1854 SN - 0893-133X PB - Nature Publishing Group DO - 10.1038/s41386-026-02483-2 UR - https://doi.org/10.1038/s41386-026-02483-2 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41386-026-02483-2", "type": "article-journal", "title": "Long-term effects of adolescent risperidone treatment on the mouse cortex", "container-title": "Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology", "author": [ { "family": "Alicea-Pauneto", "given": "Abneil D" }, { "family": "Choi", "given": "Hyowon" }, { "family": "Zhang", "given": "Wenyu" }, { "family": "Li", "given": "Xinjian" }, { "family": "Busque", "given": "Laurence A" }, { "family": "Li", "given": "Mingxuan" }, { "family": "Wu", "given": "Andrew" }, { "family": "Zhang", "given": "Evan" }, { "family": "Marc", "given": "Adam D" }, { "family": "Preall", "given": "Jonathan" }, { "family": "Wang", "given": "Kuan Hong" }, { "family": "Featherstone", "given": "Robert E" }, { "family": "Siegel", "given": "Steven J" }, { "family": "Hahn", "given": "Chang-Gyu" }, { "family": "Borgmann-Winter", "given": "Karin E" } ], "container-title-short": "Neuropsychopharmacology", "volume": "51", "issue": "10", "page": "1844-1854", "DOI": "10.1038/s41386-026-02483-2", "PMID": "42420423", "PMCID": "PMC13486667", "ISSN": "0893-133X", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41386-026-02483-2", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 8 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

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

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