OSCR

Zika virus infections of human stem cell-derived cerebral organoids reveal viral lineage-specific pathogenesis responses.

Code ↔ Paper

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

The 4 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § MATERIALS AND METHODS › Bulk RNA isolation and sequencing ↔ BulkRNA-seq_umi/nf-core_rnaseq/nf-core_rnaseq_umi.sh, the whole file · a weak match · score 0.75 · nf core, plus virus, gencode, hg38, MIT, genomes
  2. [2] § RESULTS › ZIKV-induced growth defect severity correlates with cytoarchitectural disruption ↔ scRNA-seq/020222_Gehrke.Rmd, lines 450–458 · score 0.60 · intermediate progenitor cells, neuronal progenitor cells, astrocytes, radial, Neurons
  3. [3] § RESULTS › Bulk RNA-seq reveals differential stress responses among ZIKV lineage infections ↔ scRNA-seq/020222_Gehrke.Rmd, lines 1028–1080 · score 0.59 · Log2 FC, stress response, fold change, Volcano, adj, cutoffs
  4. [4] § RESULTS › Bulk RNA-seq reveals differential stress responses among ZIKV lineage infections ↔ scRNA-seq/020222_Gehrke.Rmd, lines 1028–1080 · score 0.51 · log2 FC, stress response, overlaps, adj, folded, UG

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 authors' code

R Markdown · 1,145 lines · 44 KB · MIT · 3 matches

  1. ---
  2. title: "scRNA-seq processing"
  3. author: "Charlie Whittaker and Yann Vanrobaeys"
  4. date: "2024-07-08"
  5. output:
  6. html_document:
  7. toc: true
  8. toc_depth: 3
  9. html_notebook: default
  10. ---
  12. ```{r setup, include=FALSE}
  13. knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
  14. dev="png", message=FALSE, error=FALSE, warning=TRUE)
  15. ```
  17. ```{r}
  18. library(Seurat)
  19. library(dplyr)
  20. library(ggplot2)
  21. library(openxlsx)
  22. library(SingleR)
  23. library(ggpubr)
  24. library(readxl)
  25. library(writexl)
  27. library(openxlsx)
  28. library(tidyverse)
  29. library(reprex)
  30. library(matrixStats)
  31. library(XML)
  32. library(ggrepel)
  33. library(rtracklayer)
  34. library(DESeq2)
  35. library(apeglm)
  36. library(ComplexHeatmap)
  37. library(edgeR)
  38. library(gprofiler2)
  39. library(cluster)
  40. library(stringr)
  41. library(fgsea)
  42. ```
  44. # Setup the Seurat Object
  45. ## Load the CellRanger outputs of individual sample
  47. ```{r}
  48. d9433_d7_mock.data <- Read10X(data.dir = "CellRanger_filtered_feature_bc_matrix_dir/d9433_d7_mock/filtered_feature_bc_matrix")
  49. d9433_d7_pr.data <- Read10X(data.dir = "CellRanger_filtered_feature_bc_matrix_dir/d9434_d7_pr/filtered_feature_bc_matrix")
  50. d9433_d7_mal.data <- Read10X(data.dir = "CellRanger_filtered_feature_bc_matrix_dir/d9435_d7_mal/filtered_feature_bc_matrix")
  51. d9433_d7_ug.data <- Read10X(data.dir = "CellRanger_filtered_feature_bc_matrix_dir/d9436_d7_ug/filtered_feature_bc_matrix")
  52. ```
  54. ## Initialize the Seurat object with the raw (non-normalized data).
  56. ```{r}
  57. d9433_d7_mock <- CreateSeuratObject(counts = d9433_d7_mock.data, project = "d9433_d7_mock", min.cells = 3, min.features = 200)
  58. d9434_d7_pr <- CreateSeuratObject(counts = d9433_d7_pr.data, project = "d9433_d7_pr", min.cells = 3, min.features = 200)
  59. d9435_d7_mal <- CreateSeuratObject(counts = d9433_d7_mal.data, project = "d9433_d7_mal", min.cells = 3, min.features = 200)
  60. d9436_d7_ug <- CreateSeuratObject(counts = d9433_d7_ug.data, project = "d9433_d7_ug", min.cells = 3, min.features = 200)
  62. d9433_d7_mock[["Condition"]] <- "Mock"
  63. d9434_d7_pr[["Condition"]] <- "Puerto Rico"
  64. d9435_d7_mal[["Condition"]] <- "Malaysia"
  65. d9436_d7_ug[["Condition"]] <- "Uganda"
  66. ```
  68. ## Merging data - All samples
  70. ```{r}
  71. seurat.merge <- merge(
  72. x = d9433_d7_mock,
  73. y = list(d9434_d7_pr, d9435_d7_mal, d9436_d7_ug),
  74. add.cell.ids = c("Mock", "PR", "MAL", "UG"),
  75. project = "020222Gehrke"
  76. )
  77. ```
  78. # Standard pre-processing worflow
  80. ## Process data
  82. ```{r}
  83. seurat.merge <- NormalizeData(seurat.merge, normalization.method = "LogNormalize", scale.factor = 10000)
  84. seurat.merge <- FindVariableFeatures(seurat.merge, selection.method = "vst", nfeatures = 2000)
  85. seurat.merge <- ScaleData(seurat.merge)
  86. seurat.merge <- RunPCA(seurat.merge, features = VariableFeatures(object = seurat.merge))
  87. ElbowPlot(seurat.merge, ndims = 50)
  89. seurat.merge <- RunUMAP(seurat.merge, reduction="pca",dims=1:30)
  90. seurat.merge <- FindNeighbors(seurat.merge, dims = 1:30, verbose = FALSE)
  91. seurat.merge <- FindClusters(seurat.merge, verbose = FALSE)
  92. ```
  94. ## Dim plot to explore the unintegrated data
  96. ```{r}
  97. Idents(seurat.merge) <- 'RNA_snn_res.0.8'
  98. seurat.merge$Condition <- factor(seurat.merge$Condition, levels = c("Mock", "Puerto Rico", "Malaysia", "Uganda"))
  99. DimPlot(seurat.merge,reduction="umap",split.by="Condition",label=TRUE,repel=FALSE) + NoLegend()
  100. ```
  102. ## QCs
  104. ```{r}
  105. seurat.merge[["percent.mt"]] <- PercentageFeatureSet(seurat.merge, pattern = "^MT-")
  106. seurat.merge$log10GenesPerUMI <- log10(seurat.merge$nFeature_RNA) / log10(seurat.merge$nCount_RNA)
  108. l2.n_feature <- log2([email hidden]$nFeature_RNA)
  109. l2.n_count <- log2([email hidden]$nCount_RNA)
  110. seurat.merge <- AddMetaData(seurat.merge, l2.n_feature, "l2.n_feature")
  111. seurat.merge <- AddMetaData(seurat.merge, l2.n_count, "l2.n_count")
  113. p <- table([email hidden]$RNA_snn_res.0.8,[email hidden]$Condition)
  114. p
  116. Idents(seurat.merge) <- 'RNA_snn_res.0.8'
  117. VlnPlot(seurat.merge, features=c("percent.mt"), split.by = "Condition", ncol=1, pt.size=0)+ NoLegend()
  118. VlnPlot(seurat.merge, features=c("l2.n_feature"), split.by = "Condition", ncol=1, pt.size=0)+ NoLegend()
  119. VlnPlot(seurat.merge, features=c("l2.n_count"), split.by = "Condition", ncol=1, pt.size=0)+ NoLegend()
  121. [email hidden] %>%
  122. ggplot(aes(color=orig.ident, x=nCount_RNA, fill= orig.ident)) +
  123. geom_density(alpha = 0.2) +
  124. scale_x_log10() +
  125. theme_classic() +
  126. ylab("Cell density") +
  127. geom_vline(xintercept = 500) +
  128. ggtitle("nCount")
  130. [email hidden] %>%
  131. ggplot(aes(color=orig.ident, x=nFeature_RNA, fill= orig.ident)) +
  132. geom_density(alpha = 0.2) +
  133. scale_x_log10() +
  134. theme_classic() +
  135. ylab("Cell density") +
  136. geom_vline(xintercept = 200) +
  137. ggtitle("nFeature")
  139. [email hidden] %>%
  140. ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=percent.mt)) +
  141. geom_point() +
  142. scale_color_gradient(low = "gray90", high = "black") +
  143. stat_smooth(method=lm) +
  144. scale_x_log10() +
  145. scale_y_log10() +
  146. theme_classic() +
  147. geom_vline(xintercept = 500) +
  148. geom_hline(yintercept = 250) +
  149. facet_wrap(~orig.ident)
  151. [email hidden] %>%
  152. ggplot(aes(x=log10GenesPerUMI, color = orig.ident, fill=orig.ident)) +
  153. geom_density(alpha = 0.2) +
  154. theme_classic() +
  155. geom_vline(xintercept = 0.8)
  156. ```
  158. ## Gene per UMI filtering of the merged object
  160. ```{r}
  161. seurat.merge.filt <- subset(seurat.merge, percent.mt <= 25 & log10GenesPerUMI > 0.75)
  163. VlnPlot(seurat.merge.filt, group.by = "orig.ident",features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  165. [email hidden] %>%
  166. ggplot(aes(x=log10GenesPerUMI, color = orig.ident, fill=orig.ident)) +
  167. geom_density(alpha = 0.2) +
  168. theme_classic() +
  169. geom_vline(xintercept = 0.8)
  171. Idents(seurat.merge.filt) <- 'RNA_snn_res.0.8'
  172. DimPlot(seurat.merge.filt,reduction="umap",split.by="Condition",label=TRUE,repel=FALSE) + NoLegend()
  173. ```
  175. ## Gene level filtering of the merged object
  177. ```{r}
  178. seurat.merge.filt2 <- seurat.merge.filt
  179. # Output a logical vector for every gene on whether the more than zero counts per cell
  180. # Extract counts
  181. counts <- GetAssayData(object = seurat.merge.filt2, slot = "counts")
  183. # Output a logical vector for every gene on whether the more than zero counts per cell
  184. nonzero <- counts > 0
  186. # Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
  187. keep_genes <- Matrix::rowSums(nonzero) >= 10
  189. # Only keeping those genes expressed in more than 10 cells
  190. filtered_counts <- counts[keep_genes, ]
  192. # Reassign to filtered Seurat object
  193. seurat.merge.filt2 <- CreateSeuratObject(filtered_counts, meta.data = [email hidden])
  194. ```
  196. ## Reprocess of the filtered merged object
  198. ```{r}
  199. seurat.merge.filt2 <- NormalizeData(seurat.merge.filt2, normalization.method = "LogNormalize", scale.factor = 10000)
  200. seurat.merge.filt2 <- FindVariableFeatures(seurat.merge.filt2, selection.method = "vst", nfeatures = 2000)
  201. seurat.merge.filt2 <- ScaleData(seurat.merge.filt2)
  202. seurat.merge.filt2 <- RunPCA(seurat.merge.filt2, features = VariableFeatures(object = seurat.merge.filt2))
  203. ElbowPlot(seurat.merge.filt2, ndims = 50)
  205. seurat.merge.filt2 <- RunUMAP(seurat.merge.filt2, reduction="pca",dims=1:30)
  206. seurat.merge.filt2 <- FindNeighbors(seurat.merge.filt2, dims = 1:30, verbose = FALSE)
  207. seurat.merge.filt2 <- FindClusters(seurat.merge.filt2, verbose = FALSE)
  209. Idents(seurat.merge.filt2) <- 'RNA_snn_res.0.8'
  210. DimPlot(seurat.merge.filt2,reduction="umap",split.by="Condition",label=TRUE,repel=FALSE) + NoLegend()
  211. ```
  213. # SCT normalization, MT regression and integration
  215. ```{r}
  216. split_seurat <- SplitObject(seurat.merge.filt2, split.by = "orig.ident")
  218. for (i in 1:length(split_seurat)) {
  219. split_seurat[[i]] <- NormalizeData(split_seurat[[i]], verbose = TRUE)
  220. split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("percent.mt"))
  221. }
  223. # Select the most variable features to use for integration
  224. integ_features <- SelectIntegrationFeatures(object.list = split_seurat,
  225. nfeatures = 3000)
  227. # Prepare the SCT list object for integration
  228. split_seurat <- PrepSCTIntegration(object.list = split_seurat,
  229. anchor.features = integ_features)
  231. # Find best buddies - can take a while to run
  232. integ_anchors <- FindIntegrationAnchors(object.list = split_seurat,
  233. normalization.method = "SCT",
  234. anchor.features = integ_features)
  236. # Integrate across conditions
  237. seurat.integrated <- IntegrateData(anchorset = integ_anchors,
  238. normalization.method = "SCT")
  239. ```
  241. ## PCA view of integrated data
  243. ```{r}
  244. # Run PCA
  245. seurat.integrated <- RunPCA(object = seurat.integrated)
  247. # Plot PCA
  248. PCAPlot(seurat.integrated,split.by = "Condition")
  250. # Elbow Plot
  251. ElbowPlot(seurat.integrated, ndims = 50)
  252. ```
  254. # UMAP of integrated data
  256. ```{r}
  257. seurat.integrated <- RunUMAP(seurat.integrated, dims = 1:30, reduction = "pca")
  258. DimPlot(seurat.integrated, split.by = "Condition")
  259. ```
  261. ## Integrated Clustering
  263. ```{r}
  264. # Determine the K-nearest neighbor graph
  265. seurat.integrated <- FindNeighbors(object = seurat.integrated, dims = 1:30)
  267. # Determine the clusters for various resolutions
  268. seurat.integrated <- FindClusters(object = seurat.integrated, resolution = c(0.4, 0.8, 1.2))
  269. ```
  271. ## Dimplots of different resolutions
  273. ```{r fig.width=11, fig.height=4}
  275. Idents(seurat.integrated) <- 'integrated_snn_res.0.4'
  276. DimPlot(seurat.integrated, split.by = "Condition",label=TRUE,repel=FALSE) + NoLegend() + ggtitle("0.4 res")
  278. Idents(seurat.integrated) <- 'integrated_snn_res.0.8'
  279. DimPlot(seurat.integrated, split.by = "Condition",label=TRUE,repel=FALSE) + NoLegend() + ggtitle("0.8 res")
  281. Idents(seurat.integrated) <- 'integrated_snn_res.1.2'
  282. DimPlot(seurat.integrated, split.by = "Condition",label=TRUE,repel=FALSE) + NoLegend() + ggtitle("1.4 res")
  283. ```
  285. ## Proportion tables of integrated data
  287. ```{r}
  288. p <- table([email hidden]$integrated_snn_res.0.4,[email hidden]$orig.ident)
  289. p
  290. (round(prop.table(p,2),3)*100)
  291. ```
  293. # Investigate quality control features of the integrated clusters
  295. ```{r fig.width=5, fig.height=4}
  296. Idents(seurat.integrated,) <- 'integrated_snn_res.0.4'
  297. vln_plot_mt <- VlnPlot(seurat.integrated, features=c("percent.mt"), ncol=1, pt.size=0)+ NoLegend()
  298. vln_plot_feature <- VlnPlot(seurat.integrated, features=c("l2.n_feature"), ncol=1, pt.size=0)+ NoLegend()
  299. vln_plot_count <- VlnPlot(seurat.integrated, features=c("l2.n_count"), ncol=1, pt.size=0)+ NoLegend()
  300. UMAP_plot <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  302. ggarrange(UMAP_plot, vln_plot_mt, vln_plot_feature, vln_plot_count, ncol = 2, nrow = 2)
  303. ```
  305. # Cell type assignment of clusters
  307. ## Identiify cell type biomarkers
  309. ```{r}
  310. DefaultAssay(seurat.integrated)<-"RNA"
  311. Idents(object = seurat.integrated) <- 'integrated_snn_res.0.4'
  313. seurat.integrated.clusterMarkers <- FindAllMarkers(seurat.integrated, only.pos=TRUE, min.pct=0.25)
  315. write.xlsx(seurat.integrated.clusterMarkers, file="seurat.integrated.clusterMarkers.xlsx", rowNames=TRUE, overwrite = TRUE)
  316. ```
  318. ## Biomarkers gene expression pattern
  320. ### Biomarker TTR
  322. ```{r fig.width=11, fig.height=4}
  323. vln_plot_TTR <- VlnPlot(seurat.integrated, features = "TTR", pt.size=0) + NoLegend()
  324. UMAP_plot <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  325. feature_plot_TTR <- FeaturePlot(seurat.integrated, features = c("TTR"), min.cutoff = "q10", max.cutoff = "q90")
  327. ggarrange(vln_plot_TTR, UMAP_plot, feature_plot_TTR, ncol = 3, nrow = 1)
  328. ```
  330. ### Biomarker PCP4
  332. ```{r fig.width=11, fig.height=4}
  333. vln_plot_PCP4 <- VlnPlot(seurat.integrated, features = "PCP4", pt.size=0) + NoLegend()
  334. UMAP_plot_PCP4 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  335. feature_plot_PCP4 <- FeaturePlot(seurat.integrated, features = c("PCP4"), min.cutoff = "q10", max.cutoff = "q90")
  337. ggarrange(vln_plot_PCP4, UMAP_plot_PCP4, feature_plot_PCP4, ncol = 3, nrow = 1)
  338. ```
  340. ### Biomarker STMN2
  342. ```{r fig.width=11, fig.height=4}
  343. vln_plot_STMN2 <- VlnPlot(seurat.integrated, features = "STMN2", pt.size=0) + NoLegend()
  344. UMAP_plot_STMN2 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  345. feature_plot_STMN2 <- FeaturePlot(seurat.integrated, features = c("STMN2"), min.cutoff = "q10", max.cutoff = "q90")
  347. ggarrange(vln_plot_EOMES, UMAP_plot_EOMES, feature_plot_EOMES, ncol = 3, nrow = 1)
  348. ```
  350. ### Biomarker MAP2
  352. ```{r fig.width=11, fig.height=4}
  353. vln_plot_MAP2 <- VlnPlot(seurat.integrated, features = "MAP2", pt.size=0) + NoLegend()
  354. UMAP_plot <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  355. feature_plot_MAP2 <- FeaturePlot(seurat.integrated, features = c("MAP2"), min.cutoff = "q10", max.cutoff = "q90")
  357. ggarrange(vln_plot_MAP2, UMAP_plot, feature_plot_MAP2, ncol = 3, nrow = 1)
  358. ```
  360. ### Biomarker TOP2A
  362. ```{r fig.width=11, fig.height=4}
  363. vln_plot_TOP2A <- VlnPlot(seurat.integrated, features = "TOP2A", pt.size=0) + NoLegend()
  364. UMAP_plot <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  365. feature_plot_TOP2A <- FeaturePlot(seurat.integrated, features = c("TOP2A"), min.cutoff = "q10", max.cutoff = "q90")
  367. ggarrange(vln_plot_TOP2A, UMAP_plot, feature_plot_TOP2A, ncol = 3, nrow = 1)
  368. ```
  370. ### Biomarker CCNB2
  372. ```{r fig.width=11, fig.height=4}
  373. vln_plot_CCNB2 <- VlnPlot(seurat.integrated, features = "CCNB2", pt.size=0) + NoLegend()
  374. UMAP_plot_CCNB2 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  375. feature_plot_CCNB2 <- FeaturePlot(seurat.integrated, features = c("CCNB2"), min.cutoff = "q10", max.cutoff = "q90")
  377. ggarrange(vln_plot_TBR2, UMAP_plot_TBR2, feature_plot_TBR2, ncol = 3, nrow = 1)
  378. ```
  380. ### Biomarker NEFL
  382. ```{r fig.width=11, fig.height=4}
  383. vln_plot_NEFL <- VlnPlot(seurat.integrated, features = "NEFL", pt.size=0) + NoLegend()
  384. UMAP_plot_NEFL <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  385. feature_plot_NEFL <- FeaturePlot(seurat.integrated, features = c("NEFL"), min.cutoff = "q10", max.cutoff = "q90")
  387. ggarrange(vln_plot_NEFL, UMAP_plot_NEFL, feature_plot_NEFL, ncol = 3, nrow = 1)
  388. ```
  390. ### Biomarker MEIS2
  392. ```{r fig.width=11, fig.height=4}
  393. vln_plot_MEIS2 <- VlnPlot(seurat.integrated, features = "MEIS2", pt.size=0) + NoLegend()
  394. UMAP_plot_MEIS2 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  395. feature_plot_MEIS2 <- FeaturePlot(seurat.integrated, features = c("MEIS2"), min.cutoff = "q10", max.cutoff = "q90")
  397. ggarrange(vln_plot_MEIS2, UMAP_plot_MEIS2, feature_plot_MEIS2, ncol = 3, nrow = 1)
  398. ```
  400. ### Biomarker SOX2
  402. ```{r fig.width=11, fig.height=4}
  403. vln_plot_SOX2 <- VlnPlot(seurat.integrated, features = "SOX2", pt.size=0) + NoLegend()
  404. UMAP_plot_SOX2 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  405. feature_plot_SOX2 <- FeaturePlot(seurat.integrated, features = c("SOX2"), min.cutoff = "q10", max.cutoff = "q90")
  407. ggarrange(vln_plot_SOX2, UMAP_plot_SOX2, feature_plot_SOX2, ncol = 3, nrow = 1)
  408. ```
  410. ### Biomarker PTN
  412. ```{r fig.width=11, fig.height=4}
  413. vln_plot_PTN <- VlnPlot(seurat.integrated, features = "PTN", pt.size=0) + NoLegend()
  414. UMAP_plot_PTN <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  415. feature_plot_PTN <- FeaturePlot(seurat.integrated, features = c("PTN"), min.cutoff = "q10", max.cutoff = "q90")
  417. ggarrange(vln_plot_PTN, UMAP_plot_PTN, feature_plot_PTN, ncol = 3, nrow = 1)
  418. ```
  420. ### Biomarker S100A10
  422. ```{r fig.width=11, fig.height=4}
  423. vln_plot_S100A10 <- VlnPlot(seurat.integrated, features = "S100A10", pt.size=0) + NoLegend()
  424. UMAP_plot_S100A10 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  425. feature_plot_S100A10 <- FeaturePlot(seurat.integrated, features = c("S100A10"), min.cutoff = "q10", max.cutoff = "q90")
  427. ggarrange(vln_plot_S100A10, UMAP_plot_S100A10, feature_plot_S100A10, ncol = 3, nrow = 1)
  428. ```
  430. ### Biomarker COL1A1
  432. ```{r fig.width=11, fig.height=4}
  433. vln_plot_COL1A1 <- VlnPlot(seurat.integrated, features = "COL1A1", pt.size=0) + NoLegend()
  434. UMAP_plot_COL1A1 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  435. feature_plot_COL1A1 <- FeaturePlot(seurat.integrated, features = c("COL1A1"), min.cutoff = "q10", max.cutoff = "q90")
  437. ggarrange(vln_plot_COL1A1, UMAP_plot_COL1A1, feature_plot_COL1A1, ncol = 3, nrow = 1)
  438. ```
  440. ### Biomarker COL1A2
  442. ```{r fig.width=11, fig.height=4}
  443. vln_plot_COL1A2 <- VlnPlot(seurat.integrated, features = "COL1A2", pt.size=0) + NoLegend()
  444. UMAP_plot_COL1A2 <- DimPlot(seurat.integrated, label=TRUE,repel=FALSE) + NoLegend()
  445. feature_plot_COL1A2 <- FeaturePlot(seurat.integrated, features = c("COL1A2"), min.cutoff = "q10", max.cutoff = "q90")
  447. ggarrange(vln_plot_COL1A2, UMAP_plot_COL1A2, feature_plot_COL1A2, ncol = 3, nrow = 1)
  448. ```
  450. ## Rename cluster names
  452. ```{r fig.width=4.5, fig.height=4.5}
  453. Idents(object = seurat.integrated) <- 'integrated_snn_res.0.4'
  454. new.cluster.ids <- c("Neuronal progenitor cells", "Neuronal progenitor cells", "Intermediate progenitor cells", "Neurons", "Cell debris", "Neuronal progenitor cells", "Midbrain-Hindbrain", "Radial glial cells", "Ventral progenitors", "Ventral progenitors", "Choroid plexus", "Astrocytes", "Neuronal progenitor cells")
  455. names(new.cluster.ids) <- levels(seurat.integrated)
  456. seurat.integrated <- RenameIdents(seurat.integrated, new.cluster.ids)
  457. DimPlot(seurat.integrated, reduction = "umap", label = TRUE, pt.size = 0.5, label.size = 5, repel = TRUE) + NoLegend()
  458. ```
  460. ## Bubble plot with biomarkers
  462. ```{r fig.width=7, fig.height=6}
  463. markers.to.plot <- c("SOX2", "TOP2A", "MAP2", "NEFL", "TLE4", "EPB41L4A-AS1", "TTR", "S100A10")
  464. DotPlot(seurat.integrated, features = markers.to.plot, cols = c("black", "brown", "yellow", "green"), dot.scale = 8, split.by = "Condition") +
  465. RotatedAxis()
  466. ```
  468. # Viral expression
  470. ## Between samples and conditions
  472. ```{r}
  473. seurat.integrated$Condition <- factor(seurat.integrated$Condition, levels = c("Mock", "Malaysia", "Puerto Rico", "Uganda"))
  475. # Create individual violin plots
  476. vln_plot_1 <- VlnPlot(seurat.integrated, features = "vir-ZVMAL", group.by = "Condition", pt.size=0.3, y.max=8) + NoLegend()
  477. vln_plot_2 <- VlnPlot(seurat.integrated, features = "vir-ZVPR", group.by = "Condition", pt.size=0.3, y.max=8) + NoLegend()
  478. vln_plot_3 <- VlnPlot(seurat.integrated, features = "vir-ZVUG", group.by = "Condition", pt.size=0.3, y.max=8) + NoLegend()
  480. # Arrange the plots together
  481. ggarrange(vln_plot_1, vln_plot_2, vln_plot_3, ncol = 3, nrow = 1)
  482. ```
  484. ## Per celltype
  486. ```{r fig.width=6, fig.height=12}
  487. # Subset the Seurat object to exclude the "Cell debris" cluster
  488. seurat_no_debris <- subset(seurat.integrated, idents = "Cell debris", invert = TRUE)
  490. # Create individual violin plots
  491. vln_plot_1 <- VlnPlot(seurat_no_debris, features = "vir-ZVMAL", split.by = "Condition", pt.size=0, y.max=8)
  492. vln_plot_2 <- VlnPlot(seurat_no_debris, features = "vir-ZVPR", split.by = "Condition", pt.size=0, y.max=8)
  493. vln_plot_3 <- VlnPlot(seurat_no_debris, features = "vir-ZVUG", split.by = "Condition", pt.size=0, y.max=8)
  495. # Arrange the plots together
  496. ggarrange(vln_plot_2, vln_plot_1, vln_plot_3, ncol = 1, nrow = 3)
  497. ```
  499. ## QC per celltype per infection
  501. ```{r fig.width=10, fig.height=15}
  503. # Create individual violin plots
  504. vln_plot_mt <- VlnPlot(seurat_no_debris, features=c("percent.mt"), split.by = "Condition", ncol=1, pt.size=0)
  505. vln_plot_feature <- VlnPlot(seurat_no_debris, features=c("l2.n_feature"), split.by = "Condition", ncol=1, pt.size=0)
  506. vln_plot_count <- VlnPlot(seurat_no_debris, features=c("l2.n_count"), split.by = "Condition", ncol=1, pt.size=0)
  508. # Arrange the plots together
  509. combined_QCplot <- ggarrange(
  510. vln_plot_mt + theme(plot.margin = margin(5.5, 5.5, 5.5, 60, "pt")),
  511. vln_plot_feature + theme(plot.margin = margin(5.5, 5.5, 5.5, 60, "pt")),
  512. vln_plot_count + theme(plot.margin = margin(5.5, 5.5, 5.5, 60, "pt")),
  513. ncol = 1, nrow = 3
  514. )
  516. ggsave("combined_QCplot.pdf", plot = combined_QCplot, width = 10, height = 15)
  517. ```
  519. ## Number of cells per celltype per condition
  521. ```{r}
  522. # Create the table
  523. cell_counts <- table(Idents(seurat_no_debris), seurat_no_debris$Condition)
  525. # Convert the table to a data frame for exporting
  526. cell_counts_df <- as.data.frame.matrix(cell_counts)
  528. # Export to Excel
  529. write.xlsx(cell_counts_df, file = "cell_counts_clusters_conditions.xlsx", rowNames = TRUE)
  530. ```
  532. ## Number of infected and non infected cells
  534. ```{r}
  535. DefaultAssay(seurat_no_debris)<-"RNA"
  537. # Specify the genes of interest
  538. genes_of_interest <- c("vir-ZVPR", "vir-ZVMAL", "vir-ZVUG")
  540. # Extract expression data for the genes of interest
  541. expression_data <- FetchData(seurat_no_debris, vars = genes_of_interest)
  543. # Identify cells expressing each gene
  544. expressing_cells <- expression_data > 0
  546. # Identify cells that do not express any of the genes
  547. non_expressing_cells <- rowSums(expressing_cells) == 0
  549. # Get cell types from the active ident
  550. cell_types <- [email hidden]
  552. # Get conditions from the metadata
  553. conditions <- seurat_no_debris$Condition
  555. # Create a dataframe with cell types, conditions, and expression status
  556. expression_per_celltype_condition <- data.frame(cell_type = cell_types,
  557. condition = conditions,
  558. expressing_cells,
  559. none_expressed = non_expressing_cells)
  561. # Count the number of cells expressing each gene and not expressing any gene per cell type and condition
  562. counts_per_celltype_condition <- aggregate(. ~ cell_type + condition,
  563. data = expression_per_celltype_condition,
  564. FUN = sum)
  566. # View the results
  567. print(counts_per_celltype_condition)
  569. # Specify the output file name
  570. output_file <- "CellType_Condition_ViralExpression.xlsx"
  572. # Export the result to an Excel file
  573. write.xlsx(counts_per_celltype_condition, file = output_file, rowNames = FALSE)
  574. ```
  577. # Bubble plot interferon genes
  579. ```{r fig.width=7, fig.height=6}
  580. interferons.to.plot.bulk <- c("IFIT2", "ISG15", "IFIT3", "STAT1", "IFIT5", "IFIT1", "IFI27", "HLA-B", "IFI6", "EIF2AK2")
  582. # Create the DotPlot
  583. DotPlot(seurat.integrated, features = interferons.to.plot.bulk,
  584. cols = c("black", "brown", "yellow", "green"),
  585. dot.scale = 8, split.by = "Condition", scale.by = "size") +
  586. RotatedAxis()
  588. interferonDotPlot <- DotPlot(seurat.integrated, features = interferons.to.plot.bulk,
  589. cols = c("black", "brown", "yellow", "green"),
  590. dot.scale = 8, split.by = "Condition", scale.by = "size") +
  591. RotatedAxis()
  593. # Save the plot as a PDF
  594. ggsave(filename = "interferonDotPlot.pdf", plot = interferonDotPlot, device = "pdf", width = 10, height = 12)
  595. ```
  597. # Differentially Expressed Genes per cluster virus vs mock
  599. ```{r}
  600. seurat.integrated$Condition <- factor(seurat.integrated$Condition, levels = c("Mock", "Malaysia", "Puerto Rico", "Uganda"))
  601. seurat.integrated$ClusterNames <- [email hidden]
  602. seurat.integrated$cluster.condition <- paste(seurat.integrated$ClusterNames, seurat.integrated$Condition, sep = "_")
  604. Idents(seurat.integrated) <- "cluster.condition"
  606. table(seurat.integrated$cluster.condition)
  608. # Create a data frame from the table output
  609. df <- as.data.frame(table(seurat.integrated$cluster.condition))
  611. # Split the cluster.condition into separate columns for cluster and condition
  612. df <- tidyr::separate(df, Var1, into = c("Cluster", "Condition"), sep = "_")
  614. # Create the bar plot
  615. ggplot(df, aes(x = Cluster, y = Freq, fill = Condition)) +
  616. geom_bar(stat = "identity", position = "dodge") +
  617. labs(title = "Number of cells per cluster per condition",
  618. x = "Cluster",
  619. y = "Count",
  620. fill = "Condition") +
  621. theme_minimal()+
  622. theme(axis.text.x = element_text(angle = 45, hjust = 1))
  623. ```
  625. ```{r}
  626. # Initialize a list to store the results
  627. de_results_PR <- list()
  628. de_results_Mal <- list()
  629. de_results_Ug <- list()
  630. de_results_MalPR <- list()
  631. de_results_UgPR <- list()
  632. de_results_MalUg <- list()
  634. # Define the clusters and conditions
  635. clusters <- levels(seurat.integrated$ClusterNames)
  636. PRconditions <- c("Puerto Rico", "Mock")
  637. Malconditions <- c("Malaysia", "Mock")
  638. Ugconditions <- c("Uganda", "Mock")
  639. MalPRconditions <- c("Malaysia", "Puerto Rico")
  640. UgPRconditions <- c("Uganda", "Puerto Rico")
  641. MalUgconditions <- c("Malaysia", "Uganda")
  643. # Loop through each cluster and perform differential expression analysis for PR vs mock
  644. for (cluster in clusters) {
  645. ident.1 <- paste(cluster, PRconditions[1], sep = "_")
  646. ident.2 <- paste(cluster, PRconditions[2], sep = "_")
  648. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  649. de_results_PR[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  650. } else {
  651. message(paste("Cluster", cluster, "does not have both conditions."))
  652. }
  653. }
  655. # Loop through each cluster and perform differential expression analysis for Mal vs mock
  656. for (cluster in clusters) {
  657. ident.1 <- paste(cluster, Malconditions[1], sep = "_")
  658. ident.2 <- paste(cluster, Malconditions[2], sep = "_")
  660. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  661. de_results_Mal[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  662. } else {
  663. message(paste("Cluster", cluster, "does not have both conditions."))
  664. }
  665. }
  667. # Loop through each cluster and perform differential expression analysis for Ug vs mock
  668. for (cluster in clusters) {
  669. ident.1 <- paste(cluster, Ugconditions[1], sep = "_")
  670. ident.2 <- paste(cluster, Ugconditions[2], sep = "_")
  672. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  673. de_results_Ug[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  674. } else {
  675. message(paste("Cluster", cluster, "does not have both conditions."))
  676. }
  677. }
  679. # Loop through each cluster and perform differential expression analysis for Mal vs PR
  680. for (cluster in clusters) {
  681. ident.1 <- paste(cluster, MalPRconditions[1], sep = "_")
  682. ident.2 <- paste(cluster, MalPRconditions[2], sep = "_")
  684. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  685. de_results_MalPR[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  686. } else {
  687. message(paste("Cluster", cluster, "does not have both conditions."))
  688. }
  689. }
  691. # Loop through each cluster and perform differential expression analysis for Ug vs PR
  692. for (cluster in clusters) {
  693. ident.1 <- paste(cluster, UgPRconditions[1], sep = "_")
  694. ident.2 <- paste(cluster, UgPRconditions[2], sep = "_")
  696. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  697. de_results_UgPR[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  698. } else {
  699. message(paste("Cluster", cluster, "does not have both conditions."))
  700. }
  701. }
  703. # Loop through each cluster and perform differential expression analysis for Mal vs Ug
  704. for (cluster in clusters) {
  705. ident.1 <- paste(cluster, MalUgconditions[1], sep = "_")
  706. ident.2 <- paste(cluster, MalUgconditions[2], sep = "_")
  708. if (ident.1 %in% Idents(seurat.integrated) & ident.2 %in% Idents(seurat.integrated)) {
  709. de_results_MalUg[[paste("Cluster", cluster)]] <- FindMarkers(seurat.integrated, ident.1 = ident.1, ident.2 = ident.2, verbose = FALSE, logfc.threshold = 0)
  710. } else {
  711. message(paste("Cluster", cluster, "does not have both conditions."))
  712. }
  713. }
  714. ```
  717. # GSEA - Functional Class Sorting
  719. Preparing rank files
  720. Extract protein coding genes from assembled data files, adjust the data and write rnk file in gsea directory
  722. ## Rank files for GSEA
  724. Ranked gene lists (rnk files) were generated from differential expression results using avg_log2FC values. The process involved iterating through each cluster of differential expression results for the conditions Mal, PR, and Ug, extracting the avg_log2FC values, and renaming the columns to reflect the specific condition and cluster (e.g., Ug_v_mock.Neuronal_progenitor_cells.rnk, PR_v_mock.Neurons.rnk, etc.). These avg_log2FC values were stored in individual data frames for each celltype, which were subsequently merged into a single data frame for each condition.
  726. For each merged data frame, the gene name and corresponding avg_log2FC values were selected and saved as rnk files, ensuring only complete cases were included. These rnk files were then used for Gene Set Enrichment Analysis (GSEA) utilizing the prerank_4.3.2_v2023.1.sh script to identify enriched gene sets within the collection Gene Ontology Biological Process (c5bp).
  728. ```{r}
  729. # Initialize an empty list to store avg_log2FC dataframes
  730. avg_log2FC_list_Mal <- list()
  731. avg_log2FC_list_PR <- list()
  732. avg_log2FC_list_Ug <- list()
  733. avg_log2FC_list_MalPR <- list()
  734. avg_log2FC_list_UgPR <- list()
  735. avg_log2FC_list_MalUg <- list()
  737. # Loop through the differential expression results and extract avg_log2FC
  738. for (cluster in names(de_results_Mal)) {
  739. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  740. column_name <- paste0("mal_v_mock.", cluster_name) # Create the new column name
  742. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_Mal[[cluster]]),
  743. avg_log2FC = de_results_Mal[[cluster]]$avg_log2FC)
  744. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  745. avg_log2FC_list_Mal[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  746. }
  748. for (cluster in names(de_results_PR)) {
  749. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  750. column_name <- paste0("PR_v_mock.", cluster_name) # Create the new column name
  752. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_PR[[cluster]]),
  753. avg_log2FC = de_results_PR[[cluster]]$avg_log2FC)
  754. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  755. avg_log2FC_list_PR[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  756. }
  758. for (cluster in names(de_results_Ug)) {
  759. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  760. column_name <- paste0("Ug_v_mock.", cluster_name) # Create the new column name
  762. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_Ug[[cluster]]),
  763. avg_log2FC = de_results_Ug[[cluster]]$avg_log2FC)
  764. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  765. avg_log2FC_list_Ug[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  766. }
  768. for (cluster in names(de_results_MalPR)) {
  769. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  770. column_name <- paste0("Mal_v_PR.", cluster_name) # Create the new column name
  772. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_MalPR[[cluster]]),
  773. avg_log2FC = de_results_MalPR[[cluster]]$avg_log2FC)
  774. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  775. avg_log2FC_list_MalPR[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  776. }
  778. for (cluster in names(de_results_UgPR)) {
  779. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  780. column_name <- paste0("Ug_v_PR.", cluster_name) # Create the new column name
  782. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_UgPR[[cluster]]),
  783. avg_log2FC = de_results_UgPR[[cluster]]$avg_log2FC)
  784. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  785. avg_log2FC_list_UgPR[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  786. }
  788. for (cluster in names(de_results_MalUg)) {
  789. cluster_name <- sub("Cluster ", "", cluster) # Remove "Cluster " prefix to get the cluster name
  790. column_name <- paste0("Mal_v_Ug.", cluster_name) # Create the new column name
  792. avg_log2FC_df <- data.frame(gene_name = rownames(de_results_MalUg[[cluster]]),
  793. avg_log2FC = de_results_MalUg[[cluster]]$avg_log2FC)
  794. colnames(avg_log2FC_df)[2] <- column_name # Rename the avg_log2FC column to the new name
  795. avg_log2FC_list_MalUg[[cluster_name]] <- avg_log2FC_df # Use the cluster name without "Cluster " prefix
  796. }
  800. # Merge all dataframes in the list into one dataframe
  801. merged_df_Mal <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_Mal)
  802. merged_df_PR <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_PR)
  803. merged_df_Ug <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_Ug)
  804. merged_df_MalPR <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_MalPR)
  805. merged_df_UgPR <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_UgPR)
  806. merged_df_MalUg <- Reduce(function(x, y) merge(x, y, by = "gene_name", all = TRUE), avg_log2FC_list_MalUg)
  808. # View the final merged dataframe
  809. head(merged_df_Mal)
  810. head(merged_df_PR)
  811. head(merged_df_Ug)
  813. names <- colnames(merged_df_Mal)
  815. for (i in names[-1]){
  816. d<-merged_df_Mal %>% dplyr::select('gene_name',all_of(i))
  817. names(d)[1] <- '#Gene'
  818. d<-d[complete.cases(d), ]
  819. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  820. }
  822. names <- colnames(merged_df_PR)
  824. for (i in names[-1]){
  825. d<-merged_df_PR %>% dplyr::select('gene_name',all_of(i))
  826. names(d)[1] <- '#Gene'
  827. d<-d[complete.cases(d), ]
  828. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  829. }
  831. names <- colnames(merged_df_Ug)
  833. for (i in names[-1]){
  834. d<-merged_df_Ug %>% dplyr::select('gene_name',all_of(i))
  835. names(d)[1] <- '#Gene'
  836. d<-d[complete.cases(d), ]
  837. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  838. }
  840. names <- colnames(merged_df_MalPR)
  842. for (i in names[-1]){
  843. d<-merged_df_MalPR %>% dplyr::select('gene_name',all_of(i))
  844. names(d)[1] <- '#Gene'
  845. d<-d[complete.cases(d), ]
  846. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  847. }
  849. names <- colnames(merged_df_UgPR)
  851. for (i in names[-1]){
  852. d<-merged_df_UgPR %>% dplyr::select('gene_name',all_of(i))
  853. names(d)[1] <- '#Gene'
  854. d<-d[complete.cases(d), ]
  855. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  856. }
  858. names <- colnames(merged_df_MalUg)
  860. for (i in names[-1]){
  861. d<-merged_df_MalUg %>% dplyr::select('gene_name',all_of(i))
  862. names(d)[1] <- '#Gene'
  863. d<-d[complete.cases(d), ]
  864. write.table(d, sep='\t',file=paste0(i,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
  865. }
  867. ```
  869. ## Process luria tsv GSEA results to produce a single output file using submit_prerank_4.3.2_v2023.2_caw.sh bash script
  871. ## Identify files that need to be imported - the basedir will change according to run
  873. After running the javaGSEA version 4.3.2 bash script with msigDb version 2023.2 (Subramanian et al. 2005), the resulting data files were processed to create a comprehensive and analyzable dataset. The initial step involved identifying and importing the relevant GSEA result files. These files were located in the specified base directory and its subdirectories, matching a specific naming pattern (e.g., "^gsea_.*\.tsv$").
  875. ```{r}
  876. base_dir <- "gsea_celltypes/aug01/"
  878. # Get a list of directories in the base directory
  879. sub_dirs <- list.dirs(base_dir, recursive = FALSE)
  881. # Initialize an empty list to store file paths
  882. file_list <- list()
  884. # Loop through each sub-directory
  885. for (sub_dir in sub_dirs) {
  886. # Get a list of files in the sub-directory that match the pattern
  887. files <- list.files(path = sub_dir, pattern = "^gsea_.*\\.tsv$", full.names = TRUE)
  889. # Add the files to the file list
  890. file_list <- c(file_list, files)
  891. }
  892. ```
  894. ## Import the GSEA data into a list of lists
  896. ```{r}
  897. gsea_results.list <- lapply(file_list, read.delim)
  898. ```
  900. ## Extract relevant results names from output filenames and rename results dataframes
  902. To facilitate the analysis, the filenames were processed to extract meaningful names, indicating the specific comparison and whether the results were positive or negative. These names were used to label the corresponding data frames in the list of GSEA results.
  904. ```{r}
  905. # Define a function to process each file path
  906. process_filepath <- function(filepath) {
  907. # Remove folder info from the front of the file path
  908. short_filename <- gsub("gsea_celltypes/aug01//", "", filepath)
  909. # Replace ".rnk." with "_"
  910. short_filename <- gsub("\\.rnk\\.", "_", short_filename)
  911. # Add "_pos" if "pos" is in the remainder or "_neg" if "neg" in the remainder
  912. if (grepl("pos", short_filename)) {
  913. short_filename <- paste0(short_filename, "_pos")
  914. } else if (grepl("neg", short_filename)) {
  915. short_filename <- paste0(short_filename, "_neg")
  916. }
  917. short_filename <- gsub("\\.GseaP.*_pos$", "_pos", short_filename)
  918. short_filename <- gsub("\\.GseaP.*_neg$", "_neg", short_filename)
  919. return(short_filename)
  920. }
  922. # Apply the function to each file path in file_list
  923. short_filenames <- lapply(file_list, process_filepath)
  924. # Rename the list of gsea results
  925. names(gsea_results.list) <- short_filenames
  926. ```
  928. ## Extract relevant columns from the dataframes
  930. Next, each data frame in the list was refined to include only the relevant columns: "NAME" (gene set name), "SIZE" (number of genes in the set), "NES" (normalized enrichment score), and "FDR.q.val" (false discovery rate q-value).
  932. ```{r}
  933. # Define the columns to select
  934. columns_to_select <- c("NAME", "SIZE", "NES", "FDR.q.val")
  935. # Loop through the list of dataframes
  936. gsea_results.list <- lapply(gsea_results.list, function(df) {
  937. # Select the desired columns using dplyr::select
  938. df_selected <- dplyr::select(df, all_of(columns_to_select))
  939. # Return the modified dataframe
  940. return(df_selected)
  941. })
  942. ```
  944. ## Add a comparison column containing the name of each results list
  946. Additionally, a new column was added to each data frame to indicate the specific comparison it represents, simplifying downstream analysis.
  948. ```{r}
  949. # Loop through the list of dataframes
  950. for (comparison_name in names(gsea_results.list)) {
  951. new_name <- gsub("_pos$|_neg$", "", comparison_name)
  952. # Access the dataframe by its name
  953. df <- gsea_results.list[[comparison_name]]
  954. # Add a new column containing the comparison name
  955. df <- df %>% mutate(Comparison = new_name)
  956. # Update the dataframe in the list
  957. gsea_results.list[[comparison_name]] <- df
  958. }
  959. ```
  962. ```{r}
  963. gsea_results.list <- lapply(gsea_results.list, function(x) {
  964. x$NES <- as.numeric(x$NES)
  965. return(x)
  966. })
  967. ```
  969. ## Assemble the data into a single output file and export the results
  971. The GSEA results were then combined into a single data frame. This combined data frame included an additional column "Collection," derived from the comparison names, which identified the gene set collections (e.g., c2cp, h, c5bp).
  973. ```{r}
  974. combined_df <- bind_rows(gsea_results.list, .id = "ComparisonFull")
  975. #Use mutate and str_extract to create the new 'Collection' column
  976. combined_df <- combined_df %>%
  977. mutate(Collection = str_extract(Comparison, "[^_]+$")) %>%
  978. mutate(Comparison = str_remove(Comparison, "_c2cp|_c5bp|_h"))
  979. ```
  981. ## Pivoting dataframe
  983. To ensure unique entries and appropriate aggregation, the maximum "SIZE" value for each gene set name was determined. The data was then summarized to avoid duplicates, averaging the NES and FDR.q.val values for each gene set and comparison. This summarized data was pivoted to create a wide-format table, where NES and FDR.q.val values for different comparisons were represented as separate columns. The final data frame was enhanced by including the maximum "SIZE" values and sorting by the minimum FDR.q.val across all comparisons.
  986. ```{r}
  987. # Get the maximum SIZE value for each NAME
  988. df_size_max <- combined_df %>%
  989. group_by(NAME) %>%
  990. summarize(SIZE = max(SIZE, na.rm = TRUE), .groups = 'drop')
  992. # Ensure no duplicate entries for the same NAME and Comparison
  993. df_summary <- df %>%
  994. group_by(NAME, Comparison, Collection) %>%
  995. summarize(NES = mean(NES, na.rm = TRUE), FDR.q.val = mean(FDR.q.val, na.rm = TRUE), .groups = 'drop')
  997. # Pivot the dataframe to get NES and FDR.q.val with correct column names
  998. pivoted_df <- df_summary %>%
  999. pivot_wider(names_from = Comparison,
  1000. values_from = c(NES, FDR.q.val),
  1001. names_sep = ".") %>%
  1002. rename_with(~ sub("^(.+?)\\.NES$", "\\1.NES", .), starts_with("NES")) %>%
  1003. rename_with(~ sub("^(.+?)\\.FDR.q.val$", "\\1.FDR.q.val", .), starts_with("FDR.q.val"))
  1005. # Combine the SIZE column with the pivoted dataframe
  1006. final_df <- pivoted_df %>%
  1007. left_join(df_size_max, by = "NAME") %>%
  1008. select(NAME, SIZE, Collection, everything()) # Ensure Collection is included and in the right order
  1010. # Identify FDR.q.val columns
  1011. fdr_columns <- grep("^FDR.q.val", names(final_df), value = TRUE)
  1013. # Create MinFDR column with the minimum value of all FDR.q.val columns
  1014. final_df <- final_df %>%
  1015. rowwise() %>%
  1016. mutate(MinFDR = min(c_across(all_of(fdr_columns)), na.rm = TRUE)) %>%
  1017. ungroup()
  1019. # Sort final_df by MinFDR with smallest values at the top
  1020. final_df <- final_df %>%
  1021. arrange(MinFDR)
  1023. head(final_df)
  1025. write.xlsx(final_df, "GSEA_virus_v_mock_celltypes_c2cp_h_c5bp.xlsx")
  1026. ```
  1028. # Volcano plot
  1030. ```{r fig.width=7, fig.height=7}
  1031. library(EnhancedVolcano)
  1033. # Extract the DEG table for "Cluster Neuronal progenitor cells"
  1034. deg_table <- de_results_Ug[["Cluster Neuronal progenitor cells"]]
  1036. # Define the genes to label
  1037. genes_to_label <- c("CANX", "PDIA4", "HSP90B1", "SELENOS", "HERPUD1", "HSPA5", "SEC61B", "XBP1", "GRINA", "SEC63", "MIF", "CLU")
  1039. # Filter for significant genes in the list to label
  1040. genes_to_label_significant <- intersect(genes_to_label, rownames(deg_table)[deg_table$p_val_adj < 0.05])
  1042. # Define custom colors
  1043. keyvals <- rep('grey', nrow(deg_table))
  1044. names(keyvals) <- rownames(deg_table)
  1045. keyvals[deg_table$avg_log2FC > 0.2 & deg_table$p_val_adj < 0.05] <- 'red'
  1046. keyvals[deg_table$avg_log2FC < -0.2 & deg_table$p_val_adj < 0.05] <- 'blue'
  1048. # Create a volcano plot with custom colors
  1049. EnhancedVolcano(
  1050. deg_table,
  1051. lab = rownames(deg_table),
  1052. x = 'avg_log2FC',
  1053. y = 'p_val_adj',
  1054. xlim = c(-0.75, 0.75),
  1055. selectLab = genes_to_label_significant,
  1056. xlab = bquote(~Log[2]~ 'fold change'),
  1057. ylab = bquote(~-Log[10]~ 'p-value adjusted'),
  1058. title = "Uganda VS Mock Neuronal progenitor cells",
  1059. subtitle = bquote(italic(ER_STRESS_RESPONSE)),
  1060. pCutoff = 0.05,
  1061. FCcutoff = 0.20,
  1062. pointSize = 2.0,
  1063. labSize = 4.0,
  1064. colCustom = keyvals,
  1065. legendPosition = 'right',
  1066. legendLabSize = 10,
  1067. legendIconSize = 3.0,
  1068. drawConnectors = TRUE,
  1069. widthConnectors = 0.5,
  1070. colConnectors = 'black',
  1071. gridlines.major = FALSE,
  1072. gridlines.minor = FALSE,
  1073. border = 'full',
  1074. borderWidth = 1.0,
  1075. borderColour = 'black',
  1076. max.overlaps = Inf
  1077. )+
  1078. theme(legend.position="none")
  1080. ```
  1082. # Heatmap of specific gene sets from GSEA
  1084. ```{r}
  1085. library(ComplexHeatmap)
  1086. library(circlize)
  1088. # Convert tbl_df to data.frame
  1089. df <- as.data.frame(GSEA_LeeGeneSets_forheatmap)
  1091. # Set the first column as row names
  1092. rownames(df) <- df[[1]]
  1094. # Remove the first column from the dataframe
  1095. df <- df[ , -1]
  1097. # Filter rows where MinFDR is under 0.05
  1098. df_filtered <- df[df$MinFDR < 0.05, ]
  1100. df_filtered <- subset(df, select = -c(MinFDR))
  1102. df_filtered[is.na(df_filtered)] <- 0
  1104. # Split dataframe into groups of 3 columns
  1105. column_groups <- split(names(df), ceiling(seq_along(names(df)) / 3))
  1107. # Convert the groups into a list of matrices
  1108. list_of_matrices <- lapply(column_groups, function(cols) {
  1109. mat <- as.matrix(df[, cols])
  1110. rownames(mat) <- rownames(df)
  1111. mat
  1112. })
  1114. # Define a function to create heatmaps
  1115. create_heatmap <- function(mat) {
  1116. Heatmap(mat, name = "NES", cluster_rows = FALSE, cluster_columns = FALSE, show_row_names = TRUE, row_names_side = c("left"), row_names_gp = gpar(fontsize = 4), column_names_gp = gpar(fontsize = 4))
  1117. }
  1119. # Create heatmaps for each group
  1120. heatmaps <- lapply(list_of_matrices, create_heatmap)
  1122. # Combine heatmaps side by side using the + operator
  1123. combined_heatmap <- heatmaps[[1]]
  1124. for (i in 2:length(heatmaps)) {
  1125. combined_heatmap <- combined_heatmap + heatmaps[[i]]
  1126. }
  1128. # Draw the combined heatmap
  1129. draw(combined_heatmap, merge_legends = TRUE)
  1131. # # Save as PDF
  1132. # pdf("combined_heatmap.pdf", width = 10, height = 6) # Adjust width and height as needed
  1133. # draw(combined_heatmap, merge_legends = TRUE)
  1134. # dev.off()
  1135. ```
  1136. <<<<<<< HEAD
  1138. # write session info
  1140. Capturing information about the R session helps document your work
  1142. ```{r}
  1143. sessionInfo()
  1144. writeLines(capture.output(sessionInfo()), "220222Geh_sessionInfo.txt")
  1145. ```

020222_Gehrke.Rmd at commit c6e52b0, under MIT · at the source

Overview

Authors: Alfred T. Harding1, Yichen Zhang1, Jane-Jane Chen1, Jenna M. Antonucci1, Alexsia Richards2, Valerie Leger1, Charles A. Whittaker3, Yann S. Vanrobaeys3, Divyansh Agarwal4, Tenzin Lungjangwa2, Rudolf Jaenisch1,2, Lee Gehrke1,5
ORCID iDs: Yichen Zhang, Lee Gehrke
  1. Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
  2. Whitehead Institute for Biomedical Research, Cambridge, Massachusetts, USA
  3. Bioinformatics & Computing Core Facility of the Swanson Biotechnology Center, Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
  4. Massachusetts General Hospital, Boston, Massachusetts, USA
  5. Harvard Medical School, Boston, Massachusetts, USA
Institutions: Massachusetts Institute of Technology (United States); Whitehead Institute for Biomedical Research (United States); Massachusetts General Hospital (United States)
Journal: mBio, volume 17, issue 7, pages e00863-26
Dates: received 10 April 2026; accepted 26 May 2026; published online 15 June 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1128/mbio.00863-26 · PMID 42294680 · PMCID PMC13343971 · OpenAlex W7164856986
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), other condition (population)
Methods: Statistics, Connectivity
Keywords: virus, organoid, pathogenesis, Zika, flavivirus, stress response, tolerance
MeSH: Brain*, Organoids*, Zika Virus*, Zika Virus Infection*, Apoptosis, Host-Pathogen Interactions, Humans, Neural Stem Cells, Oxidative Stress, Reactive Oxygen Species (* major topic)
Topic: Mosquito-borne diseases and control (Public Health, Environmental and Occupational Health, Medicine), according to OpenAlex
Funding: National Institute of Allergy and Infectious Diseases (U19AI131135); National Cancer Institute (P30-CA14051)
Citations: not cited yet (Europe PMC); 87 references in the paper

