Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis.

Code ↔ Paper

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

[1] § Methods › Multiplexed error-robust in situ hybridization (MERFISH) ↔ 04_WayneEtAl_MERFISH-2.R, lines 283–338 · score 1.00 · Ap2m1, Baiap2l1, C1qa, C1qbp, C1qc, C3ar1
[2] § Methods › Multiplexed error-robust in situ hybridization (MERFISH) ↔ 01_WayneEtAl_scRNAseq.R, lines 1116–1181 · score 1.00 · Csf1r, Il1a, Il1rn, Il2rg, Tnfsf13b, Csf2ra
[3] § Methods › Quantification and statistical analysis › Analysis of scRNA-seq data ↔ 01_WayneEtAl_scRNAseq.R, lines 1116–1181 · score 0.99 · Oas1a, H2 Ab1, H2 Eb1, H2 K1, H2 D1, Cx3cr1
[4] § Methods › Quantification and statistical analysis › Analysis of scRNA-seq data ↔ 01_WayneEtAl_scRNAseq.R, lines 323–399 · score 0.99 · Atp13a5, Cd209a, Col1a1, Klrb1c, dendritic cells, Ly6g
[5] § Methods › Quantification and statistical analysis › Analysis of MERFISH data ↔ 04_WayneEtAl_MERFISH-2.R, lines 223–281 · score 0.92 · Cd3e, olfactory ensheathing cells, Cd8a, P2ry12, Cemip, Mmp9
[6] § Methods › Quantification and statistical analysis › Analysis of scRNA-seq data ↔ 03_WayneEtAl_MERFISH-1.R, lines 166–216 · score 0.76 · linear dimensional reduction, batch correction, variable features, endothelial cells, PCA, resolution
[7] § Results › Microglia upregulate inflammatory and antigen-presentation programs after GAS infections ↔ 01_WayneEtAl_scRNAseq.R, lines 978–1029 · score 0.71 · srMG2, hMG1, srMG1, srMG4, Ccl3, Ccl4
[8] § Results › GAS infection alters predicted cell–cell communication networks in the OB ↔ 05_WayneEtAl_CellChat.R, lines 52–100 · score 0.69 · ligand receptor interactions, CellChat, Inferred, networks, signaling, pathways
[9] § Results › BECs upregulate inflammatory programs and downregulate BBB-related transcripts after GAS infections ↔ 01_WayneEtAl_scRNAseq.R, lines 732–793 · score 0.68 · Mfsd2a, Itm2a, ECM, Cav1, Cavin1, junctions
[10] § Methods › Quantification and statistical analysis › Analysis of MERFISH data ↔ 03_WayneEtAl_MERFISH-1.R, lines 218–256 · score 0.61 · python script, Nearest neighbor, Endothelial cell, distance, OB, MERFISH
[11] § Results › Streptococcus-responsive microglia are enriched in the glomerular OB layer in close proximity to T cells ↔ 04_WayneEtAl_MERFISH-2.R, lines 421–486 · score 0.60 · granular layer, P2ry12, gene expression, Axl, Gpr34, Ifi30
[12] § Methods › Quantification and statistical analysis › Analysis of scRNA-seq data ↔ 01_WayneEtAl_scRNAseq.R, lines 795–860 · score 0.60 · batch correction, variable features, Elbow, Harmony, UMAP, PCA
[13] § Results › Microglia and BECs show major transcriptional shifts after multiple GAS infections ↔ 02_WayneEtAl_CrossEntropy.R, lines 1–43 · score 0.55 · scRNAseq, cross entropy, olfactory bulb, GAS infection, UMAP, mice
[14] § Results › Microglia and BECs show major transcriptional shifts after multiple GAS infections ↔ 02_WayneEtAl_CrossEntropy.R, lines 244–274 · score 0.55 · cross entropy, olfactory ensheathing cells, OECs, filtering, astrocytes, UMAP

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

R · 1,477 lines · 72 KB · no license · 6 matches

# code for scRNAseq analysis of mouse olfactory bulb after GAS infection
# analysis run with R v 4.4.1 and Seurat v 5.1.0
.rs.restartR() # restart R
rm(list=ls()) # clear environment
# load programs
library(Seurat) # version 5.1.0
library(dplyr) # version 1.1.4
library(tibble) # version 3.2.1
library(ggplot2) # version 3.5.1
library(harmony) # version 1.2.0
library(SeuratExtend) # version 1.0.0
## data processing and OB clustering ----
## importing & merging data, adding meta data & removing dead cell
# import data files & merge
ds004.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW004/raw_feature_bc_matrix")
ds005.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW005/raw_feature_bc_matrix")
ds008.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW008/raw_feature_bc_matrix")
ds009.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW009/raw_feature_bc_matrix")
ds010.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW010/raw_feature_bc_matrix")
ds011.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW011/raw_feature_bc_matrix")
ds012.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW012/raw_feature_bc_matrix")
ds013.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW013/raw_feature_bc_matrix")
ds014.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW014/raw_feature_bc_matrix")
ds015.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW015/raw_feature_bc_matrix")
ds016.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW016/raw_feature_bc_matrix")
ds017.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW017/raw_feature_bc_matrix")
ds018.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW018/raw_feature_bc_matrix")
ds021.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW021/raw_feature_bc_matrix")
ds022.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW022/raw_feature_bc_matrix")
ds023.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW023/raw_feature_bc_matrix")
ds024.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW024/raw_feature_bc_matrix")
ds025.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW025/raw_feature_bc_matrix")
ds026.data <- Read10X(data.dir = "Data/Raw data for Seurat/CW026/raw_feature_bc_matrix")
ds004 <- CreateSeuratObject(counts = ds004.data, project = "PBS", min.cells = 3, min.features = 200)
ds005 <- CreateSeuratObject(counts = ds005.data, project = "GAS", min.cells = 3, min.features = 200)
ds008 <- CreateSeuratObject(counts = ds008.data, project = "GAS", min.cells = 3, min.features = 200)
ds009 <- CreateSeuratObject(counts = ds009.data, project = "GAS", min.cells = 3, min.features = 200)
ds010 <- CreateSeuratObject(counts = ds010.data, project = "RORg", min.cells = 3, min.features = 200)
ds011 <- CreateSeuratObject(counts = ds011.data, project = "GAS", min.cells = 3, min.features = 200) # Setting orig.ident on CW011 as "GAS" but condition as "Isotype_mAb"
ds012 <- CreateSeuratObject(counts = ds012.data, project = "IL-17A_mAb", min.cells = 3, min.features = 200)
ds013 <- CreateSeuratObject(counts = ds013.data, project = "RORg", min.cells = 3, min.features = 200)
ds014 <- CreateSeuratObject(counts = ds014.data, project = "PBS", min.cells = 3, min.features = 200)
ds015 <- CreateSeuratObject(counts = ds015.data, project = "GAS", min.cells = 3, min.features = 200) # Setting orig.ident on CW015 as "GAS" but condition as "Isotype_mAb"
ds016 <- CreateSeuratObject(counts = ds016.data, project = "IL-17A_mAb", min.cells = 3, min.features = 200)
ds017 <- CreateSeuratObject(counts = ds017.data, project = "IL-17A_mAb", min.cells = 3, min.features = 200)
ds018 <- CreateSeuratObject(counts = ds018.data, project = "RORg", min.cells = 3, min.features = 200)
ds021 <- CreateSeuratObject(counts = ds021.data, project = "GAS", min.cells = 3, min.features = 200) # Setting orig.ident on CW021 as "GAS" but condition as "Csf2flox"
ds022 <- CreateSeuratObject(counts = ds022.data, project = "Csf2_CD4ko", min.cells = 3, min.features = 200)
ds023 <- CreateSeuratObject(counts = ds023.data, project = "PBS", min.cells = 3, min.features = 200)
ds024 <- CreateSeuratObject(counts = ds024.data, project = "GAS", min.cells = 3, min.features = 200)
ds025 <- CreateSeuratObject(counts = ds025.data, project = "GAS", min.cells = 3, min.features = 200) # Setting orig.ident on CW025 as "GAS" but condition as "Csf2flox"
ds026 <- CreateSeuratObject(counts = ds026.data, project = "Csf2_CD4ko", min.cells = 3, min.features = 200)
dsmerged <- merge(ds004, y = c(ds005, ds008, ds009, ds010, ds011, ds012, ds013, ds014, ds015, ds016, ds017, ds018, ds021, ds022, ds023, ds024, ds025, ds026), add.cell.ids = c("CW_004", "CW_005", "CW_008", "CW_009", "CW_010", "CW_011", "CW_012", "CW_013", "CW_014", "CW_015", "CW_016", "CW_017", "CW_018", "CW_021", "CW_022", "CW_023", "CW_024", "CW_025", "CW_026"))
rm(ds004, ds004.data, ds005, ds005.data, ds008, ds008.data, ds009, ds009.data, ds010, ds010.data, ds011, ds011.data, ds012, ds012.data, ds013, ds013.data, ds014, ds014.data, ds015, ds015.data, ds016, ds016.data, ds017, ds017.data, ds018, ds018.data, ds021, ds021.data, ds022, ds022.data, ds023, ds023.data, ds024, ds024.data, ds025, ds025.data, ds026, ds026.data)
# add meta data
# add sample to meta data
samplevec <- rep("CW004",nrow([email hidden]))
samplevec[grep("CW_005",rownames([email hidden]))] <- "CW005"
samplevec[grep("CW_008",rownames([email hidden]))] <- "CW008"
samplevec[grep("CW_009",rownames([email hidden]))] <- "CW009"
samplevec[grep("CW_010",rownames([email hidden]))] <- "CW010"
samplevec[grep("CW_011",rownames([email hidden]))] <- "CW011"
samplevec[grep("CW_012",rownames([email hidden]))] <- "CW012"
samplevec[grep("CW_013",rownames([email hidden]))] <- "CW013"
samplevec[grep("CW_014",rownames([email hidden]))] <- "CW014"
samplevec[grep("CW_015",rownames([email hidden]))] <- "CW015"
samplevec[grep("CW_016",rownames([email hidden]))] <- "CW016"
samplevec[grep("CW_017",rownames([email hidden]))] <- "CW017"
samplevec[grep("CW_018",rownames([email hidden]))] <- "CW018"
samplevec[grep("CW_021",rownames([email hidden]))] <- "CW021"
samplevec[grep("CW_022",rownames([email hidden]))] <- "CW022"
samplevec[grep("CW_023",rownames([email hidden]))] <- "CW023"
samplevec[grep("CW_024",rownames([email hidden]))] <- "CW024"
samplevec[grep("CW_025",rownames([email hidden]))] <- "CW025"
samplevec[grep("CW_026",rownames([email hidden]))] <- "CW026"
table(samplevec)
[email hidden]$sample=samplevec
# add batch to meta data
batchvec <- rep("batch A",nrow([email hidden]))
batchvec[grep("CW_008",rownames([email hidden]))] <- "batch B"
batchvec[grep("CW_009",rownames([email hidden]))] <- "batch C"
batchvec[grep("CW_010",rownames([email hidden]))] <- "batch C"
batchvec[grep("CW_011",rownames([email hidden]))] <- "batch D"
batchvec[grep("CW_012",rownames([email hidden]))] <- "batch D"
batchvec[grep("CW_013",rownames([email hidden]))] <- "batch E"
batchvec[grep("CW_014",rownames([email hidden]))] <- "batch F"
batchvec[grep("CW_015",rownames([email hidden]))] <- "batch F"
batchvec[grep("CW_016",rownames([email hidden]))] <- "batch F"
batchvec[grep("CW_017",rownames([email hidden]))] <- "batch G"
batchvec[grep("CW_021",rownames([email hidden]))] <- "batch H"
batchvec[grep("CW_022",rownames([email hidden]))] <- "batch H"
batchvec[grep("CW_023",rownames([email hidden]))] <- "batch I"
batchvec[grep("CW_024",rownames([email hidden]))] <- "batch I"
batchvec[grep("CW_018",rownames([email hidden]))] <- "batch J"
batchvec[grep("CW_025",rownames([email hidden]))] <- "batch K"
batchvec[grep("CW_026",rownames([email hidden]))] <- "batch K"
table(batchvec)
[email hidden]$batch=batchvec
# add enrichment to meta data
enrichmentvec <- rep("cd31_cd11b_FACS",nrow([email hidden]))
enrichmentvec[grep("CW_008",rownames([email hidden]))] <- "live_FACS"
enrichmentvec[grep("CW_023",rownames([email hidden]))] <- "live_FACS"
enrichmentvec[grep("CW_024",rownames([email hidden]))] <- "live_FACS"
table(enrichmentvec)
[email hidden]$enrichment=enrichmentvec
# add inoculate to meta data
inoculatevec <- rep("GAS",nrow([email hidden]))
inoculatevec[grep("CW_004",rownames([email hidden]))] <- "PBS"
inoculatevec[grep("CW_014",rownames([email hidden]))] <- "PBS"
inoculatevec[grep("CW_023",rownames([email hidden]))] <- "PBS"
table(inoculatevec)
[email hidden]$inoculate=inoculatevec
# add tissue to meta data
tissuevec <- rep("OB",nrow([email hidden]))
table(tissuevec)
[email hidden]$tissue=tissuevec
# add genotype to meta data
genotypevec <- rep("WT",nrow([email hidden]))
genotypevec[grep("CW_010",rownames([email hidden]))] <- "RORg"
genotypevec[grep("CW_013",rownames([email hidden]))] <- "RORg"
genotypevec[grep("CW_018",rownames([email hidden]))] <- "RORg"
genotypevec[grep("CW_021",rownames([email hidden]))] <- "Csf2flox"
genotypevec[grep("CW_022",rownames([email hidden]))] <- "Csf2floxCD4Cre"
genotypevec[grep("CW_025",rownames([email hidden]))] <- "Csf2flox"
genotypevec[grep("CW_026",rownames([email hidden]))] <- "Csf2floxCD4Cre"
table(genotypevec)
[email hidden]$genotype=genotypevec
# add drug to meta data
drugvec <- rep("none",nrow([email hidden]))
drugvec[grep("CW_011",rownames([email hidden]))] <- "Isotype_mAb"
drugvec[grep("CW_012",rownames([email hidden]))] <- "IL-17A_mAb"
drugvec[grep("CW_015",rownames([email hidden]))] <- "Isotype_mAb"
drugvec[grep("CW_016",rownames([email hidden]))] <- "IL-17A_mAb"
drugvec[grep("CW_017",rownames([email hidden]))] <- "IL-17A_mAb"
drugvec[grep("CW_021",rownames([email hidden]))] <- "4-OH-tamoxifen"
drugvec[grep("CW_022",rownames([email hidden]))] <- "4-OH-tamoxifen"
drugvec[grep("CW_025",rownames([email hidden]))] <- "4-OH-tamoxifen"
drugvec[grep("CW_026",rownames([email hidden]))] <- "4-OH-tamoxifen"
table(drugvec)
[email hidden]$drug=drugvec
# add condition to meta data
conditionvec <- rep("GAS",nrow([email hidden]))
conditionvec[grep("CW_004",rownames([email hidden]))] <- "PBS"
conditionvec[grep("CW_010",rownames([email hidden]))] <- "RORg"
conditionvec[grep("CW_011",rownames([email hidden]))] <- "Isotype_mAb"
conditionvec[grep("CW_012",rownames([email hidden]))] <- "IL-17A_mAb"
conditionvec[grep("CW_013",rownames([email hidden]))] <- "RORg"
conditionvec[grep("CW_014",rownames([email hidden]))] <- "PBS"
conditionvec[grep("CW_015",rownames([email hidden]))] <- "Isotype_mAb"
conditionvec[grep("CW_016",rownames([email hidden]))] <- "IL-17A_mAb"
conditionvec[grep("CW_017",rownames([email hidden]))] <- "IL-17A_mAb"
conditionvec[grep("CW_018",rownames([email hidden]))] <- "RORg"
conditionvec[grep("CW_021",rownames([email hidden]))] <- "Csf2flox"
conditionvec[grep("CW_022",rownames([email hidden]))] <- "Csf2_CD4ko"
conditionvec[grep("CW_023",rownames([email hidden]))] <- "PBS"
conditionvec[grep("CW_025",rownames([email hidden]))] <- "Csf2flox"
conditionvec[grep("CW_026",rownames([email hidden]))] <- "Csf2_CD4ko"
table(conditionvec)
[email hidden]$condition=conditionvec
# add timepoint to meta data
timepointvec <- rep("18hrs",nrow([email hidden]))
table(timepointvec)
[email hidden]$timepoint=timepointvec
# add sex to meta data
sexvec <- rep("female",nrow([email hidden]))
table(sexvec)
[email hidden]$sex=sexvec
# add strain to meta data
strainvec <- rep("C57BL6J",nrow([email hidden]))
table(strainvec)
[email hidden]$strain=strainvec
# add age to meta data
agevec <- rep("P56",nrow([email hidden])) #CW005, 008, 025 and 026 were also P56
agevec[grep("CW_004",rownames([email hidden]))] <- "P53"
agevec[grep("CW_009",rownames([email hidden]))] <- "P60"
agevec[grep("CW_010",rownames([email hidden]))] <- "P59"
agevec[grep("CW_011",rownames([email hidden]))] <- "P59"
agevec[grep("CW_012",rownames([email hidden]))] <- "P58"
agevec[grep("CW_013",rownames([email hidden]))] <- "P60"
agevec[grep("CW_014",rownames([email hidden]))] <- "P53"
agevec[grep("CW_015",rownames([email hidden]))] <- "P57"
agevec[grep("CW_016",rownames([email hidden]))] <- "P57"
agevec[grep("CW_017",rownames([email hidden]))] <- "P60"
agevec[grep("CW_018",rownames([email hidden]))] <- "P57"
agevec[grep("CW_021",rownames([email hidden]))] <- "P57"
agevec[grep("CW_022",rownames([email hidden]))] <- "P58"
agevec[grep("CW_023",rownames([email hidden]))] <- "P60"
agevec[grep("CW_024",rownames([email hidden]))] <- "P60"
table(agevec)
[email hidden]$age=agevec
# add pooled animal number to meta data
pooledvec <- rep("n=3",nrow([email hidden]))
pooledvec[grep("CW_012",rownames([email hidden]))] <- "n=2"
pooledvec[grep("CW_015",rownames([email hidden]))] <- "n=2"
pooledvec[grep("CW_016",rownames([email hidden]))] <- "n=2"
pooledvec[grep("CW_021",rownames([email hidden]))] <- "n=2"
pooledvec[grep("CW_022",rownames([email hidden]))] <- "n=2"
table(pooledvec)
[email hidden]$pooled=pooledvec
rm(samplevec, batchvec, enrichmentvec, tissuevec, genotypevec, drugvec, inoculatevec, conditionvec, timepointvec, sexvec, strainvec, agevec, pooledvec)
# remove dead cells
dsmerged[["percent.mt"]] <- PercentageFeatureSet(dsmerged, pattern = "^mt-")
Idents(object=dsmerged) <- [email hidden]$sample
VlnPlot(dsmerged, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0)
dsmerged <- subset(dsmerged, subset = nCount_RNA > 1000 & nCount_RNA < 50000 & percent.mt < 20)
VlnPlot(dsmerged, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0)
## clustering and batch correction
# normalize & generate PCA
dsmerged <- NormalizeData(dsmerged)
dsmerged <- FindVariableFeatures(dsmerged, selection.method = "vst", nfeatures = 3000)
top10 <- head(VariableFeatures(dsmerged), 10)
plot1 <- VariableFeaturePlot(dsmerged)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)
plot2
all.genes <- rownames(dsmerged)
dsmerged <- ScaleData(dsmerged, features = all.genes)
dsmerged <- RunPCA(dsmerged, features = VariableFeatures(object = dsmerged))
print(dsmerged[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(dsmerged, dims = 1:5, reduction = "pca")
ElbowPlot(dsmerged,ndims=50)
dsmerged <- FindNeighbors(dsmerged, dims = 1:39)
dsmerged <- FindClusters(dsmerged, resolution = 1)
dsmerged <- RunUMAP(dsmerged, dims = 1:39)
DimPlot(dsmerged, reduction = "umap",label=T)
DimPlot(dsmerged, reduction = "umap",group.by="orig.ident")
dsmerged <- RunTSNE(dsmerged, dims = 1:39)
DimPlot(dsmerged, reduction = "tsne",label=T)
DimPlot(dsmerged, reduction = "tsne",group.by="orig.ident")
rm(plot1, plot2, all.genes, top10)
# run harmony to reduce batch effects
pre.harmony.umap <- DimPlot(dsmerged, reduction = "umap",group.by="batch")
pre.harmony.tsne <- DimPlot(dsmerged, reduction = "tsne",group.by="batch")
dsmerged <- dsmerged %>%
RunHarmony("batch", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(dsmerged, 'harmony')
harmony_embeddings[1:5, 1:5]
dsmerged <- dsmerged %>%
RunUMAP(reduction = "harmony", dims = 1:39) %>%
FindNeighbors(reduction = "harmony", dims = 1:39) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.umap)
DimPlot(dsmerged, reduction = "umap",label=F)
DimPlot(dsmerged, reduction = "umap",label=F, split.by ="batch")
dsmerged <- dsmerged %>%
RunTSNE(reduction = "harmony", dims = 1:39) %>%
FindNeighbors(reduction = "harmony", dims = 1:39) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.tsne)
DimPlot(dsmerged, reduction = "tsne",label=T)
DimPlot(dsmerged, reduction = "tsne",group.by="condition")
saveRDS(dsmerged, file = "02_RDS files/dsmerged_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
rm(pre.harmony.umap, pre.harmony.tsne, harmony_embeddings)
## assign cell type identities
# identify cell types
genes <- VariableFeatures(dsmerged)
toplot <- CalcStats(dsmerged, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
neuroglial.markers <- c("Snap25", "Dcx", "Gfap","Aqp4", "Frzb","Cldn5", "Pecam1", "Pdgfra","Pdgfrb", "Atp13a5", "Col1a1", "Fbln1")
lymphoid.markers <- c("Ptprc", "Cd3e", "Cd4", "Cd8a", "Cd163l1", "Klrb1c", "Cd19", "Ms4a1")
myeloid.markers <- c("Ptprc", "Tmem119", "Aif1", "Itgam", "Ly6g", "Camp", "Ly6c2", "Itgax", "Xcr1", "Cd209a", "Ccr9", "Nudt17")
macrophage.vs.microglia.markers <- c("Ptprc", "Tmem119", "Cx3cr1", "P2ry12", "Itgam", "Aif1", "Hexb", "Plac8", "Mrc1", "Cd163")
DotPlot(dsmerged, features = neuroglial.markers)
# assign cluster identities
Idents(object=dsmerged) <- [email hidden]$RNA_snn_res.1
DimPlot(dsmerged, reduction = "umap",label=T)
new.cluster.ids <- c( "Macrophages", "Microglia", "Microglia", "Endothelial cells", "Olfactory ensheathing cells", "Neutrophils", "Endothelial cells", "CD4 T cells", "Microglia", "Astrocytes", "Mixed signature: Macrophages/Microglia", "Neutrophils", "Macrophages", "Dendritic cells", "Neurons", "Neurons", "Mixed signature: Macrophages/Microglia", "Pericytes", "Mixed signature: Microglia/Dendritic cells", "CD8 T cells & NK cells", "Astrocytes", "Mixed signature: Proliferating cells/CD4 T cells/Macrophages", "Dendritic cells", "Microglia/Neutrophils", "Mixed signature: Pericytes/Endothelial cells", "Macrophages", "B cells", "Olfactory ensheathing cells", "Mixed signature: Neutrophils/Macrophages", "Mixed signature: Proliferating cells/Neutrophils", "Mixed signature: Endothelial cells/Microglia", "Oligos/OPCs", "Mixed signature: Olfactory ensheathing cells/Microglia", "Fibroblasts", "Mixed signature: CD4 T cells/CD8 T cells & NK cells", "Mixed signature: Endothelial cells/Astrocytes", "Mixed signature: Astrocytes/Olfactory ensheathing cells", "Mixed signature: Astrocytes/Microglia", "Mixed signature: Endothelial cells/Neutrophils", "Mixed signature: Endothelial cells/CD4 T cells", "Mixed signature: CD4 T cells/Dendritic cells")
names(new.cluster.ids) <- levels(dsmerged)
dsmerged <- RenameIdents(dsmerged, new.cluster.ids)
DimPlot(dsmerged, reduction = "umap",label=T)+ NoLegend()
dsmerged[["cluster.names_all.cells"]] <- Idents(object = dsmerged) # stash identity
saveRDS(dsmerged, file = "02_RDS files/dsmerged_includes_mixed_signatures_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
# remove mixed signature clusters
dsmerged <- subset(x = dsmerged, idents = c("Neurons", "Astrocytes", "Olfactory ensheathing cells", "Oligos/OPCs", "Microglia", "Endothelial cells", "Pericytes", "Fibroblasts", "Macrophages", "Dendritic cells", "Neutrophils", "CD4 T cells", "CD8 T cells & NK cells", "B cells"),invert = FALSE)
DimPlot(dsmerged, reduction = "umap",label=T)
# re-order data sets for display
Idents(object=dsmerged) <- [email hidden]$cluster.names_all.cells
[email hidden]$cluster.names_all.cells <- factor(x = [email hidden]$cluster.names_all.cells, levels(dsmerged) <- c("Neurons", "Astrocytes", "Microglia", "Olfactory ensheathing cells", "Oligos/OPCs", "Endothelial cells", "Pericytes", "Fibroblasts", "Macrophages", "Neutrophils", "Dendritic cells", "CD4 T cells", "CD8 T cells & NK cells", "B cells"))
DimPlot(dsmerged, reduction = "umap",label=T, group.by = "cluster.names_all.cells") + NoLegend() ##Figure 1B##
DimPlot(dsmerged, reduction = "umap",label=T, group.by = "condition") + NoLegend() ##Figure 1C##
DimPlot(dsmerged, reduction = "umap",label=T, split.by = "enrichment") + NoLegend() ##Extended Data Figure 1D##
DimPlot(dsmerged, reduction = "umap",label=T, group.by = "batch") ##Extended Data Figure 1E##
write.csv(table(dsmerged$sample, dsmerged$cluster.names_all.cells), file = "03_quantification/CellTypeDistrBySample.csv") ##Extended Data Table 1##
saveRDS(dsmerged, file = "02_RDS files/dsmerged_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
rm(genes, toplot, neuroglial.markers, lymphoid.markers, myeloid.markers, macrophage.vs.microglia.markers, new.cluster.ids)
## PBS vs GAS plots and differential expression ----
# subset on WT PBS and GAS samples only
Idents(object=dsmerged) <- [email hidden]$condition
DimPlot(dsmerged, reduction = "umap")
dsmerged_PBS_GAS <- subset(x = dsmerged, idents = c("PBS", "GAS"),invert = FALSE)
saveRDS(dsmerged_PBS_GAS, file = "02_RDS files/dsmerged_004_005_008_009_014_023_024.rds")
Idents(object=dsmerged_PBS_GAS) <- [email hidden]$cluster.names_all.cells
[email hidden]$cluster.names_all.cells <- factor(x = [email hidden]$cluster.names_all.cells, levels(dsmerged_PBS_GAS) <- c("B cells", "CD8 T cells & NK cells", "CD4 T cells", "Dendritic cells", "Neutrophils", "Macrophages", "Fibroblasts", "Pericytes", "Endothelial cells", "Oligos/OPCs", "Olfactory ensheathing cells", "Microglia", "Astrocytes", "Neurons"))
cell.type.markers <- c("Snap25", "Dcx", "Gfap", "Aqp4", "Frzb", "Pdgfra", "Cldn5", "Pecam1", "Pdgfrb", "Atp13a5", "Col1a1", "Fbln1", "Ptprc", "Tmem119", "Aif1", "Itgam", "Ly6g", "Camp", "Itgax", "Xcr1", "Cd209a", "Ccr9", "Cd3e", "Cd4", "Cd8a", "Klrb1c", "Cd19")
DotPlot(dsmerged_PBS_GAS, features = cell.type.markers, dot.min = .05) ##Extended Data Figure 1A##
[email hidden]$cluster.names_all.cells <- factor(x = [email hidden]$cluster.names_all.cells, levels(dsmerged) <- c("Neurons", "Astrocytes", "Microglia", "Olfactory ensheathing cells", "Oligos/OPCs", "Endothelial cells", "Pericytes", "Fibroblasts", "Macrophages", "Neutrophils", "Dendritic cells", "CD4 T cells", "CD8 T cells & NK cells", "B cells"))
write.csv(table(dsmerged_PBS_GAS$cluster.names_all.cells, dsmerged_PBS_GAS$enrichment), file = "03_quantification/CellTypeDistrByEnrichment.csv") ##Extended Data Figure 1C##
## cytokine-responsive by cell type
Idents(object=dsmerged_PBS_GAS) <- [email hidden]$cluster.names_all.cells
[email hidden]$cluster.names_all.cells <- factor(x = [email hidden]$cluster.names_all.cells, levels(dsmerged_PBS_GAS) <- c("Neurons", "Astrocytes", "Microglia", "Olfactory ensheathing cells", "Oligos/OPCs", "Endothelial cells", "Pericytes", "Fibroblasts", "Macrophages", "Neutrophils", "Dendritic cells", "CD4 T cells", "CD8 T cells & NK cells", "B cells"))
cytokine.receptors <- c("Csf2rb", "Csf2ra", "Il17rc", "Il17ra", "Ifngr2", "Ifngr1")
DotPlot(dsmerged_PBS_GAS, features = cytokine.receptors, dot.min = 0.05) + ##Extended Data Figure 3D##
coord_flip() +
theme(
axis.text.x.bottom = element_text(angle = 55, hjust = 1, vjust = 0.5),
axis.text.y = element_text(hjust = 1))
rm(cell.type.markers, cytokine.receptors)
## perform differential expression by mixed effects analysis using MAST
dsmerged_PBS_GAS <- readRDS(file = "02_RDS files/dsmerged_004_005_008_009_014_023_024.rds")
library(openxlsx) # version 4.2.8
library(tibble) # version 3.3.0
cell_types <- c("Neurons", "Astrocytes", "Olfactory ensheathing cells", "Oligos/OPCs", "Microglia", "Endothelial cells", "Pericytes", "Fibroblasts", "Macrophages", "Dendritic cells", "Neutrophils", "CD4 T cells", "CD8 T cells & NK cells", "B cells")
dsmerged_PBS_GAS <- JoinLayers(dsmerged_PBS_GAS)
Idents(object = dsmerged_PBS_GAS) <- [email hidden]$cluster.names_all.cells
compare.conditions <- list(
list(ident1 = "GAS", ident2 = "PBS", result_name = "GASvsPBS")
)
run_DE <- function(object, comparisons_list, test_method = "MAST", include_sd = TRUE,
save_results = FALSE, output_file = NULL, skip_cell_subsetting = FALSE) {
results_list <- list()
if (save_results) {
if (is.null(output_file)) {
stop("Please provide a file name when save_results is TRUE.")
}
wb <- createWorkbook()
}
for (condition in comparisons_list) {
ident1 <- condition$ident1
ident2 <- condition$ident2
result_name <- condition$result_name
if (skip_cell_subsetting) {
cell_types_to_analyze <- "CellType"
} else {
cell_types_to_analyze <- cell_types
}
for (cell_type in cell_types_to_analyze) {
if (skip_cell_subsetting) { # skip subsetting - use the object directly
cell_type_object <- object
cat("Using pre-subset object for:", cell_type, "\n")
} else {
cell_type_object <- subset(x = object, idents = c(cell_type), invert = FALSE)
cat("Subsetting for cell type:", cell_type, "\n")
}
cell_type_object <- JoinLayers(cell_type_object)
Idents(object = cell_type_object) <- [email hidden]$condition
tryCatch({
cell_type_DEGs <- FindMarkers(
object = cell_type_object,
ident.1 = ident1,
ident.2 = ident2,
test.use = test_method
)
result <- rownames_to_column(cell_type_DEGs, var = "gene")
if (include_sd) { # calculate standard deviations if requested
genes_to_analyze <- result$gene
cells_ident1 <- WhichCells(cell_type_object, idents = ident1)
cells_ident2 <- WhichCells(cell_type_object, idents = ident2)
expr_data <- GetAssayData(cell_type_object, assay = "RNA", layer = "data")
mean_ident1 <- apply(expr_data[genes_to_analyze, cells_ident1, drop = FALSE], 1, mean, na.rm = TRUE)
mean_ident2 <- apply(expr_data[genes_to_analyze, cells_ident2, drop = FALSE], 1, mean, na.rm = TRUE)
sd_ident1 <- apply(expr_data[genes_to_analyze, cells_ident1, drop = FALSE], 1, sd, na.rm = TRUE)
sd_ident2 <- apply(expr_data[genes_to_analyze, cells_ident2, drop = FALSE], 1, sd, na.rm = TRUE)
result$mean_group1 <- mean_ident1[result$gene]
result$mean_group2 <- mean_ident2[result$gene]
result$sd_group1 <- sd_ident1[result$gene]
result$sd_group2 <- sd_ident2[result$gene]
result$group1_name <- ident1
result$group2_name <- ident2
}
result_key <- paste(result_name, cell_type, sep = "_")
results_list[[result_key]] <- result
if (save_results) {
sheet_name <- result_key
if (nchar(sheet_name) > 31) {
sheet_name <- substr(sheet_name, 1, 31) # excel sheet names have a 31 character limit
}
addWorksheet(wb, sheetName = sheet_name)
writeData(wb, sheet = sheet_name, x = result)
}
cat("Completed analysis for", cell_type, ":", result_name, "\n")
}, error = function(e) { # proper error handling
cat("Error in analysis for", cell_type, ":", result_name, "\n")
cat("Error message:", e$message, "\n")
})
}
}
if (save_results) {
saveWorkbook(wb, file = output_file, overwrite = TRUE)
cat("Results saved to:", output_file, "\n")
}
return(results_list)
}
DE_MAST_dsmerged <- run_DE(
object = dsmerged_PBS_GAS,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/GASvsPBS_AllCellTypes_MAST_DE.xlsx",
include_sd = FALSE,
skip_cell_subsetting = FALSE
)
## Endothelial cell subclustering and analysis ----
# open merged data file and subset on ECs
dsmerged <- readRDS(file = "02_RDS files/dsmerged_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
Idents(object=dsmerged) <- [email hidden]$cluster.names_all.cells
dsmerged_ECs <- subset(x = dsmerged, idents = c("Endothelial cells"),invert = FALSE)
## clustering, batch correction, removing mixed signature
# normalize & generate PCA
dsmerged_ECs <- NormalizeData(dsmerged_ECs)
dsmerged_ECs <- FindVariableFeatures(dsmerged_ECs, selection.method = "vst", nfeatures = 3000)
top10 <- head(VariableFeatures(dsmerged_ECs), 10)
plot1 <- VariableFeaturePlot(dsmerged_ECs)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)
plot2
all.genes <- rownames(dsmerged_ECs)
dsmerged_ECs <- ScaleData(dsmerged_ECs, features = all.genes)
dsmerged_ECs <- RunPCA(dsmerged_ECs, features = VariableFeatures(object = dsmerged_ECs))
print(dsmerged_ECs[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(dsmerged_ECs, dims = 1:5, reduction = "pca")
ElbowPlot(dsmerged_ECs,ndims=50)
dsmerged_ECs <- FindNeighbors(dsmerged_ECs, dims = 1:38)
dsmerged_ECs <- FindClusters(dsmerged_ECs, resolution = 0.4)
dsmerged_ECs <- RunUMAP(dsmerged_ECs, dims = 1:38)
DimPlot(dsmerged_ECs, reduction = "umap",label=T)
DimPlot(dsmerged_ECs, reduction = "umap",group.by="orig.ident")
dsmerged_ECs <- RunTSNE(dsmerged_ECs, dims = 1:38)
DimPlot(dsmerged_ECs, reduction = "tsne",label=T)
DimPlot(dsmerged_ECs, reduction = "tsne",group.by="orig.ident")
rm(plot1, plot2, all.genes, top10)
# run harmony to reduce batch effects
pre.harmony.umap <- DimPlot(dsmerged_ECs, reduction = "umap",group.by="batch")
pre.harmony.tsne <- DimPlot(dsmerged_ECs, reduction = "tsne",group.by="batch")
dsmerged_ECs <- dsmerged_ECs %>%
RunHarmony("batch", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(dsmerged_ECs, 'harmony')
harmony_embeddings[1:5, 1:5]
dsmerged_ECs <- dsmerged_ECs %>%
RunUMAP(reduction = "harmony", dims = 1:38) %>%
FindNeighbors(reduction = "harmony", dims = 1:38) %>%
FindClusters(resolution = 0.4) %>%
identity()
print(pre.harmony.umap)
DimPlot(dsmerged_ECs, reduction = "umap",label=F)
DimPlot(dsmerged_ECs, reduction = "umap",label=F, group.by ="batch")
dsmerged_ECs <- dsmerged_ECs %>%
RunTSNE(reduction = "harmony", dims = 1:38) %>%
FindNeighbors(reduction = "harmony", dims = 1:38) %>%
FindClusters(resolution = 0.4) %>%
identity()
print(pre.harmony.tsne)
DimPlot(dsmerged_ECs, reduction = "tsne",label=T)
DimPlot(dsmerged_ECs, reduction = "tsne",group.by="condition")
# remove mixed signature clusters
genes <- VariableFeatures(dsmerged_ECs)
toplot <- CalcStats(dsmerged_ECs, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
neuroglial.markers <- c("Snap25", "Dcx", "Gfap","Aqp4", "Frzb","Cldn5", "Pecam1", "Pdgfra","Pdgfrb", "Atp13a5", "Col1a1", "Fbln1")
lymphoid.markers <- c("Ptprc", "Cd3e", "Cd4", "Cd8a", "Cd163l1", "Klrb1c", "Cd19", "Ms4a1")
myeloid.markers <- c("Ptprc", "Tmem119", "Aif1", "Itgam", "Ly6g", "Camp", "Ly6c2", "Itgax", "Xcr1", "Cd209a", "Ccr9", "Nudt17")
DotPlot(dsmerged_ECs, features = neuroglial.markers)
FeaturePlot(dsmerged_ECs, c("nCount_RNA", "nFeature_RNA", "percent.mt")) # cluster 4 is low viability
FeaturePlot(dsmerged_ECs, c("Aqp4", "Gfap", "Cspg5", "Olig1")) # cluster 6 has mixed signature: astrocytes/oligos/endothelial cells
FeaturePlot(dsmerged_ECs, c("Pdgfrb", "Atp13a5")) # cluster 7 has mixed signature: pericytes/endothelial cells
FeaturePlot(dsmerged_ECs, c("Cx3cr1", "Ptprc", "Cd3e")) # cluster 8 has mixed signature: microglia/endothelial cells; also a couple T cells
FeaturePlot(dsmerged_ECs, c("Frzb", "Plp1", "Col1a1", "Fbln1", "Pdgfra")) # cluster 9 has mixed signature: olfactory ensheathing cells/fibroblasts/endothelial cells
FeaturePlot(dsmerged_ECs, c("Top2a")) # cluster 9 has mixed signature: olfactory ensheathing cells/fibroblasts/endothelial cells
Idents(object=dsmerged_ECs) <- [email hidden]$RNA_snn_res.0.4
dsmerged_ECs <- subset(x = dsmerged_ECs, idents = c(4, 6, 7, 8, 9),invert = TRUE)
DimPlot(dsmerged_ECs, reduction = "umap",label=T)
DimPlot(dsmerged_ECs, reduction = "umap", split.by = "inoculate")
saveRDS(dsmerged_ECs, file = "02_RDS files/dsmerged_ECs_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
## subset on PBS and GAS ECs only
Idents(object=dsmerged_ECs) <- [email hidden]$condition
DimPlot(dsmerged_ECs, reduction = "umap")
dsmerged_ECs_PBS_GAS <- subset(x = dsmerged_ECs, idents = c("PBS", "GAS"),invert = FALSE)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap")
# recluster
dsmerged_ECs_PBS_GAS <- NormalizeData(dsmerged_ECs_PBS_GAS)
dsmerged_ECs_PBS_GAS <- FindVariableFeatures(dsmerged_ECs_PBS_GAS, selection.method = "vst", nfeatures = 3000)
top10 <- head(VariableFeatures(dsmerged_ECs_PBS_GAS), 10)
plot1 <- VariableFeaturePlot(dsmerged_ECs_PBS_GAS)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)
plot2
all.genes <- rownames(dsmerged_ECs_PBS_GAS)
dsmerged_ECs_PBS_GAS <- ScaleData(dsmerged_ECs_PBS_GAS, features = all.genes)
dsmerged_ECs_PBS_GAS <- RunPCA(dsmerged_ECs_PBS_GAS, features = VariableFeatures(object = dsmerged_ECs_PBS_GAS))
print(dsmerged_ECs_PBS_GAS[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(dsmerged_ECs_PBS_GAS, dims = 1:5, reduction = "pca")
ElbowPlot(dsmerged_ECs_PBS_GAS,ndims=50)
dsmerged_ECs_PBS_GAS <- FindNeighbors(dsmerged_ECs_PBS_GAS, dims = 1:32)
dsmerged_ECs_PBS_GAS <- FindClusters(dsmerged_ECs_PBS_GAS, resolution = 0.4)
dsmerged_ECs_PBS_GAS <- RunUMAP(dsmerged_ECs_PBS_GAS, dims = 1:32)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap",label=T)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap", group.by="condition") ##Extended Data Figure 1F##
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap", group.by="batch") ##Extended Data Figure 1G##
dsmerged_ECs_PBS_GAS <- RunTSNE(dsmerged_ECs_PBS_GAS, dims = 1:32)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "tsne",label=T)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "tsne",group.by="condition")
rm(plot1, plot2, all.genes, top10)
# run harmony
pre.harmony.umap <- DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap",group.by="batch", label = F)
pre.harmony.tsne <- DimPlot(dsmerged_ECs_PBS_GAS, reduction = "tsne",group.by="batch", label = F)
dsmerged_ECs_PBS_GAS <- dsmerged_ECs_PBS_GAS %>%
RunHarmony("batch", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(dsmerged_ECs_PBS_GAS, 'harmony')
harmony_embeddings[1:5, 1:5]
dsmerged_ECs_PBS_GAS <- dsmerged_ECs_PBS_GAS %>%
RunUMAP(reduction = "harmony", dims = 1:32) %>%
FindNeighbors(reduction = "harmony", dims = 1:32) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.umap)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap",label=T)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "umap",label=F, group.by="inoculate")
dsmerged_ECs_PBS_GAS <- dsmerged_ECs_PBS_GAS %>%
RunTSNE(reduction = "harmony", dims = 1:32) %>%
FindNeighbors(reduction = "harmony", dims = 1:32) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.tsne)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "tsne",label=T)
DimPlot(dsmerged_ECs_PBS_GAS, reduction = "tsne",label=F, group.by="inoculate")
saveRDS(dsmerged_ECs_PBS_GAS, file = "02_RDS files/dsmerged_ECs_004_005_008_009_014_023_024.rds")
## Differential expression on endothelial cells by MAST
dsmerged_ECs <- readRDS(file = "02_RDS files/dsmerged_ECs_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
cell_types <- c("Endothelial cells")
dsmerged_ECs <- JoinLayers(dsmerged_ECs)
Idents(object = dsmerged_ECs) <- [email hidden]$condition
compare.conditions <- list(
list(ident1 = "GAS", ident2 = "PBS", result_name = "GASvsPBS"),
list(ident1 = "RORg", ident2 = "GAS", result_name = "RORgvsGAS"),
list(ident1 = "IL-17A_mAb", ident2 = "Isotype_mAb", result_name = "IL17A.mAb.vs.Isotype"),
list(ident1 = "Csf2_CD4ko", ident2 = "Csf2flox", result_name = "Csf2.CD4ko.vs.Csf2fl"))
DE_MAST_dsmerged_ECs <- run_DE(
object = dsmerged_ECs,
comparisons_list = compare.conditions,
test_method = "MAST",
output_file = "03_quantification/ECs_AllConditions_MAST_DE.xlsx",
include_sd = FALSE,
save_results = TRUE,
skip_cell_subsetting = TRUE)
## EC PBS vs GAS aggregate heat maps
dsmerged_ECs_PBS_GAS <- readRDS(file = "02_RDS files/dsmerged_ECs_004_005_008_009_014_023_024.rds")
Idents(object=dsmerged_ECs_PBS_GAS) <- [email hidden]$condition
generate_heatmap <- function(seurat_object,
marker_set,
sample_levels = c("CW004", "CW014", "CW023", "CW005", "CW008", "CW009", "CW024"),
output_path = "04_plots/Figure 1",
plot_width = 3.6115,
plot_height = 4.66,
color_min = -2,
color_max = 2) {
seurat_object_name <- deparse(substitute(seurat_object))
marker_set_name <- deparse(substitute(marker_set))
heatmap_filename_pdf <- paste0("Heatmap_", seurat_object_name, "_", marker_set_name, ".pdf")
Idents(object = seurat_object) <- [email hidden]$sample
sample_aggregates <- AggregateExpression(seurat_object, return.seurat = TRUE)
levels(x = sample_aggregates) <- sample_levels
# Create heatmap with labels and save as PDF
heatmap_pdf <- DoHeatmap(sample_aggregates, features = marker_set, draw.lines = FALSE) +
scale_fill_gradient2(low = "blue", high = "red", limits = c(color_min, color_max))
ggsave(
filename = heatmap_filename_pdf,
plot = heatmap_pdf,
device = "pdf",
path = output_path,
scale = 1,
width = plot_width,
height = plot_height,
units = "in",
dpi = 300,
bg = "transparent"
)
}
sample.levels <- c("CW004", "CW014", "CW023", "CW005", "CW008", "CW009", "CW024")
type_I_interferon <- c("Irf7", "Isg15", "Isg20", "Ifit1", "Ifit2", "Ifit3", "Ifitm1", "Ifitm2", "Ifitm3", "Ifi44", "Ifi47", "Ifi204")
antigen_presentation <- c("B2m", "Cd74", "H2-Aa", "H2-Ab1", "H2-Eb1", "H2-K1", "H2-D1", "H2-Q4", "H2-Q7", "H2-T22", "Tap1", "Tap2")
EC_response_to_LPS <- c("Lcn2", "Ch25h", "Tmem252", "Apod", "Litaf", "Ctla2a", "Pltp", "Ptn", "Tspan13", "Cxcl12", "Utrn", "Fcgrt")
ECM_markers <- c("Hmcn1", "Ltbp4", "Cst3", "Ptn", "Itga4", "Mmp15", "Itga1", "Adamts10", "Adamts5", "Lamc3", "Col4a3", "Col14a1")
generate_heatmap(dsmerged_ECs_PBS_GAS, type_I_interferon, ##Figure 1F##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -1.75,
color_max = 2.2)
generate_heatmap(dsmerged_ECs_PBS_GAS, antigen_presentation, ##Figure 1G##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -1.75,
color_max = 2.2)
generate_heatmap(dsmerged_ECs_PBS_GAS, EC_response_to_LPS, ##Figure 1H##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -1.75,
color_max = 2.2)
generate_heatmap(dsmerged_ECs_PBS_GAS, ECM_markers, ##Figure 1I##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -1.75,
color_max = 2.2)
junctional_adhesion_markers <- c("Amot", "Amotl1", "Amotl2", "Arhgap17", "Ash1l", "Cdh5", "Cgnl1", "Cldn12", "Cldn5", "Ctnna1", "Ctnnb1", "Dlg1", "Dlgap4", "F11r", "Jam2", "Lin7c", "Magi3", "Marveld2", "Ocln", "Pard3", "Pard6g", "Scrib", "Tjap1", "Tjp1", "Tjp2", "Wnk1", "Icam1", "Vcam1")
transporters <- c("Abcb1a", "Abcc4", "Abcg2", "Atp10a", "Slc16a1", "Slc16a2", "Slc16a4", "Slc19a3", "Slc1a1", "Slc22a8", "Slc25a20", "Slc25a32", "Slc25a33", "Slc30a1", "Slc31a1", "Slc35f2", "Slc38a3", "Slc38a5", "Slc39a10", "Slc40a1", "Slc46a3", "Slc6a17", "Slc6a6", "Slc7a1", "Slc7a3", "Slco1a4", "Slco1c1", "Slco2b1")
RMT_vs_bulk_trancytosis <- c("Bsg", "Igf1r", "Insr", "Ldlr", "Ldlrad3", "Lrp10", "Lrp8", "Scarb1", "Slc2a1", "Slc3a2", "Slc7a5", "Tfrc", "Ap2a1", "Ap2a2", "Ap2b1", "Cav1", "Cav2", "Cavin1", "Cavin2", "Cavin3", "Cltb", "Cltc", "Epn1", "Eps15", "Fcho2", "Hip1", "Picalm", "Rin3")
other_CNS_EC_markers <- c("Apcdd1", "Car4", "Clic5", "Cobll1", "Ddc", "Degs2", "Ebf1", "Foxq1", "Itih5", "Itm2a", "Jag2", "Kitl", "Mfsd2a", "Nebl", "Pcdh19", "Pglyrp1", "Prom1", "Rad54b", "Sema3c", "Sema7a", "Sgpp2", "Spock2", "Tbx1", "Tbx3", "Tnfrsf19", "Tpd52l1", "Vwa1", "Zic3")
generate_heatmap(dsmerged_ECs_PBS_GAS, junctional_adhesion_markers, ##Extended Figure 1H##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -2.1,
color_max = 2.5)
generate_heatmap(dsmerged_ECs_PBS_GAS, transporters, ##Extended Figure 1I##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -2.1,
color_max = 2.5)
generate_heatmap(dsmerged_ECs_PBS_GAS, RMT_vs_bulk_trancytosis, ##Extended Figure 1J##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -2.1,
color_max = 2.5)
generate_heatmap(dsmerged_ECs_PBS_GAS, other_CNS_EC_markers, ##Extended Figure 1K##
sample_levels = sample.levels,
output_path = "04_plots/Figure 1",
color_min = -2.1,
color_max = 2.5)
rm(sample.levels, type_I_interferon, antigen_presentation, EC_response_to_LPS, ECM_markers, junctional_adhesion_markers, transporters, RMT_vs_bulk_trancytosis, other_CNS_EC_markers)
## Microglia clustering and analysis ----
dsmerged <- readRDS(file = "02_RDS files/dsmerged_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
Idents(object=dsmerged) <- [email hidden]$cluster.names_all.cells
dsmerged_microglia <- subset(x = dsmerged, idents = c("Microglia"),invert = FALSE)
## clustering, batch correction, removing mixed signature
# normalize & generate PCA
dsmerged_microglia <- NormalizeData(dsmerged_microglia)
dsmerged_microglia <- FindVariableFeatures(dsmerged_microglia, selection.method = "vst", nfeatures = 3000)
top10 <- head(VariableFeatures(dsmerged_microglia), 10)
plot1 <- VariableFeaturePlot(dsmerged_microglia)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)
plot2
all.genes <- rownames(dsmerged_microglia)
dsmerged_microglia <- ScaleData(dsmerged_microglia, features = all.genes)
dsmerged_microglia <- RunPCA(dsmerged_microglia, features = VariableFeatures(object = dsmerged_microglia))
print(dsmerged_microglia[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(dsmerged_microglia, dims = 1:5, reduction = "pca")
ElbowPlot(dsmerged_microglia,ndims=50)
dsmerged_microglia <- FindNeighbors(dsmerged_microglia, dims = 1:42)
dsmerged_microglia <- FindClusters(dsmerged_microglia, resolution = 0.4)
dsmerged_microglia <- RunUMAP(dsmerged_microglia, dims = 1:42)
DimPlot(dsmerged_microglia, reduction = "umap",label=T)
DimPlot(dsmerged_microglia, reduction = "umap",group.by="orig.ident")
dsmerged_microglia <- RunTSNE(dsmerged_microglia, dims = 1:42)
DimPlot(dsmerged_microglia, reduction = "tsne",label=T)
DimPlot(dsmerged_microglia, reduction = "tsne",group.by="orig.ident")
rm(plot1, plot2, all.genes, top10)
# run harmony to reduce batch effects
pre.harmony.umap <- DimPlot(dsmerged_microglia, reduction = "umap",group.by="batch")
pre.harmony.tsne <- DimPlot(dsmerged_microglia, reduction = "tsne",group.by="batch")
dsmerged_microglia <- dsmerged_microglia %>%
RunHarmony("batch", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(dsmerged_microglia, 'harmony')
harmony_embeddings[1:5, 1:5]
dsmerged_microglia <- dsmerged_microglia %>%
RunUMAP(reduction = "harmony", dims = 1:42) %>%
FindNeighbors(reduction = "harmony", dims = 1:42) %>%
FindClusters(resolution = 0.4) %>%
identity()
print(pre.harmony.umap)
DimPlot(dsmerged_microglia, reduction = "umap",label=F)
DimPlot(dsmerged_microglia, reduction = "umap",label=F, group.by = "batch")
dsmerged_microglia <- dsmerged_microglia %>%
RunTSNE(reduction = "harmony", dims = 1:42) %>%
FindNeighbors(reduction = "harmony", dims = 1:42) %>%
FindClusters(resolution = 0.4) %>%
identity()
print(pre.harmony.tsne)
DimPlot(dsmerged_microglia, reduction = "tsne",label=T)
DimPlot(dsmerged_microglia, reduction = "tsne", group.by = "batch")
rm(pre.harmony.umap, pre.harmony.tsne, harmony_embeddings)
# remove mixed signature clusters
genes <- VariableFeatures(dsmerged_microglia)
toplot <- CalcStats(dsmerged_microglia, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
neuroglial.markers <- c("Snap25", "Dcx", "Gfap","Aqp4", "Frzb","Cldn5", "Pecam1", "Pdgfra","Pdgfrb", "Atp13a5", "Col1a1", "Fbln1")
lymphoid.markers <- c("Ptprc", "Cd3e", "Cd4", "Cd8a", "Cd163l1", "Klrb1c", "Cd19", "Ms4a1")
myeloid.markers <- c("Ptprc", "Tmem119", "Aif1", "Itgam", "Ly6g", "Camp", "Ly6c2", "Itgax", "Xcr1", "Cd209a", "Ccr9", "Nudt17")
macrophage.vs.microglia.markers <- c("Ptprc", "Tmem119", "Cx3cr1", "P2ry12", "Itgam", "Aif1", "Hexb", "Plac8", "Mrc1", "Cd163")
DotPlot(dsmerged_microglia, features = neuroglial.markers)
FeaturePlot(dsmerged_microglia, features = c("Snap25", "Dcx")) # cluster 6 has mixed signature: microglia/neurons
FeaturePlot(dsmerged_microglia, features = c("Gfap", "Aqp4", "Pecam1")) # cluster 7 has mixed signature: microglia/astrocytes/neutrophils
Idents(object=dsmerged_microglia) <- [email hidden]$RNA_snn_res.0.4
dsmerged_microglia <- subset(x = dsmerged_microglia, idents = c(6, 7),invert = TRUE)
DimPlot(dsmerged_microglia, reduction = "umap",label=T)
DimPlot(dsmerged_microglia, reduction = "umap", split.by = "condition")
dsmerged_microglia[["microglia.cluster"]] <- Idents(object = dsmerged_microglia) # stash identity
genes <- VariableFeatures(dsmerged_microglia)
toplot <- CalcStats(dsmerged_microglia, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
# assign cluster identities
genes <- VariableFeatures(dsmerged_microglia)
toplot <- CalcStats(dsmerged_microglia, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
Idents(object=dsmerged_microglia) <- [email hidden]$microglia.cluster
DimPlot(dsmerged_microglia, reduction = "umap", label=T)
new.cluster.ids <- c("hMG1", "srMG1", "srMG2", "srMG3", "srMG4", "hMG2")
names(new.cluster.ids) <- levels(dsmerged_microglia)
dsmerged_microglia <- RenameIdents(dsmerged_microglia, new.cluster.ids)
DimPlot(dsmerged_microglia, reduction = "umap",label=T)+ NoLegend()
dsmerged_microglia[["microglia.type"]] <- Idents(object = dsmerged_microglia) # stash identity
Idents(object=dsmerged_microglia) <- [email hidden]$microglia.type
[email hidden]$microglia.type <- factor(x = [email hidden]$microglia.type, levels(dsmerged_microglia) <- c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"))
DimPlot(dsmerged_microglia, reduction = "umap",label=T, split.by = "microglia.type")
rm(genes, toplot, lymphoid.markers, macrophage.vs.microglia.markers, myeloid.markers, neuroglial.markers)
saveRDS(dsmerged_microglia, file = "02_RDS files/dsmerged_microglia_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
## subset on PBS and GAS microglia only
Idents(object=dsmerged_microglia) <- [email hidden]$condition
DimPlot(dsmerged_microglia, reduction = "umap")
dsmerged_microglia_PBS_GAS <- subset(x = dsmerged_microglia, idents = c("PBS", "GAS"),invert = FALSE)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap")
# re-cluster and run harmony
dsmerged_microglia_PBS_GAS <- NormalizeData(dsmerged_microglia_PBS_GAS)
dsmerged_microglia_PBS_GAS <- FindVariableFeatures(dsmerged_microglia_PBS_GAS, selection.method = "vst", nfeatures = 3000)
top10 <- head(VariableFeatures(dsmerged_microglia_PBS_GAS), 10)
plot1 <- VariableFeaturePlot(dsmerged_microglia_PBS_GAS)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)
plot2
all.genes <- rownames(dsmerged_microglia_PBS_GAS)
dsmerged_microglia_PBS_GAS <- ScaleData(dsmerged_microglia_PBS_GAS, features = all.genes)
dsmerged_microglia_PBS_GAS <- RunPCA(dsmerged_microglia_PBS_GAS, features = VariableFeatures(object = dsmerged_microglia_PBS_GAS))
print(dsmerged_microglia_PBS_GAS[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(dsmerged_microglia_PBS_GAS, dims = 1:5, reduction = "pca")
ElbowPlot(dsmerged_microglia_PBS_GAS,ndims=50)
dsmerged_microglia_PBS_GAS <- FindNeighbors(dsmerged_microglia_PBS_GAS, dims = 1:28)
dsmerged_microglia_PBS_GAS <- FindClusters(dsmerged_microglia_PBS_GAS, resolution = 1)
dsmerged_microglia_PBS_GAS <- RunUMAP(dsmerged_microglia_PBS_GAS, dims = 1:28)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=T)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",group.by="orig.ident")
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",group.by="batch")
dsmerged_microglia_PBS_GAS <- RunTSNE(dsmerged_microglia_PBS_GAS, dims = 1:28)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "tsne",label=T)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "tsne",group.by="orig.ident")
rm(plot1, plot2, all.genes, top10)
pre.harmony.umap <- DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap", group.by="batch", label = F)
pre.harmony.tsne <- DimPlot(dsmerged_microglia_PBS_GAS, reduction = "tsne", group.by="batch", label = F)
dsmerged_microglia_PBS_GAS <- dsmerged_microglia_PBS_GAS %>%
RunHarmony("batch", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(dsmerged_microglia_PBS_GAS, 'harmony')
harmony_embeddings[1:5, 1:5]
dsmerged_microglia_PBS_GAS <- dsmerged_microglia_PBS_GAS %>%
RunUMAP(reduction = "harmony", dims = 1:28) %>%
FindNeighbors(reduction = "harmony", dims = 1:28) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.umap)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=T)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=F, group.by = "batch") ##Extended Data Figure S2A##
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=F, group.by = "inoculate") ##Figure 2A##
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=F, split.by = "microglia.type") ##Figure 2A##
dsmerged_microglia_PBS_GAS <- dsmerged_microglia_PBS_GAS %>%
RunTSNE(reduction = "harmony", dims = 1:28) %>%
FindNeighbors(reduction = "harmony", dims = 1:28) %>%
FindClusters(resolution = 1) %>%
identity()
print(pre.harmony.tsne)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "tsne",label=T)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "tsne",label=F, group.by="clust")
dsmerged_microglia_PBS_GAS[["clusters.res0.4_PBS.GAS.microglia"]] <- Idents(object = dsmerged_microglia_PBS_GAS) # stash identity
rm(pre.harmony.umap, pre.harmony.tsne, harmony_embeddings)
# assign cluster identities
genes <- VariableFeatures(dsmerged_microglia_PBS_GAS)
toplot <- CalcStats(dsmerged_microglia_PBS_GAS, features = genes, method = "zscore", order = "p", n = 5)
Heatmap(toplot, lab_fill = "zscore")
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$clusters.res0.4_PBS.GAS.microglia
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=T)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=F, split.by="microglia.type")
new.cluster.ids <- c("hMG1", "hMG1", "hMG1", "hMG1", "srMG1", "srMG2", "srMG3", "srMG4", "hMG2", "srMG1", "srMG3")
names(new.cluster.ids) <- levels(dsmerged_microglia_PBS_GAS)
dsmerged_microglia_PBS_GAS <- RenameIdents(dsmerged_microglia_PBS_GAS, new.cluster.ids)
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=T)+ NoLegend()
dsmerged_microglia_PBS_GAS[["microglia.type_PBS.GAS"]] <- Idents(object = dsmerged_microglia_PBS_GAS) # stash identity
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$microglia.type_PBS.GAS
[email hidden]$microglia.type_PBS.GAS <- factor(x = [email hidden]$microglia.type_PBS.GAS, levels(dsmerged_microglia_PBS_GAS) <- c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"))
DimPlot(dsmerged_microglia_PBS_GAS, reduction = "umap",label=T, split.by = "microglia.type_PBS.GAS")
table(dsmerged_microglia_PBS_GAS$microglia.type, dsmerged_microglia_PBS_GAS$inoculate) ##Extended Data Figure 2B##
exAM.score.list <- c("Fos", "Jun", "Dusp1", "Hspa1a", "Hist1h1d", "Hist1h2ac", "Nfkbid", "Nfkbiz")
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$condition
Dotplot.exAM <- DotPlot(dsmerged_microglia_PBS_GAS, exAM.score.list) + coord_flip() ##Extended Data Figure 2D##
ggsave("04_plots/Figure 2/Dotplot_exAM_dsmerged_microglia.pdf",
plot = Dotplot.exAM,
width = 4, height = 6,
units = "in")
dsmerged_microglia_PBS_GAS <- AddModuleScore(dsmerged_microglia_PBS_GAS,
features = list(exAM.score.list),
name = "exAM.score")
[email hidden]$exAM.score <- [email hidden]$exAM.score1
[email hidden]$exAM.score1 <- NULL
head(x=dsmerged_microglia_PBS_GAS)
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$condition
p1 <- FeatureScatter(dsmerged_microglia_PBS_GAS, feature1 = "Ccl3", feature2 = "Ccl4") + scale_color_manual(values = c('gray40','skyblue'))
p2 <- FeatureScatter(dsmerged_microglia_PBS_GAS, feature1 = "Ccl3", feature2 = "exAM.score") + scale_color_manual(values = c('gray40','skyblue'))
p3 <- FeatureScatter(dsmerged_microglia_PBS_GAS, feature1 = "Ccl4", feature2 = "exAM.score") + scale_color_manual(values = c('gray40','skyblue'))
library(cowplot) # version 1.1.3
exAM_plots <- plot_grid(p1, p2, p3, nrow = 3)
ggsave("04_plots/Figure 2/Scatter_exAM_plots.pdf",
plot = exAM_plots,
width = 4, height = 10,
units = "in")
rm(genes, toplot, lymphoid.markers, macrophage.vs.microglia.markers, myeloid.markers, neuroglial.markers, p1, p2, p3, exAM.score.list, exAM_plots)
saveRDS(dsmerged_microglia_PBS_GAS, file = "02_RDS files/dsmerged_microglia_004_005_008_009_014_023_024.rds")
## Differential expression on microglia by MAST
dsmerged_microglia <- readRDS(file = "02_RDS files/dsmerged_microglia_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
cell_types <- c("Microglia")
dsmerged_microglia <- JoinLayers(dsmerged_microglia)
Idents(object = dsmerged_microglia) <- [email hidden]$condition
compare.conditions <- list(
list(ident1 = "GAS", ident2 = "PBS", result_name = "GASvsPBS"),
list(ident1 = "RORg", ident2 = "GAS", result_name = "RORgvsGAS"),
list(ident1 = "IL-17A_mAb", ident2 = "Isotype_mAb", result_name = "IL17A.mAb.vs.Isotype"),
list(ident1 = "Csf2_CD4ko", ident2 = "Csf2flox", result_name = "Csf2.CD4ko.vs.Csf2fl"))
DE_MAST_dsmerged_microglia <- run_DE(
object = dsmerged_microglia,
comparisons_list = compare.conditions,
test_method = "MAST",
output_file = "03_quantification/Microglia_AllConditions_MAST_DE.xlsx",
include_sd = FALSE,
save_results = TRUE,
skip_cell_subsetting = TRUE)
## Microglia PBS vs GAS aggregate heat maps
dsmerged_microglia_PBS_GAS <- readRDS(file = "02_RDS files/dsmerged_microglia_004_005_008_009_014_023_024.rds")
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$condition
sample.levels <- c("CW004", "CW014", "CW023", "CW005", "CW008", "CW009", "CW024")
microglial.genes <- c("B2m", "H2-D1", "H2-K1", "Tap1", "Cd74", "H2-Aa", "H2-Ab1", "H2-Eb1", "Ccl2", "Ccl3", "Ccl4", "Ccl5", "Ccl12", "Cxcl10", "Il1b", "Tnf", "Stat1", "Stat2", "Irf1", "Irf7", "Isg15", "Ifi30", "Ifi204", "Ifi211", "Ifitm3", "Oas1a")
generate_heatmap(dsmerged_microglia_PBS_GAS, microglial.genes, ##Figure 2C##
sample_levels = sample.levels,
output_path = "04_plots/Figure 2",
color_min = -2,
color_max = 2)
## Microglia PBS vs GAS heatmaps by subcluster
generate_heatmap_Z <- function(seurat_object,
marker_set,
cluster_levels = NULL,
output_path = "04_plots/Figure 1",
plot_width = 3.6115,
plot_height = 4.66,
color_min = -2,
color_max = 2) {
seurat_object_name <- deparse(substitute(seurat_object))
marker_set_name <- deparse(substitute(marker_set))
heatmap_filename_pdf <- paste0("Heatmap_", seurat_object_name, "_", marker_set_name, "_zscore.pdf")
cluster_aggregates <- AggregateExpression(seurat_object, return.seurat = TRUE)
cluster_aggregates <- ScaleData(cluster_aggregates, features = marker_set)
heatmap_pdf <- DoHeatmap(cluster_aggregates,
features = marker_set,
slot = "scale.data", # Use z-scored data
draw.lines = FALSE,
raster = FALSE) +
scale_fill_gradient2(low = "blue", high = "red", limits = c(color_min, color_max)) +
labs(fill = "Z-score") +
theme(
axis.text.y.left = element_blank(),
axis.text.y.right = element_text(size = 8),
axis.text.x.top = element_blank(),
axis.text.x.bottom = element_text(angle = 45, hjust = 1, size = 8)
) +
scale_y_discrete(position = "right") +
scale_x_discrete(position = "bottom")
ggsave(
filename = heatmap_filename_pdf,
plot = heatmap_pdf,
device = "pdf",
path = output_path,
scale = 1,
width = plot_width,
height = plot_height,
units = "in",
dpi = 300,
bg = "transparent"
)
}
Idents(object=dsmerged_microglia_PBS_GAS) <- [email hidden]$microglia.type
cytokine.signaling <- c("Csf2ra", "Ccl5", "Cxcl9", "Ccl7", "Il1rn", "Cxcl10", "Cxcl16", "Ccl12", "Ccl2", "Cxcl13", "Il2rg", "Il1b", "Mif", "Tnfsf13b", "Il1a", "Tnf", "Il4ra", "Csf2rb2", "Ccl4", "Csf2rb", "Ccl3", "Osm", "Csf1")
homeostatic <- c("Serinc3", "Tgfrb1", "P2ry12", "Sall1", "Siglech", "Mef2a", "Cd33", "Smad3", "Tgfb1", "Gpr34", "Adgrg1", "Csf1r", "Fcrls", "Ccr5", "Bin1", "Tmem119", "Cx3cr1", "Olfml3", "Glul", "Jun", "Mafb", "Egr1", "Spi1")
disease.associated <- c("Ctsl", "Lpl", "Ctsb", "Cd63", "Lyz2", "Apoe", "Spp1", "Msr1", "Axl", "Il1b", "Itgax", "Cst7", "Fth1", "Nos2", "Cybb", "Lgals3", "Grn", "Bhlhe40", "Vegfa", "Fabp5", "Csf1", "Timp2", "Gnas", "Cd9", "Mertk", "Ccl6", "Ctsd", "Trem2", "Tyrobp")
IFN <- c("Stat2", "Oasl2", "Iigp1", "Ifi206", "Irf7", "Igtp", "Ifi211", "Ifit1", "Oasl1", "Oas3", "Ifi47", "Ifi205", "Ifi209", "Ifit3", "Ifit2", "Ifi213", "Isg20", "Stat1", "Cxcl9", "Ifitm3", "Ifi35", "Ifi204", "Ifi207", "Isg15", "Cxcl10", "Oas1a", "Ifitm2", "Ifitm1", "Irf1")
generate_heatmap_Z(dsmerged_microglia_PBS_GAS, cytokine.signaling, ##Figure 2K##
cluster_levels = c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"),
output_path = "04_plots/Figure 2",
color_min = -2.2,
color_max = 2.2)
generate_heatmap_Z(dsmerged_microglia_PBS_GAS, homeostatic, ##Figure 2L##
cluster_levels = c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"),
output_path = "04_plots/Figure 2",
color_min = -2.2,
color_max = 2.2)
generate_heatmap_Z(dsmerged_microglia_PBS_GAS, disease.associated, ##Figure 2M##
cluster_levels = c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"),
output_path = "04_plots/Figure 2",
color_min = -2.2,
color_max = 2.2)
generate_heatmap_Z(dsmerged_microglia_PBS_GAS, IFN, ##Figure 2N##
cluster_levels = c("hMG1", "hMG2", "srMG1", "srMG2", "srMG3", "srMG4"),
output_path = "04_plots/Figure 2",
color_min = -2.2,
color_max = 2.2)
rm(generate_heatmap_Z, cytokine.signaling, homeostatic, disease.associated, IFN)
## RORgt heatmaps and GO analysis ----
dsmerged_ECs <- readRDS(file = "02_RDS files/dsmerged_ECs_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
Idents(object=dsmerged_ECs) <- [email hidden]$condition
dsmerged_ECs_RORg_WT_GAS <- subset(x = dsmerged_ECs, idents = c("GAS", "RORg"),invert = FALSE)
dsmerged_ECs_RORg_WT_GAS <- JoinLayers(dsmerged_ECs_RORg_WT_GAS)
Idents(object=dsmerged_ECs_RORg_WT_GAS) <- [email hidden]$sample
sample.levels <- c("CW005", "CW008", "CW009", "CW024", "CW010", "CW013", "CW018")
genelist.ECs <-c("Itih5", "Tjp1", "Slc16a1", "Mfsd2a", "Bsg", "Ch25h", "Icam1", "Lcn2", "Lrg1", "B2m", "Cd74", "H2-D1", "H2-K1", "Tap1", "Ifi208", "Ifitm2", "Ifitm3", "Ifi47", "Iigp1", "Igtp", "Irgm1", "Gbp2")
generate_heatmap(dsmerged_ECs_RORg_WT_GAS, genelist.ECs, ##Figure 4B##
sample_levels = sample.levels,
output_path = "04_plots/Figure 4",
color_min = -2,
color_max = 2.2)
dsmerged_microglia <- readRDS(file = "02_RDS files/dsmerged_microglia_004_005_008_009_010_011_012_013_014_015_016_017_018_021_022_023_024_025_026.rds")
Idents(object=dsmerged_microglia) <- [email hidden]$condition
dsmerged_microglia_RORg_WT_GAS <- subset(x = dsmerged_microglia, idents = c("GAS", "RORg"),invert = FALSE)
dsmerged_microglia_RORg_WT_GAS <- JoinLayers(dsmerged_microglia_RORg_WT_GAS)
Idents(object=dsmerged_microglia_RORg_WT_GAS) <- [email hidden]$sample
sample.levels <- c("CW005", "CW008", "CW009", "CW024", "CW010", "CW013", "CW018")
genelist.microglia <-c("B2m", "Cd74", "H2-D1", "H2-K1", "Tap1", "Tap2", "H2-Ab1", "H2-Eb1", "Cx3cr1", "Gpr34", "P2ry12", "Bhlhe40", "Cd63", "Cst7", "Spp1", "Ccl2", "Ccl3", "Ccl4", "Ccl5", "Cxcl10", "Il1b", "Tnf")
generate_heatmap(dsmerged_microglia_RORg_WT_GAS, RORvsWT_microglia, ##Figure 4E##
sample_levels = sample.levels,
output_path = "04_plots/Figure 4",
color_min = -2,
color_max = 2.2)
## Run DE with standard deviation for GSEA Signal2Noise parameter
compare.conditions <- list(
list(ident1 = "RORg", ident2 = "GAS", result_name = "RORgvsGAS")
)
DE.MAST.RORg.ECs <- run_DE(
object = dsmerged_ECs_RORg_WT_GAS,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/RORgvsGAS_ECs_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
DE.MAST.RORg.microglia <- run_DE(
object = dsmerged_microglia_RORg_WT_GAS,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/RORgvsGAS_Microglia_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
## GSEA using GO pathway lists
library(readxl) # version 1.4.5
library(clusterProfiler) # version 4.12.6
library(enrichplot) # version 1.24.4
library(org.Mm.eg.db) # version 3.19.1
run_gsea_analysis <- function(cell_type,
input_file,
compare_conditions,
sheet_name = "RORgvsGAS",
p_cutoff = 0.05,
show_categories = 5,
manual_go_terms = NULL,
plot_width = 5,
plot_height = 6) {
cat("Running GSEA analysis for", cell_type, "\n")
de_data <- read_excel(input_file, sheet = sheet_name)
excluded_patterns <- c("^mt-", "^Rps", "^Rpl", "^Mrps", "^Mrpl")
clean_data <- de_data[!is.na(de_data$gene) &
!is.na(de_data$Signal2Noise) &
!grepl(paste(excluded_patterns, collapse = "|"), de_data$gene), ]
gene_list <- setNames(clean_data$Signal2Noise, clean_data$gene)
gene_list <- sort(gene_list, decreasing = TRUE)
gse_results <- gseGO(geneList = gene_list,
OrgDb = org.Mm.eg.db,
ont = "BP",
keyType = "SYMBOL",
minGSSize = 10,
maxGSSize = 800,
pvalueCutoff = p_cutoff)
table_output <- paste0("03_quantification/GSEA_", cell_type, "_", compare_conditions, ".csv")
plot_output <- paste0("04_plots/Figure 4/GSEA_dotplot_", cell_type, "_", compare_conditions, ".pdf")
write.csv(as.data.frame(gse_results), table_output, row.names = FALSE)
gse_plot <- gse_results # prepare data for plotting
if (!is.null(manual_go_terms)) { # filter for manual go terms if provided
gse_plot@result <- gse_plot@result[
gse_plot@result$ID %in% manual_go_terms |
gse_plot@result$Description %in% manual_go_terms, ]
cat("Plotting", nrow(gse_plot@result), "manually selected GO terms\n")
}
plot_categories <- ifelse(!is.null(manual_go_terms), nrow(gse_plot@result), show_categories)
gse_plot@result$qvalue <- -gse_plot@result$enrichmentScore # replace qvalue after selection
p <- dotplot(gse_plot,
color = "qvalue",
x = "NES",
showCategory = plot_categories,
split = ".sign") +
scale_color_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, name = "enrichmentScore") +
xlim(-3.5, 3.5) +
ggtitle(paste("GSEA:", cell_type, compare_conditions))
ggsave(plot_output, plot = p, width = plot_width, height = plot_height, dpi = 300)
print(p)
cat("Complete! Files saved:", table_output, "and", plot_output, "\n")
return(gse_results)
}
# Note: manually added "Signal2Noise" column in Excel and removed "_CellType" from sheet name
EC_results <- run_gsea_analysis( ##Figure 4B##
cell_type = "ECs",
manual_go_terms = c("GO:0035458", "GO:0001568", "GO:0019882", "GO:0034341", "GO:0016477", "GO:0007155", "GO:0001568", "GO:0050671", "GO:0030334"),
input_file = "03_quantification/RORgvsGAS_ECs_MAST_DE_withSD.xlsx",
compare_conditions = "RORg_vs_GAS"
)
microglia_results <- run_gsea_analysis( ##Figure 4E##
cell_type = "Microglia",
manual_go_terms = c("GO:0019882", "GO:0006119", "GO:0034341", "GO:0016477", "GO:0001817", "GO:0042129", "GO:0035458", "GO:0097396"),
input_file = "03_quantification/RORgvsGAS_Microglia_MAST_DE_withSD.xlsx",
compare_conditions = "RORg_vs_GAS"
)
## Csf2-CD4KO heatmaps ----
Idents(object=dsmerged_ECs) <- [email hidden]$condition
dsmerged_ECs_Csf2CD4 <- subset(x = dsmerged_ECs, idents = c("Csf2_CD4ko", "Csf2flox"),invert = FALSE)
dsmerged_ECs_Csf2CD4 <- JoinLayers(dsmerged_ECs_Csf2CD4)
Idents(object=dsmerged_ECs_Csf2CD4) <- [email hidden]$sample
sample.levels <- c("CW021", "CW025", "CW022", "CW026")
generate_heatmap(dsmerged_ECs_Csf2CD4, genelist.ECs, ##Extended Data Figure 4I##
sample_levels = sample.levels,
output_path = "04_plots/Figure 4",
color_min = -2,
color_max = 2.2)
Idents(object=dsmerged_microglia) <- [email hidden]$condition
dsmerged_microglia_Csf2CD4 <- subset(x = dsmerged_microglia, idents = c("Csf2_CD4ko", "Csf2flox"),invert = FALSE)
dsmerged_microglia_Csf2CD4 <- JoinLayers(dsmerged_microglia_Csf2CD4)
Idents(object=dsmerged_microglia_Csf2CD4) <- [email hidden]$sample
generate_heatmap(dsmerged_microglia_Csf2CD4, genelist.microglia, ##Extended Data Figure 4J##
sample_levels = sample.levels,
output_path = "04_plots/Figure 4",
color_min = -2,
color_max = 2.2)
compare.conditions <- list(
list(ident1 = "Csf2_CD4ko", ident2 = "Csf2flox", result_name = "Csf2_CD4kovsFlox")
)
DE.MAST.Csf2CD4ko.ECs <- run_DE(
object = dsmerged_ECs_Csf2CD4,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/Csf2_CD4kovsFlox_ECs_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
DE.MAST.Csf2CD4ko.microglia <- run_DE(
object = dsmerged_microglia_Csf2CD4,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/Csf2_CD4kovsFlox_microglia_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
EC_results <- run_gsea_analysis( ##Extended Data Figure 4J##
cell_type = "ECs",
sheet_name = "Csf2_CD4kovsFlox",
manual_go_terms = c("GO:0035456", "GO:0009617", "GO:0019883", "GO:0034341", "GO:1990868"),
input_file = "03_quantification/Csf2_CD4kovsFlox_ECs_MAST_DE_withSD.xlsx",
compare_conditions = "Csf2_CD4kovsFlox"
)
microglia_results <- run_gsea_analysis( ##Extended Data Figure 4J##
cell_type = "Microglia",
sheet_name = "Csf2_CD4kovsFlox",
manual_go_terms = c("GO:0034097", "GO:0034341", "GO:0034097"),
input_file = "03_quantification/Csf2_CD4kovsFlox_microglia_MAST_DE_withSD.xlsx",
compare_conditions = "Csf2_CD4kovsFlox"
)
## IL-17A mAb heatmaps and ridge plots ----
Idents(object=dsmerged_ECs) <- [email hidden]$condition
dsmerged_ECs_IL17A.mAb.Isotype <- subset(x = dsmerged_ECs, idents = c("IL-17A_mAb", "Isotype_mAb"),invert = FALSE)
dsmerged_ECs_IL17A.mAb.Isotype <- JoinLayers(dsmerged_ECs_IL17A.mAb.Isotype)
Idents(object=dsmerged_ECs_IL17A.mAb.Isotype) <- [email hidden]$sample
sample.levels <- c("CW011", "CW015", "CW012", "CW016", "CW017")
generate_heatmap(dsmerged_ECs_IL17A.mAb.Isotype, genelist.ECs, ##Figure 5B##
sample_levels = sample.levels,
output_path = "04_plots/Figure 5",
color_min = -2,
color_max = 2.2)
Idents(object=dsmerged_microglia) <- [email hidden]$condition
dsmerged_microglia_IL17A.mAb.Isotype <- subset(x = dsmerged_microglia, idents = c("IL-17A_mAb", "Isotype_mAb"),invert = FALSE)
dsmerged_microglia_IL17A.mAb.Isotype <- JoinLayers(dsmerged_microglia_IL17A.mAb.Isotype)
Idents(object=dsmerged_microglia_IL17A.mAb.Isotype) <- [email hidden]$sample
generate_heatmap(dsmerged_microglia_IL17A.mAb.Isotype, genelist.microglia, ##Extended Data Figure 5D##
sample_levels = sample.levels,
output_path = "04_plots/Figure 5",
color_min = -2,
color_max = 2.2)
compare.conditions <- list(
list(ident1 = "IL-17A_mAb", ident2 = "Isotype_mAb", result_name = "IL17AmAbvsIsotype")
)
DE.MAST.IL17A_mAb.ECs <- run_DE(
object = dsmerged_ECs_IL17A.mAb.Isotype,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/IL17AmAbvsIsotype_ECs_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
DE.MAST.IL17A_mAb.microglia <- run_DE(
object = dsmerged_microglia_IL17A.mAb.Isotype,
comparisons_list = compare.conditions,
save_results = TRUE,
test_method = "MAST",
output_file = "03_quantification/IL17AmAbvsIsotype_microglia_MAST_DE_withSD.xlsx",
include_sd = TRUE,
skip_cell_subsetting = TRUE
)
EC_results <- run_gsea_analysis( ##Figure 5C##
cell_type = "ECs",
sheet_name = "IL17AmAbvsIsotype",
manual_go_terms = c("GO:0002396", "GO:0001959", "GO:0019221", "GO:0032496", "GO:0006119", "GO:0045333", "GO:0000902", "GO:0009615", "GO:0009617"),
input_file = "03_quantification/IL17AmAbvsIsotype_ECs_MAST_DE_withSD.xlsx",
compare_conditions = "IL-17A mAb vs Isotype"
)
microglia_results <- run_gsea_analysis( ##Extended Data Figure 5E##
cell_type = "Microglia",
sheet_name = "IL17AmAbvsIsotype",
manual_go_terms = c("GO:0006119", "GO:0019882", "GO:0001816", "GO:0034341", "GO:0032728"),
input_file = "03_quantification/IL17AmAbvsIsotype_Microglia_MAST_DE_withSD.xlsx",
compare_conditions = "IL-17A mAb vs Isotype"
)
# Ridge plots for MHC I antigen presentation
Idents(object=dsmerged_ECs) <- [email hidden]$condition
dsmerged_ECs_antigen <- subset(x = dsmerged_ECs, idents = c("Csf2_CD4ko", "Csf2flox"), invert = TRUE)
Idents(object=dsmerged_microglia) <- [email hidden]$condition
dsmerged_microglia_antigen <- subset(x = dsmerged_microglia, idents = c("Csf2_CD4ko", "Csf2flox"), invert = TRUE)
Idents(object=dsmerged) <- [email hidden]$condition
dsmerged_antigen <- subset(x = dsmerged, idents = c("Csf2_CD4ko", "Csf2flox"), invert = TRUE)
Idents(object=dsmerged_antigen) <- [email hidden]$cluster.names_all.cells
dsmerged_astrocytes_antigen <- subset(x = dsmerged_antigen, idents = c("Astrocytes"), invert = FALSE)
Idents(object=dsmerged_astrocytes_antigen) <- [email hidden]$condition
dsmerged_OECs_antigen <- subset(x = dsmerged_antigen, idents = c("Olfactory ensheathing cells"), invert = FALSE)
Idents(object=dsmerged_OECs_antigen) <- [email hidden]$condition
make_ridge_plot <- function(data, feature) {
RidgePlot(data, features = feature) +
scale_fill_manual(values = c("gray80", "skyblue", "mediumseagreen", "wheat", "goldenrod1")) +
xlim(-0.45, 5) +
theme_minimal() +
NoLegend() +
theme(
plot.title = element_text(size = 14, face = "italic", hjust = 0.5),
axis.text.y = element_blank(),
axis.title = element_blank(),
panel.grid.major = element_line(color = "gray90", size = 0.5),
panel.grid.minor = element_blank()
)
}
# Now create your plots with one line each:
p1 <- make_ridge_plot(dsmerged_astrocytes_antigen, "B2m")
p2 <- make_ridge_plot(dsmerged_astrocytes_antigen, "H2-D1")
p3 <- make_ridge_plot(dsmerged_astrocytes_antigen, "H2-K1")
p4 <- make_ridge_plot(dsmerged_OECs_antigen, "B2m")
p5 <- make_ridge_plot(dsmerged_OECs_antigen, "H2-D1")
p6 <- make_ridge_plot(dsmerged_OECs_antigen, "H2-K1")
p7 <- make_ridge_plot(dsmerged_ECs_antigen, "B2m")
p8 <- make_ridge_plot(dsmerged_ECs_antigen, "H2-D1")
p9 <- make_ridge_plot(dsmerged_ECs_antigen, "H2-K1")
p10 <- make_ridge_plot(dsmerged_microglia_antigen, "B2m")
p11 <- make_ridge_plot(dsmerged_microglia_antigen, "H2-D1")
p12 <- make_ridge_plot(dsmerged_microglia_antigen, "H2-K1")
combined_ridge_plot <- plot_grid(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12,
nrow = 4)
ggsave("04_plots/Figure 5/AntigenPresentation_Ridgeplots.pdf", ##Extended Data Figure 5F##
plot = combined_ridge_plot,
width = 8, height = 10,
units = "in")

01_WayneEtAl_scRNAseq.R at commit 28b8a8e, no license · at the source

Overview

Authors: Charlotte R. Wayne¹, Uğur Akcan¹, Travis E. Faust², Violeta Durán-Laforet², Danny Jamoul¹, Luca Bremner¹, Nicole Ampatey¹, Büşra T. Akcan¹, Sarah J. Ho¹, Bogoljub Ciric³, Shannon L. Delaney⁴, Wendy S. Vargas¹, Susan Swedo⁵, Vilas Menon¹, Dorothy P. Schafer², Tyler Cutforth¹, Dritan Agalliu^1,6

ORCID iDs: Charlotte R. Wayne, Travis E. Faust, Vilas Menon, Dorothy P. Schafer, Dritan Agalliu

Department of Neurology, Columbia University Irving Medical Center,New York, NY USA
Department of Neurobiology, Brudnick Neuropsychiatric Research Institute, University of Massachusetts Chan Medical School,Worcester, MA USA
Department of Neurology, Thomas Jefferson University,Philadelphia, PA USA
Lyme and Tick-Borne Diseases Research Center, Department of Psychiatry, Columbia University Irving Medical Center,New York, NY USA
Section of Behavioral Pediatrics, National Institute of Mental Health,Bethesda, MD USA
Department of Pathology and Cell Biology, Columbia University Irving Medical Center,New York, NY USA

Institutions: Columbia University Irving Medical Center (United States); University of Massachusetts Chan Medical School (United States); Thomas Jefferson University (United States); National Institute of Mental Health (United States)

Journal: Nature communications, volume 17, issue 1, article 9658

Dates: received 13 April 2023; accepted 13 July 2026; published online 11 August 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41467-026-76232-w · PMID 42717209 · PMCID PMC13558752 · OpenAlex W7202203061

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), mouse (organism), other condition (population), cellular / molecular (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions, Evoked potentials, Machine learning, fMRI & imaging

Keywords: Neuroimmunology, Infection, Paediatric neurological disorders

MeSH: Cytokines*, Encephalitis*, Endothelial Cells*, Microglia*, Streptococcal Infections*, Th17 Cells*, Animals, Blood-Brain Barrier, Brain, Chemokines, Disease Models, Animal, Female, Granulocyte-Macrophage Colony-Stimulating Factor, Humans, Interleukin-17, Mice, Mice, Inbred C57BL, Mice, Knockout, Receptors, Interleukin-17, Signal Transduction (* major topic)

Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex

Funding: National Institute of Mental Health (R01MH112849, 56MH109987, R01MH113743); HHS | NIH | National Heart, Lung, and Blood Institute (NHLBI) (R61HL159949, R33HL159949); U.S. Department of Health & Human Services | NIH | National Eye Institute (NEI) (R01EY033994); National Institute on Aging (RF1AG078352, RF1AG068281, R01AG068281); International OCD Foundation (International Obsessive-Compulsive Disorder Foundation); NINDS (2T32NS064928, R01NS117533); BrightFocus Foundation (BrightFocus) (A2022006F); Alzheimer’s Association (AARF-22-923219); Massachusetts Life Sciences Center

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

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

AgalliuLab/s-pyogenes-scRNAseq

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 28b8a8e8ac4fdf7fb9e078ef138eb37f4cb856c6, 27 October 2025

Languages: R (5)

Size: 5 files, 5 scripts

Software Heritage: not archived

Found in: “Code availability”

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

Tools: Seurat (5 files), tidyverse (5 files), ggplot2 (4 files), patchwork (2 files), circlize (1 file), clusterProfiler (1 file), ComplexHeatmap (1 file), cowplot (1 file), ggpubr (1 file), Harmony (1 file), reshape2 (1 file), reticulate (1 file), rstatix (1 file), SingleCellExperiment (1 file), UMAP (1 file)

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

27 September 2026: the link answers

5 files

01_WayneEtAl_scRNAseq.R, R, 1,477 lines, 6 matches
02_WayneEtAl_CrossEntropy.R, R, 274 lines, 2 matches
03_WayneEtAl_MERFISH-1.R, R, 373 lines, 2 matches
04_WayneEtAl_MERFISH-2.R, R, 1,509 lines, 3 matches
05_WayneEtAl_CellChat.R, R, 366 lines, 1 match

At the source: github.com/AgalliuLab/s-pyogenes-scRNAseq

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: AgalliuLab/s-pyogenes-scRNAseq
it says that the code is available on request

Read it in the paper: doi.org/10.1038/s41467-026-76232-w.

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

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

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: NCBI GEO GSE221106
it points to the authors' code: AgalliuLab/s-pyogenes-scRNAseq
it says that the code is available on request

Read it in the paper: doi.org/10.1038/s41467-026-76232-w.

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 17 authors, 3 keywords, 20 MeSH terms, 9 funders, 103 references.

Cite

This paper

Wayne, C. R., Akcan, U., Faust, T. E., Durán-Laforet, V., Jamoul, D., Bremner, L., Ampatey, N., Akcan, B. T., Ho, S. J., Ciric, B., Delaney, S. L., Vargas, W. S., Swedo, S., Menon, V., Schafer, D. P., Cutforth, T., & Agalliu, D. (2026). Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis. Nature communications, 17(1), 9658. https://doi.org/10.1038/s41467-026-76232-w

BibTeX

@article{wayne2026th17, author = {Wayne, Charlotte R. and Akcan, Uğur and Faust, Travis E. and Durán-Laforet, Violeta and Jamoul, Danny and Bremner, Luca and Ampatey, Nicole and Akcan, Büşra T. and Ho, Sarah J. and Ciric, Bogoljub and Delaney, Shannon L. and Vargas, Wendy S. and Swedo, Susan and Menon, Vilas and Schafer, Dorothy P. and Cutforth, Tyler and Agalliu, Dritan}, title = {{Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis}}, journal = {Nature communications}, year = {2026}, month = aug, volume = {17}, number = {1}, pages = {9658}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-76232-w}, url = {https://doi.org/10.1038/s41467-026-76232-w}, pmid = {42717209}, pmcid = {PMC13558752} }

RIS

TY - JOUR AU - Wayne, Charlotte R. AU - Akcan, Uğur AU - Faust, Travis E. AU - Durán-Laforet, Violeta AU - Jamoul, Danny AU - Bremner, Luca AU - Ampatey, Nicole AU - Akcan, Büşra T. AU - Ho, Sarah J. AU - Ciric, Bogoljub AU - Delaney, Shannon L. AU - Vargas, Wendy S. AU - Swedo, Susan AU - Menon, Vilas AU - Schafer, Dorothy P. AU - Cutforth, Tyler AU - Agalliu, Dritan TI - Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/08/11 VL - 17 IS - 1 SP - 9658 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-76232-w UR - https://doi.org/10.1038/s41467-026-76232-w LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-76232-w", "type": "article-journal", "title": "Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis", "container-title": "Nature communications", "author": [ { "family": "Wayne", "given": "Charlotte R." }, { "family": "Akcan", "given": "Uğur" }, { "family": "Faust", "given": "Travis E." }, { "family": "Durán-Laforet", "given": "Violeta" }, { "family": "Jamoul", "given": "Danny" }, { "family": "Bremner", "given": "Luca" }, { "family": "Ampatey", "given": "Nicole" }, { "family": "Akcan", "given": "Büşra T." }, { "family": "Ho", "given": "Sarah J." }, { "family": "Ciric", "given": "Bogoljub" }, { "family": "Delaney", "given": "Shannon L." }, { "family": "Vargas", "given": "Wendy S." }, { "family": "Swedo", "given": "Susan" }, { "family": "Menon", "given": "Vilas" }, { "family": "Schafer", "given": "Dorothy P." }, { "family": "Cutforth", "given": "Tyler" }, { "family": "Agalliu", "given": "Dritan" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "9658", "DOI": "10.1038/s41467-026-76232-w", "PMID": "42717209", "PMCID": "PMC13558752", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-76232-w", "language": "en", "issued": { "date-parts": [ [ 2026, 8, 11 ] ] } }

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.1038/s41593-026-02367-0 [code]: A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, SingleCellExperiment, reticulate, 10 other tools, cellular / molecular, 4 references
[2] 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: Harmony, SingleCellExperiment, UMAP, 10 other tools, genetics / omics, other condition, cellular / molecular, 2 references
[3] doi:10.1038/s41586-026-10629-x [code]: Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: Harmony, reticulate, rstatix, 10 other tools, genetics / omics, mouse, cellular / molecular, 1 reference
[4] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, reticulate, 10 other tools, genetics / omics, other condition, cellular / molecular, 1 reference
[5] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, SingleCellExperiment, rstatix, 10 other tools, genetics / omics, mouse, cellular / molecular, 1 reference
[6] doi:10.1038/s41586-026-10214-2 [code]: Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, SingleCellExperiment, reticulate, 10 other tools, genetics / omics, other condition, mouse, 1 reference
[7] doi:10.1038/s41398-026-04200-5 [code]: Postmortem brain single-nucleus and bulk gene expression analyses identify shared and distinct abnormalities in bipolar disorder and major depressive disorder.
Journal: Translational psychiatry
In common: Harmony, SingleCellExperiment, rstatix, 9 other tools, genetics / omics, cellular / molecular, 2 references
[8] doi:10.1016/j.celrep.2026.117500 [code]: Spatio-molecular gene expression reflects dorsal anterior cingulate cortex structure and function in the human brain.
Journal: Cell reports
In common: Harmony, SingleCellExperiment, reticulate, 9 other tools, genetics / omics, cellular / molecular, 2 references
[9] doi:10.1126/sciadv.aeb4265 [code]: Single-nucleus profiling reveals a core disease signature and cell type-specific vulnerabilities in early Rett syndrome.
Journal: Science advances
In common: Harmony, SingleCellExperiment, rstatix, 8 other tools, genetics / omics, other condition, mouse, 1 other category, 2 references
[10] doi:10.1038/s42003-026-10034-0 [code]: Region- and cell type-specific changes in gene expression in the cerebellum after classical fear conditioning.
Journal: Communications biology
In common: SingleCellExperiment, reticulate, rstatix, 9 other tools, genetics / omics, mouse, cellular / molecular, 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.