Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.

Code ↔ Paper

13 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 13 matches

[1] § Methods › Bioinformatics and data analysis › SHARE-seq data processing and analysis ↔ Astro_GR_KO_snRNA-final.ipynb, lines 66–74 · score 0.90 · FindVariableFeatures, NormalizeData, ScaleData, RunPCA, RunUMAP, doublets
[2] § Astrocyte GR promotes neural circuit maturation ↔ Astro_GR_KO_snRNA-final.ipynb, lines 794–811 · score 0.81 · L5NP, L5PT, L6CT, L6IT, L5IT, excitatory neuron
[3] § Atlas of experience-dependent V1 development ↔ Astro_GR_KO_snRNA-final.ipynb, lines 399–454 · score 0.79 · L5NP, L5PT, L6CT, L5IT, CR, Endo
[4] § Methods › Bioinformatics and data analysis › Human brain multiome analysis ↔ Human_brain_development_multiome_analysis-final.ipynb, lines 432–497 · score 0.79 · L4IT, Human brain, L6IT, L5IT, multiome, L2
[5] § Methods › Bioinformatics and data analysis › Mouse developmental astrocyte snRNA-seq analysis ↔ Astrocyte_mouse_atlas_snRNAseq-final.ipynb, lines 90–97 · score 0.75 · AddModuleScore, GR repressed, GR activated, Seurat, Mouse, astrocyte
[6] § Atlas of experience-dependent V1 development ↔ Human_brain_development_multiome_analysis-final.ipynb, lines 432–497 · score 0.74 · L5PT, L6CT, L5IT, vascular, Fosl2, Endo
[7] § Methods › Bioinformatics and data analysis › CUT&RUN data processing and analysis ↔ GR_NFIA_CUTnRUN_analysis-final.ipynb, lines 221–226 · score 0.73 · ChIP, BEDtools intersect, GR CUT, NFIA, edgeR, overlapped
[8] § Methods › Bioinformatics and data analysis › Human brain multiome analysis ↔ Human_brain_development_multiome_analysis-final.ipynb, lines 1142–1166 · score 0.64 · FindMotifs, enriched motifs, human, multiome, clustered, brain
[9] § Methods › Bioinformatics and data analysis › SHARE-seq data processing and analysis ↔ Human_brain_development_multiome_analysis-final.ipynb, lines 365–381 · score 0.57 · FindMarkers, chromVAR score, MAST, FC, row, motifs
[10] § Methods › Bioinformatics and data analysis › Bulk nuclear RNA-seq data processing and analysis ↔ Share_seqV2_example.sh, lines 851–903 · score 0.56 · featureCounts, RNA seq, exons, log, mm10
[11] § Methods › Bioinformatics and data analysis › MERFISH data processing and analysis ↔ Astro_GR_KO_snRNA-final.ipynb, lines 643–696 · score 0.54 · glmLRT, edgeR, subset, matrix, DR, NR
[12] § GR regulates postnatal astrocyte maturation ↔ Human_brain_development_multiome_analysis-final.ipynb, lines 562–614 · score 0.52 · chromVAR score, CTX, infant, Adoles, Tri, age
[13] § Methods › Bioinformatics and data analysis › CUT&RUN data processing and analysis ↔ Share_seqV2_example.sh, lines 605–655 · score 0.51 · mapping quality, samtools, duplicate, BAM, CUT, filtering

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

Jupyter notebook · 2,255 lines · 102 KB · no license · 5 matches