Abstract

Zika virus (ZIKV) has multiple lineages and strains that cause a range of disease severity, underscoring the need to elucidate differential neuropathogenesis mechanisms. Here, we performed systematic, side-by-side comparisons of African, Asian, and American ZIKV lineage infections using cerebral organoids derived from human embryonic stem cells, a relevant human model experimental system. African lineage ZIKV strains, as well as the ancestral Asian Malaysia strain, persistently infected neural progenitor cells, causing apoptosis and severe disruption of ventricular cytoarchitecture. In contrast, contemporary Asian and American lineage viruses were cleared from ventricles, coinciding with low apoptosis and reduced neuropathology. Single-cell RNA sequencing demonstrated upregulated cell-type-specific antiviral signaling during American lineage infections, coinciding with viral clearance from ventricular progenitor cells. Conversely, pathogenic African lineage infections were associated with apoptosis, reduced STAT2 and IFIT1 protein levels, and enhanced activation of stress pathways. African lineage and ancestral Malaysian strain infections induced mitochondrial oxidative stress. Scavenging the reactive oxygen species improved ventricular cytoarchitecture and progenitor survival, but without reducing viral titers. Together, these findings suggest that lineage- and strain-specific host stress responses, rather than viral burden alone, contribute to ZIKV-induced neurodevelopmental damage. An implication of this study is that host-directed therapeutic strategies, used to improve host tolerance to viral infection, may benefit clinical outcomes.

IMPORTANCE: This study provides a systematic comparison of Zika virus (ZIKV) lineage infections in a relevant human model system to correlate mechanisms of host cellular responses with neuropathogenesis. By analyzing African, Asian, and American ZIKV infections side by side in cerebral organoids derived from human embryonic stem cells, we link persistent infection of neural progenitor cells to elevated cellular stress responses and structural disruption of organoid ventricles. Structural disruption was reduced by adding a hydroxyl radical scavenger, but without lowering viral titer. These data are important in strongly suggesting that host responses to viral infection can be of equal or greater importance than viral burden in determining pathogenesis. This study provides mechanistic insight into how closely related viral lineage infections result in divergent outcomes in the developing brain. The data have added importance for suggesting the potential value of host-directed therapeutics to improve ZIKV tolerance without addressing viral titer.

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

Repository

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

KochInstitute-Bioinformatics/Gehrke_d7_infection

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c6e52b0ceece280e5e84e9dc544ee3d855cb5537, 18 February 2026
Languages: R (9), Shell (6)
Size: 238 files, 15 scripts
Software Heritage: not archived
Found in: “DATA AVAILABILITY”
Holds: README, license file, 2 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: DESeq2 (3 files), tidyverse (3 files), ComplexHeatmap (2 files), edgeR (2 files), circlize (1 file), ggplot2 (1 file), ggpubr (1 file), Nextflow (1 file), Seurat (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files
At the source: github.com/KochInstitute-Bioinformatics/Gehrke_d7_infection.git

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;
  • 15 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

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

Data

Datasets cited

Data availability

The GitHub repository for the code used for bulk and single-cell RNAseq analyses can be found at https://github.com/KochInstitute-Bioinformatics/Gehrke_d7_infection.git. Bulk and single-cell RNA-Seq data are available from the Gene Expression Omnibus under accession number GSE29779 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE297797).

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 7 keywords, 10 MeSH terms, 2 funders, 83 references.

Cite

This paper

Harding, A. T., Zhang, Y., Chen, J.-J., Antonucci, J. M., Richards, A., Leger, V., Whittaker, C. A., Vanrobaeys, Y. S., Agarwal, D., Lungjangwa, T., Jaenisch, R., & Gehrke, L. (2026). Zika virus infections of human stem cell-derived cerebral organoids reveal viral lineage-specific pathogenesis responses. mBio, 17(7), e00863-26. https://doi.org/10.1128/mbio.00863-26

BibTeX

@article{harding2026zika,
author = {Harding, Alfred T. and Zhang, Yichen and Chen, Jane-Jane and Antonucci, Jenna M. and Richards, Alexsia and Leger, Valerie and Whittaker, Charles A. and Vanrobaeys, Yann S. and Agarwal, Divyansh and Lungjangwa, Tenzin and Jaenisch, Rudolf and Gehrke, Lee},
title = {{Zika virus infections of human stem cell-derived cerebral organoids reveal viral lineage-specific pathogenesis responses}},
journal = {mBio},
year = {2026},
month = jun,
volume = {17},
number = {7},
pages = {e00863--26},
publisher = {American Society for Microbiology (ASM)},
issn = {2150-7511},
doi = {10.1128/mbio.00863-26},
url = {https://doi.org/10.1128/mbio.00863-26},
pmid = {42294680},
pmcid = {PMC13343971}
}

RIS

TY - JOUR
AU - Harding, Alfred T.
AU - Zhang, Yichen
AU - Chen, Jane-Jane
AU - Antonucci, Jenna M.
AU - Richards, Alexsia
AU - Leger, Valerie
AU - Whittaker, Charles A.
AU - Vanrobaeys, Yann S.
AU - Agarwal, Divyansh
AU - Lungjangwa, Tenzin
AU - Jaenisch, Rudolf
AU - Gehrke, Lee
TI - Zika virus infections of human stem cell-derived cerebral organoids reveal viral lineage-specific pathogenesis responses
T2 - mBio
J2 - mBio
PY - 2026
DA - 2026/06/15
VL - 17
IS - 7
SP - e00863
EP - 26
SN - 2150-7511
PB - American Society for Microbiology (ASM)
DO - 10.1128/mbio.00863-26
UR - https://doi.org/10.1128/mbio.00863-26
LA - en
ER -

CSL-JSON

{
"id": "10.1128/mbio.00863-26",
"type": "article-journal",
"title": "Zika virus infections of human stem cell-derived cerebral organoids reveal viral lineage-specific pathogenesis responses",
"container-title": "mBio",
"author": [
{
"family": "Harding",
"given": "Alfred T."
},
{
"family": "Zhang",
"given": "Yichen"
},
{
"family": "Chen",
"given": "Jane-Jane"
},
{
"family": "Antonucci",
"given": "Jenna M."
},
{
"family": "Richards",
"given": "Alexsia"
},
{
"family": "Leger",
"given": "Valerie"
},
{
"family": "Whittaker",
"given": "Charles A."
},
{
"family": "Vanrobaeys",
"given": "Yann S."
},
{
"family": "Agarwal",
"given": "Divyansh"
},
{
"family": "Lungjangwa",
"given": "Tenzin"
},
{
"family": "Jaenisch",
"given": "Rudolf"
},
{
"family": "Gehrke",
"given": "Lee"
}
],
"container-title-short": "mBio",
"volume": "17",
"issue": "7",
"page": "e00863-26",
"DOI": "10.1128/mbio.00863-26",
"PMID": "42294680",
"PMCID": "PMC13343971",
"ISSN": "2150-7511",
"publisher": "American Society for Microbiology (ASM)",
"URL": "https://doi.org/10.1128/mbio.00863-26",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
15
]
]
}
}

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: edgeR, DESeq2, circlize, 5 other tools, 4 references
[2] doi:10.1038/s41593-026-02300-5 [code]
Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: edgeR, DESeq2, circlize, 5 other tools, other condition
[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: edgeR, DESeq2, circlize, 5 other tools, other condition
[4] doi:10.1016/j.isci.2026.115657 [code]
Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: edgeR, DESeq2, circlize, 5 other tools, other condition
[5] 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: Nextflow, edgeR, circlize, 4 other tools, other condition
[6] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: edgeR, DESeq2, circlize, 5 other tools
[7] doi:10.1093/braincomms/fcag239 [code]
Extracellular matrix remodelling in degenerative cervical myelopathy.
Journal: Brain communications
In common: edgeR, DESeq2, circlize, 5 other tools
[8] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: edgeR, DESeq2, circlize, 5 other tools
[9] doi:10.1016/j.isci.2026.115573 [code]
Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: edgeR, DESeq2, circlize, 5 other tools
[10] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: edgeR, DESeq2, circlize, 5 other tools

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.