A scalable human-zebrafish xenotransplantation model reveals gastrosome-mediated processing of dying neurons by human microglia.

Code ↔ Paper

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

[1] § Materials and methods › Data analysis › RNA sequencing ↔ 03_custom/01_custom_analysis_postmeeting.Rmd, lines 104–171 · score 1.00 · online pathway databases, Multi dimensional scaling, leading fold change, GRCh38 assembly, domain expert, unsupervised manner
[2] § Materials and methods › Maintenance of iPSC and generation of fluorescent reporter lines ↔ 00_get_data/00_get_data.sh, lines 53–101 · score 0.70 · pluripotent stem cell, iPSC, Nat, WTC, days, human

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

R Markdown · 2,286 lines · 78 KB · GPL-3.0 · 1 match

---
title: "Ambra Villani's microglia RNA-seq"
author: "Izaskun Mallona, Mark D. Robinson's lab, UZH"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
html_document:
toc: true
toc_float: true
code_folding: hide
code_download: true
number_sections: true
df_print: kable
theme: lumen
params:
seed: 665
---
```{r warning = FALSE, message = FALSE}
## .libPaths('/home/imallona/R/R4_bioc314')
library(ggplot2)
## library(EnhancedVolcano)
## library(iSEE)
## library(knitr)
## library(DRIMSeq)
## library(EnsDb.Hsapiens.v86)
library(SingleCellExperiment)
library(corrplot)
#library(GGally)
#library(scuttle)
#library(knitr)
#library(greybox) # cramer's distance
#library(viridis)
#library(regioneR)
#library(GenomicRanges)
#library(rtracklayer)
## camera start
## library(dplyr)
## library(msigdbr)
## library(CAMERA)
#require(topGO)
#require(org.Hs.eg.db)
## library("ggVennDiagram")
library(recount3)
library(tibble)
library(dplyr)
## library(tidyr)
library(limma)
library(edgeR)
library(sva)
library(DT)
## library(reshape2)
```
```{r edgeR-load-pkg}
suppressPackageStartupMessages({
library(dplyr)
library(tximport)
library(tximeta)
library(SingleCellExperiment)
library(edgeR)
library(ggplot2)
library(msigdbr)
library(EnhancedVolcano)
})
```
```{r}
## render on error
knitr::knit_hooks$set(error = function(x, options) {
knitr::knit_exit()
})
```
```{r}
knitr::opts_chunk$set(fig.width = 5,
fig.height = 5,
cache = TRUE,
include = TRUE,
fig.path = "plots_post/",
dev = c("png", "svg"),
cache.lazy = FALSE,
warning = TRUE,
message = TRUE)
## plan("multiprocess", workers = NTHREADS)
## options(future.globals.maxSize= 1.4e9)
## options(bitmapType='cairo')
```
# Aim
> Compare to published datasets (see block1, 2 and 3). The highest priority,m is to compare to those in block 1 (from Abud et al.).
Methods: We processed our RNA-seq, McQuade's TREM raw reads with ARMOR (PMID 31088905). Briefly, reads were aligned and counted to the human genome (GRCh38 assembly and Gencode release 43) with salmon v1.4.0 and with STAR 2.7.7a. Differential expression analysis was carried out with edgeR v3.36.0 using salmon outputs.
We downloaded Abud's and McQuade's iPSCs data from the Sequence Read Archive (SRA) accessions SRP092075 and SRP155574, respectively, using recount3 (PMID 34844637) and GRCh38 as reference genome and Gencode's annotation.
Results: We detect noticeable age effects (within Villani's WTs), relatively weaker mutation effects (within Villani's mutant), and batch effects when comparing to Abud's, McQuade's data. The MDS plot after batch correction of Villani's WT + Villani's mutants + Abud's looks the best; the same MDS without Villani's mutants not so.
See `Comparison of Villani's, Abud's, McQuade's data`.
> Do they express typical microglial genes (I have a list)
Methods: We generated multi-dimensional scaling (MDS) plots based on normalized count data using edgeR v3.36.0. Briefly, MDS plots depict similarities between samples (including replicates and batches) in an unsupervised manner, depicting the two leading fold-change dimensions which explain the largest proportion of variation in gene expression across samples.
Results: Please use the iSEE deployments:
- http://imlspenticton.uzh.ch:3740/ambra_villani_microglia_only_2023/
- http://imlspenticton.uzh.ch:3740/ambra_villani_microglia_plus_abud_2023/
- http://imlspenticton.uzh.ch:3740/ambra_villani_microglia_plus_abud_and_mcquade_2023/
Mind the iSEE deployments show 'counts', 'CPMs' and 'logCPMs'. My recommendation is to check 'CPMs' or 'logCPMs', so readouts are already normalized/don't depend on sequencing depth. CPMs and logCPMs were calculated with `edgeR::cpm` and with `prior.count = 2`.
> They express genes relevant in neurodegenerative disorders (I have a list)
Results: See the iSEE deployments above.
> 1. The different mutants still becomes microglia (add to the MDS?)
> a. Slc37a2 I have 2 different mutants (clones) same gene (but different mutation)
> b. TREM2 I just have 1 mutant
>2. They are indeed mutants
> a. RNA level reduced compared to WT. As WT for comparison use the 2 from av165 (same round of differentiation)
> b. Alignment to WT sequence?
Results: please browse the STAR BAM files.
> 3. The TREM2 mutant has the typical signature of microglia missing this gene according to literature
> a. Compare to McQuade TREM2 KO (highlighted in yellow in the spreadsheet)
> b. Lipid metabolism affected (expected)
Methods: We ran a geneset enrichment analysis using `camera` from limma v3.50.3 (PMID PMC3458527) and [mSigDB](http://software.broadinstitute.org/gsea/msigdb) annotations on differentially expressed genes. Selected genesets included:
- C2, curated gene sets from online pathway databases, publications in PubMed, and knowledge of domain experts
- C5, Gene Ontology
- C8, cell type signatures
Results: please check sections named `Geneset analysis`.
# Notes
Plots in PNG and SVG format <a href="plots_post/">can be browsed and downloaded here</a>.
# Comparison of Villani's, Abud's, McQuade's data
## Metadata
```{r}
options(ggplot2.discrete.fill = function() scale_fill_brewer(palette = "Set3"))
options(ggplot2.discrete.color = function() scale_color_brewer(palette = "Set3"))
```
```{r}
meta <- read.csv(text="treatment,names,age
iPSC_MG_(Abud),SRR4450428,unknown
iPSC_MG_(Abud),SRR4450429,unknown
iPSC_MG_(Abud),SRR4450430,unknown
iPSC_MG_(Abud),SRR4450431,unknown
iPSC_MG_(Abud),SRR4450432,unknown
iPSC_MG_(Abud),SRR4450433,unknown
iHPC_(Abud),SRR4450434,unknown
iHPC_(Abud),SRR4450435,unknown
iHPC_(Abud),SRR4450436,unknown
iPSC_(Abud),SRR4450437,unknown
iPSC_(Abud),SRR4450438,unknown
iPSC_(Abud),SRR4450439,unknown
iPSC_(Abud),SRR4450440,unknown
CD14_M_(Abud),SRR4450441,unknown
CD14_M_(Abud),SRR4450442,unknown
CD14_M_(Abud),SRR4450443,unknown
CD14_M_(Abud),SRR4450444,unknown
CD14_M_(Abud),SRR4450445,unknown
CD16_M_(Abud),SRR4450446,unknown
CD16_M_(Abud),SRR4450447,unknown
CD16_M_(Abud),SRR4450448,unknown
CD16_M_(Abud),SRR4450449,unknown
Fetal_MG_(Abud),SRR4450450,unknown
Fetal_MG_(Abud),SRR4450451,unknown
Fetal_MG_(Abud),SRR4450452,unknown
Adult_MG_(Abud),SRR4450453,unknown
Adult_MG_(Abud),SRR4450454,unknown
Adult_MG_(Abud),SRR4450455,unknown
Blood_DC_(Abud),SRR4450456,unknown
Blood_DC_(Abud),SRR4450457,unknown
Blood_DC_(Abud),SRR4450458,unknown
iPSC_MG_WT_(McQuade),SRR12608405,unknown
iPSC_MG_WT_(McQuade),SRR12608406,unknown
iPSC_MG_WT_(McQuade),SRR12608407,unknown
iPSC_MG_WT_(McQuade),SRR12608408,unknown
iPSC_MG_TREM2_KO_(McQuade),SRR12608409,unknown
iPSC_MG_TREM2_KO_(McQuade),SRR12608410,unknown
iPSC_MG_TREM2_KO_(McQuade),SRR12608411,unknown
iPSC_MG_TREM2_KO_(McQuade),SRR12608412,unknown
iPSC_MG_2.0_(McQuade),SRR7613924,unknown
iPSC_MG_2.0_(McQuade),SRR7613925,unknown
iPSC_MG_2.0_(McQuade),SRR7613926,unknown
iHPC_2.0_(McQuade),SRR8180288,unknown
iHPC_2.0_(McQuade),SRR8180289,unknown
iHPC_2.0_(McQuade),SRR8180290,unknown
iPSC_MG,20220209.A-av111_WT_001_R1,day_34
iPSC_MG,20220209.A-av111_WT_002_R1,day_34
iPSC_MG,20220209.A-av111_WT_003_R1,day_34
iPSC_MG,20220209.A-tw029_WT_2_R1,day_38
iPSC_MG,20220209.A-tw029_WT_3_R1,day_38
iPSC_MG,20220504.B-WTC_1_R1,day_28
iPSC_MG,20220504.B-WTC_2_R1,day_28
iPSC_MG_Slc37a2_mut1,20220504.B-SLC-EX2B1_1_R1,day_28
iPSC_MG_Slc37a2_mut1,20220504.B-SLC-EX2B1_2_R1,day_28
iPSC_MG_Slc37a2_mut2,20220504.B-SLC-EX6C5_1_R1,day_28
iPSC_MG_Slc37a2_mut2,20220504.B-SLC-EX6C5_2_R1,day_28
iPSC_MG_TREM2_mut,20220504.B-TREM2-A1_1_R1,day_28
iPSC_MG_TREM2_mut,20220504.B-TREM2-A1_2_R1,day_28")
```
```{r}
meta$origin <- 'Villani'
meta$origin[grepl('Abud', meta$treatment)] <- 'Abud'
meta$origin[grepl('McQuade', meta$treatment)] <- 'McQuade iPSC'
meta$origin[grepl('SRR126', meta$names)] <- 'McQuade TREM'
```
Data included in this analysis:
```{r, results = 'asis'}
(DT::datatable(meta, rownames = FALSE,
extensions = 'Buttons',
filter = "none",
options = list(pageLength = 5, autowidth = TRUE,
dom = 'Blftip',
buttons = c('copy', 'csv', 'excel'))))
```
## Overview of Villani's data {.tabset .tabset-pills}
```{r}
d <- readRDS('/home/imallona/avillani_microglia/ambra_only_output/outputR/shiny_sce.rds')$sce_gene
```
```{r}
logcpms <- assay(d, "logcpm")
mds <- limma::plotMDS(logcpms, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE,
var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = colnames(logcpms),
MDS1 = mds$x,
MDS2 = mds$y)
}
mds <- mds %>%
dplyr::full_join(data.frame(colData(d)), by = "names")
```
### MDS {.tabset .tabset-pills}
Multidimensional scaling (MDS) plots samples on a two-dimensional scatterplot so that distances on the plot approximate the typical log2 fold changes between the samples.
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
for (item in c('treatment', 'timepoint', 'names')) {
cat('#### ', item, '\n\n')
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = get(item)), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS", item)) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
labs(color = item) +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
cat('\n\n')
}
```
```{r}
assay(d, 'cpms') <- cpm(d)
saveRDS(d, file = 'ambra_villani_only_sce.rds')
```
## Villani's wildtypes and mutants and Abud's data {.tabset .tabset-pills}
```{r, message = FALSE, warning = FALSE}
srp <- "SRP092075"
abud <- recount3::create_rse_manual(
project = srp,
project_home = "data_sources/sra",
organism = "human",
annotation = "gencode_v29",
type = "gene")
assay(abud, "counts") <- transform_counts(abud)
```
```{r, eval = TRUE}
abud <- subset(abud,
select = colData(abud)$external_id %in% meta$names)
## abud_annot <- merge(abud_annot, colData(abud), by.x = 'srr', by.y = 'external_id')
```
```{r}
abudlcpm <- edgeR::cpm(assay(abud, 'counts'),
log = TRUE,
prior.count = 2)
```
```{r}
## table(rowData(d)$gene_id %in% rowData(abud)$gene_id)
rowData(d)$unversioned_gene_id <- sapply(strsplit(rowData(d)$gene_id, '\\.'), function(x) return(x[[1]]))
rowData(abud)$unversioned_gene_id <- sapply(strsplit(rowData(abud)$gene_id, '\\.'), function(x) return(x[[1]]))
## table(rowData(d)$unversioned_gene_id %in% rowData(abud)$unversioned_gene_id)
shared <- intersect(rowData(d)$unversioned_gene_id, rowData(abud)$unversioned_gene_id)
dict <- data.frame(unversioned_gene_id = shared,
armor = rownames(rowData(d)[rowData(d)$unversioned_gene_id %in% shared,]))
dict <- merge(dict, rowData(abud)[rowData(abud)$unversioned_gene_id %in% shared,
c('unversioned_gene_id', 'gene_id')],
by = 'unversioned_gene_id')
colnames(dict)[3] <- 'recount'
dict <- dict[grep('PAR_Y', dict$recount, invert = TRUE),]
```
```{r}
## remove the PAR_Y recount3 which introduce duplicates, i.e.
## dim(dict)
## print(grep('PAR_Y', dict$recount, value = TRUE))
## dim(dict)
dict <- dict[grep('PAR_Y', dict$recount, invert = TRUE),]
```
### MDS before batch correction
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(cbind(logcpms[dict$armor,],
abudlcpm[dict$recount,]), top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
rownames(mds) <- mds$names
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Abud', no = 'Ambra')
mds <- merge(mds, meta, by = 'names')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
```
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS before ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
```
### Batch correction
```{r}
# parametric adjustment
combat_abud <- sva::ComBat(dat=cbind(logcpms[dict$armor,],
abudlcpm[dict$recount,]),
batch = c(rep('ambra', ncol(logcpms)),
rep('abud', ncol(abudlcpm))),
mod = NULL, par.prior = TRUE, prior.plots = TRUE)
```
### MDS after batch correction {.tabset .tabset-pills}
#### Batches
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(combat_abud, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Abud', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
```
#### Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(fill=guide_legend(ncol=2)) +
theme(aspect.ratio=1) +
scale_color_brewer(palette = "Set3"))
```
```{r}
tmp <- list(logcpms = cbind(logcpms[dict$armor,],
abudlcpm[dict$recount,]),
counts = cbind(counts(d)[dict$armor,],
assay(abud, 'raw_counts')[dict$recount,]))
rownames(meta) <- meta$names
tmp2 <- meta[colnames(tmp$counts),]
tmp$cpms <- edgeR::cpm(tmp$counts, log = FALSE, prior.count = 2)
rownames(mds) <- mds[,1]
mds <- mds[colnames(tmp$counts),]
tmp3 <- SingleCellExperiment(assays = tmp,
colData = data.frame(tmp2),
rowData = rowData(d)[dict$armor,],
metadata = list(study = "Villani, Abud"),
reducedDims = list(mds = mds[,2:3]))
saveRDS(object = tmp3,
file = 'villani_plus_abud_sce.rds')
rm(tmp, tmp2)
```
```{r}
rm(tmp3)
```
## Villani's wildtypes (no mutants) and Abud's data {.tabset .tabset-pills}
```{r}
ambra_wts <- meta[meta$treatment == 'iPSC_MG', 'names']
```
### MDS before batch correction {.tabset .tabset-pills}
```{r fids, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(cbind(logcpms[dict$armor,ambra_wts],
abudlcpm[dict$recount,]), top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Abud', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
for (item in c('treatment', 'age', 'names', 'origin')) {
cat('#### ', item, '\n\n')
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = get(item)), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS before ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
labs(color = item) +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
cat('\n\n')
}
```
### Batch correction
```{r}
# parametric adjustment; explore the parametric too
combat_abud <- sva::ComBat(dat=cbind(logcpms[dict$armor,ambra_wts],
abudlcpm[dict$recount,]),
batch = c(rep('ambra', length(ambra_wts)),
rep('abud', ncol(abudlcpm))),
mod = NULL, par.prior = TRUE, prior.plots = FALSE)
```
### MDS after batch correction {.tabset .tabset-pills}
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(combat_abud, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Abud', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
for (item in c('origin', 'treatment', 'names', 'age')) {
cat('#### ', item, '\n\n')
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = get(item)), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
cat('\n\n')
}
```
#### Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = age), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(fill=guide_legend(ncol=2)) +
theme(aspect.ratio=1) +
scale_color_brewer(palette = "Set3"))
```
## Villani's wildtypes and mutants and McQuade TREM data {.tabset .tabset-pills}
```{r}
trem <- readRDS('/home/imallona/avillani_microglia/mcquade_trem_only_output/outputR/shiny_sce.rds')$sce_gene
trem <- subset(trem,
select = colData(trem)$names %in% meta$names)
tremlcpm <- assay(trem, "logcpm")
```
### MDS before batch correction
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(cbind(logcpms[dict$armor,],
tremlcpm[dict$armor,]), top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'McQuade TREM', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS before ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
```
### Batch correction
```{r}
# parametric adjustment
combat_trem <- sva::ComBat(dat=cbind(logcpms[dict$armor,],
tremlcpm[dict$armor,]),
batch = c(rep('ambra', ncol(logcpms)),
rep('trem', ncol(tremlcpm))),
mod = NULL, par.prior = TRUE, prior.plots = TRUE)
```
### MDS after batch correction {.tabset .tabset-pills}
#### Batch
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(combat_trem, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
# mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Trem', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
```
#### Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%",vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(fill=guide_legend(ncol=2)) +
theme(aspect.ratio=1) +
scale_color_brewer(palette = "Set3"))
```
## Villani's wildtypes and mutants and McQuade's iPSCs data {.tabset .tabset-pills}
```{r, message = FALSE, warning = FALSE}
srp <- 'SRP155574'
mc <- recount3::create_rse_manual(
project = srp,
project_home = "data_sources/sra",
organism = "human",
annotation = "gencode_v29",
type = "gene")
assay(mc, "counts") <- transform_counts(mc)
```
```{r}
mc <- subset(mc,
select = colData(mc)$external_id %in% meta$names)
assay(mc, "counts") <- transform_counts(mc)
```
```{r}
mclcpm <- edgeR::cpm(assay(mc, 'counts'),
log = TRUE,
prior.count = 2)
```
```{r}
## table(rowData(d)$gene_id %in% rowData(mc)$gene_id)
rowData(d)$unversioned_gene_id <- sapply(strsplit(rowData(d)$gene_id, '\\.'), function(x) return(x[[1]]))
rowData(mc)$unversioned_gene_id <- sapply(strsplit(rowData(mc)$gene_id, '\\.'), function(x) return(x[[1]]))
## table(rowData(d)$unversioned_gene_id %in% rowData(mc)$unversioned_gene_id)
shared <- intersect(rowData(d)$unversioned_gene_id, rowData(mc)$unversioned_gene_id)
dict <- data.frame(unversioned_gene_id = shared,
armor = rownames(rowData(d)[rowData(d)$unversioned_gene_id %in% shared,]))
dict <- merge(dict, rowData(mc)[rowData(mc)$unversioned_gene_id %in% shared, c('unversioned_gene_id', 'gene_id')],
by = 'unversioned_gene_id')
colnames(dict)[3] <- 'recount'
dict <- dict[grep('PAR_Y', dict$recount, invert = TRUE),]
```
### MDS before batch correction
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(cbind(logcpms[dict$armor,],
mclcpm[dict$recount,]), top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
#mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Mc', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=0.5) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS before ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio = 1) +
theme(legend.position='right'))
```
### Batch correction
```{r}
# nonparametric adjustment
combat_mc <- sva::ComBat(dat=cbind(logcpms[dict$armor,],
mclcpm[dict$recount,]),
batch = c(rep('ambra', ncol(logcpms)),
rep('mc', ncol(mclcpm))),
mod = NULL, par.prior = TRUE, prior.plots = TRUE)
```
### MDS after batch correction {.tabset .tabset-pills}
#### Batch
```{r, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(combat_mc, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
#mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Mc', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right')) #, legend.box='vertical', legend.text=element_text(size=9), legend.title=element_text(size=11),
#axis.text=element_text(size=10), axis.title=element_text(size=12), title=element_text(size=12)))
```
#### Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%",vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(fill=guide_legend(ncol=2)) +
theme(aspect.ratio=1) +
scale_color_brewer(palette = "Set3"))
```
# Joint embedding of Villani's plus public data, combined {.tabset .tabset-pills}
```{r}
pre <- cbind(logcpms[dict$armor,],
abudlcpm[dict$recount,],
tremlcpm[dict$armor,],
mclcpm[dict$recount,])
```
```{r, fig.width = 30, fig.height = 30}
## dim(cbind(logcpms[dict$armor,], mclcpm[dict$recount,]))
cnames <- colnames(pre)
colnames(pre) <- meta[meta$names == colnames(pre), 'treatment']
## corrplot(cor(pre, method= 'spearman'),
## method = 'circle', type = 'lower', insig = 'blank',
## number.cex = 0.8, order = 'hclust', diag = TRUE , addCoef.col = 'white')
```
```{r}
## take the original colnames back
colnames(pre) <- cnames
```
## MDS before batch correction {.tabset .tabset-pills}
### Batch
```{r firstmds, fig.width = 8, fig.height = 5}
mds <- limma::plotMDS(pre, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y, var.explained = mds$var.explained)
}
## mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'public', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 10, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=1) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS before ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio = 1) +
theme(legend.position='right'))
```
## Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = origin), alpha=0.5) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%",vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(color=guide_legend(ncol=2)) +
theme(aspect.ratio=1))
```
## Batch correction
```{r}
pre <- pre[,meta$names]
stopifnot(all(colnames(pre) == meta$names))
# nonparametric adjustment is as bad
post <- sva::ComBat(dat=pre,
batch = as.factor(meta$origin),
mod = NULL, par.prior = TRUE, prior.plots = TRUE)
```
## MDS after batch correction {.tabset .tabset-pills}
### Batch
```{r, fig.width = 8, fig.height = 5}
rm(mds)
mds <- limma::plotMDS(post, top = 500, labels = NULL, pch = NULL,
cex = 1, dim.plot = c(1, 2), ndim = 7,
gene.selection = "common",
xlab = NULL, ylab = NULL, plot = FALSE, var.explained = TRUE)
vexp1 <- round(mds$var.explained[1]*100)
vexp2 <- round(mds$var.explained[2]*100)
if (!is.null(mds$cmdscale.out)) {
## Bioc 3.12 and earlier
mds <- mds$cmdscale.out
colnames(mds) <- paste0("MDS", seq_len(ncol(mds)))
mds <- as.data.frame(mds) %>% tibble::rownames_to_column(var = "names")
} else {
mds <- data.frame(names = rownames(mds$distance.matrix.squared),
MDS1 = mds$x,
MDS2 = mds$y)
}
# mds$origin <- ifelse(grepl('SRR', mds$names), yes = 'Mc', no = 'Ambra')
## mds <- na.omit( mds %>% dplyr::full_join(data.frame(meta), by = "names"))
mds <- merge(mds, meta, by = 'names')
rownames(mds) <- mds$names
```
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = origin), alpha=0.5) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%", vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),
title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(aspect.ratio=1) +
theme(legend.position='right'))
```
### Annotation
```{r, fig.width = 8, fig.height = 4, results = 'asis'}
print(ggplot(mds, aes(x=MDS1, y=MDS2, color = treatment, shape = origin), alpha=0.5) +
geom_point(size=3) +
labs(x = sprintf("Dim. 1, var. exp. %s%%",vexp1),
y = sprintf("Dim. 2, var. exp. %s%%", vexp2),title= paste("MDS after ComBat")) +
coord_cartesian(xlim=range(mds$MDS1, na.rm=TRUE),
ylim=range(mds$MDS2, na.rm=TRUE), expand=TRUE) +
theme_bw() +
theme(legend.position='right') +
guides(color=guide_legend(ncol=2)) +
theme(aspect.ratio=1))
```
```{r}
pre_counts <- cbind(assay(d, 'counts')[dict$armor,],
assay(abud, 'counts')[dict$recount,],
assay(trem, 'counts')[dict$armor,],
assay(mc, 'counts')[dict$recount,])
pre_counts <- pre_counts[,colnames(pre)]
pre_cpms <- edgeR::cpm(pre_counts,
log = FALSE,
prior.count = 2)
rownames(mds) <- mds[,1]
mds <- mds[colnames(pre_counts),]
saveRDS(object = SingleCellExperiment(list(logcpms = pre,
counts = pre_counts,
cpms = pre_cpms),
colData = data.frame(meta, row.names = meta$names),
rowData = rowData(d)[dict$armor,],
metadata = list(study = "Villani, Abud, McQuade"),
reducedDims = list(mds = data.frame(mds[,2:3], row.names = mds[,1]))),
file = 'all_data_sce.rds')
```
# Differential expression: mutants vs WTs (regardless of the WTs age)
Original ARMOR run, duplicated here to consolidate outputs.
Here, we perform differential gene expression analysis with edgeR
[@Robinson2010edgeR] followed by gene set analysis with camera [@Wu2012camera],
based on abundance estimates from Salmon. For more detailed information of each
step, please refer to the
[edgeR user guide](https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf).
```{r edgeR-print-se}
## List of SummarizedExperiment objects (gene/transcript level)
se <- readRDS('/home/imallona/avillani_microglia/ambra_only_output/outputR/tximeta_se.rds')
## Get gene-level SummarizedExperiment object
sg <- se$sg
metadata <- colData(sg)
sg
```
## Plot total number of reads per sample
```{r edgeR-plot-totalcount}
ggplot(data.frame(totCount = colSums(assay(sg, "counts")),
sample = colnames(assay(sg, "counts")),
stringsAsFactors = FALSE),
aes(x = sample, y = totCount)) + geom_bar(stat = "identity") +
theme_bw() + xlab("") + ylab("Total read count") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
```
## Create DGEList and include average transcript length offsets
A `DGEList` is the main object `edgeR` requires to perform the DGE analysis. It
is designed to store read counts and associated information. After creating this
object, we add offsets, which are average transcript length correction terms
[@Soneson2016tximport],
and scale them so they are consistent with library sizes (sequencing depth for
each sample).
Then we calculate normalization factors to scale the raw library sizes and
minimize the log-fold changes between the samples for most genes. Here the
trimmed mean of M-values between each pair of samples (TMM) is used by default
[@Robinson2010TMM].
Finally we add gene annotation information.
```{r edgeR-dge-generate}
dge0 <- tximeta::makeDGEList(sg)
dge0$genes <- as.data.frame(rowRanges(sg))
```
## Calculate logCPMs and add as an assay
We calculate log-counts per million (CPMs) because they are useful descriptive
measures for the expression level of a gene. Note, however, that the normalized
values are not used for the differential expression analysis. By default, the
normalized library sizes are used in the computation.
We add the logCPMs to one of the fields (or assay) of the first gene-level
`SummarizedExperiment` object `sg`. At the end of the analysis, we will use this
object again to export the results of all the genes we started with.
```{r edgeR-add-logcpm}
logcpms <- edgeR::cpm(dge0, offset = dge0$offset, log = TRUE,
prior.count = 2)
dimnames(logcpms) <- dimnames(dge0$counts)
stopifnot(all(rownames(logcpms) == rownames(sg)),
all(colnames(logcpms) == colnames(sg)))
assay(sg, "logcpm") <- logcpms
```
Next, we specify the design matrix of the experiment, defining which sample
annotations will be taken into account in the statistical modeling.
```{r edgeR-define-design}
design <- '~ 0 + treatment'
metadata <- read.table('/home/imallona/avillani_microglia/ARMOR/metadata_ambra_only.tsv', header = TRUE)
stopifnot(all(colnames(dge0) == metadata$names))
des <- model.matrix(as.formula(design), data = metadata)
```
## Filter out lowly expressed genes
Next we determine which genes have sufficiently large counts to be retained in
the statistical analysis, and remove the rest. After removing genes, we
recalculate the normalization factors.
```{r edgeR-filter-genes}
dim(dge0)
keep <- edgeR::filterByExpr(dge0, design = des)
dge <- dge0[keep, ]
dim(dge)
```
## Estimate dispersion and fit QL model
We model the count data using a quasi-likelihood (QL) negative binomial (NB)
generalized log-linear model, which accounts for gene-specific variability from
both biological and technical sources. Before fitting the model, we estimate
the NB dispersion (overall biological variability across all genes), and the QL
dispersion (gene-specific) using the `estimateDisp()` function.
It is also good practice to look at the relationship between the biological
coefficient of variation (NB dispersion) and the gene abundance (in logCPMs).
```{r edgeR-estimate-disp}
## Estimate dispersion and fit model
dge <- estimateDisp(dge, design = des)
qlfit <- glmQLFit(dge, design = des)
## Plot dispersions
plotBCV(dge)
```
## Define contrasts
Before testing for differences in gene expression, we define the contrasts
we wish to test for. Here we represent the constrasts as a numeric matrix:
```{r edgeR-define-contrasts}
contrast <- c('treatmentiPSC_MG_Slc37a2_mut1-treatmentiPSC_MG',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG',
'treatmentiPSC_MG_TREM2_mut-treatmentiPSC_MG',
'treatmentiPSC_MG_Slc37a2_mut1-treatmentiPSC_MG_TREM2_mut',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG_TREM2_mut',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG_Slc37a2_mut1')
(contrasts <- as.data.frame(makeContrasts(contrasts = contrast, levels = des)))
```
```{r edgeR-perform-tests}
signif3 <- function(x) signif(x, digits = 3)
edgeR_res <- lapply(contrasts, function(cm) {
qlf <- glmQLFTest(qlfit, contrast = cm)
tt <- topTags(qlf, n = Inf, sort.by = "none")$table
tt %>%
dplyr::mutate(mlog10PValue = -log10(PValue)) %>%
dplyr::mutate_at(vars(one_of(c("logFC", "logCPM", "F",
"PValue", "FDR", "mlog10PValue"))),
list(signif3))
})
```
## DGE tests and MA plots {.tabset .tabset-pills}
We can visualize the test results by plotting the logCPM (average) vs the logFC,
and coloring genes with an adjusted p-value below 0.05 (or another specificed
FDR threshold). A plot is drawn for every contrast.
```{r, results = 'asis'}
if (is(edgeR_res, "data.frame")) {
print(ggplot(edgeR_res, aes(x = logCPM, y = logFC, color = FDR <= 0.05)) +
geom_point() + theme_bw() +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black")))
} else {
for (nm in names(edgeR_res)) {
cat('### ', nm, '\n\n')
print(ggplot(edgeR_res[[nm]], aes(x = logCPM, y = logFC, color = FDR <= 0.05)) +
geom_point() + theme_bw() +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
ggtitle(nm))
cat('\n\n')
}
}
```
## Explore DGE results {.tabset .tabset-pills}
We export the results into text files that can be opened using any text editor.
The volcanoes plot the uncorrected p-values, but the coloring (significance) is based on FDR-adjusted p-values.
Independently of the FDRs, genes with an absolute log2 fold-change greater than 0.5 are also reported (coloring and vertical lines will at the negative and positive values of the FC cut-off).
```{r edgeR-save-results, results = 'asis', fig.width = 8, fig.height = 8}
## Write results to text files and make MA plots
if (is(edgeR_res, "data.frame")) {
write.table(edgeR_res %>% dplyr::arrange(PValue) %>%
dplyr::select(-dplyr::any_of("tx_ids")),
file = "edgeR_dge_results.txt",
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
} else {
for (nm in names(edgeR_res)) {
cat('### ', nm, '\n\n')
edgeR_res[[nm]]$simple_gene_name <- edgeR_res[[nm]]$gene_name
edgeR_res[[nm]]$simple_gene_name[is.na(edgeR_res[[nm]]$simple_gene_name)] <- edgeR_res[[nm]]$gene_id[is.na(edgeR_res[[nm]]$simple_gene_name)]
print(EnhancedVolcano(edgeR_res[[nm]],
lab = edgeR_res[[nm]]$simple_gene_name,
x = 'logFC',
y = 'PValue',
pCutoffCol = 'FDR',
title = 'differential gene expression',
subtitle = nm,
pCutoff = 0.05,
FCcutoff = 0.5,
pointSize = 3.0,
drawConnectors = TRUE,
labSize = 6.0))
fn <- paste0("edgeR_dge_results_", nm, "_all_ages.txt")
write.table(edgeR_res[[nm]] %>% dplyr::arrange(PValue) %>%
dplyr::select(-dplyr::any_of("tx_ids")),
file = fn,
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
cat(sprintf('\n\n<a href="%s">Download results</a>\n\n', fn))
## (DT::datatable(edgeR_res[[nm]], rownames = FALSE,
## extensions = 'Buttons',
## filter = "none",
## options = list(pageLength = 5, autowidth = TRUE,
## dom = 'Blftip',
## buttons = c('copy', 'csv', 'excel'))))
cat('\n\n')
## write.table(edgeR_res[[nm]] %>% dplyr::arrange(PValue) %>%
## dplyr::select(-dplyr::any_of("tx_ids")),
## file = paste0("edgeR_dge_results_", nm, ".txt"),
## sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
}
}
```
## Geneset analysis {.tabset .tabset-pills}
We will use `camera` to perform an enrichment analysis for a collection of
gene sets from the [mSigDB](http://software.broadinstitute.org/gsea/msigdb),
packaged in the `msigdbr` R package. Here, we load the gene set definitions
and select which ones to include in the analysis. Genesets included:
- C2, curated gene sets from online pathway databases, publications in PubMed, and knowledge of domain experts
- C5, Gene Ontology
- C8, cell type signatures
```{r camera-load-genesets, eval = !is.null(genesets), include = !is.null(genesets)}
## genesets <- 'H,C1,C2,C3,C5,C8'
genesets <- 'C2,C5,C8'
genesets <- strsplit(gsub(" ","",genesets), ",")[[1]]
organism <- 'human'
## Retrieve gene sets and combine in a tibble
m_df <- bind_rows(lapply(genesets,
function(x) msigdbr(species = organism, category = x)))
```
```{r camera-text2, echo= FALSE, results = 'asis', eval = !is.null(genesets)}
cat("We consider only gene sets where the
number of genes shared with the data set is not too small and not too large.
`camera` is a competitive gene set test that accounts for correlations among
the genes within a gene set.")
```
```{r camera-filter-gene-sets, eval = !is.null(genesets), include = !is.null(genesets)}
minSize <- 3
maxSize <- 500
## Get index for genes in each gene set in the DGEList
indexList <- limma::ids2indices(
gene.sets = lapply(split(m_df, f = m_df$gs_name), function(w) w$gene_symbol),
identifiers = dge$genes$gene_name,
remove.empty = TRUE
)
## Filter out too small or too large gene sets
gsSizes <- vapply(indexList, length, 0)
indexList <- indexList[gsSizes >= minSize & gsSizes <= maxSize]
```
```{r camera-check-indexList-length, eval = !is.null(genesets), include = FALSE}
## Check if the index list is empty after filtering
if (length(indexList) == 0){
genesets <- NULL
empty <- TRUE
} else {
empty <- FALSE
}
```
```{r, echo = FALSE, results = 'asis', eval = !is.null(genesets) && empty}
cat("**NOTE:**
The index list is empty after filtering and `camera` cannot be run. Either try
different gene categories, try different filtering parameters or disable the
gene set analysis in the `config.yaml` file by setting `run_camera: False`.")
```
```{r, eval = !is.null(genesets), include = !is.null(genesets)}
camera_res <- lapply(contrasts, function(cm) {
camera(dge, index = indexList, design = des, contrast = cm,
inter.gene.cor = NA)
})
```
```{r camera-save-results, results = 'asis', eval = !is.null(genesets), include = !is.null(genesets)}
## Write results to text files
if (is(camera_res, "data.frame")) {
## write.table(camera_res %>% tibble::rownames_to_column("GeneSet") %>%
## dplyr::arrange(PValue),
## file = "camera_dge_results.txt",
## sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
} else {
for (nm in names(camera_res)) {
cat('### ', nm, '\n\n')
cat('Significantly (before multiple testing correction, pvalue 0.05) enriched/depleted genesets for this contrast:')
print(sum(camera_res[[nm]]$PValue < 0.05))
cat('\n\n')
cat('Significantly (FDR < 0.05) enriched/depleted genesets for this contrast:')
print(sum(camera_res[[nm]]$FDR < 0.05))
cat('\n\n')
fn <- paste0("camera_dge_results_", nm, "_all_ages.txt")
write.table(camera_res[[nm]] %>%
tibble::rownames_to_column("GeneSet") %>%
dplyr::arrange(PValue),
file = fn,
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
cat(sprintf('\n\n<a href="%s">Download geneset analysis results</a>\n\n', fn))
cat('\n\n')
}
}
```
```{r camera-save-se, eval = !is.null(genesets), include = !is.null(genesets)}
geneSets <- lapply(indexList, function(i) dge$genes$gene_name[i])
saveRDS(list(cameraRes = camera_res,
geneSets = geneSets), file = "camera_gsa_all_ages.rds")
```
# Differential expression: mutants vs 28-days-old WTs
```{r}
sg <- sg[,colData(sg)$timepoint == 28]
metadata <- colData(sg)
```
## Plot total number of reads per sample
```{r fig.width = 4, fig.height = 4}
ggplot(data.frame(totCount = colSums(assay(sg, "counts")),
sample = colnames(assay(sg, "counts")),
stringsAsFactors = FALSE),
aes(x = sample, y = totCount)) + geom_bar(stat = "identity") +
theme_bw() + xlab("") + ylab("Total read count") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
```
## Create DGEList and include average transcript length offsets
A `DGEList` is the main object `edgeR` requires to perform the DGE analysis. It
is designed to store read counts and associated information. After creating this
object, we add offsets, which are average transcript length correction terms
[@Soneson2016tximport],
and scale them so they are consistent with library sizes (sequencing depth for
each sample).
Then we calculate normalization factors to scale the raw library sizes and
minimize the log-fold changes between the samples for most genes. Here the
trimmed mean of M-values between each pair of samples (TMM) is used by default
[@Robinson2010TMM].
Finally we add gene annotation information.
```{r}
dge0 <- tximeta::makeDGEList(sg)
dge0$genes <- as.data.frame(rowRanges(sg))
```
## Calculate logCPMs and add as an assay
We calculate log-counts per million (CPMs) because they are useful descriptive
measures for the expression level of a gene. Note, however, that the normalized
values are not used for the differential expression analysis. By default, the
normalized library sizes are used in the computation.
We add the logCPMs to one of the fields (or assay) of the first gene-level
`SummarizedExperiment` object `sg`. At the end of the analysis, we will use this
object again to export the results of all the genes we started with.
```{r}
logcpms <- edgeR::cpm(dge0, offset = dge0$offset, log = TRUE,
prior.count = 2)
dimnames(logcpms) <- dimnames(dge0$counts)
stopifnot(all(rownames(logcpms) == rownames(sg)),
all(colnames(logcpms) == colnames(sg)))
assay(sg, "logcpm") <- logcpms
```
Next, we specify the design matrix of the experiment, defining which sample
annotations will be taken into account in the statistical modeling.
```{r}
design <- '~ 0 + treatment'
## metadata <- read.table('/home/imallona/avillani_microglia/ARMOR/metadata_ambra_only.tsv', header = TRUE)
stopifnot(all(colnames(dge0) == metadata$names))
des <- model.matrix(as.formula(design), data = metadata)
```
## Filter out lowly expressed genes
Next we determine which genes have sufficiently large counts to be retained in
the statistical analysis, and remove the rest. After removing genes, we
recalculate the normalization factors.
```{r}
dim(dge0)
keep <- edgeR::filterByExpr(dge0, design = des)
dge <- dge0[keep, ]
dim(dge)
```
## Estimate dispersion and fit QL model
We model the count data using a quasi-likelihood (QL) negative binomial (NB)
generalized log-linear model, which accounts for gene-specific variability from
both biological and technical sources. Before fitting the model, we estimate
the NB dispersion (overall biological variability across all genes), and the QL
dispersion (gene-specific) using the `estimateDisp()` function.
It is also good practice to look at the relationship between the biological
coefficient of variation (NB dispersion) and the gene abundance (in logCPMs).
```{r}
## Estimate dispersion and fit model
dge <- estimateDisp(dge, design = des)
qlfit <- glmQLFit(dge, design = des)
## Plot dispersions
plotBCV(dge)
```
## Define contrasts
Before testing for differences in gene expression, we define the contrasts
we wish to test for. Here we represent the constrasts as a numeric matrix:
```{r}
contrast <- c('treatmentiPSC_MG_Slc37a2_mut1-treatmentiPSC_MG',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG',
'treatmentiPSC_MG_TREM2_mut-treatmentiPSC_MG',
'treatmentiPSC_MG_Slc37a2_mut1-treatmentiPSC_MG_TREM2_mut',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG_TREM2_mut',
'treatmentiPSC_MG_Slc37a2_mut2-treatmentiPSC_MG_Slc37a2_mut1')
(contrasts <- as.data.frame(makeContrasts(contrasts = contrast, levels = des)))
```
```{r}
signif3 <- function(x) signif(x, digits = 3)
edgeR_res <- lapply(contrasts, function(cm) {
qlf <- glmQLFTest(qlfit, contrast = cm)
tt <- topTags(qlf, n = Inf, sort.by = "none")$table
tt %>%
dplyr::mutate(mlog10PValue = -log10(PValue)) %>%
dplyr::mutate_at(vars(one_of(c("logFC", "logCPM", "F",
"PValue", "FDR", "mlog10PValue"))),
list(signif3))
})
```
## DGE tests and MA plots {.tabset .tabset-pills}
We can visualize the test results by plotting the logCPM (average) vs the logFC,
and coloring genes with an adjusted p-value below 0.05 (or another specificed
FDR threshold). A plot is drawn for every contrast.
```{r, results = 'asis'}
if (is(edgeR_res, "data.frame")) {
print(ggplot(edgeR_res, aes(x = logCPM, y = logFC, color = FDR <= 0.05)) +
geom_point() + theme_bw() +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black")))
} else {
for (nm in names(edgeR_res)) {
cat('### ', nm, '\n\n')
print(ggplot(edgeR_res[[nm]], aes(x = logCPM, y = logFC, color = FDR <= 0.05)) +
geom_point() + theme_bw() +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
ggtitle(nm))
cat('\n\n')
}
}
```
## Explore DGE results {.tabset .tabset-pills}
We export the results into text files that can be opened using any text editor.
The volcanoes plot the uncorrected p-values, but the coloring (significance) is based on FDR-adjusted p-values.
Independently of the FDRs, genes with an absolute log2 fold-change greater than 0.5 are also reported (coloring and vertical lines will at the negative and positive values of the FC cut-off).
```{r, results = 'asis', fig.width = 8, fig.height = 8}
## Write results to text files and make MA plots
if (is(edgeR_res, "data.frame")) {
write.table(edgeR_res %>% dplyr::arrange(PValue) %>%
dplyr::select(-dplyr::any_of("tx_ids")),
file = "edgeR_dge_results.txt",
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
} else {
for (nm in names(edgeR_res)) {
cat('### ', nm, '\n\n')
edgeR_res[[nm]]$simple_gene_name <- edgeR_res[[nm]]$gene_name
edgeR_res[[nm]]$simple_gene_name[is.na(edgeR_res[[nm]]$simple_gene_name)] <- edgeR_res[[nm]]$gene_id[is.na(edgeR_res[[nm]]$simple_gene_name)]
print(EnhancedVolcano(edgeR_res[[nm]],
lab = edgeR_res[[nm]]$simple_gene_name,
x = 'logFC',
y = 'PValue',
pCutoffCol = 'FDR',
title = 'differential gene expression',
subtitle = nm,
pCutoff = 0.05,
FCcutoff = 0.5,
drawConnectors = TRUE,
pointSize = 3.0,
labSize = 6.0))
fn <- paste0("edgeR_dge_results_", nm, "_28_days_only.txt")
write.table(edgeR_res[[nm]] %>% dplyr::arrange(PValue) %>%
dplyr::select(-dplyr::any_of("tx_ids")),
file = fn,
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
cat(sprintf('\n\n<a href="%s">Download results</a>\n\n', fn))
## (DT::datatable(edgeR_res[[nm]], rownames = FALSE,
## extensions = 'Buttons',
## filter = "none",
## options = list(pageLength = 5, autowidth = TRUE,
## dom = 'Blftip',
## buttons = c('copy', 'csv', 'excel'))))
cat('\n\n')
## write.table(edgeR_res[[nm]] %>% dplyr::arrange(PValue) %>%
## dplyr::select(-dplyr::any_of("tx_ids")),
## file = paste0("edgeR_dge_results_", nm, ".txt"),
## sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
}
}
```
## Geneset analysis {.tabset .tabset-pills}
```{r, eval = !is.null(genesets), include = !is.null(genesets)}
## genesets <- 'H,C1,C2,C3,C5,C8'
genesets <- 'C2,C5,C8'
genesets <- strsplit(gsub(" ","",genesets), ",")[[1]]
organism <- 'human'
## Retrieve gene sets and combine in a tibble
m_df <- bind_rows(lapply(genesets,
function(x) msigdbr(species = organism, category = x)))
minSize <- 3
maxSize <- 500
## Get index for genes in each gene set in the DGEList
indexList <- limma::ids2indices(
gene.sets = lapply(split(m_df, f = m_df$gs_name), function(w) w$gene_symbol),
identifiers = dge$genes$gene_name,
remove.empty = TRUE
)
## Filter out too small or too large gene sets
gsSizes <- vapply(indexList, length, 0)
indexList <- indexList[gsSizes >= minSize & gsSizes <= maxSize]
```
```{r, eval = !is.null(genesets), include = FALSE}
## Check if the index list is empty after filtering
if (length(indexList) == 0){
genesets <- NULL
empty <- TRUE
} else {
empty <- FALSE
}
```
```{r, eval = !is.null(genesets), include = !is.null(genesets)}
geneSets <- lapply(indexList, function(i) dge$genes$gene_name[i])
saveRDS(list(cameraRes = camera_res,
geneSets = geneSets), file = "camera_gsa_28_days_only.rds")
```
```{r camera-perform-tests, cache = FALSE, eval = !is.null(genesets), include = !is.null(genesets)}
camera_res <- lapply(contrasts, function(cm) {
camera(dge, index = indexList, design = des, contrast = cm,
inter.gene.cor = NA)
})
```
```{r, results = 'asis', eval = !is.null(genesets), include = !is.null(genesets)}
## Write results to text files
if (is(camera_res, "data.frame")) {
## write.table(camera_res %>% tibble::rownames_to_column("GeneSet") %>%
## dplyr::arrange(PValue),
## file = "camera_dge_results.txt",
## sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
} else {
for (nm in names(camera_res)) {
cat('### ', nm, '\n\n')
cat('Significantly (before multiple testing correction, pvalue 0.05) enriched/depleted genesets for this contrast:')
print(sum(camera_res[[nm]]$PValue < 0.05))
cat('\n\n')
cat('Significantly (FDR < 0.05) enriched/depleted genesets for this contrast:')
print(sum(camera_res[[nm]]$FDR < 0.05))
cat('\n\n')
fn <- paste0("camera_dge_results_", nm, "_28_days_only.txt")
write.table(camera_res[[nm]] %>%
tibble::rownames_to_column("GeneSet") %>%
dplyr::arrange(PValue),
file = fn,
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
cat(sprintf('\n\n<a href="%s">Download geneset analysis results</a>\n\n', fn))
cat('\n\n')
}
}
```
# Consistency check
```{r}
a0 <- readRDS('ambra_villani_only_sce.rds')
a1 <- readRDS('villani_plus_abud_sce.rds')
a2 <- readRDS('all_data_sce.rds')
set.seed(54)
genes <- sample(rownames(a0), 10, replace = FALSE)
id <- '20220504.B-TREM2-A1_2_R1'
stopifnot(all(assay(a0, 'logcpm')[genes,id] == assay(a1, 'logcpms')[genes,id]) )
stopifnot(all(assay(a0, 'logcpm')[genes,id] == assay(a2, 'logcpms')[genes,id]) )
stopifnot(all(assay(a0, 'counts')[genes,id] == assay(a1, 'counts')[genes,id]) )
stopifnot(all(assay(a0, 'counts')[genes,id] == assay(a2, 'counts')[genes,id]) )
## stopifnot(all(assay(a0, 'cpms')[genes,id] == assay(a1, 'cpms')[genes,id]) )
## stopifnot(all(assay(a0, 'cpms')[genes,id] == assay(a2, 'cpms')[genes,id]))
print('passed')
```
# Plot custom logCPMs
# Timestamp
```{r sessionInfo2, cache = FALSE}
date()
sessionInfo()
## devtools::session_info()
```

01_custom_analysis_postmeeting.Rmd, under GPL-3.0 · at the source

Overview

Authors: Ambra Villani¹, Jana Wittmann¹, Tamara Wyss¹, Izaskun Mallona^1,2, Irene Santisteban Ortiz^3,4, Nathalie Tichy¹, Corinna Maria Biermeier¹, Monique Pena^3,4, Ayush Aditya Pal¹, Darren Gilmour¹, Simon T Schafer^3,4, Francesca Peri¹

ORCID iDs: Ambra Villani, Jana Wittmann, Izaskun Mallona, Corinna Maria Biermeier, Ayush Aditya Pal, Darren Gilmour, Simon T Schafer, Francesca Peri

Department of Molecular Life Sciences, University of Zurich, Zurich, Switzerland
SIB, Swiss Institute of Bioinformatics, Zurich, Switzerland
Department of Psychiatry and Psychotherapy, School of Medicine and Health, Technical University of Munich, Munich, Germany
Center for Organoid Systems, Munich Institute for Biomedical Engineering, Technical University of Munich, Garching, Germany

Institutions: University of Zurich (Switzerland); SIB Swiss Institute of Bioinformatics (Switzerland); Technical University of Munich (Germany)

Journal: Communications biology, volume 9, issue 1, article 785

Dates: received 15 August 2025; accepted 17 March 2026; published online 9 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s42003-026-09948-6 · PMID 41957412 · PMCID PMC13250125 · OpenAlex W4413855248

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.08.27.672578

Status: code verified

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

Methods: Statistics, Evoked potentials, fMRI & imaging

Keywords: Cell death and immune response, Apoptosis

MeSH: Microglia*, Neurons*, Zebrafish*, Animals, Apoptosis, Brain, Efferocytosis, Humans, Induced Pluripotent Stem Cells, Transplantation, Heterologous (* major topic)

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

Funding: Swiss National Science Foundation (212794, 310030_204834, 310030, 310030_212794)

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

Research resources: RRID:SCR_007370

Abstract

Microglia engulf dying neurons through efferocytosis, a critical function in both development and disease. How microglia process the engulfed neuronal material—especially lipids—remains poorly understood, despite its central role in neurodegeneration. Thus, we developed HuZIBRA, a scalable in vivo xenotransplantation model in which human iPSC-derived microglia-like cells (iMGLs) are introduced into the developing zebrafish brain (zf-hiMG), a system characterized by high levels of neuronal cell death and amenable to precise genetic and pharmacological manipulation. We show that human microglia-like cells recognize and engulf apoptotic zebrafish neurons, indicating conserved efferocytic mechanisms. In these cells, engulfed neuronal material accumulates into a distinct, lipid-rich intracellular compartment, the gastrosome, which we also observed in iMGLs placed in a human brain-like environment. The size of the human gastrosome dynamically reflects neuronal cell death levels and is regulated by key genes, including TREM2 and SLC37A2. Pharmacological inhibition of the cholesterol transporter NPC1 induces gastrosome expansion and lipid accumulation, recapitulating pathological features of Niemann-Pick disease type C. Thus, HuZIBRA provides a powerful in vivo platform to uncover cell-autonomous adaptive responses of human microglia to apoptotic and metabolic stress, with the gastrosome emerging as a key integrator of neuronal debris processing and disease-relevant lipid metabolism.

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

Repository

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

Zenodo 18997850

License: GPL-3.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: “Data availability”

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

Tools: ggplot2 (9 files), SingleCellExperiment (9 files), ComplexHeatmap (6 files), circlize (3 files), limma (3 files), tidyverse (3 files), edgeR (2 files), Snakemake (1 file)

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

29 September 2026: the link answers (HTTP 200)

16 files

00_get_data/00_get_data.sh, Shell, 179 lines, 1 match
00_get_data/01_get_reference_data.sh, Shell, 17 lines
01_config/00_write_armor_config_files.sh, Shell, 390 lines
02_run_armor/00_run_armor.sh, Shell, 28 lines
03_custom/00_custom_analysis_premeeting.Rmd, R, 950 lines
03_custom/01_custom_analysis_postmeeting.Rmd, R, 2,286 lines, 1 match
03_custom/02_minor_plots.Rmd, R, 166 lines
04_further_custom/00_minor_plots.Rmd, R, 340 lines
04_further_custom/01_custom_heatmap.Rmd, R, 683 lines
04_further_custom/02_other_plots.Rmd, R, 230 lines
04_further_custom/03_browsing_old_edgeR_results.sh, Shell, 50 lines
05_without_tmem119/01_custom_heatmap_without_tmem.Rmd, R, 693 lines
05_without_tmem119/02_other_plots_without_tmem.Rmd, R, 239 lines
06_abud_vs_villani_distances/01_mds_and_corrplot.Rmd, R, 361 lines
LICENSE, License, 674 lines
README.md, Text, 68 lines

At the source: zenodo.org/records/18997850

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

Tracing map

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

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
14 scripts, each with its path and the digest of its content;
2 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:GSE324751, at NCBI GEO; found in “Data availability”

Data availability

Numerical data related to this manuscript are provided in Supplementary Data 1. The imaging data generated in this study are available from the corresponding author on reasonable request. Raw and processed RNA-seq data generated in this study are available at GEO GSE324751 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE324751). Reanalyzed data were downloaded from SRA accessions SRP092075 and SRP155574. Code to analyze the RNA-seq data is available at https://zenodo.org/records/18997850 under the GPLv3 terms.

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

Versions

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

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 2 keywords, 10 MeSH terms, 1 funder, 70 references, 1 RRID.

Cite

This paper

Villani, A., Wittmann, J., Wyss, T., Mallona, I., Santisteban Ortiz, I., Tichy, N., Biermeier, C. M., Pena, M., Pal, A. A., Gilmour, D., Schafer, S. T., & Peri, F. (2026). A scalable human-zebrafish xenotransplantation model reveals gastrosome-mediated processing of dying neurons by human microglia. Communications biology, 9(1), 785. https://doi.org/10.1038/s42003-026-09948-6

BibTeX

@article{villani2026scalable, author = {Villani, Ambra and Wittmann, Jana and Wyss, Tamara and Mallona, Izaskun and Santisteban Ortiz, Irene and Tichy, Nathalie and Biermeier, Corinna Maria and Pena, Monique and Pal, Ayush Aditya and Gilmour, Darren and Schafer, Simon T and Peri, Francesca}, title = {{A scalable human-zebrafish xenotransplantation model reveals gastrosome-mediated processing of dying neurons by human microglia}}, journal = {Communications biology}, year = {2026}, month = apr, volume = {9}, number = {1}, pages = {785}, publisher = {Nature Publishing Group}, issn = {2399-3642}, doi = {10.1038/s42003-026-09948-6}, url = {https://doi.org/10.1038/s42003-026-09948-6}, pmid = {41957412}, pmcid = {PMC13250125} }

RIS

TY - JOUR AU - Villani, Ambra AU - Wittmann, Jana AU - Wyss, Tamara AU - Mallona, Izaskun AU - Santisteban Ortiz, Irene AU - Tichy, Nathalie AU - Biermeier, Corinna Maria AU - Pena, Monique AU - Pal, Ayush Aditya AU - Gilmour, Darren AU - Schafer, Simon T AU - Peri, Francesca TI - A scalable human-zebrafish xenotransplantation model reveals gastrosome-mediated processing of dying neurons by human microglia T2 - Communications biology J2 - Commun Biol PY - 2026 DA - 2026/04/09 VL - 9 IS - 1 SP - 785 SN - 2399-3642 PB - Nature Publishing Group DO - 10.1038/s42003-026-09948-6 UR - https://doi.org/10.1038/s42003-026-09948-6 LA - en ER -

CSL-JSON

{ "id": "10.1038/s42003-026-09948-6", "type": "article-journal", "title": "A scalable human-zebrafish xenotransplantation model reveals gastrosome-mediated processing of dying neurons by human microglia", "container-title": "Communications biology", "author": [ { "family": "Villani", "given": "Ambra" }, { "family": "Wittmann", "given": "Jana" }, { "family": "Wyss", "given": "Tamara" }, { "family": "Mallona", "given": "Izaskun" }, { "family": "Santisteban Ortiz", "given": "Irene" }, { "family": "Tichy", "given": "Nathalie" }, { "family": "Biermeier", "given": "Corinna Maria" }, { "family": "Pena", "given": "Monique" }, { "family": "Pal", "given": "Ayush Aditya" }, { "family": "Gilmour", "given": "Darren" }, { "family": "Schafer", "given": "Simon T" }, { "family": "Peri", "given": "Francesca" } ], "container-title-short": "Commun Biol", "volume": "9", "issue": "1", "page": "785", "DOI": "10.1038/s42003-026-09948-6", "PMID": "41957412", "PMCID": "PMC13250125", "ISSN": "2399-3642", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s42003-026-09948-6", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 9 ] ] } }

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: SingleCellExperiment, edgeR, limma, 4 other tools, cellular / molecular, 7 references
[2] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: SingleCellExperiment, edgeR, limma, 4 other tools, cellular / molecular, 3 references
[3] 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: SingleCellExperiment, edgeR, limma, 4 other tools, cellular / molecular, 1 reference
[4] doi:10.1016/j.neuron.2026.07.007 [code]: Human-specific SRGAP2 paralogs synchronize neotenic microglial maturation and synaptic development.
Journal: Neuron
In common: ggplot2, tidyverse, 7 references
[5] doi:10.1038/s41386-026-02406-1 [code]: Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators.
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
In common: limma, circlize, ComplexHeatmap, 2 other tools, cellular / molecular, 3 references
[6] doi:10.1038/s41467-026-76341-6 [code]: Neonatal inflammation disrupts a temporally restricted postnatal Numb-enriched microglial state in mice.
Journal: Nature communications
In common: SingleCellExperiment, circlize, ComplexHeatmap, 2 other tools, 3 references
[7] doi:10.1038/s41467-026-71790-5 [code]: Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification.
Journal: Nature communications
In common: Snakemake, edgeR, circlize, 3 other tools, cellular / molecular, 1 reference
[8] doi:10.1038/s41467-026-70989-w [code]: MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos.
Journal: Nature communications
In common: SingleCellExperiment, edgeR, circlize, 3 other tools, zebrafish, 1 reference
[9] doi:10.1111/adb.70179 [code]: Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala.
Journal: Addiction biology
In common: SingleCellExperiment, edgeR, limma, 4 other tools, cellular / molecular
[10] doi:10.1016/j.isci.2026.115573 [code]: Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: SingleCellExperiment, edgeR, limma, 4 other tools, cellular / molecular

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.