# %% [markdown]
# # Load libraries
# %%
# Load in libraries
library(Seurat)
library(ggplot2)
library(data.table)
library(stringr)
library(Signac)
library(tradeSeq)
library(BSgenome.Hsapiens.UCSC.hg38)
library(EnsDb.Hsapiens.v86)
library(monaLisa)
library(edgeR)
library(DESeq2)
library(RColorBrewer)
library(pheatmap)
library(paletteer)
library(JASPAR2020)
library(TFBSTools)
library(BiocParallel)
library(paletteer)
library(ghibli)
library(Rmisc)
library(ggbreak)
library(colorspace)
library(pheatmap)
library(RColorBrewer)
library(viridis)
library(ComplexHeatmap)
library(circlize)
library(clusterProfiler)
library(ChIPseeker)
# %%
# Create 'regionize' function for conveniency
regionize <- function(x) {
y <- gsub(',', '', gsub(':', '-', x))
paste0(unlist(strsplit(y,'-'))[1],'-',
as.numeric(unlist(strsplit(y,'-'))[2]) - 5000,'-',
as.numeric(unlist(strsplit(y,'-'))[3]) + 5000)
}
# %%
#library(future)
#plan("multicore", workers = 2)
#options(future.globals.maxSize = 100 * 1024 ^ 3) # for 180 Gb RAM
# %% [markdown]
# # Initial processing
# %%
# Load in data
human_multiome <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/snMultiome_atlas_Seurat_object.rds')
# %%
DefaultAssay(human_multiome) <- 'ATAC'
# %%
frag_list <- list('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-18-PFC_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-19-V1_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-1-PFC_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-20-CTX_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-32-V1_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-37-CTX_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-41-PFC-2_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-41-V1_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-43-PFC_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-45-CTX_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/ARKFrozen-8-V1_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/GW27-2-7-18-PFC_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/HDBR-14584-CTX_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/HDBR-14831-CTX_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/HDBR-14834-TC_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/HDBR-15020-FB_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-1325-BA10-2_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-1325-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-1671-BA10-2_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-1671-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4267-BA10-2_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4341-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4341-BA9_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4373-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4373-BA9-3_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4392-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4392-BA9_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4458-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-4458-BA9-3_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5162-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5162-BA9_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5376-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5376-BA9_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5554-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5554-BA9_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-5900-BA17_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-M1154-BA10-2_atac_fragments.tsv.gz',
'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/NIH-M2837-BA10-2_atac_fragments.tsv.gz')
# %%
frags <- Fragments(human_multiome) # get list of fragment objects
Fragments(human_multiome) <- NULL # remove fragment information from assay
frags[[1]] <- UpdatePath(frags[[1]], new.path = frag_list[[1]])
frags[[2]] <- UpdatePath(frags[[2]], new.path = frag_list[[2]])
frags[[3]] <- UpdatePath(frags[[3]], new.path = frag_list[[3]])
frags[[4]] <- UpdatePath(frags[[4]], new.path = frag_list[[4]])
frags[[5]] <- UpdatePath(frags[[5]], new.path = frag_list[[5]])
frags[[6]] <- UpdatePath(frags[[6]], new.path = frag_list[[6]])
frags[[7]] <- UpdatePath(frags[[7]], new.path = frag_list[[7]])
frags[[8]] <- UpdatePath(frags[[8]], new.path = frag_list[[8]])
frags[[9]] <- UpdatePath(frags[[9]], new.path = frag_list[[9]])
frags[[10]] <- UpdatePath(frags[[10]], new.path = frag_list[[10]])
frags[[11]] <- UpdatePath(frags[[11]], new.path = frag_list[[11]])
frags[[12]] <- UpdatePath(frags[[12]], new.path = frag_list[[12]])
frags[[13]] <- UpdatePath(frags[[13]], new.path = frag_list[[13]])
frags[[14]] <- UpdatePath(frags[[14]], new.path = frag_list[[14]])
frags[[15]] <- UpdatePath(frags[[15]], new.path = frag_list[[15]])
frags[[16]] <- UpdatePath(frags[[16]], new.path = frag_list[[16]])
frags[[17]] <- UpdatePath(frags[[17]], new.path = frag_list[[17]])
frags[[18]] <- UpdatePath(frags[[18]], new.path = frag_list[[18]])
frags[[19]] <- UpdatePath(frags[[19]], new.path = frag_list[[19]])
frags[[20]] <- UpdatePath(frags[[20]], new.path = frag_list[[20]])
frags[[21]] <- UpdatePath(frags[[21]], new.path = frag_list[[21]])
frags[[22]] <- UpdatePath(frags[[22]], new.path = frag_list[[22]])
frags[[23]] <- UpdatePath(frags[[23]], new.path = frag_list[[23]])
frags[[24]] <- UpdatePath(frags[[24]], new.path = frag_list[[24]])
frags[[25]] <- UpdatePath(frags[[25]], new.path = frag_list[[25]])
frags[[26]] <- UpdatePath(frags[[26]], new.path = frag_list[[26]])
frags[[27]] <- UpdatePath(frags[[27]], new.path = frag_list[[27]])
frags[[28]] <- UpdatePath(frags[[28]], new.path = frag_list[[28]])
frags[[29]] <- UpdatePath(frags[[29]], new.path = frag_list[[29]])
frags[[30]] <- UpdatePath(frags[[30]], new.path = frag_list[[30]])
frags[[31]] <- UpdatePath(frags[[31]], new.path = frag_list[[31]])
frags[[32]] <- UpdatePath(frags[[32]], new.path = frag_list[[32]])
frags[[33]] <- UpdatePath(frags[[33]], new.path = frag_list[[33]])
frags[[34]] <- UpdatePath(frags[[34]], new.path = frag_list[[34]])
frags[[35]] <- UpdatePath(frags[[35]], new.path = frag_list[[35]])
frags[[36]] <- UpdatePath(frags[[36]], new.path = frag_list[[36]])
frags[[37]] <- UpdatePath(frags[[37]], new.path = frag_list[[37]])
frags[[38]] <- UpdatePath(frags[[38]], new.path = frag_list[[38]])
Fragments(human_multiome) <- frags # assign update list of fragment objects back to the assay
# %%
# Save RDS object (w/fragments file added)
saveRDS(human_multiome, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/human_multiome_w_frags.rds')
# %% [markdown]
# # Call peaks & runChromVAR on whole dataset
# %%
# Read RDS object (w/fragments file added)
human_multiome <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/human_multiome_w_frags.rds')
# %%
DefaultAssay(human_multiome) <- 'ATAC'
peaks_group <- CallPeaks(
object = human_multiome,
group.by = 'Group',
outdir = '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte',
fragment.tempdir = '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte',
macs2.path = '/n/groups/neuroduo/Bruno/jupytervenv_gcc92_gdal314/bin/macs2',
)
# %%
# Save object
saveRDS(peaks_group, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/peaks_group.rds')
# %%
# Quantify reads in peaks
macs2_counts <- FeatureMatrix(
fragments = Fragments(human_multiome),
features = peaks_group,
process_n = 5000,
cells = colnames(human_multiome)
)
# %%
# Create new chromatin assay with peaks
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
human_multiome[["peaks_group"]] <- CreateChromatinAssay(
counts = macs2_counts,
fragments = Fragments(human_multiome),
annotation = annotation,
genome = "hg38"
)
# %%
# Keep peaks only on standard chromosomes
peaks.keep <- seqnames(granges(human_multiome)) %in% standardChromosomes(granges(human_multiome))
# %%
human_multiome[["peaks_group_filt"]] <- CreateChromatinAssay(
counts = human_multiome[['peaks_group']]@data[as.vector(peaks.keep), ],
fragments = Fragments(human_multiome),
annotation = annotation,
genome = "hg38"
)
# %%
# Save object after peaks as chromatin assay
saveRDS(human_multiome, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250422_human_multiome_V2_peaks_group.rds')
# %%
# Run region stats
DefaultAssay(human_multiome) <- 'peaks_group_filt'
human_multiome <- RegionStats(human_multiome, genome = BSgenome.Hsapiens.UCSC.hg38)
# %%
# Get JASPAR database
jaspar_pfm <- getMatrixSet(
x = JASPAR2020,
opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)
# %%
# Add jaspar motifs
human_multiome <- AddMotifs(
object = human_multiome,
genome = BSgenome.Hsapiens.UCSC.hg38,
pfm = jaspar_pfm,
assay = 'peaks_group_filt'
)
# %%
# Run chromVAR for human dataset
human_multiome <- RunChromVAR(
object = human_multiome,
new.assay.name = "chromvar_peaks_group",
genome = BSgenome.Hsapiens.UCSC.hg38,
assay = 'peaks_group_filt'
)
# %%
saveRDS(human_multiome, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250423_human_multiome_jaspar_V2_peaks_group_chromvar.rds')
# %% [markdown]
# # Find differential motifs across all cell types & plot heatmap for Fig. 3
# %%
human_multiome_group_chromvar <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250423_human_multiome_jaspar_V2_peaks_group_chromvar.rds')
# %%
DefaultAssay(human_multiome_group_chromvar) <- 'chromvar_peaks_group'
# %%
# For each terminal-differentiated EN type: combine RG-vRG, RG-tRG, RG-oRG, IPC-EN, EN-Newborn, either EN-IT-Immature (l23-IT,L4-IT, L5-IT, L6-IT)
# or EN-Non-IT-Immature (L5-ET, L56_NP, L6-CT, L6b)
# For each terminal-differentiated IN type: combine IN-CGE-Immature w/ IN-CGE-VIP or IN-CGE-SNCG
# combine IN-MGE-Immature w/ IN-MGE-SST or IN-MGE-PV
# Combine Oligodendrocye-Immature w/Oligodendrocyte
# Combine Astrocyte-Immature w/Astrocyte-protoplasmic, Astrocyte-fibrous
# %%
# keep only chromvar assay for analysis (to prevent memory issue)
human_multiome_group_chromvar_only <- human_multiome_group_chromvar[['chromvar_peaks_group']]
# %%
# Convert to seurat object
human_multiome_group_chromvar_only_seu <- CreateSeuratObject(human_multiome_group_chromvar_only)
[email hidden] <- [email hidden]
# %%
# Generate list of downsampled Seurat objects for each mature type
down_obj_list <- list()
# IT neurons
EN_IT_list <- c('EN-L2_3-IT','EN-L4-IT','EN-L5-IT','EN-L6-IT')
for (i in EN_IT_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,
'RG-vRG',
'RG-tRG',
'RG-oRG',
'IPC-EN',
'EN-Newborn',
'EN-IT-Immature'))
if(dim(seurat_sub)[2] > 4000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =4000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 4000) {
down_obj_list[[i]] <- seurat_sub}
}
# ET neurons
EN_nonIT_list <- c('EN-L5-ET','EN-L5_6-NP','EN-L6-CT','EN-L6b')
for (i in EN_nonIT_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,'RG-vRG',
'RG-tRG',
'RG-oRG',
'IPC-EN',
'EN-Newborn',
'EN-Non-IT-Immature'))
if(dim(seurat_sub)[2] > 10000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =10000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 10000) {
down_obj_list[[i]] <- seurat_sub}
}
# CGE neurons
IN_CGE_list <- c('IN-CGE-VIP','IN-CGE-SNCG')
for (i in IN_CGE_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,'IN-CGE-Immature'))
if(dim(seurat_sub)[2] > 10000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =10000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 10000) {
down_obj_list[[i]] <- seurat_sub}
}
# MGE neurons
IN_MGE_list <- c('IN-MGE-SST','IN-MGE-PV')
for (i in IN_MGE_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,'IN-MGE-Immature'))
if(dim(seurat_sub)[2] > 6000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =6000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 6000) {
down_obj_list[[i]] <- seurat_sub}
}
# Oligodendrocyte
oligo_list <- c('Oligodendrocyte')
for (i in oligo_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,'Oligodendrocyte-Immature'))
if(dim(seurat_sub)[2] > 10000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =10000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 10000) {
down_obj_list[[i]] <- seurat_sub}
}
# Other cells
other_list <- c('Microglia','Vascular','OPC','IN-Mix-LAMP5')
for (i in other_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated == i)
if(dim(seurat_sub)[2] > 10000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =10000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 10000) {
down_obj_list[[i]] <- seurat_sub}
}
# Astrocytes
astro_list <- c('Astrocyte-Protoplasmic')
for (i in astro_list) {
seurat_sub <- subset(human_multiome_group_chromvar_only_seu, subset = type_updated %in% c(i,
'Astrocyte-Fibrous',
'Astrocyte-Immature'))
if(dim(seurat_sub)[2] > 10000) {
seurat_down <- seurat_sub[,sample(colnames(seurat_sub), size =10000, replace=F)]
down_obj_list[[i]] <- seurat_down}
else if (dim(seurat_sub)[2] <= 10000) {
down_obj_list[[i]] <- seurat_sub}
}
# %%
# Run FindMarkers on chromvar score for each cell type
DefaultAssay(human_multiome_group_chromvar) <- 'peaks_group_filt'
diff_chromvar_results <- lapply(down_obj_list, function(x){
Idents(x) <- 'Group'
test_res <- FindMarkers(
object = x,
ident.1 = 'Adolescence',
ident.2 = 'First_trimester',
mean.fxn = rowMeans,
test.use='MAST',
fc.name = "avg_diff")
test_res$TF <- ConvertMotifID(Motifs(human_multiome_group_chromvar),id=rownames(test_res))
return(test_res)
})
# %%
# Convert -log10(padj) for each type & calculate scaled padj for induced TFs
diff_chromvar_results_induced_TF <- lapply(diff_chromvar_results, function(x){
x$log10padj <- -log10(x$p_val_adj)
x_up <- x[x$avg_diff>0,]
x_up$scaled_padj <- x_up$log10padj / x_up[1,'log10padj']
return(x_up)
})
# Set cell type column
for (i in names(diff_chromvar_results_induced_TF)){
diff_chromvar_results_induced_TF[[i]]$celltype <- i
}
# %%
# Save results
saveRDS(diff_chromvar_results_induced_TF, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250425_diff_chromvar_results_induced_TF.rds')
# %%
# Read results
diff_chromvar_results_induced_TF <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250425_diff_chromvar_results_induced_TF.rds')
# %%
# Combine into dataframe
df <- do.call(rbind, diff_chromvar_results_induced_TF)
# %%
# Collect top 2 motifs per cell type
top_motifs <- list()
top_motifs <- lapply(diff_chromvar_results_induced_TF, function(x){
top2 <- x[1:2,'TF']
})
# Combine into single df
top_motif_df <- reshape2::melt(do.call(cbind, top_motifs))
# Collect unique motifs for plotting
uniq_motifs <- unique(top_motif_df$value)
# %%
# Subset results for unique motifs
df_motif_ploting <- df[df$TF %in% uniq_motifs,
c('TF','scaled_padj','celltype')]
# %%
# Cast dataframe to wide form
plot_df_cast <- reshape2::acast(setDT(df_motif_ploting), celltype~TF,value.var='scaled_padj')
plot_df_cast[is.na(plot_df_cast)] <- 0
# %%
# Order cell types by class
plot_df_cast_reorder <- plot_df_cast[c('EN-L2_3-IT',
'EN-L4-IT',
'EN-L5-IT',
'EN-L5-ET',
'EN-L5_6-NP',
'EN-L6-IT',
'EN-L6-CT',
'EN-L6b',
'IN-MGE-PV',
'IN-MGE-SST',
'IN-CGE-SNCG',
'IN-CGE-VIP',
'IN-Mix-LAMP5',
'Vascular',
'Microglia',
'Oligodendrocyte',
'OPC',
'Astrocyte-Protoplasmic'),]
rownames(plot_df_cast_reorder)[1] <- 'L2/3'
rownames(plot_df_cast_reorder)[2] <- 'L4IT'
rownames(plot_df_cast_reorder)[3] <- 'L5IT'
rownames(plot_df_cast_reorder)[4] <- 'L5PT'
rownames(plot_df_cast_reorder)[5] <- 'L5/6NP'
rownames(plot_df_cast_reorder)[6] <- 'L6IT'
rownames(plot_df_cast_reorder)[7] <- 'L6CT'
rownames(plot_df_cast_reorder)[8] <- 'L6b'
rownames(plot_df_cast_reorder)[9] <- 'PV+'
rownames(plot_df_cast_reorder)[10] <- 'SST+'
rownames(plot_df_cast_reorder)[11] <- 'SNCG+'
rownames(plot_df_cast_reorder)[12] <- 'VIP+'
rownames(plot_df_cast_reorder)[13] <- 'LAMP5+'
rownames(plot_df_cast_reorder)[14] <- 'Endo'
rownames(plot_df_cast_reorder)[15] <- 'Micro'
rownames(plot_df_cast_reorder)[16] <- 'Oligo'
rownames(plot_df_cast_reorder)[18] <- 'Astro'
# Order motifs
plot_df_cast_reorder <- plot_df_cast_reorder[,c('EGR1',
'Wt1',
'EGR3',
'ZNF75D',
'MEF2A',
'MEF2B',
'MEF2C',
'FOS',
'FOS::JUN',
'FOS::JUND',
'FOSL2',
'JUN(var.2)',
'Smad2::Smad3',
'BATF3',
'DBP',
'TEF',
'GABPA',
'ELF1',
'ELF3',
'SPI1',
'NR4A2',
'NR2C1',
'ASCL1',
'NFIC',
'NR3C1',
'NR3C2')]
# %%
# Re-label Smad2::Smad3 for space
colnames(plot_df_cast_reorder)[colnames(plot_df_cast_reorder) == 'Smad2::Smad3'] <- 'Smad2/3'
colnames(plot_df_cast_reorder)[colnames(plot_df_cast_reorder) == 'NR3C1'] <- 'GR'
colnames(plot_df_cast_reorder)[colnames(plot_df_cast_reorder) == 'NR3C2'] <- 'MR'
# %%
options(repr.plot.width = 8.7, repr.plot.height = 6.2,repr.plot.res=500)
col_fun = colorRamp2(seq(0, 1,length.out=9), rev(brewer.pal(9, "RdYlBu")))
p <- ComplexHeatmap::pheatmap(plot_df_cast_reorder,
color = col_fun,
cluster_rows = F,
cluster_cols = F,
fontsize = 20,
column_split = c(rep("Group1", 24), rep("Group2", 2)),
column_gap = unit(2, "mm"),
column_title = NULL,
row_names_side = "left",
border_color = 'grey30',
legend=F,
show_rownames = T)
custom_lgd <- Legend(
col_fun = col_fun,
title = "-log10\n(padj)",
legend_height = unit(0.1, "cm"),
legend_width = unit(.5, "cm"),
at = c(0,1),
title_position = "topleft",
labels = c('ns','max'),
title_gp = gpar(fontsize = 18),
labels_gp = gpar(fontsize = 16),
title_gap = unit(4, "mm") # <-- This WILL work here
)
draw(p, heatmap_legend_side="right", heatmap_legend_list = list(custom_lgd))
# %% [markdown]
# # Generate astrocyte GR chromVAR line plots, separately for V1 & PFC
# %%
# Load in dataset
human_multiome_group_chromvar <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250423_human_multiome_jaspar_V2_peaks_group_chromvar.rds')
# %%
DefaultAssay(human_multiome_group_chromvar) <- 'chromvar_peaks_group'
# %%
# keep only chromvar assay for analysis (to prevent memory issue)
human_multiome_group_chromvar_only <- human_multiome_group_chromvar[['chromvar_peaks_group']]
# %%
# Convert to seurat object
human_multiome_group_chromvar_only_seu <- CreateSeuratObject(human_multiome_group_chromvar_only)
[email hidden] <- [email hidden]
# %%
# Subset to astrocytes
human_multiome_group_chromvar_only_seu_astro <- subset(human_multiome_group_chromvar_only_seu, subset = subclass == 'Astrocyte')
# %%
options(repr.plot.width=3.9, repr.plot.height=4.5,repr.plot.res=500)
# Create dataframe of chromvar score
astro_nr3c1_df <- data.frame(counts = GetAssayData(human_multiome_group_chromvar_only_seu_astro, slot = "data")['MA0113.3',],
Region = human_multiome_group_chromvar_only_seu_astro$region_summary,
group = human_multiome_group_chromvar_only_seu_astro$Group,
pcd = human_multiome_group_chromvar_only_seu_astro$Estimated_postconceptional_age_in_days)
# Add time variable for plotting
astro_nr3c1_df$Time <- ifelse(astro_nr3c1_df$group=='First_trimester',1,
ifelse(astro_nr3c1_df$group=='Second_trimester',2,
ifelse(astro_nr3c1_df$group=='Third_trimester',3,
ifelse(astro_nr3c1_df$group=='Infancy',4,5))))
astro_nr3c1_df_summary <- summarySE(astro_nr3c1_df, measurevar="counts", conf.interval = 0.95, groupvars=c("Region","Time"))
# Replace General with CTX for plotting
astro_nr3c1_df_summary$Region <- gsub('General','CTX',astro_nr3c1_df_summary$Region)
scaleFUN <- function(x) sprintf("%.1f", x)
ggplot(astro_nr3c1_df_summary, aes(x=Time, y=counts, colour=Region)) +
geom_ribbon(aes(x=as.numeric(Time),ymax = counts + ci, ymin = counts - ci,fill=Region),
alpha = 0.3,
linetype=0) +
geom_line(aes(group=Region),linewidth=2) +
geom_point(size=4) +
geom_hline(yintercept=0,lwd=1,linetype=1,alpha=0.5) + theme_classic() +
ylab('GR chromVAR score') +
xlab('Stage') +
scale_y_continuous(labels=scaleFUN) +
theme(legend.position='right',
axis.text=element_text(size=16,color='black'),
legend.text=element_text(size=16,color='black'),
legend.margin=margin(0,8,0,0),
axis.text.x=element_text(angle=90,vjust=0.5,margin=margin(5,0,0,0)),
axis.title=element_text(size=18,color='black'),
legend.title=element_text(size=18,color='black'),
axis.title.y = element_text(margin = margin(0,5,0,0)),
axis.title.x = element_text(margin = margin(12,0,0,0))) + coord_cartesian(ylim=c(-1.5,1.5)) +
scale_x_continuous(breaks=seq(1,5,1),
labels = c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.')) +
scale_color_manual(name = "Human\nbrain region", values=c('grey70',
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[1],
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[8]))+
scale_fill_manual(name = "Human\nbrain region", values=c('grey70',
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[1],
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[8]))
# %%
options(repr.plot.width = 10, repr.plot.height = 3.5, repr.plot.res = 800)
VlnPlot(human_multiome_group_chromvar_only_seu_astro, group.by='Estimated_postconceptional_age_in_days', features= 'MA0113.3',pt.size=0) +
ggtitle('') +
theme(legend.position = 'none',
axis.ticks.x = element_blank(),
axis.title = element_text(size=14,color='black'),
axis.title.x = element_text(margin=margin(10,0,0,0)),
axis.text = element_text(size=12,color='black'),
axis.text.x = element_text(color='black',angle=90,hjust=1,vjust=0.5),
axis.line.x=element_blank()) +
geom_boxplot(fill='grey90',outlier.size=0,width=0.4,lwd=0.3,color='black') +
geom_hline(yintercept=-.9,alpha=0.5) +
ylab('Human astrocyte\nGR chromVAR score') +
xlab('Estimated days post-conception') +
scale_fill_manual(values = c(rep(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],3),
rep(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],7),
rep(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],5),
rep(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],5),
rep(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23],4)))
# %% [markdown]
# # Plot chromVAR results for astrocytes as separate volcano plots for V1 & PFC
# %%
# Run chromVAR separately on V1 & PFC astrocytes
human_multiome_group_chromvar_only_seu_astro <- subset(human_multiome_group_chromvar_only_seu, subset = subclass == 'Astrocyte')
# %%
# Subset to V1 and PFC cells
V1_astro <- human_multiome_group_chromvar_only_seu_astro[,human_multiome_group_chromvar_only_seu_astro$region_summary %in% c('General','V1')]
PFC_astro <- human_multiome_group_chromvar_only_seu_astro[,human_multiome_group_chromvar_only_seu_astro$region_summary %in% c('General','PFC')]
# %%
# Test in V1 astrocytes
Idents(V1_astro) <- 'Group'
test_V1_astro <- FindMarkers(
object = V1_astro,
ident.1 = 'Adolescence',
ident.2 = 'First_trimester',
mean.fxn = rowMeans,
test.use='MAST',
fc.name = "avg_diff")
DefaultAssay(human_multiome_group_chromvar) <-'peaks_group_filt'
test_V1_astro$TF <- ConvertMotifID(Motifs(human_multiome_group_chromvar),id=rownames(test_V1_astro))
# Test in PFC astrocytes
Idents(PFC_astro) <- 'Group'
test_PFC_astro <- FindMarkers(
object = PFC_astro,
ident.1 = 'Adolescence',
ident.2 = 'First_trimester',
mean.fxn = rowMeans,
test.use='MAST',
fc.name = "avg_diff")
test_PFC_astro$TF <- ConvertMotifID(Motifs(human_multiome_group_chromvar),id=rownames(test_PFC_astro))
# %%
# Make Volcano plot
options(repr.plot.width=3, repr.plot.height=4.3,repr.plot.res=900)
# Add a column of NAs
test_V1_astro$significant <- "NO"
# Set differential motifs to "induced or repressed"
test_V1_astro$significant[ (test_V1_astro$avg_diff > 1) & (test_V1_astro$p_val_adj < 0.01)] <- "Up"
test_V1_astro$significant[ (test_V1_astro$avg_diff < -1) & (test_V1_astro$p_val_adj < 0.01)] <- "Down"
test_V1_astro$significant <- factor(test_V1_astro$significant, levels = c("NO","Up","Down"))
test_pos <- test_V1_astro[test_V1_astro$avg_diff>0,]
p8 <- ggplot(data = test_pos %>% arrange(significant), aes(x = avg_diff, y = -log10(p_val), col = significant)) + geom_point() +
theme_classic() +
scale_color_manual(values = c("grey70",paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[8])) +
theme(axis.text = element_text(size = 20, color = 'black'),
axis.title = element_text(size = 20),
plot.margin=margin(0,40,0,0),
plot.title = element_text(size = 20, hjust = 0.5),
legend.position = 'None') +
xlab('Adoles. vs. 1st Tri.\nchromVAR FC') +
ylab('-log10(p-value)') + ylim(c(0,160)) + xlim(c(0,3))
p8
# %%
# Make Volcano plot
options(repr.plot.width=3, repr.plot.height=4.3,repr.plot.res=900)
# Add a column of NAs
test_PFC_astro$significant <- "NO"
# Set differential motifs to "induced or repressed"
test_PFC_astro$significant[ (test_PFC_astro$avg_diff > 1) & (test_PFC_astro$p_val_adj < 0.01)] <- "Up"
test_PFC_astro$significant[ (test_PFC_astro$avg_diff < -1) & (test_PFC_astro$p_val_adj < 0.01)] <- "Down"
test_PFC_astro$significant <- factor(test_PFC_astro$significant, levels = c("NO","Up","Down"))
test_pos <- test_PFC_astro[test_PFC_astro$avg_diff>0,]
p8 <- ggplot(data = test_pos %>% arrange(significant), aes(x = avg_diff, y = -log10(p_val), col = significant)) + geom_point() +
theme_classic() +
scale_color_manual(values = c("grey70",paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[1])) +
theme(axis.text = element_text(size = 20, color = 'black'),
axis.title = element_text(size = 20),
plot.margin=margin(0,40,0,0),
plot.title = element_text(size = 20, hjust = 0.5),
legend.position = 'None') +
xlab('Adoles. vs. 1st Tri.\nchromVAR FC') +
ylab('-log10(p-value)') + ylim(c(0,160)) + xlim(c(0,3))
p8
# %% [markdown]
# # Call peaks only on astrocytes
# %%
# Read RDS object (w/fragments file added)
human_multiome <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/human_multiome_w_frags.rds')
# %%
# Subset astrocytes
human_multiome_astro <- subset(human_multiome, subset = subclass == 'Astrocyte')
# %%
# Call peaks on astrocytes
DefaultAssay(human_multiome_astro) <- 'ATAC'
peaks <- CallPeaks(
object = human_multiome_astro,
macs2.path = '/n/groups/neuroduo/Bruno/jupytervenv_gcc92_gdal314/bin/macs2',
group.by = 'Group'
)
# %%
# Save object
saveRDS(peaks,'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/peaks_240121.rds')
# %%
# Quantify reads in peaks
macs2_counts <- FeatureMatrix(
fragments = Fragments(human_multiome_astro),
features = peaks,
cells = colnames(human_multiome_astro)
)
# %%
# Save object
saveRDS(macs2_counts,'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/macs2_counts_240121.rds')
# %%
# Create new chromatin assay with peaks
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
human_multiome_astro[["peaks"]] <- CreateChromatinAssay(
counts = macs2_counts,
fragments = Fragments(human_multiome_astro),
annotation = annotation,
genome = "hg38"
)
# %%
# Save object after peaks as chromatin assay
saveRDS(human_multiome_astro,'/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %% [markdown]
# # Compute astrocyte developmental DEGs (devDEGs) for Fig. 3
# %%
# Read RDS file
human_multiome_astro_v2 <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Aggregate counts
human_multiome_astro_agg <- AggregateExpression(human_multiome_astro_v2,
assays = 'RNA',
group.by = 'Ident',
slot='counts',
return.seurat = TRUE)
# %%
# Create metadata
meta_df <- unique(data.frame([email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,'Ident'],
[email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,c('Group','region_summary')]))
metadata_df <- data.frame(Sample = meta_df[,1],
Group = meta_df[,2],
Region = meta_df[,3],
Group_region = paste0(meta_df[,2], '_', meta_df[,3]))
# %%
# Run DESeq2 analysis (use this to match DESeq analysis used for bulk astrocyte RNA-seq)
deseq_obj <- DESeqDataSetFromMatrix(countData = GetAssayData(human_multiome_astro_agg,slot='counts'),
colData = metadata_df[match(colnames(GetAssayData(human_multiome_astro_agg,slot='counts')),metadata_df$Sample),],
design = ~ Group)
# %%
# Run Deseq2
deseq_obj <- DESeq(deseq_obj)
resdf <- na.omit(as.data.frame(results(deseq_obj, contrast=c("Group","Adolescence","First_trimester"),lfcThreshold=1)))
# %%
# Normalize data
ntd <- normTransform(deseq_obj)
# %%
# Compute mean expression per group for plotting
df <- data.frame(First_trimester = rowMeans(assay(ntd)[,c(13:16)]),
Second_trimester = rowMeans(assay(ntd)[,c(1:11)]),
Third_trimester = rowMeans(assay(ntd)[,c(12,21,36:38)]),
Infant = rowMeans(assay(ntd)[,c(17:20,24:29)]),
Adolescence = rowMeans(assay(ntd)[,c(22,23,30:35)]))
names(df) <- c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.')
# %%
# Set colors for plotting
my_palette <- colorRampPalette(rev(brewer.pal(9,'RdBu')))(1000) # 1000 colors
my_breaks <- seq(-2, 2, length.out = 1001)
# %%
# Generate initial plot
options(repr.plot.width = 2, repr.plot.height = 4,repr.plot.res=500)
t <- df[rownames(resdf[resdf$padj<0.01,]),]
p <- ComplexHeatmap::pheatmap(t(scale(t(t))),
cluster_cols = FALSE,
color = my_palette,
cluster_rows = T,
clustering_method = 'ward.D2',
cutree_rows=4,
treeheight_row = unit(6, "mm"),
show_rownames = FALSE,
heatmap_legend_param = list(title = "Scaled RNA",at = c(-2,0,2),
labels = c('-2','0','2'),
legend_height = unit(1, "cm"),
legend_width = unit(1.5, "cm"),
title_position = "topcenter",
legend_direction = "horizontal",
title_gp = gpar(fontsize = 10),
labels_gp = gpar(fontsize = 10)),
border_color=NA,
border=TRUE)
draw(p, heatmap_legend_side="bottom")
# %%
# Extract clustering results
p2 <- pheatmap::pheatmap(t,
cluster_cols = FALSE,
cluster_rows = T,
clustering_method = 'ward.D2',
cutree_rows=4,
scale = 'row',
silent=T)
t.clust <- as.data.frame(cbind(t,
cluster = cutree(p2$tree_row,
k = 4)))
table(t.clust$cluster)
# %%
# Suppose you want the slices to appear in order 3, 1, 2 & without dendrogram
new_order <- c(1, 3, 2, 4)
# Convert to factor with desired order
row_split <- factor(cutree(p2$tree_row, k = 4), levels = new_order)
row_order <- order.dendrogram(as.dendrogram(p2$tree_row))
options(repr.plot.width = 1.8, repr.plot.height = 4,repr.plot.res=900)
t <- df[rownames(resdf[resdf$padj<0.01,]),]
p <- ComplexHeatmap::pheatmap(t(scale(t(t))),
cluster_cols = FALSE,
color = my_palette,
cluster_rows = F,
clustering_method = 'ward.D2',
row_split = row_split,
row_order = row_order,
row_title =NULL,
treeheight_row = unit(6, "mm"),
show_rownames = FALSE,
border = TRUE,
heatmap_legend_param = list(title = "Scaled RNA",at = c(-2,0,2),
labels = c('-2','0','2'),
legend_height = unit(1, "cm"),
legend_width = unit(1.5, "cm"),
title_position = "topcenter",
legend_direction = "horizontal",
title_gp = gpar(fontsize = 10),
labels_gp = gpar(fontsize = 10)),
border_color=NA)
draw(p, heatmap_legend_side="bottom")
# %%
# Process data
sample_df <- data.frame(t(scale(t(data.frame(assay(ntd)[rownames(t.clust),])))))
sample_df$gene <- rownames(sample_df)
sample_df_melt <- reshape2::melt(sample_df, id='gene')
sample_df_melt$Group <- rep(metadata_df[match(colnames(GetAssayData(human_multiome_astro_agg,slot='counts')),metadata_df$Sample),c('Group')],
each=dim(t.clust)[1])
sample_df_melt$cluster <- as.character(rep(t.clust$cluster, times = (dim(sample_df)[2]-1)))
# %%
# Add time variable for plotting
sample_df_melt$Time <- ifelse(sample_df_melt$Group=='First_trimester',1,
ifelse(sample_df_melt$Group=='Second_trimester',2,
ifelse(sample_df_melt$Group=='Third_trimester',3,
ifelse(sample_df_melt$Group=='Infancy',4,5))))
sample_df_melt_summary <- summarySE(sample_df_melt, measurevar="value", conf.interval = 0.95, groupvars=c('Time','cluster'))
# %%
# Generate scaled expression for gene clusters as line plot
options(repr.plot.width=3.4, repr.plot.height=4.5,repr.plot.res=500)
scaleFUN <- function(x) sprintf("%.1f", x)
ggplot(sample_df_melt_summary, aes(x=Time, y=value, group=cluster, color=cluster)) +
geom_ribbon(aes(x=as.numeric(Time),ymax = value + se, ymin = value - se,fill=cluster,color=cluster),
alpha = 0.3,
linetype=0) +
geom_line(aes(group=cluster),linewidth=2) +
geom_point(size=4) +
geom_hline(yintercept=0,lwd=1,linetype=1,alpha=0.5) + theme_classic() +
ylab('Scaled expression') +
xlab('Stage') +
scale_y_continuous(labels=scaleFUN) +
theme(legend.position='right',
axis.text=element_text(size=16,color='black'),
legend.text=element_text(size=16,color='black'),
legend.margin=margin(0,8,0,0),
axis.text.x=element_text(angle=90,vjust=0.5,margin=margin(5,0,0,0)),
axis.title=element_text(size=18,color='black'),
legend.title=element_text(size=18,color='black'),
axis.title.y = element_text(margin = margin(0,5,0,0)),
axis.title.x = element_text(margin = margin(12,0,0,0))) +
scale_x_continuous(breaks=seq(1,5,1),
labels = c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.')) +
scale_color_manual(name = 'Cluster', values = as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[1:4])) +
scale_fill_manual(name = 'Cluster', values = as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[1:4]))
# %% [markdown]
# # Compute overlap b/w devDEG clusters & mouse GR targets for Fig. 3
# %%
# Convert mouse genes to human
library(dplyr)
mouse_human_genes = read.csv("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt",sep="\t")
convert_mouse_to_human <- function(gene_list){
output = c()
for(gene in gene_list){
class_key = (mouse_human_genes %>% filter(Symbol == gene & Common.Organism.Name=="mouse, laboratory"))[['DB.Class.Key']]
if(!identical(class_key, integer(0)) ){
human_genes = (mouse_human_genes %>% filter(DB.Class.Key == class_key & Common.Organism.Name=="human"))[,"Symbol"]
for(human_gene in human_genes){
output = append(output,human_gene)
}
}
}
return (output)
}
# %%
# Read in GR-target gene data
P21_DEG_relax <- readRDS('/n/groups/neuroduo/Bruno/Astro_NR_DR_RNA_analysis/P21_GR_dep_res_relax.rds')
# %%
# Convert mouse to human gene symbols
P21_DEG_relax_sig_human <- convert_mouse_to_human(rownames(P21_DEG_relax[P21_DEG_relax$padj<0.05,]))
# %%
# Compute overlap with each cluster
overlap_df <- data.frame(overlap = c(table(rownames(t.clust[t.clust$cluster == 1,]) %in%
P21_DEG_relax_sig_human)[2] / length(rownames(t.clust[t.clust$cluster == 1,])),
table(rownames(t.clust[t.clust$cluster == 2,]) %in%
P21_DEG_relax_sig_human)[2] / length(rownames(t.clust[t.clust$cluster == 2,])),
table(rownames(t.clust[t.clust$cluster == 3,]) %in%
P21_DEG_relax_sig_human)[2] / length(rownames(t.clust[t.clust$cluster == 3,])),
table(rownames(t.clust[t.clust$cluster == 4,]) %in%
P21_DEG_relax_sig_human)[2] / length(rownames(t.clust[t.clust$cluster == 4,]))),
cluster = c('1','2','3','4'))
overlap_df$overlap <- overlap_df$overlap * 100
overlap_df$x <- '1'
# %%
# Calculate overlaps for Fisher test
ov1 <- rbind(table(rownames(t.clust[t.clust$cluster == 1,]) %in%
P21_DEG_relax_sig_human),
table(rownames(resdf[rownames(resdf) %notin% rownames(t.clust[t.clust$cluster ==1,]),]) %in%
P21_DEG_relax_sig_human))
ov2 <- rbind(table(rownames(t.clust[t.clust$cluster == 2,]) %in%
P21_DEG_relax_sig_human),
table(rownames(resdf[rownames(resdf) %notin% rownames(t.clust[t.clust$cluster ==2,]),]) %in%
P21_DEG_relax_sig_human))
ov3 <- rbind(table(rownames(t.clust[t.clust$cluster == 3,]) %in%
P21_DEG_relax_sig_human),
table(rownames(resdf[rownames(resdf) %notin% rownames(t.clust[t.clust$cluster ==3,]),]) %in%
P21_DEG_relax_sig_human))
ov4 <- rbind(table(rownames(t.clust[t.clust$cluster == 4,]) %in%
P21_DEG_relax_sig_human),
table(rownames(resdf[rownames(resdf) %notin% rownames(t.clust[t.clust$cluster ==4,]),]) %in%
P21_DEG_relax_sig_human))
# %%
# Add Fisher test results to object
overlap_df$pval <- c(fisher.test(ov1)$p.value,
fisher.test(ov2)$p.value,
fisher.test(ov3)$p.value,
fisher.test(ov4)$p.value)
overlap_df$logpval <- -log10(overlap_df$pval)
# %%
# Set order for plotting
overlap_df$cluster <- factor(overlap_df$cluster, levels = rev(c('1','3','2','4')))
# %%
# Generate dot plot of overlap and Fisher test p-value
options(repr.plot.width=1.6, repr.plot.height=2,repr.plot.res=1200)
ggplot(overlap_df, aes(x=x, y=cluster, size = overlap, fill=logpval)) +
geom_point(shape=21,alpha=1,stroke=0.4) + theme_bw() + scale_radius(breaks = c(6,18)) + theme(axis.title = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_line(color = "grey98"),
axis.text.y = element_blank(),
legend.title=element_text(hjust=0,size=10),
legend.text=element_text(size=10),
legend.margin=margin(0,5,0,0),
legend.key.height = unit(1,'mm'),
legend.key.width = unit(4,'mm'),
legend.position='right',
axis.ticks.x = element_blank()) + labs(size = "% overlap",fill='-log10\n(p-val)') +
scale_fill_viridis(option='rocket',breaks=c(2,25),labels=c('ns','25'))
# %% [markdown]
# # Link peaks to astrocyte devDEGs for motif analysis
# %%
library(future)
plan("multicore", workers = 8)
options(future.globals.maxSize = 100 * 1024 ^ 3) # for 100 Gb RAM
# %%
# Read RDS file
human_multiome_astro_v2 <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Run region stats
DefaultAssay(human_multiome_astro_v2) <- 'peaks'
human_multiome_astro_v2 <- RegionStats(human_multiome_astro_v2, genome = BSgenome.Hsapiens.UCSC.hg38)
# %%
# Get JASPAR database
jaspar_pfm <- getMatrixSet(
x = JASPAR2020,
opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)
# %%
# Add jaspar motifs
human_multiome_astro_v2 <- AddMotifs(
object = human_multiome_astro_v2,
genome = BSgenome.Hsapiens.UCSC.hg38,
pfm = jaspar_pfm,
assay = 'peaks'
)
# %%
# Save RDS file w/motif object for astrocyte-called peaks
saveRDS(human_multiome_astro_v2, '/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Run LinkPeaks on devDEGs
DefaultAssay(human_multiome_astro_v2) <- 'RNA'
c1_links <- LinkPeaks(human_multiome_astro_v2,
expression.assay = 'RNA',
peak.assay = 'peaks',
min.distance = 5000,
distance = 300000,
genes.use = rownames(t.clust[t.clust$cluster == 1,]))
c2_links <- LinkPeaks(human_multiome_astro_v2,
expression.assay = 'RNA',
peak.assay = 'peaks',
min.distance = 5000,
distance = 300000,
genes.use = rownames(t.clust[t.clust$cluster == 2,]))
c3_links <- LinkPeaks(human_multiome_astro_v2,
expression.assay = 'RNA',
peak.assay = 'peaks',
min.distance = 5000,
distance = 300000,
genes.use = rownames(t.clust[t.clust$cluster == 3,]))
c4_links <- LinkPeaks(human_multiome_astro_v2,
expression.assay = 'RNA',
peak.assay = 'peaks',
min.distance = 5000,
distance = 300000,
genes.use = rownames(t.clust[t.clust$cluster == 4,]))
# %%
# Convert results to dataframe
DefaultAssay(c1_links) <- 'peaks'
c1_links_df <- as.data.frame(Links(c1_links))
DefaultAssay(c2_links) <- 'peaks'
c2_links_df <- as.data.frame(Links(c2_links))
DefaultAssay(c3_links) <- 'peaks'
c3_links_df <- as.data.frame(Links(c3_links))
DefaultAssay(c4_links) <- 'peaks'
c4_links_df <- as.data.frame(Links(c4_links))
# %%
# Find enriched motifs in top-ranked peaks for each cluster
enriched.motifs1 <- FindMotifs(
object = c1_links,
c1_links_df[(c1_links_df$pvalue <0.001) &
(c1_links_df$score>0),'peak']
)
enriched.motifs2 <- FindMotifs(
object = c2_links,
c2_links_df[(c2_links_df$pvalue <0.001) &
(c2_links_df$score>0),'peak']
)
enriched.motifs3 <- FindMotifs(
object = c3_links,
c3_links_df[(c3_links_df$pvalue <0.001) &
(c3_links_df$score>0),'peak']
)
enriched.motifs4 <- FindMotifs(
object = c4_links,
c4_links_df[(c4_links_df$pvalue <0.001) &
(c4_links_df$score>0),'peak']
)
# %%
# Generate plot for cluster 1
options(repr.plot.width=3, repr.plot.height=4,repr.plot.res=900)
enriched.motifs1$motif.name <- gsub('FOSL1::JUN\$var\\.2\$', 'FOSL1::\nJUN(var.2)', enriched.motifs1$motif.name)
enriched.motifs1$motif.name <- factor(enriched.motifs1$motif.name, levels = enriched.motifs1$motif.name)
p3 <- ggplot(enriched.motifs1[1:5,], aes(x=motif.name, y=-log10(p.adjust))) +
geom_bar(stat = "identity",width=0.7,fill=as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[1])) + theme_classic() +
xlab('TF motif') +
ylab('-log10(padj)') +
theme(axis.text = element_text(color='black',size=18),
axis.text.x = element_text(angle=45,hjust=1),
axis.title = element_text(color='black',size=18),
plot.margin=margin(5,8,5,0),
axis.title.x=element_blank())
# %%
# Generate plot for cluster 2
options(repr.plot.width=3, repr.plot.height=4,repr.plot.res=900)
enriched.motifs2$motif.name <- gsub('NR3C1', 'GR', enriched.motifs2$motif.name)
enriched.motifs2$motif.name <- gsub('NR3C2', 'MR', enriched.motifs2$motif.name)
enriched.motifs2$motif.name <- factor(enriched.motifs2$motif.name, levels = enriched.motifs2$motif.name)
p2 <- ggplot(enriched.motifs2[1:5,], aes(x=motif.name, y=-log10(p.adjust))) +
geom_bar(stat = "identity",width=0.7,fill=as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[2])) + theme_classic() +
xlab('TF motif') +
ylab('-log10(padj)') +
theme(axis.text = element_text(color='black',size=18),
axis.text.x = element_text(angle=45,hjust=1),
plot.margin=margin(40,8,-10,0),
axis.title = element_text(color='black',size=18),
axis.title.x=element_blank())
# %%
# Generate plot for cluster 3
options(repr.plot.width=3, repr.plot.height=4,repr.plot.res=900)
enriched.motifs3$motif.name <- factor(enriched.motifs3$motif.name, levels = enriched.motifs3$motif.name)
p4 <- ggplot(enriched.motifs3[1:5,], aes(x=motif.name, y=-log10(p.adjust))) +
geom_bar(stat = "identity",width=0.7,fill=as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[3])) + theme_classic() +
xlab('TF motif') +
ylab('-log10(padj)') +
theme(axis.text = element_text(color='black',size=18),
axis.text.x = element_text(angle=45,hjust=1),
plot.margin=margin(5,0,5,0),
axis.title = element_text(color='black',size=18),
axis.title.x=element_blank())
# %%
# Generate plot for cluster 4
options(repr.plot.width=3, repr.plot.height=4,repr.plot.res=900)
enriched.motifs4$motif.name <- factor(enriched.motifs4$motif.name, levels = enriched.motifs4$motif.name)
p1 <- ggplot(enriched.motifs4[1:5,], aes(x=motif.name, y=-log10(p.adjust))) +
geom_bar(stat = "identity",width=0.7,fill=as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[4])) + theme_classic() +
xlab('TF motif') +
ylab('-log10(padj)') +
theme(axis.text = element_text(color='black',size=18),
axis.text.x = element_text(angle=45,hjust=1),
plot.margin=margin(40,0,-10,0),
axis.title = element_text(color='black',size=18),
axis.title.x=element_blank())
# %%
# Assemble plots
options(repr.plot.width=5.2, repr.plot.height=6.8,repr.plot.res=1200)
cowplot::plot_grid(p3,p4,p2,p1,ncol=2,align='hv')
# %% [markdown]
# # Identify astrocyte developmental DARs (devDARs) using DESeq2
# %%
# Read RDS file w/motif object for astrocyte-called peaks
human_multiome_astro_v2 <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Aggregate ATAC peak counts per sample
human_multiome_astro_agg <- AggregateExpression(human_multiome_astro_v2,
assays = 'peaks',
group.by = 'Ident',
slot = 'counts',
return.seurat = TRUE)
# %%
# Create metadata
meta_df <- unique(data.frame([email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,'Ident'],
[email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,c('Group','region_summary')]))
metadata_df <- data.frame(Sample = meta_df[,1],
Group = meta_df[,2],
Region = meta_df[,3],
Group_region = paste0(meta_df[,2], '_', meta_df[,3]))
# %%
# Run DESeq2 analysis on ATAC peaks
deseq_obj <- DESeqDataSetFromMatrix(countData = GetAssayData(human_multiome_astro_agg,slot='counts'),
colData = metadata_df[match(colnames(GetAssayData(human_multiome_astro_agg,slot='counts')),metadata_df$Sample),],
design = ~ Group)
# Run Deseq2
deseq_obj <- DESeq(deseq_obj)
resdf <- na.omit(as.data.frame(results(deseq_obj,contrast=c("Group","Adolescence","First_trimester"),lfcThreshold=1)))
# %%
# Normalize data
ntd <- normTransform(deseq_obj)
# %%
# Compute mean expression per group for plotting
df <- data.frame(First_trimester = rowMeans(assay(ntd)[,c(13:16)]),
Second_trimester = rowMeans(assay(ntd)[,c(1:11)]),
Third_trimester = rowMeans(assay(ntd)[,c(12,21,36:38)]),
Infant = rowMeans(assay(ntd)[,c(17:20,24:29)]),
Adolescence = rowMeans(assay(ntd)[,c(22,23,30:35)]))
names(df) <- c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.')
# %%
# Create color palette
my_palette <- colorRampPalette(rev(brewer.pal(9,'RdBu')))(1000) # 1000 colors
my_breaks <- seq(-2, 2, length.out = 1001)
# %%
# Generate initial heatmap
options(repr.plot.width = 2, repr.plot.height = 4,repr.plot.res=500)
t <- df[rownames(resdf[resdf$padj<0.05,]),]
p <- ComplexHeatmap::pheatmap(t(scale(t(t))),
cluster_cols = FALSE,
color = my_palette,
cluster_rows = T,
clustering_method = 'ward.D2',
cutree_rows=4,
treeheight_row = unit(6, "mm"),
show_rownames = FALSE,
heatmap_legend_param = list(title = "Scaled ATAC",at = c(-2,0,2),
labels = c('-2','0','2'),
legend_height = unit(1, "cm"),
legend_width = unit(1.5, "cm"),
title_position = "topcenter",
legend_direction = "horizontal",
title_gp = gpar(fontsize = 10),
labels_gp = gpar(fontsize = 10)),
border_color=NA,
border=TRUE)
draw(p, heatmap_legend_side="bottom")
# %%
# Extract clustering results
t <- df[rownames(resdf[resdf$padj<0.05,]),]
p2 <- pheatmap::pheatmap(t,
cluster_cols = FALSE,
cluster_rows = T,
clustering_method = 'ward.D2',
cutree_rows=4,
scale = 'row',
silent=T)
t.clust <- as.data.frame(cbind(t,
cluster = cutree(p2$tree_row,
k = 4)))
table(t.clust$cluster)
# %%
# Suppose you want the slices to appear in order & without dendrogram
new_order <- c(3, 1, 4, 2)
# Convert to factor with desired order
row_split <- factor(cutree(p2$tree_row, k = 4), levels = new_order)
row_order <- order.dendrogram(as.dendrogram(p2$tree_row))
options(repr.plot.width = 1.8, repr.plot.height = 4,repr.plot.res=900)
t <- df[rownames(resdf[resdf$padj<0.05,]),]
p <- ComplexHeatmap::pheatmap(t(scale(t(t))),
cluster_cols = FALSE,
color = my_palette,
cluster_rows = F,
clustering_method = 'ward.D2',
row_split = row_split,
row_order = row_order,
row_title =NULL,
treeheight_row = unit(6, "mm"),
show_rownames = FALSE,
border = TRUE,
heatmap_legend_param = list(title = "Scaled ATAC",at = c(-2,0,2),
labels = c('-2','0','2'),
legend_height = unit(1, "cm"),
legend_width = unit(1.5, "cm"),
title_position = "topcenter",
legend_direction = "horizontal",
title_gp = gpar(fontsize = 10),
labels_gp = gpar(fontsize = 10)),
border_color=NA)
draw(p, heatmap_legend_side="bottom")
# %%
# Process data
sample_df <- data.frame(t(scale(t(data.frame(assay(ntd)[rownames(t.clust),])))))
sample_df$peak <- rownames(sample_df)
sample_df_melt <- reshape2::melt(sample_df, id='peak')
sample_df_melt$Group <- rep(metadata_df[match(colnames(GetAssayData(human_multiome_astro_agg,slot='counts')),metadata_df$Sample),c('Group')],
each=dim(t.clust)[1])
sample_df_melt$cluster <- as.character(rep(t.clust$cluster, times = (dim(sample_df)[2]-1)))
# %%
# Add time variable for plotting
sample_df_melt$Time <- ifelse(sample_df_melt$Group=='First_trimester',1,
ifelse(sample_df_melt$Group=='Second_trimester',2,
ifelse(sample_df_melt$Group=='Third_trimester',3,
ifelse(sample_df_melt$Group=='Infancy',4,5))))
sample_df_melt_summary <- summarySE(sample_df_melt, measurevar="value", conf.interval = 0.95, groupvars=c('Time','cluster'))
# %%
# Generate line plots of scaled accessibility within peak clusters
options(repr.plot.width=3.4, repr.plot.height=4.5,repr.plot.res=500)
scaleFUN <- function(x) sprintf("%.1f", x)
ggplot(sample_df_melt_summary, aes(x=Time, y=value, group=cluster, color=cluster)) +
geom_ribbon(aes(x=as.numeric(Time),ymax = value + se, ymin = value - se,fill=cluster,color=cluster),
alpha = 0.3,
linetype=0) +
geom_line(aes(group=cluster),linewidth=2) +
geom_point(size=4) +
geom_hline(yintercept=0,lwd=1,linetype=1,alpha=0.5) + theme_classic() +
ylab('Scaled accessibility') +
xlab('Stage') +
scale_y_continuous(labels=scaleFUN) +
theme(legend.position='right',
axis.text=element_text(size=16,color='black'),
legend.text=element_text(size=16,color='black'),
legend.margin=margin(0,8,0,0),
axis.text.x=element_text(angle=90,vjust=0.5,margin=margin(5,0,0,0)),
axis.title=element_text(size=18,color='black'),
legend.title=element_text(size=18,color='black'),
axis.title.y = element_text(margin = margin(0,5,0,0)),
axis.title.x = element_text(margin = margin(12,0,0,0))) +
scale_x_continuous(breaks=seq(1,5,1),
labels = c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.')) +
scale_color_manual(name = 'Cluster', values = as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[c(3,4,1,2)])) +
scale_fill_manual(name = 'Cluster', values = as.vector(paletteer::paletteer_d("fishualize::Pseudocheilinus_tetrataenia")[c(3,4,1,2)]))
# %% [markdown]
# # Compare to astrocyte devDARs to liftover mouse GR binding sites
# %%
# Read in GR sites (lifted over already to hg38 w/UCSC)
GR_sites <- rtracklayer::import.bed('/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/CR01_count_edger_report_str_FC25.col.human.bed')
# %%
# Collect ATAC peaks in each cluster
c1 <- t.clust[t.clust$cluster==1,]
c2 <- t.clust[t.clust$cluster==2,]
c3 <- t.clust[t.clust$cluster==3,]
c4 <- t.clust[t.clust$cluster==4,]
# %%
# Convert human ATAC sites to granges objects
c1_gr <- GRanges(
seqnames = do.call(rbind,strsplit(rownames(c1),'-'))[,1],
ranges = IRanges(start = as.numeric(do.call(rbind,strsplit(rownames(c1),'-'))[,2]),
end = as.numeric(do.call(rbind,strsplit(rownames(c1),'-'))[,3])))
c2_gr <- GRanges(
seqnames = do.call(rbind,strsplit(rownames(c2),'-'))[,1],
ranges = IRanges(start = as.numeric(do.call(rbind,strsplit(rownames(c2),'-'))[,2]),
end = as.numeric(do.call(rbind,strsplit(rownames(c2),'-'))[,3])))
c3_gr <- GRanges(
seqnames = do.call(rbind,strsplit(rownames(c3),'-'))[,1],
ranges = IRanges(start = as.numeric(do.call(rbind,strsplit(rownames(c3),'-'))[,2]),
end = as.numeric(do.call(rbind,strsplit(rownames(c3),'-'))[,3])))
c4_gr <- GRanges(
seqnames = do.call(rbind,strsplit(rownames(c4),'-'))[,1],
ranges = IRanges(start = as.numeric(do.call(rbind,strsplit(rownames(c4),'-'))[,2]),
end = as.numeric(do.call(rbind,strsplit(rownames(c4),'-'))[,3])))
# Generate total ATAC sites as background for Fisher's test
DefaultAssay(human_multiome_astro_v2) <- 'peaks'
total_gr <- GRanges(
seqnames = do.call(rbind,strsplit(rownames(human_multiome_astro_v2),'-'))[,1],
ranges = IRanges(start = as.numeric(do.call(rbind,strsplit(rownames(human_multiome_astro_v2),'-'))[,2]),
end = as.numeric(do.call(rbind,strsplit(rownames(human_multiome_astro_v2),'-'))[,3])))
# %%
# Compute overlap with each cluster
overlap_df <- data.frame(overlap = c(length(findOverlaps(GR_sites, c1_gr)) / length(c1_gr),
length(findOverlaps(GR_sites, c2_gr)) / length(c2_gr),
length(findOverlaps(GR_sites, c3_gr)) / length(c3_gr),
length(findOverlaps(GR_sites, c4_gr)) / length(c4_gr)),
oldcluster = c('1','2','3','4'),
newcluster = c('3','1','4','2'))
overlap_df$overlap <- overlap_df$overlap * 100
overlap_df$x <- '1'
# %%
# Calculate overlaps for Fisher test
ov1 <- rbind(data.frame(neg = length(c1_gr) - length(findOverlaps(c1_gr,GR_sites)),
pos = length(findOverlaps(c1_gr,GR_sites))),
data.frame(neg = length(total_gr[-subjectHits(findOverlaps(c1_gr,total_gr))]) - length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c1_gr,total_gr))])),
pos = length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c1_gr,total_gr))]))))
ov2 <- rbind(data.frame(neg = length(c2_gr) - length(findOverlaps(c2_gr, GR_sites)),
pos = length(findOverlaps(c2_gr,GR_sites))),
data.frame(neg = length(total_gr[-subjectHits(findOverlaps(c2_gr,total_gr))]) - length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c2_gr,total_gr))])),
pos = length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c2_gr,total_gr))]))))
ov3 <- rbind(data.frame(neg = length(c3_gr) - length(findOverlaps(c3_gr,GR_sites)),
pos = length(findOverlaps(c3_gr,GR_sites))),
data.frame(neg = length(total_gr[-subjectHits(findOverlaps(c3_gr,total_gr))]) - length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c3_gr,total_gr))])),
pos = length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c3_gr,total_gr))]))))
ov4 <- rbind(data.frame(neg = length(c4_gr) - length(findOverlaps(c4_gr,GR_sites)),
pos = length(findOverlaps(c4_gr,GR_sites))),
data.frame(neg = length(total_gr[-subjectHits(findOverlaps(c4_gr,total_gr))]) - length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c4_gr,total_gr))])),
pos = length(findOverlaps(GR_sites,total_gr[-subjectHits(findOverlaps(c4_gr,total_gr))]))))
# %%
# Add Fisher test results to object
overlap_df$pval <- c(fisher.test(ov1,alternative='less')$p.value,
fisher.test(ov2,alternative='less')$p.value,
fisher.test(ov3,alternative='less')$p.value,
fisher.test(ov4,alternative='less')$p.value)
overlap_df$logpval <- -log10(overlap_df$pval)
# %%
# Set order for plotting
overlap_df$oldcluster <- factor(overlap_df$oldcluster, levels = rev(c('3','1','4','2')))
# %%
# Print results
overlap_df
# %%
# Generate dot plot of % overlap and Fisher's test p-value
options(repr.plot.width=1.6, repr.plot.height=2.0,repr.plot.res=1200)
ggplot(overlap_df, aes(x=x, y=oldcluster, size = overlap, fill=logpval)) +
geom_point(shape=21,alpha=1,stroke=0.4) + theme_bw() + scale_radius(breaks = c(4,14)) + theme(axis.title = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_line(color = "grey98"),
axis.text.y = element_blank(),
legend.title=element_text(hjust=0,size=10),
legend.text=element_text(size=10),
legend.margin=margin(0,5,0,0),
legend.key.height = unit(1,'mm'),
legend.key.width = unit(4,'mm'),
legend.position='right',
axis.ticks.x = element_blank()) + labs(size = "% overlap",fill='-log10\n(p-val)') +
scale_fill_viridis(option='rocket',breaks=c(2,50),labels=c('ns','50'))
# %% [markdown]
# # Generate human astrocyte UMAPs
# %%
# Read RDS file w/motif object for astrocyte-called peaks
human_multiome_astro_v2 <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Generate plots
options(repr.plot.width = 4.5*1.5, repr.plot.height = 3*1.5, repr.plot.res=500)
human_multiome_astro_v2$Group <- factor(human_multiome_astro_v2$Group,levels=
c('First_trimester','Second_trimester','Third_trimester','Infancy','Adolescence'))
p1 <- DimPlot(human_multiome_astro_v2, group.by = 'Group',shuffle=T) + theme_void() + theme(plot.title=element_blank(),
plot.margin=margin(0,10,0,0),
legend.position='none') +
scale_color_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23]))
human_multiome_astro_v2$region_summary <- factor(human_multiome_astro_v2$region_summary,levels = c('General',
'PFC',
'V1'))
p2 <- DimPlot(human_multiome_astro_v2, group.by = 'region_summary',shuffle=T) + theme_void() + theme(plot.title=element_blank(),
plot.margin=margin(0,0,0,10),
legend.position='none') +
scale_color_manual(values=c('grey70',
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[1],
paletteer::paletteer_d("colorBlindness::Blue2Orange8Steps")[8]))
p1 + p2
# %% [markdown]
# # Generate human and mouse astrocyte RNA & track plots
# %%
# Load in dataset
human_multiome_astro_v2 <- readRDS('/n/groups/neuroduo/Bruno/Wang_2025_human_astrocyte/250427_human_multiome_astro_v2.rds')
# %%
# Aggregate RNA counts by sample
human_multiome_astro_agg <- AggregateExpression(human_multiome_astro_v2,
assays = 'RNA',
group.by = 'Ident',
return.seurat = TRUE)
# %%
# Create metadata
meta_df <- unique(data.frame([email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,'Ident'],
[email hidden][human_multiome_astro_v2$Ident %in% human_multiome_astro_v2$Ident,c('Group','region_summary')]))
metadata_df <- data.frame(Sample = meta_df[,1],
Group = meta_df[,2],
Region = meta_df[,3],
Group_region = paste0(meta_df[,2], '_', meta_df[,3]))
# %%
# Create deseq2 object
deseq_obj <- DESeqDataSetFromMatrix(countData = GetAssayData(human_multiome_astro_agg,slot='counts'),
colData = metadata_df[match(colnames(GetAssayData(human_multiome_astro_agg,slot='counts')),metadata_df$Sample),],
design = ~ Group)
# %%
# Generate boxplots
options(repr.plot.width=2.4, repr.plot.height=4,repr.plot.res=500)
goi <- 'ETNPPL'
data <- plotCounts(deseq_obj, gene = goi, intgroup = c('Group'),returnData = TRUE)
data$Group <- gsub('First_trimester', '1st Tri.',data$Group)
data$Group <- gsub('Second_trimester', '2nd Tri.',data$Group)
data$Group <- gsub('Third_trimester', '3rd Tri.',data$Group)
data$Group <- gsub('Infancy', 'Infant', data$Group)
data$Group <- gsub('Adolescence', 'Adoles.', data$Group)
data$Group <- factor(data$Group, levels = c('1st Tri.','2nd Tri.',
'3rd Tri.','Infant','Adoles.'))
ggplot(data, aes(x=Group, y=count, color=Group,fill=Group)) +
geom_boxplot(alpha=0.2,outlier.size=0)+
geom_point(size=1.5, shape=21,position = position_jitterdodge(jitter.width=1)) + theme_classic() +
theme(axis.title = element_text(size=16,color='black'),
axis.text = element_text(size=14,color='black'),
axis.text.x = element_text(angle=90,vjust=0.5,hjust=1),
legend.text = element_text(size=14,color='black'),
axis.title.y = element_text(margin = margin(0,2,0,0)),
axis.title.x = element_text(margin = margin(10,0,0,0)),
plot.margin=margin(5,22,0,1),
legend.box.margin=margin(0,0,-10,-10),
legend.position='none') + #ylim(c(0,max(data$count) * 1.1)) +
scale_color_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
scale_fill_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
guides(fill = "none") + labs(color="") + labs(x = 'Stage', y = paste0(goi, ' expression'))
# %%
# Generate boxplots
options(repr.plot.width=2.4, repr.plot.height=4,repr.plot.res=500)
goi <- 'GJB6'
data <- plotCounts(deseq_obj, gene = goi, intgroup = c('Group'),returnData = TRUE)
data$Group <- gsub('First_trimester', '1st Tri.',data$Group)
data$Group <- gsub('Second_trimester', '2nd Tri.',data$Group)
data$Group <- gsub('Third_trimester', '3rd Tri.',data$Group)
data$Group <- gsub('Infancy', 'Infant', data$Group)
data$Group <- gsub('Adolescence', 'Adoles.', data$Group)
data$Group <- factor(data$Group, levels = c('1st Tri.','2nd Tri.',
'3rd Tri.','Infant','Adoles.'))
ggplot(data, aes(x=Group, y=count, color=Group,fill=Group)) +
geom_boxplot(alpha=0.2,outlier.size=0)+
geom_point(size=1.5, shape=21,position = position_jitterdodge(jitter.width=1)) + theme_classic() +
theme(axis.title = element_text(size=16,color='black'),
axis.text = element_text(size=14,color='black'),
axis.text.x = element_text(angle=90,vjust=0.5,hjust=1),
legend.text = element_text(size=14,color='black'),
axis.title.y = element_text(margin = margin(0,2,0,0)),
axis.title.x = element_text(margin = margin(10,0,0,0)),
plot.margin=margin(5,22,0,1),
legend.box.margin=margin(0,0,-10,-10),
legend.position='none') + #ylim(c(0,max(data$count) * 1.1)) +
scale_color_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
scale_fill_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
guides(fill = "none") + labs(color="") + labs(x = 'Stage', y = paste0(goi, ' expression'))
# %%
# Generate boxplots
options(repr.plot.width=2.4, repr.plot.height=4,repr.plot.res=500)
goi <- 'HTRA1'
data <- plotCounts(deseq_obj, gene = goi, intgroup = c('Group'),returnData = TRUE)
data$Group <- gsub('First_trimester', '1st Tri.',data$Group)
data$Group <- gsub('Second_trimester', '2nd Tri.',data$Group)
data$Group <- gsub('Third_trimester', '3rd Tri.',data$Group)
data$Group <- gsub('Infancy', 'Infant', data$Group)
data$Group <- gsub('Adolescence', 'Adoles.', data$Group)
data$Group <- factor(data$Group, levels = c('1st Tri.','2nd Tri.',
'3rd Tri.','Infant','Adoles.'))
ggplot(data, aes(x=Group, y=count, color=Group,fill=Group)) +
geom_boxplot(alpha=0.2,outlier.size=0)+
geom_point(size=1.5, shape=21,position = position_jitterdodge(jitter.width=1)) + theme_classic() +
theme(axis.title = element_text(size=16,color='black'),
axis.text = element_text(size=14,color='black'),
axis.text.x = element_text(angle=90,vjust=0.5,hjust=1),
legend.text = element_text(size=14,color='black'),
axis.title.y = element_text(margin = margin(0,2,0,0)),
axis.title.x = element_text(margin = margin(10,0,0,0)),
plot.margin=margin(5,22,0,1),
legend.box.margin=margin(0,0,-10,-10),
legend.position='none') + #ylim(c(0,max(data$count) * 1.1)) +
scale_color_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
scale_fill_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
guides(fill = "none") + labs(color="") + labs(x = 'Stage', y = paste0(goi, ' expression'))
# %%
# Generate boxplots
options(repr.plot.width=2.4, repr.plot.height=4,repr.plot.res=500)
goi <- 'PHYHD1'
data <- plotCounts(deseq_obj, gene = goi, intgroup = c('Group'),returnData = TRUE)
data$Group <- gsub('First_trimester', '1st Tri.',data$Group)
data$Group <- gsub('Second_trimester', '2nd Tri.',data$Group)
data$Group <- gsub('Third_trimester', '3rd Tri.',data$Group)
data$Group <- gsub('Infancy', 'Infant', data$Group)
data$Group <- gsub('Adolescence', 'Adoles.', data$Group)
data$Group <- factor(data$Group, levels = c('1st Tri.','2nd Tri.',
'3rd Tri.','Infant','Adoles.'))
ggplot(data, aes(x=Group, y=count, color=Group,fill=Group)) +
geom_boxplot(alpha=0.2,outlier.size=0)+
geom_point(size=1.5, shape=21,position = position_jitterdodge(jitter.width=1)) + theme_classic() +
theme(axis.title = element_text(size=16,color='black'),
axis.text = element_text(size=14,color='black'),
axis.text.x = element_text(angle=90,vjust=0.5,hjust=1),
legend.text = element_text(size=14,color='black'),
axis.title.y = element_text(margin = margin(0,2,0,0)),
axis.title.x = element_text(margin = margin(10,0,0,0)),
plot.margin=margin(5,22,0,1),
legend.box.margin=margin(0,0,-10,-10),
legend.position='none') + #ylim(c(0,max(data$count) * 1.1)) +
scale_color_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
scale_fill_manual(values = c(colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[1],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[8],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[12],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[17],
colorRampPalette(paletteer_d("PNWColors::Sunset2"))(24)[23])) +
guides(fill = "none") + labs(color="") + labs(x = 'Stage', y = paste0(goi, ' expression'))
# %%
# Change group labels for plotting
human_multiome_astro_v2$Group_plotting <- ifelse(human_multiome_astro_v2$Group == 'First_trimester','1st Tri.',
ifelse(human_multiome_astro_v2$Group == 'Second_trimester','2nd Tri.',
ifelse(human_multiome_astro_v2$Group == 'Third_trimester','3rd Tri.',
ifelse(human_multiome_astro_v2$Group == 'Infancy','Infant','Adoles.'))))
human_multiome_astro_v2$Group_plotting <- factor(human_multiome_astro_v2$Group_plotting,
levels = c('1st Tri.',
'2nd Tri.',
'3rd Tri.',
'Infant',
'Adoles.'))
# %%
# Generate ETNPPL track
options(repr.plot.width = 3, repr.plot.height = 4.4, repr.plot.res = 700)
roi <- 'chr4-108739887-108745637'
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
goi <- 'ETNPPL'
human_multiome_astro_v2 <- LinkPeaks(human_multiome_astro_v2, peak.assay = 'peaks', expression.assay = 'RNA', genes.use = goi,
score_cutoff = 0.01)
cov_plot <- CoveragePlot(
object = human_multiome_astro_v2,
region = roi,
assay='peaks',
annotation = FALSE,
peaks = FALSE,
links =FALSE,
group.by = 'Group_plotting'
) + theme(axis.line.y = element_blank(),
text = element_text(size=18),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.label.y = element_text(hjust=1),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
DefaultAssay(human_multiome_astro_v2) <-'peaks'
gene_plot <- AnnotationPlot(
object = human_multiome_astro_v2,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
peak_plot <- PeakPlot(
object = human_multiome_astro_v2,
region = roi
) + theme(axis.text.x = element_text(color = 'black', size = 7),
axis.line.y = element_blank())
df <- as.data.frame(Links(human_multiome_astro_v2))
poi <- paste0(peak_plot$data$seqnames,'-',
peak_plot$data$start,'-',
peak_plot$data$end)
peak_plot$data$PCC <- df[df$peak %in% poi,'score']
peak_plot <- peak_plot + aes(color = peak_plot$data$PCC) +
theme(legend.position='left',
axis.line.y = element_blank(),
legend.text=element_text(size=12),
legend.title=element_text(size=14,margin=margin(0,0,8,0)),
legend.key.height = unit(0.1, "cm"),
legend.margin = margin(20,-110,0,0),
legend.background = element_rect(fill = "transparent", color = NA),
legend.box.background = element_rect(fill = "transparent", color = NA),
legend.key.width = unit(0.2, "cm"),
plot.margin=margin(0,0,0,0)) + labs(color = "PCC") + ylab('') +
scale_color_gradientn(colors = brewer.pal(1000, 'Reds'),
limits=c(0,ceiling(max(peak_plot$data$PCC)*100)/100 ),
breaks = seq(0, ceiling(max(peak_plot$data$PCC)*100)/100,
by = ceiling((max(peak_plot$data$PCC))*100)/100)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
CombineTracks(
plotlist = list(cov_plot, peak_plot,gene_plot),
heights = c(16, 0.5,3),
widths = c(10, 2)) &
scale_fill_manual(values = paletteer_d("PNWColors::Sunset2")) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate GJB6 track for supplement
options(repr.plot.width = 3, repr.plot.height = 4.4, repr.plot.res = 700)
roi <- regionize('chr13:20228936-20233410')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
goi <- 'GJB6'
human_multiome_astro_v2 <- LinkPeaks(human_multiome_astro_v2, peak.assay = 'peaks', expression.assay = 'RNA', genes.use = goi,
score_cutoff = 0.01)
cov_plot <- CoveragePlot(
object = human_multiome_astro_v2,
region = roi,
assay='peaks',
annotation = FALSE,
peaks = FALSE,
links =FALSE,
group.by = 'Group_plotting'
) + theme(axis.line.y = element_blank(),
text = element_text(size=18),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.label.y = element_text(hjust=1),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
DefaultAssay(human_multiome_astro_v2) <-'peaks'
gene_plot <- AnnotationPlot(
object = human_multiome_astro_v2,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
peak_plot <- PeakPlot(
object = human_multiome_astro_v2,
region = roi
) + theme(axis.text.x = element_text(color = 'black', size = 7),
axis.line.y = element_blank())
df <- as.data.frame(Links(human_multiome_astro_v2))
poi <- paste0(peak_plot$data$seqnames,'-',
peak_plot$data$start,'-',
peak_plot$data$end)
peak_plot$data$PCC <- df[df$peak %in% poi,'score']
peak_plot <- peak_plot + aes(color = peak_plot$data$PCC) +
theme(legend.position='left',
axis.line.y = element_blank(),
legend.text=element_text(size=12),
legend.title=element_text(size=14,margin=margin(0,0,8,0)),
legend.key.height = unit(0.1, "cm"),
legend.margin = margin(20,-110,0,0),
legend.background = element_rect(fill = "transparent", color = NA),
legend.box.background = element_rect(fill = "transparent", color = NA),
legend.key.width = unit(0.2, "cm"),
plot.margin=margin(0,0,0,0)) + labs(color = "PCC") +ylab('') +
scale_color_gradientn(colors = brewer.pal(1000, 'Reds'),
limits=c(0,ceiling(max(peak_plot$data$PCC)*100)/100 ),
breaks = seq(0, ceiling(max(peak_plot$data$PCC)*100)/100,
by = ceiling((max(peak_plot$data$PCC))*100)/100)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
CombineTracks(
plotlist = list(cov_plot, peak_plot,gene_plot),
heights = c(16, 0.5,3),
widths = c(10, 2)) &
scale_fill_manual(values = paletteer_d("PNWColors::Sunset2")) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate GJB6 track for supplement
options(repr.plot.width = 3, repr.plot.height = 4.4, repr.plot.res = 700)
roi <- regionize('chr9:128917629-128922835')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
goi <- 'LRRC8A'
human_multiome_astro_v2 <- LinkPeaks(human_multiome_astro_v2, peak.assay = 'peaks', expression.assay = 'RNA', genes.use = goi,
score_cutoff = 0,pvalue_cutoff = 1)
# Collect peaks for plotting
roi_gr <- GRanges(seqnames = str_split(roi,'-')[[1]][1],
IRanges(start = as.numeric(str_split(roi,'-')[[1]][2]),
end = as.numeric(str_split(roi,'-')[[1]][3])))
links_gr <- GRanges(seqnames = do.call(rbind,str_split(Links(human_multiome_astro_v2)$peak,'-'))[,1],
IRanges(start = as.numeric(do.call(rbind,str_split(Links(human_multiome_astro_v2)$peak,'-'))[,2]),
end = as.numeric(do.call(rbind,str_split(Links(human_multiome_astro_v2)$peak,'-'))[,3])))
cov_plot <- CoveragePlot(
object = human_multiome_astro_v2,
region = roi,
assay='peaks',
annotation = FALSE,
peaks = FALSE,
links =FALSE,
group.by = 'Group_plotting'
) + theme(axis.line.y = element_blank(),
text = element_text(size=18),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.label.y = element_text(hjust=1),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
DefaultAssay(human_multiome_astro_v2) <-'peaks'
gene_plot <- AnnotationPlot(
object = human_multiome_astro_v2,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 4 # Increase gene font size
peak_plot <- PeakPlot(
object = human_multiome_astro_v2,
region = roi
) + theme(axis.text.x = element_text(color = 'black', size = 7),
axis.line.y = element_blank())
df <- as.data.frame(Links(human_multiome_astro_v2))
poi <- paste0(peak_plot$data$seqnames,'-',
peak_plot$data$start,'-',
peak_plot$data$end)
peak_plot$data$PCC <- df[queryHits(findOverlaps(links_gr, roi_gr)),'score']
peak_plot <- peak_plot + aes(color = peak_plot$data$PCC) +
theme(legend.position='left',
axis.line.y = element_blank(),
legend.text=element_text(size=12),
legend.title=element_text(size=14,margin=margin(0,0,8,0)),
legend.key.height = unit(0.1, "cm"),
legend.margin = margin(20,-110,0,0),
legend.background = element_rect(fill = "transparent", color = NA),
legend.box.background = element_rect(fill = "transparent", color = NA),
legend.key.width = unit(0.2, "cm"),
plot.margin=margin(0,0,0,0)) + labs(color = "PCC") +ylab('') +
scale_color_gradientn(colors = brewer.pal(1000, 'Reds'),
limits=c(floor(min(peak_plot$data$PCC)*100)/100,ceiling(max(peak_plot$data$PCC)*100)/100 ),
breaks = seq(floor(min(peak_plot$data$PCC)*100)/100, ceiling(max(peak_plot$data$PCC)*100)/100,
by = ceiling(max(peak_plot$data$PCC)*100)/100 - floor(min(peak_plot$data$PCC)*100)/100)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
CombineTracks(
plotlist = list(cov_plot, peak_plot,gene_plot),
heights = c(16, 0.5,3),
widths = c(10, 2)) &
scale_fill_manual(values = paletteer_d("PNWColors::Sunset2")) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate CSRP1 track
options(repr.plot.width = 3, repr.plot.height = 4.4, repr.plot.res = 700)
roi <- regionize('chr1:201498321-201503221')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
goi <- 'CSRP1'
human_multiome_astro_v2 <- LinkPeaks(human_multiome_astro_v2, peak.assay = 'peaks', expression.assay = 'RNA', genes.use = goi,
score_cutoff = 0,pvalue_cutoff = 1)
cov_plot <- CoveragePlot(
object = human_multiome_astro_v2,
region = roi,
assay='peaks',
annotation = FALSE,
peaks = FALSE,
links =FALSE,
group.by = 'Group_plotting'
) + theme(axis.line.y = element_blank(),
text = element_text(size=18),
axis.ticks.y = element_blank(),
#axis.title.y = element_blank(),
axis.label.y = element_text(hjust=1),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
DefaultAssay(human_multiome_astro_v2) <-'peaks'
gene_plot <- AnnotationPlot(
object = human_multiome_astro_v2,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
peak_plot <- PeakPlot(
object = human_multiome_astro_v2,
region = roi
) + theme(axis.text.x = element_text(color = 'black', size = 7),
axis.line.y = element_blank())
df <- as.data.frame(Links(human_multiome_astro_v2))
poi <- paste0(peak_plot$data$seqnames,'-',
peak_plot$data$start,'-',
peak_plot$data$end)
peak_plot$data$PCC <- df[df$peak %in% poi,'score']
peak_plot <- peak_plot + aes(color = peak_plot$data$PCC) +
theme(legend.position='left',
axis.line.y = element_blank(),
legend.text=element_text(size=12),
legend.title=element_text(size=14,margin=margin(0,0,8,0)),
legend.key.height = unit(0.1, "cm"),
legend.margin = margin(20,-110,0,0),
legend.background = element_rect(fill = "transparent", color = NA),
legend.box.background = element_rect(fill = "transparent", color = NA),
legend.key.width = unit(0.2, "cm"),
plot.margin=margin(0,0,0,0)) + labs(color = "PCC") +ylab('') +
scale_color_gradientn(colors = brewer.pal(1000, 'Reds'),
limits=c(floor(min(peak_plot$data$PCC)*100)/100,ceiling(max(peak_plot$data$PCC)*100)/100 ),
breaks = seq(floor(min(peak_plot$data$PCC)*100)/100, ceiling(max(peak_plot$data$PCC)*100)/100,
by = ceiling(max(peak_plot$data$PCC)*100)/100 - floor(min(peak_plot$data$PCC)*100)/100)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
CombineTracks(
plotlist = list(cov_plot, peak_plot,gene_plot),
heights = c(16, 0.5,3),
widths = c(10, 2)) &
scale_fill_manual(values = paletteer_d("PNWColors::Sunset2")) & theme(plot.margin = margin(0,15,0,0))
# %%
# Load in mouse multiome data
final_obj_motif_chromVAR <- readRDS('/n/groups/neuroduo/Bruno/RDS_files/230810_final_obj_motif_chromVAR.rds')
# %%
# Generate mouse Etnppl track
options(repr.plot.width = 2.7, repr.plot.height = 4.4, repr.plot.res = 500)
roi <- 'chr3-130632373-130635599'
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
cov_plot <- BigwigTrack(
bigwig = list('DR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_DR_GR_merge.bw',
'NR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_NR_GR_merge.bw',
'GR\nCon' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/delCre_pos_GR_merge.bw',
'GR\nKO' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/Cre_pos_GR_merge.bw',
'ATAC\nGR-Con' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_deltaCre_GR.merge.bw',
'ATAC\nGR-KO' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_Cre_GR.merge.bw'),
region = roi,
smooth=5) + theme(legend.position='none',
axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
text = element_text(size=18),
axis.label.y = element_text(hjust=1),
axis.line.y = element_blank(),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) )) +
scale_fill_manual(values = c("#675478","#D49E78",
rep(c('#EBAE00','#019CB5'),2)))
gene_plot <- AnnotationPlot(
object = final_obj_motif_chromVAR,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
CombineTracks(
plotlist = list(cov_plot,gene_plot),
heights = c(16, 3)) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate mouse Gjb6 track
options(repr.plot.width = 2.7, repr.plot.height = 4.4, repr.plot.res = 900)
roi <- regionize('chr14:57130226-57134463')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
cov_plot <- BigwigTrack(
bigwig = list('DR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_DR_GR_merge.bw',
'NR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_NR_GR_merge.bw',
'GR\nCon' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/delCre_pos_GR_merge.bw',
'GR\nKO' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/Cre_pos_GR_merge.bw',
'ATAC\nGR-Con' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_deltaCre_GR.merge.bw',
'ATAC\nGR-KO' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_Cre_GR.merge.bw'),
region = roi,
smooth=10) + theme(legend.position='none',
axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
text = element_text(size=18),
axis.label.y = element_text(hjust=1),
axis.line.y = element_blank(),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) )) +
scale_fill_manual(values = c("#675478","#D49E78",
rep(c('#EBAE00','#019CB5'),2)))
gene_plot <- AnnotationPlot(
object = final_obj_motif_chromVAR,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
CombineTracks(
plotlist = list(cov_plot,gene_plot),
heights = c(16, 3)) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate mouse Csrp1 track
options(repr.plot.width = 2.7, repr.plot.height = 4.4, repr.plot.res = 900)
roi <- regionize('chr1:135732451-135737114')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
cov_plot <- BigwigTrack(
bigwig = list('DR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_DR_GR_merge.bw',
'NR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_NR_GR_merge.bw',
'GR\nCon' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/delCre_pos_GR_merge.bw',
'GR\nKO' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/Cre_pos_GR_merge.bw',
'ATAC\nGR-Con' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_deltaCre_GR.merge.bw',
'ATAC\nGR-KO' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_Cre_GR.merge.bw'),
region = roi,
smooth=10) + theme(legend.position='none',
axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
text = element_text(size=18),
axis.label.y = element_text(hjust=1),
axis.line.y = element_blank(),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) )) +
scale_fill_manual(values = c("#675478","#D49E78",
rep(c('#EBAE00','#019CB5'),2)))
gene_plot <- AnnotationPlot(
object = final_obj_motif_chromVAR,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 5 # Increase gene font size
CombineTracks(
plotlist = list(cov_plot,gene_plot),
heights = c(16, 3)) & theme(plot.margin = margin(0,15,0,0))
# %%
# Generate mouse Lrrc8a/Phyhd1 tracks
options(repr.plot.width = 2.7, repr.plot.height = 4.4, repr.plot.res = 900)
roi <- regionize('chr2:30,263,378-30,268,744')
min <- floor(as.numeric(str_split(roi,'-')[[1]][2]) / 100) * 100
max <- floor(as.numeric(str_split(roi,'-')[[1]][3]) / 100) * 100
cov_plot <- BigwigTrack(
bigwig = list('DR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_DR_GR_merge.bw',
'NR\nGR' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_DR_NR/Astro_NR_GR_merge.bw',
'GR\nCon' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/delCre_pos_GR_merge.bw',
'GR\nKO' = '/n/groups/neuroduo/Bruno/Cell_type_GR_CnR_v2/Astro_cre_delta/Cre_pos_GR_merge.bw',
'ATAC\nGR-Con' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_deltaCre_GR.merge.bw',
'ATAC\nGR-KO' = '/n/groups/neuroduo/Bruno/Astro_GR_KO_ATAC/BAMs/ATAC_Cre_GR.merge.bw'),
region = roi,
smooth=2) + theme(legend.position='none',
axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
text = element_text(size=18),
axis.label.y = element_text(hjust=1),
axis.line.y = element_blank(),
plot.title = element_text(hjust=1)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) )) +
scale_fill_manual(values = c("#675478","#D49E78",
rep(c('#EBAE00','#019CB5'),2)))
Annotation(final_obj_motif_chromVAR)[Annotation(final_obj_motif_chromVAR)$gene_name == 'Gm28035'] <- NULL
gene_plot <- AnnotationPlot(
object = final_obj_motif_chromVAR,
region = roi
) + scale_color_manual(values='black') + ylab('') +
theme(
axis.text = element_text(size=12,color='black'),
axis.title = element_text(size=14,color='black'),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)) +
scale_x_continuous(limits = c(min,max), breaks = seq(min, max, by = (max - min) ))
gene_plot$layers[[4]]$aes_params$size <- 4.1 # Increase gene font size
CombineTracks(
plotlist = list(cov_plot,gene_plot),
heights = c(16, 3)) & theme(plot.margin = margin(0,15,0,0))

Human_brain_development_multiome_analysis-final.ipynb at commit 2536cae, no license · at the source

Overview

Authors: Bruno Gegenhuber¹, Takuma Sonoda^1,2, Lisa Traunmüller¹, Christopher P. Davis¹, Shon A. Koren¹, Eric C. Griffith¹, Chinfei Chen^1,2, Michael E. Greenberg¹

ORCID iDs: Lisa Traunmüller, Eric C. Griffith, Chinfei Chen, Michael E. Greenberg

Institutions: Harvard University (United States); Boston Children's Hospital (United States)

Journal: Nature, volume 655, issue 8125, pages 1233-1241

Dates: received 9 June 2025; accepted 8 April 2026; published online 20 May 2026; in print 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41586-026-10512-9 · PMID 42162428 · PMCID PMC13421323 · OpenAlex W7161786132

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), mouse (organism), cellular / molecular (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions, Evoked potentials, Connectivity, fMRI & imaging, Single-unit activity, calcium imaging

Keywords: Molecular neuroscience, Astrocyte, Glial development, Epigenetics and plasticity, Visual system

MeSH: Astrocytes*, Neuronal Plasticity*, Receptors, Glucocorticoid*, Signal Transduction*, Animals, Chromatin, Female, Gene Expression Regulation, Developmental, Glucocorticoids, Humans, Male, Mice, Visual Cortex (* major topic)

Topic: Neurogenesis and neuroplasticity mechanisms (Developmental Neuroscience, Neuroscience), according to OpenAlex

Funding: NINDS NIH HHS (F32 NS134623, T32 NS007473, R35 NS143029, F32 NS112455)

Citations: cited by 2 papers (Europe PMC); 97 references in the paper

Abstract

Sensory experience refines neural circuits during critical periods of postnatal development1–3. Although neuronal activity is known to orchestrate the circuit wiring that underlies this process4,5, the environmental cues that restrain developmental plasticity as animals mature are less clear. Here we examine the experience-dependent maturation of the mouse primary visual cortex across postnatal development using paired single-cell transcriptomic and chromatin accessibility sequencing. In addition to identifying the activity-dependent gene programs that emerge within each cortical cell type, we find that light exposure drives astrocyte maturation through cell-type-specific recruitment of the glucocorticoid receptor (encoded by Nr3c1) to chromatin. Astrocyte glucocorticoid receptor signalling activates an extensive gene regulatory program that is partially conserved in human brain development and promotes maturation processes that may regulate critical period closure. Collectively, these findings reveal that astrocyte glucocorticoid receptor signalling restricts neuronal plasticity. Glucocorticoid regulation of astrocyte maturation may also contribute to the effects of early-life stress across the brain, and the disruption of this process may increase susceptibility to neuropsychiatric disease.

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

Repositories

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

masai1116/SHARE-seq-alignmentV2

License: GPL-3.0

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 68f80379f4bf64eff9b4f1962b864a9330456a34, 20 April 2023

Languages: R (6), Python (4), Shell (1)

Size: 55 files, 11 scripts

Software Heritage: not archived

Found in: the text, “SHARE-seq data processing and analysis”

Holds: README, license file

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

Tools: Biopython (2 files), NumPy (2 files), tidyverse (2 files), BEDTools (1 file), data.table (1 file), Matplotlib (1 file), pandas (1 file), pysam (1 file), reshape2 (1 file), SAMtools (1 file), SciPy (1 file), STAR (1 file), Subread (featureCounts) (1 file)

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

28 September 2026: the link answers

13 files

Read_distribution.R, R, 30 lines
Share_seqV2_example.sh, Shell, 903 lines, 2 matches
UMI_gene_perCell_plot_v3.R, R, 146 lines
fastq.process.py3.v0.8.py, Python, 1,219 lines
fastq.process.py3.v0.9.py, Python, 1,227 lines
lib_size_sc_V3_bulk.R, R, 40 lines
lib_size_sc_V5_single_species.R, R, 106 lines
lib_size_sc_V5_species_mixing.R, R, 197 lines
pySinglesGetAlnStats_sai.py, Python, 206 lines
rm_dup_barcode_UMI_v3.py, Python, 123 lines
sum_reads_v2.R, R, 47 lines
LICENSE.md, License, 675 lines
README.md, Text, 112 lines

At the source: github.com/masai1116/SHARE-seq-alignmentV2

broadinstitute.github.io/picard

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: the text, “CUT&RUN data processing and analysis”

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

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

At the source: broadinstitute.github.io/picard/

brunogegenhuber/GR_gene_reg

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 2536caed625a91bb02fb1a630bd5107543286489, 20 March 2026

Languages: Jupyter (8)

Size: 9 files, 8 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 8 notebooks

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

Tools: DESeq2 (6 files), ggplot2 (6 files), edgeR (5 files), tidyverse (5 files), cowplot (4 files), reshape2 (4 files), circlize (3 files), clusterProfiler (3 files), ComplexHeatmap (3 files), Seurat (3 files), data.table (2 files), ggpubr (2 files), pheatmap (2 files)

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

28 September 2026: the link answers

6 files

Astro_GR_KO_ATAC-final.ipynb, Jupyter, 125 lines
Astro_GR_KO_bulk_RNA_analysis-final.ipynb, Jupyter, 919 lines
Astro_GR_KO_snRNA-final.ipynb, Jupyter, 1,271 lines, 4 matches
Astrocyte_mouse_atlas_snRNAseq-final.ipynb, Jupyter, 335 lines, 1 match
GR_NFIA_CUTnRUN_analysis-final.ipynb, Jupyter, 1,453 lines, 1 match
Human_brain_development_multiome_analysis-final.ipynb, Jupyter, 2,255 lines, 5 matches
repository limit reached (2,000 files or 30 MB): the rest is at the source (3 files)

At the source: github.com/brunogegenhuber/GR_gene_reg

Code availability

Custom scripts can be found at https://github.com/brunogegenhuber/GR_gene_reg.

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

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:

3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
17 scripts, each with its path and the digest of its content;
13 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

doi:10.5061/dryad.2280gb612, at Dryad; found in “Data availability”
figshare:29971903, at figshare; found in “Data availability”
geo:GSE306265, at NCBI GEO; found in “Data availability”

Data availability

All sequencing data generated in this study have been deposited in GEO (GSE306265 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE306265)). MERFISH data are available at 10.6084/m9.figshare.29971903 (ref. 95). The following publicly available datasets were also analysed: mouse kidney GR ChIP–seq (GSE115368 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE115368)), mouse liver GR ChIP–seq (GSE72087 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72087)), mouse primary brown preadipocyte GR ChIP–seq (GSE76619 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76619)), mouse embryonic fibroblast GR ChIP–seq (GSE69947 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69947)), mouse mammary gland GR ChIP–seq (GSE74826 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74826)), mouse primary bone marrow-derived macrophage GR ChIP–seq (GSE99887 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99887)), human brain development single-nucleus multiome-seq (10.5061/dryad.2280gb612 (ref. 96)), and mouse astrocyte development snRNA-seq (https://singlecell.broadinstitute.org/single_cell/study/SCP2719/a-multi-region-transcriptomic-atlas-of-developmental-cell-type-diversity-in-mouse-brain). The following publicly available databases were used for data processing: NCBI mm10 genome and RefSeq gene annotation (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001635.20), NCBI hg38 genome and RefSeq gene annotation (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.26), and UCSC mm10 refGene annotation (http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/genes). Source data are provided with this paper.

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 2, 28 September 2026

Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 5 keywords, 13 MeSH terms, 1 funder, 96 references.

Cite

This paper

Gegenhuber, B., Sonoda, T., Traunmüller, L., Davis, C. P., Koren, S. A., Griffith, E. C., Chen, C., & Greenberg, M. E. (2026). Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity. Nature, 655(8125), 1233-1241. https://doi.org/10.1038/s41586-026-10512-9

BibTeX

@article{gegenhuber2026astrocyte, author = {Gegenhuber, Bruno and Sonoda, Takuma and Traunmüller, Lisa and Davis, Christopher P. and Koren, Shon A. and Griffith, Eric C. and Chen, Chinfei and Greenberg, Michael E.}, title = {{Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity}}, journal = {Nature}, year = {2026}, month = may, volume = {655}, number = {8125}, pages = {1233--1241}, publisher = {Nature Portfolio}, issn = {0028-0836}, doi = {10.1038/s41586-026-10512-9}, url = {https://doi.org/10.1038/s41586-026-10512-9}, pmid = {42162428}, pmcid = {PMC13421323} }

RIS

TY - JOUR AU - Gegenhuber, Bruno AU - Sonoda, Takuma AU - Traunmüller, Lisa AU - Davis, Christopher P. AU - Koren, Shon A. AU - Griffith, Eric C. AU - Chen, Chinfei AU - Greenberg, Michael E. TI - Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity T2 - Nature J2 - Nature PY - 2026 DA - 2026/05/20 VL - 655 IS - 8125 SP - 1233 EP - 1241 SN - 0028-0836 PB - Nature Portfolio DO - 10.1038/s41586-026-10512-9 UR - https://doi.org/10.1038/s41586-026-10512-9 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41586-026-10512-9", "type": "article-journal", "title": "Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity", "container-title": "Nature", "author": [ { "family": "Gegenhuber", "given": "Bruno" }, { "family": "Sonoda", "given": "Takuma" }, { "family": "Traunmüller", "given": "Lisa" }, { "family": "Davis", "given": "Christopher P." }, { "family": "Koren", "given": "Shon A." }, { "family": "Griffith", "given": "Eric C." }, { "family": "Chen", "given": "Chinfei" }, { "family": "Greenberg", "given": "Michael E." } ], "container-title-short": "Nature", "volume": "655", "issue": "8125", "page": "1233-1241", "DOI": "10.1038/s41586-026-10512-9", "PMID": "42162428", "PMCID": "PMC13421323", "ISSN": "0028-0836", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41586-026-10512-9", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 20 ] ] } }

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.celrep.2026.117073 [code]: Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Subread (featureCounts), pysam, BEDTools, 16 other tools, mouse, cellular / molecular, 6 references
[2] doi:10.21203/rs.3.rs-9927928/v1 [code]: Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: Subread (featureCounts), STAR, pysam, 15 other tools, mouse, 6 references
[3] doi:10.1126/sciadv.aed2952 [code]: Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: Subread (featureCounts), STAR, pysam, 16 other tools, cellular / molecular, 3 references
[4] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pysam, edgeR, circlize, 15 other tools, cellular / molecular, 6 references
[5] doi:10.7554/elife.107393 [code]: Chromosome-scale genome assembly of the European common cuttlefish <i>Sepia officinalis</i>.
Journal: eLife
In common: Subread (featureCounts), STAR, pysam, 15 other tools, cellular / molecular, 4 references
[6] doi:10.1186/s13059-026-04177-w [code]: Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: pysam, BEDTools, SAMtools, 13 other tools, mouse, cellular / molecular, 7 references
[7] doi:10.1038/s41467-026-69944-6 [code]: Multi-modal dissection of cell-type specific TDP-43 pathology in the motor cortex.
Journal: Nature communications
In common: pysam, BEDTools, SAMtools, 15 other tools, 5 references
[8] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: Subread (featureCounts), STAR, SAMtools, 16 other tools, 3 references
[9] doi:10.1038/s41467-026-71803-3 [code]: Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: BEDTools, SAMtools, DESeq2, 11 other tools, mouse, cellular / molecular, 9 references
[10] doi:10.1038/s41467-026-76675-1 [code]: Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: STAR, pysam, Biopython, 16 other tools, 1 reference

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.

Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

masai1116/SHARE-seq-alignmentV2

broadinstitute.github.io/picard

brunogegenhuber/GR_gene_reg

Code availability

Tracing map

Data

Datasets cited

Data availability

Versions

Version 2, 28 September 2026

Version 1, 28 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity