OSCR

NvashA function reveals temporal differences in neural subtype generation in cnidarians.

Code ↔ Paper

5 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 5 matches
  1. [1] § Materials and methods › Single-cell RNA sequencing and analysis ↔ 05_NvashA_subsetting_trajectory.Rmd, lines 51–67 · score 0.98 · Nvcnido fos1, Pharyngeal Ectoderm, NvfoxA, NvsnailA, Neural Progenitor Cells, likeE
  2. [2] § Results › Single-cell atlas of NvashA-expressing cells identifies embryonic and larval-born neuronal subtypes ↔ 05_NvashA_subsetting_trajectory.Rmd, lines 538–584 · score 0.90 · NvTBA1, n2 g1, n1 l2, n1 l3, N1.g.early, NGD.1
  3. [3] § Materials and methods › Trajectory inference and pseudotime ↔ 05_NvashA_subsetting_trajectory.Rmd, lines 1033–1153 · score 0.90 · SeuratWrappers, cluster_cells, learn_graph, clicking, interactively, nodes
  4. [4] § Results › Larval neurons spatially overlap with previously identified neuronal domains along the oral-aboral axis ↔ 05_NvashA_subsetting_trajectory.Rmd, lines 538–584 · score 0.78 · n1 l1, n2 g1, NGD.1, Nv2.1205, Nv2.2620, axis
  5. [5] § Results › Single-cell atlas of NvashA-expressing cells identifies embryonic and larval-born neuronal subtypes ↔ 05_NvashA_subsetting_trajectory.Rmd, lines 317–331 · score 0.56 · neuron.gast, neuron.pl, trajectory, gland, cnidocyte, NvashA

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,242 lines · 63 KB · no license · 5 matches

  1. ---
  2. title: "ashA figures for pub"
  3. output: html_document
  4. date: "2024-10-07"
  5. ---
  7. ```{r}
  8. library(dplyr)
  9. library(Seurat)
  10. library(patchwork)
  11. library(ggplot2)
  12. library(janitor)
  13. library(sctransform)
  14. library(DoubletFinder)
  15. library(glmGamPoi)
  16. library(BiocManager)
  17. library(multtest)
  18. library(metap)
  19. library(scCustomize)
  20. library(qs)
  21. library(igraph)
  22. library(tidyr)
  23. library(magrittr)
  24. ```
  27. # Supplemental Figure 1
  28. #gastrula stage
  29. ```{r}
  30. Stella_combined_filt_gastrula <- FindNeighbors(Stella_combined_filt20dim, dims = 1:20)
  31. Stella_combined_filt_gastrula <- FindClusters(Stella_combined_filt_gastrula, resolution = 1)
  32. Stella_combined_filt_gastrula <- RunUMAP(Stella_combined_filt_gastrula, dims = 1:20, verbose = FALSE)
  34. DimPlot(Stella_combined_filt_gastrula, reduction = "umap", label = FALSE, pt.size = 0.005, label.size = 5)
  35. ```
  37. ```{r}
  38. library(magrittr)
  39. [email hidden] %<>%
  40. mutate(CombinedCluster = case_when(seurat_clusters %in% c(13, 14, 27) ~ "Gland",
  41. seurat_clusters %in% c(25, 16, 20, 29, 30, 23, 10, 19) ~ "Cnidocyte",
  42. seurat_clusters %in% c(11) ~ "Neural Progenitor Cell",
  43. seurat_clusters %in% c(17, 18, 28) ~ "Neuron",
  44. seurat_clusters %in% c(5) ~ "Pharyngeal Ectoderm",
  45. seurat_clusters %in% c(3, 7, 26) ~ "Mesendoderm",
  46. seurat_clusters %in% c(1, 2, 24, 6, 9, 22, 8, 21, 15, 0, 4, 12) ~ "Trunk/Aboral Ectoderm"))
  47. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  48. Stella_combined_filt_gastrula <- SetIdent(Stella_combined_filt_gastrula, value = [email hidden]$CombinedCluster)
  49. ```
  51. #Need to set the active ident. CombinedCluster is the combined and named clusters from above. seurat_clusters is the original cluster numbers
  52. ```{r}
  53. p5<-DimPlot(Stella_combined_filt_gastrula, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 6)
  54. p5
  56. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  57. Stella_combined_filt_gastrula <- SetIdent(Stella_combined_filt_gastrula, value = [email hidden]$CombinedCluster)
  59. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1", "NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g009665000.1", "NV2g000252000.1", "NV2g011441000.1", "NV2g000472000.1", "NV2g003168000.1", "NV2g012726000.1")
  61. library(patchwork)
  63. plt <- DotPlot(Stella_combined_filt_gastrula, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  64. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("NvsoxC"), italic("NvsoxB(2)"), italic("Nvath-like"), italic("NvashA"), italic("Nvelav1"), italic("NvfoxA"), italic("NvsnailA"), italic("Nvwnt2"), italic("Nvsix3/6")))
  65. plt <- plt + ggtitle('Marker genes gastrula annotation')
  66. plt
  67. ```
  69. ```{r, fig.width = 5, fig.height = 4}
  70. # identify your favorite gene in the UMAP plot
  71. goi <- "NV2g009665000.1"
  73. p5 <- FeaturePlot_scCustom(seurat_object = Stella_combined_filt_gastrula, features = goi, pt.size = 0.0005)
  74. p5 <- p5 + ggtitle("NvashA in gastrula")
  75. p5
  76. ```
  79. #early planula
  80. ```{r}
  81. Stella_early_pl
  82. dim(Stella_early_pl)
  83. ncol(Stella_early_pl)
  84. ```
  86. ```{r}
  87. Stella_early_pl <- FindNeighbors(Stella_early_pl_raw, dims = 1:20)
  88. Stella_early_pl <- FindClusters(Stella_early_pl, resolution = 1)
  89. Stella_early_pl <- RunUMAP(Stella_early_pl, dims = 1:20, verbose = FALSE)
  91. DimPlot(Stella_early_pl, reduction = "umap")
  92. ```
  94. ```{r}
  95. library(magrittr)
  96. [email hidden] %<>%
  97. mutate(CombinedCluster = case_when(seurat_clusters %in% c(7, 16) ~ "Gland",
  98. seurat_clusters %in% c(18, 19, 20, 11) ~ "Cnidocyte",
  99. seurat_clusters %in% c(14, 12) ~ "Neural Progenitor Cell",
  100. seurat_clusters %in% c(15, 22) ~ "Neuron",
  101. seurat_clusters %in% c(8, 24, 21) ~ "Pharyngeal Ectoderm",
  102. seurat_clusters %in% c(6) ~ "Mesendoderm",
  103. seurat_clusters %in% c(0, 1, 2, 3, 4, 5, 9, 10, 13, 17, 23) ~ "Trunk/Aboral Ectoderm"))
  104. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  105. Stella_early_pl <- SetIdent(Stella_early_pl, value = [email hidden]$CombinedCluster)
  106. ```
  108. ```{r}
  109. p5<-DimPlot(Stella_early_pl, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 6)
  110. p5
  112. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  113. Stella_early_pl <- SetIdent(Stella_early_pl, value = [email hidden]$CombinedCluster)
  115. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1", "NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g009665000.1", "NV2g000252000.1", "NV2g011441000.1", "NV2g000472000.1", "NV2g003168000.1", "NV2g012726000.1")
  117. plt <- DotPlot(Stella_early_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  118. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("NvsoxC"), italic("NvsoxB(2)"), italic("Nvath-like"), italic("NvashA"), italic("Nvelav1"), italic("NvfoxA"), italic("NvsnailA"), italic("Nvwnt2"), italic("Nvsix3/6")))
  119. plt <- plt + ggtitle('Marker genes early planula annotation')
  120. plt
  121. ```
  123. ```{r}
  124. goi <- "NV2g009665000.1" #
  125. p5 <- FeaturePlot_scCustom(seurat_object = Stella_early_pl, features = goi, pt.size = 0.0005)
  126. p5 <- p5 + ggtitle("NvashA in early planula")
  127. p5
  128. ```
  131. #mid planula
  132. ```{r}
  133. [email hidden] %<>%
  134. mutate(CombinedCluster = case_when(seurat_clusters %in% c(10) ~ "Gland",
  135. seurat_clusters %in% c(17, 15, 9) ~ "Cnidocyte",
  136. seurat_clusters %in% c(6, 16) ~ "Neural Progenitor Cell",
  137. seurat_clusters %in% c(12, 20) ~ "Neuron",
  138. seurat_clusters %in% c(5, 19, 18) ~ "Pharyngeal Ectoderm",
  139. seurat_clusters %in% c(13) ~ "Mesendoderm",
  140. seurat_clusters %in% c(0, 1, 2, 3, 4, 7, 8, 11, 14) ~ "Trunk/Aboral Ectoderm"))
  141. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  142. Stella_mid_pl <- SetIdent(Stella_mid_pl, value = [email hidden]$CombinedCluster)
  144. ```
  146. ```{r}
  147. p5<-DimPlot(Stella_mid_pl, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 6)
  148. p5
  150. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  151. Stella_mid_pl <- SetIdent(Stella_mid_pl, value = [email hidden]$CombinedCluster)
  153. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1", "NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g009665000.1", "NV2g000252000.1", "NV2g011441000.1", "NV2g000472000.1", "NV2g003168000.1", "NV2g012726000.1")
  155. plt <- DotPlot(Stella_mid_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  156. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("NvsoxC"), italic("NvsoxB(2)"), italic("Nvath-like"), italic("NvashA"), italic("Nvelav1"), italic("NvfoxA"), italic("NvsnailA"), italic("Nvwnt2"), italic("Nvsix3/6")))
  157. plt <- plt + ggtitle('Marker genes mid planula annotation')
  158. plt
  159. ```
  161. ```{r}
  162. goi <- "NV2g009665000.1"
  163. p5 <- FeaturePlot_scCustom(seurat_object = Stella_mid_pl, features = goi, pt.size = 0.0005)
  164. p5 <- p5 + ggtitle("NvashA in mid planula")
  165. p5
  166. ```
  168. #late planula
  169. #grouping clusters and cluster identification
  170. ```{r}
  171. library(magrittr)
  172. [email hidden] %<>%
  173. mutate(CombinedCluster = case_when(seurat_clusters %in% c(10, 31) ~ "Gland",
  174. seurat_clusters %in% c(14, 21, 20, 22, 30, 18) ~ "Cnidocyte",
  175. seurat_clusters %in% c(13, 19) ~ "Neural Progenitor Cell",
  176. seurat_clusters %in% c(3, 32, 28, 27, 23) ~ "Neuron",
  177. seurat_clusters %in% c(8, 12, 25, 17, 15) ~ "Pharyngeal Ectoderm",
  178. seurat_clusters %in% c(7) ~ "Mesendoderm",
  179. seurat_clusters %in% c(0, 1, 2, 4, 5, 6, 9, 11, 16, 24, 26, 29) ~ "Trunk/Aboral Ectoderm"))
  180. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  181. Stella_late_pl <- SetIdent(Stella_late_pl, value = [email hidden]$CombinedCluster)
  183. ```
  185. #late planula- annotation figure for pub
  186. ```{r}
  187. p5<-DimPlot(Stella_late_pl, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 6)
  188. p5
  190. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  191. Stella_late_pl <- SetIdent(Stella_late_pl, value = [email hidden]$CombinedCluster)
  193. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1", "NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g009665000.1", "NV2g000252000.1", "NV2g011441000.1", "NV2g000472000.1", "NV2g003168000.1", "NV2g012726000.1")
  195. plt <- DotPlot(Stella_late_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  196. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("NvsoxC"), italic("NvsoxB(2)"), italic("Nvath-like"), italic("NvashA"), italic("Nvelav1"), italic("NvfoxA"), italic("NvsnailA"), italic("Nvwnt2"), italic("Nvsix3/6")))
  197. plt <- plt + ggtitle('Marker genes late planula annotation')
  198. plt
  199. ```
  201. ```{r}
  202. goi <- "NV2g009665000.1"
  203. p5 <- FeaturePlot_scCustom(seurat_object = Stella_late_pl, features = goi, pt.size = 0.0005)
  204. p5 <- p5 + ggtitle("NvashA in late planula")
  205. p5
  206. ```
  209. # Figure 1
  210. # new ashA analysis on subsetted objects
  211. ```{r}
  212. DATA_DIR <- "ashA_combined/data"
  213. RESULTS_DIR <- "ashA_combined/results"
  214. ```
  217. # Subsetting NvashA + cells from each timepoint
  218. ```{r}
  219. #gastrula
  221. R1 <- RidgePlot(object = Stella_combined_filt_gastrula, features = 'NV2g009665000.1')
  222. R1 <- R1 + ggtitle('NvashA expression levels across each cluster at gastrula')
  223. R1
  225. ashA_gastrula <- subset(Stella_combined_filt_gastrula, subset = `NV2g009665000.1` > 0.65)
  226. dim(ashA_gastrula); dim(Stella_combined_filt_gastrula) # cells express NvashA
  227. save(ashA_gastrula, file = file.path(DATA_DIR, "ashA_gastrula.Rda"))
  228. ```
  230. ```{r}
  231. #early planula
  233. R1 <- RidgePlot(object = Stella_early_pl, features = 'NV2g009665000.1')
  234. R1 <- R1 + ggtitle('NvashA expression levels across each cluster at early planula')
  235. R1
  237. ashA_early_pl <- subset(Stella_early_pl, subset = `NV2g009665000.1` > 0.65)
  238. dim(ashA_early_pl); dim(Stella_early_pl) # cells express NvashA
  239. save(ashA_early_pl, file = file.path(DATA_DIR, "ashA_early_pl.Rda"))
  240. ```
  242. ```{r}
  243. #mid planula
  245. R1 <- RidgePlot(object = Stella_mid_pl, features = 'NV2g009665000.1')
  246. R1 <- R1 + ggtitle('NvashA expression levels across each cluster at mid planula')
  247. R1
  249. ashA_mid_pl <- subset(Stella_mid_pl, subset = `NV2g009665000.1` > 0.65)
  250. dim(ashA_mid_pl); dim(Stella_mid_pl) # cells express NvashA
  251. save(ashA_mid_pl, file = file.path(DATA_DIR, "ashA_mid_pl.Rda"))
  252. ```
  254. ```{r}
  255. #late planula
  257. R1 <- RidgePlot(object = Stella_late_pl, features = 'NV2g009665000.1')
  258. R1 <- R1 + ggtitle('NvashA expression levels across each cluster at late planula')
  259. R1
  261. ashA_late_pl <- subset(Stella_late_pl, subset = `NV2g009665000.1` > 0.65)
  262. dim(ashA_late_pl); dim(Stella_late_pl) # cells express NvashA
  263. save(ashA_late_pl, file = file.path(DATA_DIR, "ashA_late_pl.Rda"))
  264. ```
  266. # merge all files and make "timepoint" identity
  267. ```{r}
  268. ashA_combinedv4 <-merge(x = ashA_gastrula, y = list(ashA_early_pl, ashA_mid_pl, ashA_late_pl), add.cell.ids = c("gastrula", "early planula", "mid planula", "late planula"))
  269. ```
  271. ```{r}
  272. ids <- rownames([email hidden])
  273. table(stringr::str_split_fixed(string = ids, pattern = "_", n = 3)[,1])
  274. timepoint <- stringr::str_split_fixed(string = ids, pattern = "_", n = 3)[,1]
  275. [email hidden]$timepoint <- timepoint
  277. Idents(ashA_combinedv4) <- "timepoint"
  278. Idents(ashA_combinedv4)
  279. ```
  281. #analysis of merged NvashA subset data
  282. ```{r}
  283. ashA_combinedv4_20dims <- FindNeighbors(ashA_combinedv4, dims = 1:20)
  284. ashA_combinedv4_20dims <- FindClusters(ashA_combinedv4_20dims, resolution = 1)
  285. head(Idents(ashA_combinedv4_20dims), 10)
  286. ashA_combinedv4_20dims <- RunUMAP(ashA_combinedv4_20dims, dims = 1:20, verbose = FALSE)
  288. p1.ashA <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "timepoint", label = FALSE, pt.size = 0.1, label.size = 5)
  289. p2.ashA <- DimPlot(ashA_combinedv4_20dims, reduction = "umap",group.by = "orig.ident" , pt.size = 0.1)
  290. p3.ashA <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", label = TRUE, pt.size = 0.1, label.size = 8)
  291. p4.ashA <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 8)
  293. p1.ashA
  294. p2.ashA
  295. p3.ashA
  296. p4.ashA
  297. ```
  299. # figure for annotation of NvashA clusters - DotPlot
  300. ```{r fig.width = 3, fig.height = 3}
  302. ashA_combinedv4_20dims <- SetIdent(ashA_combinedv4_20dims, value = [email hidden]$seurat_clusters)
  303. [email hidden] <- factor([email hidden],
  304. levels = c("16", "4", "3", "10", "12", "1", "15", "20", "6", "9", "11", "13", "5",
  305. "2", "0", "7",
  306. "14", "18", "17", "22", "21", "19", "8"))
  308. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1", "NV2g010927000.1", "NV2g018311000.1", "NV2g007706000.1", "NV2g001852000.1", "NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g000252000.1","NV2g002719000.1", "NV2g015023000.1")
  311. plt <- DotPlot(ashA_combinedv4_20dims, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  312. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("Nvshk1"), italic("Nv2.18311"), italic("Nvnh2l1-like-1"), italic("Nvdmrt-E"), italic("NvsoxC"), italic("NvsoxB2"), italic("Nvath-like"), italic("Nvelav1"), italic("NvTBA1-like7"), italic("Nve429")))
  313. plt <- plt + ggtitle('Marker genes across each cluster in NvashA subset')
  314. plt
  315. ```
  317. #Grouping clusters - NvashA subset
  318. #Need to set the active ident. CombinedCluster is the combined and named clusters from above. seurat_clusters is the original cluster numbers
  319. ```{r}
  320. library(magrittr)
  321. [email hidden] %<>%
  322. mutate(CombinedCluster = case_when(seurat_clusters %in% c(12, 10, 3, 4) ~ "cnidocyte",
  323. seurat_clusters %in% c(16) ~ "gland",
  324. seurat_clusters %in% c(1, 15, 20) ~ "cni.gast",
  325. seurat_clusters %in% c(6, 9, 13, 11, 5) ~ "neuron.gast-pl",
  326. seurat_clusters %in% c(2, 0, 7, 14, 17, 18, 19, 21, 22) ~ "neuron.pl",
  327. #seurat_clusters %in% c(8) ~ "unknown"
  328. ))
  329. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("gland", "cnidocyte", "cni.gast", "neuron.gast-pl", "neuron.pl"))
  330. ashA_combinedv4_20dims <- SetIdent(ashA_combinedv4_20dims, value = [email hidden]$CombinedCluster)
  331. ```
  333. ```{r}
  334. p5.ashA_20<-DimPlot(ashA_combinedv4_20dims, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 6)
  335. p5.ashA_20
  336. ```
  339. # Supplemental Figure 2
  340. # make UMAPs for all genes at once
  341. ```{r}
  342. list_of_genes <- c("NV2g006608000.1", "NV2g004477000.1", "NV2g016402000.1", "NV2g000252000.1","NV2g002719000.1", "NV2g015023000.1")
  343. list_of_names <- c( "Nvath-like", "NvsoxB(2)", "NvsoxC", "Nvelav1", "Nvtba1-like7", "NVE490")
  345. if(length(list_of_names) == length(list_of_genes)){
  346. for(i in 1:length(list_of_genes)){
  347. plt <- FeaturePlot_scCustom(ashA_combinedv4_20dims, features = list_of_genes[i], pt.size = 0.005)
  348. plt <- plt + ggtitle(list_of_names[i])
  349. print(plt)
  350. ggsave(plt, file = paste0("/ashA_combined/umaps/", list_of_names[i], "-umap.jpg"), device = "jpg", width = 5, height = 4)
  351. }
  352. } else {
  353. message("Error! - fix your names/ids!")
  354. }
  355. ```
  357. # Extract colors for timepoints
  358. ```{r}
  360. p1.ashA <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "timepoint", label = FALSE, pt.size = 0.1, label.size = 5)
  362. # extract colors
  363. p <- ggplot_build(p1.ashA)
  364. colors <- unique(p$data[[1]]$colour)
  365. library(magrittr)
  366. [email hidden] %<>%
  367. mutate(colorG = ifelse(timepoint == "gastrula", colors[1], "grey50"),
  368. colorEP = ifelse(timepoint == "early planula", colors[2], "grey50"),
  369. colorMP = ifelse(timepoint == "mid planula", colors[3], "grey50"),
  370. colorLP = ifelse(timepoint == "late planula", colors[4], "grey50"))
  372. p1.ashA_G <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "colorG",
  373. label = FALSE, pt.size = 0.1, label.size = 5, cols = c(colors[1], "grey80"))
  374. p1.ashA_G <- p1.ashA_G + NoLegend() + ggtitle("gastrula")
  375. p1.ashA_G
  376. ggsave(filename = "p1ashA_G.jpg", plot = p1.ashA_24, height = 4, width = 5)
  379. p1.ashA_EP <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "colorEP",
  380. label = FALSE, pt.size = 0.1, label.size = 5, cols = c(colors[2], "grey80"))
  381. p1.ashA_EP <- p1.ashA_EP + NoLegend() + ggtitle("early planula")
  382. p1.ashA_EP
  383. ggsave(filename = "p1ashA_EP.jpg", plot = p1.ashA_EP, height = 4, width = 5)
  386. p1.ashA_MP <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "colorMP",
  387. label = FALSE, pt.size = 0.1, label.size = 5, cols = c(colors[3], "grey80"))
  388. p1.ashA_MP <- p1.ashA_MP + NoLegend() + ggtitle("mid planula")
  389. p1.ashA_MP
  390. ggsave(filename = "p1ashA_MP.jpg", plot = p1.ashA_MP, height = 4, width = 5)
  393. p1.ashA_LP <- DimPlot(ashA_combinedv4_20dims, reduction = "umap", group.by = "colorLP",
  394. label = FALSE, pt.size = 0.1, label.size = 5, cols = c(colors[4], "grey80"))
  395. p1.ashA_LP <- p1.ashA_LP + NoLegend() + ggtitle("late planula")
  396. p1.ashA_LP
  397. ggsave(filename = "p1ashA_LP.jpg", plot = p1.ashA_LP, height = 4, width = 5)
  398. ```
  401. # Supplemental Figure 3
  402. # Find DEGs
  403. ```{r}
  404. markers_ashA_combinedv4_20dims <- FindAllMarkers(ashA_combinedv4_20dims)
  405. library(rio)
  406. library(dplyr)
  407. markers_ashA_combinedv4_20dims %>%
  408. tibble::rownames_to_column("gene_id") %>%
  409. rio::export(., file = "/markers_ashA_combinedv4_20dims.csv")
  410. ```
  413. # Determine if the markers are unique markers. Genes = all neural from ashA FindAllMarkers top genes
  414. ```{r fig.width = 8, fig.height = 3}
  416. ashA_combinedv4_20dims <- SetIdent(ashA_combinedv4_20dims, value = [email hidden]$seurat_clusters)
  417. [email hidden] <- factor([email hidden],
  418. levels = c("16", "4", "3", "10", "12", "1", "15", "20", "6", "9", "11", "13", "5",
  419. "2", "0", "7",
  420. "14", "18", "17", "22", "19", "21", "8"))
  422. list_of_genes <- c("NV2g016299000.1", "NV2g014774000.1", "NV2g012706000.1", "NV2g005875000.1", "NV2g014348000.1", #6
  423. "NV2g021230000.1", "NV2g014384000.1", "NV2g003558000.1", "NV2g023233000.1", "NV2g018400000.1", #9
  424. "NV2g011343000.1", "NV2g024863000.1", "NV2g025171000.1", "NV2g015162000.1", "NV2g024857000.1", #11
  425. "NV2g011772000.1", "NV2g014153000.1", "NV2g025073000.1", "NV2g004915000.1", "NV2g021248000.1", #13
  426. "NV2g021302000.1", "NV2g021840000.1", "NV2g024129000.1", "NV2g023280000.1", "NV2g001446000.1", #5
  427. "NV2g014792000.1", "NV2g003097000.1", "NV2g016402000.1", "NV2g013863000.1", "NV2g004360000.1", #2
  428. "NV2g007811000.1", "NV2g023729000.1", "NV2g010520000.1", "NV2g005103000.1", "NV2g003747000.1", #0
  429. "NV2g024838000.1", "NV2g017174000.1", "NV2g008406000.1", "NV2g023195000.1", "NV2g016035000.1", #7
  430. "NV2g007753000.1", "NV2g001757000.1", "NV2g004704000.1", "NV2g004096000.1", "NV2g007867000.1", #14
  431. "NV2g015930000.1", "NV2g002620000.1", "NV2g002627000.1", "NV2g002623000.1", "NV2g020225000.1", #18
  432. "NV2g001999000.1", "NV2g011608000.1", "NV2g015217000.1", "NV2g015222000.1", "NV2g017107000.1", #17
  433. "NV2g021058000.1", "NV2g013128000.1", "NV2g016346000.1", "NV2g014663000.1", "NV2g009262000.1", #22
  434. "NV2g018588000.1", "NV2g003875000.1", "NV2g015362000.1", "NV2g015156000.1", "NV2g015982000.1", #19
  435. "NV2g001204000.1", "NV2g001205000.1", "NV2g007735000.1", "NV2g018362000.1", "NV2g020402000.1", #21
  436. "NV2g014532000.1", "NV2g002896000.1", "NV2g013603000.1", "NV2g004414000.1", "NV2g009834000.1" #8
  437. )
  439. plt4 <- DotPlot(ashA_combinedv4_20dims, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  441. plt4 <- plt4 + scale_x_discrete(labels = expression(
  442. italic("Nvprga"), italic("Nv2.14774"), italic("Nv2.12706"), italic("Nv2.5875"), italic("Nvpapss-like-1"), #6
  443. italic("Nvats6-like-7"), italic("Nv2.14384"), italic("Nv2.3558"), italic("Nv2.23233"), italic("Nvdd3-like-2"), #9
  444. italic("Nvaa2BR-like-2"), italic("Nvmuc1-like-2"), italic("Nvegfl8-like-1"), italic("Nv2.15162"), italic("Nvvwa7-like-2"), #11
  445. italic("Nvblml2-like-1"), italic("Nv2.14153"), italic("Nvaa3R-like-7"), italic("Nvtaar6-like-2"), italic("Nv2.21248"), #13
  446. italic("Nvprga-R.199"), italic("Nv2.21840"), italic("Nvshb-like-6"), italic("Nvgr101-like-10"), italic("Nv2.1446"), #5
  447. italic("Nvrl12-like-1"), italic("Nvrs3A-like-1"), italic("NvsoxC"), italic("Nv2.13863"), italic("Nvrl13-like-1"), #2
  448. italic("Nv2.7811"), italic("Nv2.23729"), italic("Nv2.10520"), italic("Nv2.5103"), italic("Nv2.3747"), #0
  449. italic("Nv2.24838"), italic("Nvaebp1-like-1"), italic("Nvats7-like-1"), italic("Nvlox5-like-3"), italic("Nvnas4-like-2"), #7
  450. italic("Nvviaat-like-40"), italic("Nvfez"), italic("NvCNTP1-like-14"), italic("NvSC6A1-like-9"), italic("Nvbtn1-like-6"), #14
  451. italic("Nv2.15930"), italic("Nv2.2620"), italic("Nv2.2627"), italic("Nv2.2623"), italic("Nv2.20225"), #18
  452. italic("Nvviaat-like-31"), italic("Nvbha15-like-1"), italic("Nv2.15217"), italic("Nv2.15222"), italic("Nvnr1BA-like-1"), #17
  453. italic("Nv2.21058"), italic("Nvaa2AR-like-3"), italic("Nvopn4-like-13"), italic("Nv2.14663"), italic("Nv2.9262"), #22
  454. italic("Nvach10-like-11"), italic("Nvqrfpr-like-31"), italic("Nvhce2-like-1"), italic("Nv.15156"), italic("Nv2.15982"), #19
  455. italic("Nv2.1204"), italic("Nv2.1205"), italic("Nvopn4B-like-6"), italic("Nvadrb2-like-40"), italic("Nvcckar-like-2"), #21
  456. italic("Nv2.14532"), italic("Nvtda6-like-1"), italic("Nvgp2-like-44"), italic("Nv2.4414"), italic("Nv2.9834") #8
  457. ))
  458. plt4 <- plt4 + ggtitle('top 5 DEGs from each cluster on NvashA subset')
  459. plt4
  460. ggsave(plt4, file = "/top_markers_ashA_DEGs.jpg", width = 20, height =8)
  461. ```
  463. # top 5 NvashA markers (DEGs) on late planula dataset
  464. ```{r}
  465. library(magrittr)
  466. [email hidden] %<>%
  467. mutate(CombinedCluster = case_when(seurat_clusters %in% c(10, 31) ~ "Gland",
  468. seurat_clusters %in% c(14, 21, 20, 22, 30, 18) ~ "Cnidocyte",
  469. seurat_clusters %in% c(13, 19) ~ "Neural Progenitor Cell",
  470. seurat_clusters %in% c(3, 32, 28, 27, 23) ~ "Neuron",
  471. seurat_clusters %in% c(8, 12, 25, 17, 15) ~ "Pharyngeal Ectoderm",
  472. seurat_clusters %in% c(7) ~ "Mesendoderm",
  473. seurat_clusters %in% c(0, 1, 2, 4, 5, 6, 9, 11, 16, 24, 26, 29) ~ "Trunk/Aboral Ectoderm"))
  474. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  475. Stella_late_pl <- SetIdent(Stella_late_pl, value = [email hidden]$CombinedCluster)
  477. ```
  479. #Need to set the active ident. CombinedCluster is the combined and named clusters. seurat_clusters is the original cluster numbers ***used in paper
  480. ```{r, fig.width = 20, fig.height = 5}
  482. [email hidden]$CombinedCluster <- factor([email hidden]$CombinedCluster, levels = c("Gland", "Cnidocyte", "Neural Progenitor Cell", "Neuron", "Pharyngeal Ectoderm", "Mesendoderm", "Trunk/Aboral Ectoderm"))
  483. Stella_late_pl <- SetIdent(Stella_late_pl, value = [email hidden]$CombinedCluster)
  485. list_of_genes <- c("NV2g016299000.1", "NV2g014774000.1", "NV2g012706000.1", "NV2g005875000.1", "NV2g014348000.1", #6
  486. "NV2g021230000.1", "NV2g014384000.1", "NV2g003558000.1", "NV2g023233000.1", "NV2g018400000.1", #9
  487. "NV2g011343000.1", "NV2g024863000.1", "NV2g025171000.1", "NV2g015162000.1", "NV2g024857000.1", #11
  488. "NV2g011772000.1", "NV2g014153000.1", "NV2g025073000.1", "NV2g004915000.1", "NV2g021248000.1", #13
  489. "NV2g021302000.1", "NV2g021840000.1", "NV2g024129000.1", "NV2g023280000.1", "NV2g001446000.1", #5
  490. "NV2g014792000.1", "NV2g003097000.1", "NV2g016402000.1", "NV2g013863000.1", "NV2g004360000.1", #2
  491. "NV2g007811000.1", "NV2g023729000.1", "NV2g010520000.1", "NV2g005103000.1", "NV2g003747000.1", #0
  492. "NV2g024838000.1", "NV2g017174000.1", "NV2g008406000.1", "NV2g023195000.1", "NV2g016035000.1", #7
  493. "NV2g007753000.1", "NV2g001757000.1", "NV2g004704000.1", "NV2g004096000.1", "NV2g007867000.1", #14
  494. "NV2g015930000.1", "NV2g002620000.1", "NV2g002627000.1", "NV2g002623000.1", "NV2g020225000.1", #18
  495. "NV2g001999000.1", "NV2g011608000.1", "NV2g015217000.1", "NV2g015222000.1", "NV2g017107000.1", #17
  496. "NV2g021058000.1", "NV2g013128000.1", "NV2g016346000.1", "NV2g014663000.1", "NV2g009262000.1", #22
  497. "NV2g018588000.1", "NV2g003875000.1", "NV2g015362000.1", "NV2g015156000.1", "NV2g015982000.1", #19
  498. "NV2g001204000.1", "NV2g001205000.1", "NV2g007735000.1", "NV2g018362000.1", "NV2g020402000.1", #21
  499. "NV2g014532000.1", "NV2g002896000.1", "NV2g013603000.1", "NV2g004414000.1", "NV2g009834000.1" #8
  500. )
  503. plt <- DotPlot(Stella_late_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  504. plt <- plt + scale_x_discrete(labels = expression(
  505. italic("Nvprga"), italic("Nv2.14774"), italic("Nv2.12706"), italic("Nv2.5875"), italic("Nvpapss-like-1"), #6
  506. italic("Nvats6-like-7"), italic("Nv2.14384"), italic("Nv2.3558"), italic("Nv2.23233"), italic("Nvdd3-like-2"), #9
  507. italic("Nvaa2BR-like-2"), italic("Nvmuc1-like-2"), italic("Nvegfl8-like-1"), italic("Nv2.15162"), italic("Nvvwa7-like-2"), #11
  508. italic("Nvblml2-like-1"), italic("Nv2.14153"), italic("Nvaa3R-like-7"), italic("Nvtaar6-like-2"), italic("Nv2.21248"), #13
  509. italic("Nvprga-R.199"), italic("Nv2.21840"), italic("Nvshb-like-6"), italic("Nvgr101-like-10"), italic("Nv2.1446"), #5
  510. italic("Nvrl12-like-1"), italic("Nvrs3A-like-1"), italic("NvsoxC"), italic("Nv2.13863"), italic("Nvrl13-like-1"), #2
  511. italic("Nv2.7811"), italic("Nv2.23729"), italic("Nv2.10520"), italic("Nv2.5103"), italic("Nv2.3747"), #0
  512. italic("Nv2.24838"), italic("Nvaebp1-like-1"), italic("Nvats7-like-1"), italic("Nvlox5-like-3"), italic("Nvnas4-like-2"), #7
  513. italic("Nvviaat-like-40"), italic("Nvfez"), italic("NvCNTP1-like-14"), italic("NvSC6A1-like-9"), italic("Nvbtn1-like-6"), #14
  514. italic("Nv2.15930"), italic("Nv2.2620"), italic("Nv2.2627"), italic("Nv2.2623"), italic("Nv2.20225"), #18
  515. italic("Nvviaat-like-31"), italic("Nvbha15-like-1"), italic("Nv2.15217"), italic("Nv2.15222"), italic("Nvnr1BA-like-1"), #17
  516. italic("Nv2.21058"), italic("Nvaa2AR-like-3"), italic("Nvopn4-like-13"), italic("Nv2.14663"), italic("Nv2.9262"), #22
  517. italic("Nvach10-like-11"), italic("Nvqrfpr-like-31"), italic("Nvhce2-like-1"), italic("Nv.15156"), italic("Nv2.15982"), #19
  518. italic("Nv2.1204"), italic("Nv2.1205"), italic("Nvopn4B-like-6"), italic("Nvadrb2-like-40"), italic("Nvcckar-like-2"), #21
  519. italic("Nv2.14532"), italic("Nvtda6-like-1"), italic("Nvgp2-like-44"), italic("Nv2.4414"), italic("Nv2.9834") #8
  520. ))
  521. plt <- plt + ggtitle('top 5 subset markers on late planula')
  522. plt
  523. ggsave(plt, file = "/top_markers_LP_dotplot.jpg", width = 20, height = 5)
  524. ```
  527. # Supplemental Figure 4
  528. # Find DEGs from Cole et al. 2024's neuroglandular dataset
  529. ```{r}
  530. markers_neurogland.all <- FindAllMarkers(neurogland.all)
  531. library(rio)
  532. library(dplyr)
  533. markers_neurogland.all %>%
  534. tibble::rownames_to_column("gene_id") %>%
  535. rio::export(., file = "/markers_neurogland.all.csv")
  536. ```
  538. # top genes from Alison Cole's (Cole et al., 2024) neugoglandular datawset on NvashA subset
  539. ```{r fig.width = 9, fig.height = 3}
  541. ashA_combinedv4_20dims <- SetIdent(ashA_combinedv4_20dims, value = [email hidden]$seurat_clusters)
  542. [email hidden] <- factor([email hidden],
  543. levels = c("16", "4", "3", "10", "12", "1", "15", "20", "6", "9", "11", "13", "5",
  544. "2", "0", "7",
  545. "14", "18", "17", "22", "19", "21", "8"))
  547. list_of_genes <- c("NV2g011600000.1", "NV2g016299000.1", "NV2g000651000.1", "NV2g016823000.1", "NV2g011343000.1", #N1.L2
  548. "NV2g014153000.1", "NV2g025073000.1", "NV2g011772000.1", "NV2g004915000.1", "NV2g009042000.1", #N1.L1
  549. "NV2g012348000.1", "NV2g021230000.1","NV2g003558000.1", "NV2g021840000.1", "NV2g024129000.1", #N1.L3
  550. "NV2g015787000.1", "NV2g016037000.1","NV2g011524000.1", "NV2g007282000.1", "NV2g023984000.1", #N2.1
  551. "NV2g020120000.1", "NV2g008437000.1","NV2g007753000.1", "NV2g010940000.1", "NV2g017198000.1", #N2.2
  552. "NV2g018186000.1", "NV2g001757000.1", "NV2g023678000.1", "NV2g002044000.1","NV2g000371000.1", #N2.3
  553. "NV2g016206000.1", "NV2g010682000.1","NV2g025643000.1", "NV2g006888000.1", "NV2g004704000.1", #N2.4
  554. "NV2g002620000.1", "NV2g002627000.1", "NV2g002623000.1", "NV2g000952000.1", "NV2g005984000.1", #N1.g.early
  555. "NV2g006171000.1", "NV2g015217000.1", "NV2g001449000.1", "NV2g014600000.1", "NV2g008266000.1", #N2.g1
  556. "NV2g021058000.1", "NV2g016346000.1", "NV2g014663000.1", "NV2g011930000.1", "NV2g014773000.1", #NGD.1
  557. "NV2g018872000.1", "NV2g001204000.1", "NV2g001205000.1", "NV2g001202000.1", "NV2g020328000.1", #N1.2
  558. "NV2g019212000.1", "NV2g008558000.1", "NV2g008740000.1", "NV2g019209000.1", "NV2g006433000.1",
  559. "NV2g017383000.1", "NV2g002154000.1", "NV2g022717000.1", "NV2g024546000.1", "NV2g013148000.1"
  560. )
  563. plt9 <- DotPlot(ashA_combinedv4_20dims, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  564. plt9 <- plt9 + ggtitle("subtype markers on ashA subset")
  565. plt9 <- plt9 + scale_x_discrete(labels = expression(
  566. italic("Nv2.11600"), italic("NvPRGa"), italic("NvMFD6B-like-1"), italic("NvAQP2-like-1"), italic("NvAA2BR-like-2"),
  567. italic("Nv2.14153"), italic("NvAA3R-like-7"), italic("NvBLML2-like-1"), italic("NvTAAR6-like-2"), italic("NvTLL-like2"),
  568. italic("NvPR1-like-1"), italic("NvATS6-like-7"), italic("Nv2.3558"), italic("Nv2.21840"), italic("NvSHB-like-6"),
  569. italic("Nv2.15787"), italic("Nv2.16037"), italic("NvTBA3-like-2"), italic("NvFibrinogen-like"), italic("NvEDIL3-like-40"),
  570. italic("NvTMPS6-like-2"), italic("NvQGRFa"), italic("NvVIAAT-like-40"), italic("NvGABR1-like-2"), italic("NvFBCD1-like-22"),
  571. italic("NvGBRG3-like-1"), italic("Nvfez"), italic("NvFXD3A-like-1"), italic("Nv2.2044"), italic("NvGBRR1-like-5"), #N2.3
  572. italic("NvECT-like-2"), italic("NvNLGNY-like-2"), italic("NvTBX5A-like-4"), italic("NvACHA3-like-9"), italic("NvCNTP1-like-14"), #N2.4
  573. italic("Nv2.2620"), italic("Nv2.2627"), italic("Nv2.2623"), italic("NvCP17A-like-28"), italic("NvPGCA-like-20"), #N1.g.early
  574. italic("NvTBA1-like-1"), italic("Nv2.15217"), italic("NvPRDM14d"), italic("NvCALM-like-34"), italic("Nv2.8266"), #N2.g1
  575. italic("Nv2.21058"), italic("NvOPN4-like-13"), italic("Nv2.14663"), italic("NvKCNJ5-like-1"), italic("Nv2.14773"), #NGD.1
  576. italic("NvEmx3"), italic("Nv2.1204"), italic("Nv2.1205"), italic("Nv2.1202"), italic("NvTM11D-like-4"), #N1.2
  577. italic("NvFBP1-like-40"), italic("NvVITRN-like-3"), italic("NvCO7A1-like-1"), italic("NvFBP1-like-50"), italic("NvSAAR1-like-1"),
  578. italic("NvCD63-like-12"), italic("Nv2.2154"), italic("Nv2.22717"), italic("NvDMBT1-like-28"), italic("NvSON-like-3")
  579. ))
  581. plt9
  582. ggsave(plt9, file = "/subtpemarker_dotplot.jpg", width = 20, height =8)
  584. ```
  586. #Nvasha markers on Cole's neuroglandular dataset
  587. ```{r, fig.width = 8, fig.height = 8}
  589. genes <- c(
  590. # "PRGa", "NV2.14774", "NV2.12706", "NV2.5875", "PAPSS-like-1", #6
  591. # "ATS6-like-7", "NV2.14384", "NV2.3558", "NV2.23233", "DD3-like-2", #9
  592. "AA2BR-like-2", "MUC1-like-2", "EGFL8-like-1", "NV2.15162", "VWA7-like-2", #11
  593. "BLML2-like-1", "NV2.14153", "AA3R-like-7", "TAAR6-like-2", "NV2.21248", #13
  594. "PRGa-R.199", "NV2.21840", "SHB-like-6", "GR101-like-10", "NV2.1446", #5
  595. # "RL12-like-1", "RS3A-like-1", "SoxC", "NV2.13863", "RL13-like-1", #2
  596. # "NV2.7811", "NV2.23729", "NV2.10520", "NV2.5103", "NV2.3747", #0
  597. "NV2.24838", "AEBP1-like-1", "ATS7-like-1", "LOX5-like-3", "NAS4-like-2", #7
  598. "VIAAT-like-40", "fez", "CNTP1-like-14", "SC6A1-like-9", "BTN1-like-6", #14
  599. "NV2.15930", "NV2.2620", "NV2.2627", "NV2.2623", "NV2.20225", #18
  600. "VIAAT-like-31", "BHA15-like-1", "NV2.15217", "NV2.15222", "NR1BA-like-1", #17
  601. "NV2.21058", "AA2AR-like-3", "OPN4-like-13", "NV2.14663", "NV2.9262", #22
  602. "ACH10-like-11", "QRFPR-like-31", "HCE2-like-1", "NV2.15156", "NV2.15982", #19
  603. "NV2.1204", "NV2.1205", "OPN4B-like-6", "ADRB2-like-40", "CCKAR-like-2", #21
  604. "NV2.14532", "TDA6-like-1", "GP2-like-44", "NV2.4414", "NV2.9834") #8
  606. plt9 <- DotPlot(neurogland.all, features = genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  607. plt9 <- plt9 + ggtitle("NvashA cluster markers/DEGs on neuroglandular dataset")
  608. plt9 <- plt9 + scale_x_discrete(labels = expression(
  609. #6, #9
  610. italic("Nvaa2BR-like-2"), italic("Nvmuc1-like-2"), italic("Nvegfl8-like-1"), italic("Nv2.15162"), italic("Nvvwa7-like-2"), #11
  611. italic("Nvblml2-like-1"), italic("Nv2.14153"), italic("Nvaa3R-like-7"), italic("Nvtaar6-like-2"), italic("Nv2.21248"), #13
  612. italic("Nvprga-R.199"), italic("Nv2.21840"), italic("Nvshb-like-6"), italic("Nvgr101-like-10"), italic("Nv2.1446"), #5
  613. #2, #0
  614. italic("Nv2.24838"), italic("Nvaebp1-like-1"), italic("Nvats7-like-1"), italic("Nvlox5-like-3"), italic("Nvnas4-like-2"), #7
  615. italic("Nvviaat-like-40"), italic("Nvfez"), italic("NvCNTP1-like-14"), italic("NvSC6A1-like-9"), italic("Nvbtn1-like-6"), #14
  616. italic("Nv2.15930"), italic("Nv2.2620"), italic("Nv2.2627"), italic("Nv2.2623"), italic("Nv2.20225"), #18
  617. italic("Nvviaat-like-31"), italic("Nvbha15-like-1"), italic("Nv2.15217"), italic("Nv2.15222"), italic("Nvnr1BA-like-1"), #17
  618. italic("Nv2.21058"), italic("Nvaa2AR-like-3"), italic("Nvopn4-like-13"), italic("Nv2.14663"), italic("Nv2.9262"), #22
  619. italic("Nvach10-like-11"), italic("Nvqrfpr-like-31"), italic("Nvhce2-like-1"), italic("Nv.15156"), italic("Nv2.15982"), #19
  620. italic("Nv2.1204"), italic("Nv2.1205"), italic("Nvopn4B-like-6"), italic("Nvadrb2-like-40"), italic("Nvcckar-like-2"), #21
  621. italic("Nv2.14532"), italic("Nvtda6-like-1"), italic("Nvgp2-like-44"), italic("Nv2.4414"), italic("Nv2.9834") #8
  622. ))
  623. plt9
  624. ggsave(plt9, file = "/ashAmarkersNGdata.jpg", width = 15, height = 10)
  625. ```
  628. #Supplemental Figure 5
  629. #cnidocyte subtype dotplot - full clusters
  630. ```{r fig.width = 8, fig.height = 3}}
  632. ashA_combinedv4_20dims <- SetIdent(ashA_combinedv4_20dims, value = [email hidden]$seurat_clusters)
  633. [email hidden] <- factor([email hidden],
  634. levels = c("16", "4", "3", "10", "12", "1", "15", "20", "6", "9", "11", "13", "5",
  635. "2", "0", "7",
  636. "14", "18", "17", "22", "21", "19", "8"))
  638. features_CnidoSub <- c("NV2g000363000.1", "NV2g007706000.1", "NV2g007069000.1", "NV2g001852000.1",
  639. "NV2g014792000.1", "NV2g003097000.1", "NV2g005947000.1", "NV2g015684000.1",
  640. "NV2g010927000.1", "NV2g018311000.1", "NV2g009989000.1", "NV2g003202000.1",
  641. "NV2g007207000.1", "NV2g002780000.1", "NV2g007209000.1", "NV2g019217000.1",
  642. "NV2g025705000.1", "NV2g012610000.1", "NV2g009823000.1", "NV2g004514000.1",
  643. "NV2g021967000.1", "NV2g014593000.1", "NV2g014594000.1", "NV2g019221000.1",
  644. "NV2g014815000.1", "NV2g007709000.1",
  645. "NV2g009503000.1", "NV2g008284000.1", "NV2g012900000.1", "NV2g007818000.1",
  646. "NV2g011818000.1", "NV2g012689000.1", "NV2g010115000.1", "NV2g011819000.1",
  647. "NV2g009120000.1", "NV2g000411000.1", "NV2g002044000.1", "NV2g004676000.1",
  648. "NV2g014822000.1", "NV2g018300000.1", "NV2g014159000.1", "NV2g018114000.1",
  649. "NV2g009028000.1", "NV2g015716000.1", "NV2g003814000.1", "NV2g020507000.1")
  652. plt <- DotPlot(ashA_combinedv4_20dims, features = features_CnidoSub) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  653. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmarco-like-2"), italic("Nvnh2l1-like-1"), italic("Nvh3-like-17"), italic("Nvdmrt-E"),
  654. italic("Nvrl12-like-1"), italic("Nvrs3a-like-1"), italic("Nvrs8-like-1",), italic("Nvrs12-like-1"),
  655. italic("Nvshk1"), italic("Nv2.18311"), italic("Nvco6-like-4"), italic("Nvhmcn1-like-87"),
  656. italic("Nvvmp6"), italic("Nvhmcn1-like-89"), italic("Nv2.7209"), italic("Nv2.19217"),
  657. italic("Nv2.25705"), italic("Nvmac1"), italic("Nv2.9823"), italic("Nvcop3-like-1"),
  658. italic("Nv2.21967"), italic("Nv2.14593"), italic("Nv2.14594"), italic("Nv2.19221"),
  659. italic("Nv2.14815"), italic("Nv2.7709"),
  660. italic("Nv2.9503"), italic("Nv2.8284"), italic("Nvcapa-like-1"), italic("Nvggha-like-11"),
  661. italic("Nvco6a-like-6"), italic("Nv2.12689"), italic("Nv2.10115"), italic("NV2.11819"),
  662. italic("Nvnupr1-like-1"), italic("Nvak1a1-like-1"), italic("Nv2.2044"), italic("Nvtraf4-like-19"),
  663. italic("Nv2.14822"), italic("Nv2.18300"), italic("Nvlama2-like-2.1"), italic("Nvisp2-like-2"),
  664. italic("Nvgabr2-like-21"), italic("Nvfbp1-like-2"), italic("Nvnpff2-like-90"), italic("Nvpk1l2-like-31")))
  665. plt <- plt + ggtitle('Cnidocyte subtype markers in NvashA subset')
  666. plt
  667. ggsave(plt, file = "/CnidoSubtypeDotplotv2.jpg", width = 12, height = 8)
  670. #removed NV2g015045000.1,NV2g025210000.1, NV2g018310000.1, NV2g006718000.1, NV2g007385000.1, NV2g016293000.1, NV2g025609000.1, NV2g008572000.1 NV2g014765000.1, NV2g007708000.1, NV2g019873000.1 = not found
  671. ```
  673. #umaps for cnidocyte genes on NvashA subset
  674. ```{r}
  675. list_of_genes <- c("NV2g019749000.1", "NV2g010686000.1", "NV2g010927000.1", "NV2g018311000.1", "NV2g007706000.1", "NV2g001852000.1")
  676. list_of_names <- c("Nvcnido-fos1", "Nvncol3", "Nvshk1", "Nv2.18311", "Nvnh2l1-like-1", "Nvdmrt-E")
  678. if(length(list_of_names) == length(list_of_genes)){
  679. for(i in 1:length(list_of_genes)){
  680. plt <- FeaturePlot_scCustom(ashA_combinedv4_20dims, features = list_of_genes[i], pt.size = 0.005)
  681. plt <- plt + ggtitle(list_of_names[i])
  682. print(plt)
  683. ggsave(plt, file = paste0("/ashA_combined/umaps/", list_of_names[i], "-umap.jpg"), device = "jpg", width = 5, height = 4)
  684. }
  685. } else {
  686. message("Error! - fix your names/ids!")
  687. }
  688. ```
  691. #Supplemental Figure 6
  692. #UMAPs for ashA target genes on ashA subset
  693. ```{r}
  694. list_of_genes <- c("NV2g023192000.1", "NV2g024331000.1", "NV2g005875000.1", "NV2g012875000.1")
  695. list_of_names <- c( "Nvlwamide-like", "NvglraA3-like", "Nvserum amyloid A-like", "Nvabcc4-like")
  697. if(length(list_of_names) == length(list_of_genes)){
  698. for(i in 1:length(list_of_genes)){
  699. plt <- FeaturePlot_scCustom(ashA_combinedv4_20dims, features = list_of_genes[i], pt.size = 0.005)
  700. plt <- plt + ggtitle(list_of_names[i])
  701. print(plt)
  702. }
  703. } else {
  704. message("Error! graphs won't genereate - fix your names/ids!")
  705. }
  706. ```
  708. # NvashA known target genes on gastrula dataset
  709. ```{r}
  710. Stella_combined_filt_gastrula <- SetIdent(Stella_combined_filt_gastrula, value = [email hidden]$seurat_clusters)
  711. [email hidden] <- factor([email hidden],
  712. levels = c("13", "14", "27", "25", "16", "20", "29", "30",
  713. "23", "10", "19", "11",
  714. "17", "18", "28", "5", "3", "7", "26", "1", "2", "24", "6", "9", "22", "8", "21", "15", "0", "4", "12"))
  716. list_of_genes <- c("NV2g009665000.1", "NV2g023192000.1", "NV2g024331000.1", "NV2g005875000.1", "NV2g012875000.1")
  718. plt3 <- DotPlot(Stella_combined_filt_gastrula, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  719. plt3 <- plt3 + scale_x_discrete(labels = expression(italic("NvashA"), italic("Nvlwamide-like"), italic("NvglraA3-like"), italic("Nvserum amyloid A-like"), italic("Nvabcc4-like")))
  720. plt3 <- plt3 + ggtitle('NvashA targets on gastrula dataset')
  721. plt3
  722. ```
  724. # NvashA known target genes on mid planula dataset
  725. ```{r}
  726. Stella_mid_pl <- SetIdent(Stella_mid_pl, value = [email hidden]$seurat_clusters)
  727. [email hidden] <- factor([email hidden],
  728. levels = c("10", "17", "15", "9", "6", "16", "12", "20",
  729. "5", "19", "18", "13",
  730. "0", "1", "2", "3", "4", "7", "8", "11", "14"))
  732. list_of_genes <- c("NV2g009665000.1", "NV2g023192000.1", "NV2g024331000.1", "NV2g005875000.1", "NV2g012875000.1")
  734. plt3 <- DotPlot(Stella_mid_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  735. plt3 <- plt3 + scale_x_discrete(labels = expression(italic("NvashA"), italic("Nvlwamide-like"), italic("NvglraA3-like"), italic("Nvserum amyloid A-like"), italic("Nvabcc4-like")))
  736. plt3 <- plt3 + ggtitle('NvashA targets on mid planula dataset')
  737. plt3
  738. ```
  740. # NvashA known targets on late planula dataset
  741. ```{r}
  742. Stella_late_pl <- SetIdent(Stella_late_pl, value = [email hidden]$seurat_clusters)
  743. [email hidden] <- factor([email hidden],
  744. levels = c("10", "31", "14", "21", "20", "22", "30", "18",
  745. "13", "19", "3", "32", "28", "27", "23",
  746. "8", "12", "25", "17", "15", "7",
  747. "0", "1", "2", "4", "5", "6", "9", "11", "16","24", "26", "29"))
  749. list_of_genes <- c("NV2g009665000.1", "NV2g023192000.1", "NV2g024331000.1", "NV2g005875000.1", "NV2g012875000.1")
  751. plt3 <- DotPlot(Stella_late_pl, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  752. plt3 <- plt3 + scale_x_discrete(labels = expression(italic("NvashA"), italic("Nvlwamide-like"), italic("NvglraA3-like"), italic("Nvserum amyloid A-like"), italic("Nvabcc4-like")))
  753. plt3 <- plt3 + ggtitle('NvashA targets on late planula dataset')
  754. plt3
  755. ```
  758. #FIGURE 2
  759. #Subsetting neual only cells from NvashA+ subset
  760. ```{r}
  761. ashA_neural_sub3 <- subset(x = ashA_combinedv4_20dims, idents = c("6", "9", "11", "13", "5",
  762. "2", "0", "7",
  763. "14", "18", "17", "22", "19", "21"))
  765. save(ashA_neural_sub3,
  766. file = file.path(DATA_DIR, "ashA_neural_sub3.Rda"))
  767. ```
  769. ```{r}
  770. # keep original cluster IDs (the ones it was subset by)
  771. ashA_neural_sub3$cluster_orig <- Idents(ashA_neural_sub3)
  773. # preserve level order from the parent object
  774. ashA_neural_sub3$cluster_orig <- factor(
  775. ashA_neural_sub3$cluster_orig,
  776. levels = levels(ashA_combinedv4_20dims) # keeps original cluster order
  777. )
  779. # (optional) also stash numeric seurat_clusters if present
  780. if ("seurat_clusters" %in% colnames([email hidden])) {
  781. ashA_neural_sub3$seurat_clusters_orig <- ashA_neural_sub3$seurat_clusters
  782. }
  784. table(ashA_neural_sub3$cluster_orig)
  786. ```
  788. ```{r}
  789. # Recalculate within the subset
  790. ashA_neural_sub3 <- NormalizeData(ashA_neural_sub3)
  791. ashA_neural_sub3 <- FindVariableFeatures(ashA_neural_sub3, selection.method = "vst", nfeatures = 2000)
  792. ashA_neural_sub3 <- ScaleData(ashA_neural_sub3)
  793. ashA_neural_sub3 <- RunPCA(ashA_neural_sub3, npcs = 30)
  794. ```
  796. 15 dims
  797. ```{r}
  798. ashA_neural_sub3_15dims <- FindNeighbors(ashA_neural_sub3, dims = 1:15)
  799. ashA_neural_sub3_15dims <- FindClusters(ashA_neural_sub3_15dims, resolution = 1)
  800. head(Idents(ashA_neural_sub3_15dims), 10)
  801. ashA_neural_sub3_15dims <- RunUMAP(ashA_neural_sub3_15dims, dims = 1:15, verbose = FALSE)
  803. p1 <- DimPlot(ashA_neural_sub3_15dims, reduction = "umap", group.by = "timepoint", label = FALSE, pt.size = 0.1, label.size = 5)
  804. p2 <- DimPlot(ashA_neural_sub3_15dims, reduction = "umap",group.by = "orig.ident" , pt.size = 0.1)
  805. p3 <- DimPlot(ashA_neural_sub3_15dims, reduction = "umap", label = TRUE, pt.size = 0.1, label.size = 8)
  806. p4 <- DimPlot(ashA_neural_sub3_15dims, reduction = "umap", label = FALSE, pt.size = 0.1, label.size = 8)
  808. ashA_neural_sub3_15dims$cluster_new <- Idents(ashA_neural_sub3_15dims)
  810. p1
  811. p2
  812. p3
  813. p4
  815. save(ashA_neural_sub3_15dims3,
  816. file = file.path(DATA_DIR, "ashA_neural_sub3_15dims.Rda"))
  818. ```
  820. ```{r}
  821. # plot originals without changing Idents()
  822. p_orig <- DimPlot(
  823. ashA_neural_sub3_15dims,
  824. reduction = "umap",
  825. group.by = "cluster_orig",
  826. label = TRUE, repel = TRUE, pt.size = 0.1
  827. ) + ggtitle("UMAP — Original cluster IDs")
  829. # New clusters (current Idents) for comparison
  830. p_new <- DimPlot(
  831. ashA_neural_sub3_15dims,
  832. reduction = "umap",
  833. label = TRUE, repel = TRUE, pt.size = 0.1
  834. ) + ggtitle("UMAP — New clusters")
  836. p_new; p_orig
  838. ```
  840. #Supplemental Figure 7
  841. #annotation_final_for paper
  842. ```{r fig.width = 9, fig.height = 7}
  843. list_of_genes <- c("NV2g012902000.1", "NV2g018581000.1", "NV2g019749000.1", "NV2g010686000.1","NV2g016402000.1", "NV2g004477000.1", "NV2g006608000.1", "NV2g009665000.1", "NV2g023192000.1", "NV2g000252000.1"
  844. )
  846. plt <- DotPlot(ashA_neural_sub3_15dims, features = list_of_genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  847. plt <- plt + scale_x_discrete(labels = expression(italic("Nvmucin"), italic("Nvnot-likeE"), italic("Nvcnido-fos1"), italic("Nvncol3"), italic("NvsoxC"), italic("NvsoxB(2)"), italic("Nvath-like"), italic("NvashA"), italic("NvLWamide-like"), italic("Nvelav1")
  848. ))
  849. plt <- plt + ggtitle('NvashA neural subset')
  850. plt
  851. ```
  853. #used this code to find all markers
  854. ```{r}
  855. markers_ashA_neural_sub3_15dims <- FindAllMarkers(ashA_neural_sub3_15dims)
  856. library(rio)
  857. library(dplyr)
  858. markers_ashA_neural_sub3_15dims %>%
  859. tibble::rownames_to_column("gene_id") %>%
  860. rio::export(., file = "/markers_ashA_neural_sub3_15dims.csv")
  861. ```
  864. ```{r}
  865. library(dplyr)
  866. library(readr)
  867. library(readxl)
  869. ### 1. Read marker file (from Seurat FindAllMarkers output)
  870. markers <- read_csv("markers_ashA_neural_sub3_15dims.csv")
  872. ### 2. Read annotation file
  873. annotations <- read_csv("manifest_mappings.csv")
  875. ### 3. Extract top 5 genes per cluster *in file order*
  876. top_markers <- markers %>%
  877. group_by(cluster) %>%
  878. mutate(row_in_cluster = row_number()) %>%
  879. filter(row_in_cluster <= 5) %>%
  880. ungroup() %>%
  881. arrange(cluster, row_in_cluster)
  883. ### 4. Join annotations
  884. # top_markers$gene is Nv2gID
  885. # annotations$mapped_geneID is the matching column
  886. top_markers_annot <- top_markers %>%
  887. left_join(annotations, by = c("gene" = "mapped_geneID"))
  889. ### 5. Construct feature vector (Nv2g IDs) and label vector (gene names)
  890. list_of_genes <- top_markers_annot$gene # Seurat features
  891. labels_vec <- top_markers_annot$gene.name # Human-readable names
  893. ### 6. Replace missing names with gene IDs
  894. labels_vec[is.na(labels_vec)] <- list_of_genes[is.na(labels_vec)]
  896. ```
  898. ```{r fig.width = 20, fig.height = 10}
  899. plt5 <- DotPlot(
  900. ashA_neural_sub3_15dims,
  901. features = list_of_genes,
  902. scale = FALSE
  903. ) +
  904. RotatedAxis() +
  905. theme(axis.text.x = element_text(angle = 90)) +
  906. FontSize(10)
  908. plt5 <- plt5 + scale_x_discrete(labels = labels_vec)
  910. plt5 <- plt5 + ggtitle("Top 5 subset markers on NvashA neural subset")
  911. plt5
  912. ```
  914. #used this code to easily see the top 5 genes
  915. ```{r}
  916. # This will PRINT a block:
  917. clusters <- sort(unique(top_markers$cluster))
  919. cat("list_of_genes <- c(\n")
  920. for (cl in clusters) {
  921. genes_cl <- top_markers$gene[top_markers$cluster == cl]
  922. cat(" ", paste0('"', genes_cl, '"', collapse = ", "), ", #", cl, "\n", sep = "")
  923. }
  924. cat(")\n")
  926. ```
  928. ```{r}
  929. # Print human-readable gene names grouped by cluster
  930. clusters <- sort(unique(top_markers_annot$cluster))
  932. cat("labels_vec <- c(\n")
  933. for (cl in clusters) {
  934. names_cl <- top_markers_annot$gene.name[top_markers_annot$cluster == cl]
  935. # Replace any NA names with the Nv2g ID so nothing is blank
  936. names_cl[is.na(names_cl)] <- top_markers_annot$gene[top_markers_annot$cluster == cl][is.na(names_cl)]
  938. cat(" ", paste0('"', names_cl, '"', collapse = ", "), ", #", cl, "\n", sep = "")
  939. }
  940. cat(")\n")
  942. ```
  944. ```{r fig.width = 20, fig.height = 10}
  947. list_of_genes <- c("NV2g021840000.1", "NV2g021302000.1", "NV2g024129000.1", "NV2g001446000.1", "NV2g012348000.1", #0
  948. "NV2g007811000.1", "NV2g012793000.1", "NV2g018662000.1", "NV2g014260000.1", "NV2g010809000.1", #1
  949. "NV2g014348000.1", "NV2g014774000.1", "NV2g001200000.1", "NV2g016299000.1", "NV2g019361000.1", #2
  950. "NV2g018827000.1", "NV2g018825000.1", "NV2g002299000.1", "NV2g002596000.1", "NV2g020327000.1", #3
  951. "NV2g011343000.1", "NV2g023614000.1", "NV2g012875000.1", "NV2g024863000.1", "NV2g003856000.1", #4
  952. "NV2g011441000.1", "NV2g008165000.1", "NV2g011806000.1", "NV2g011442000.1", "NV2g000333000.1", #5
  953. "NV2g011772000.1", "NV2g014153000.1", "NV2g025073000.1", "NV2g004915000.1", "NV2g021248000.1", #6
  954. "NV2g006813000.1", "NV2g006810000.1", "NV2g008789000.1", "NV2g017754000.1", "NV2g006809000.1", #7
  955. "NV2g007753000.1", "NV2g001757000.1", "NV2g004704000.1", "NV2g020120000.1", "NV2g004096000.1", #8
  956. "NV2g013863000.1", "NV2g002898000.1", "NV2g014792000.1", "NV2g006333000.1", "NV2g004360000.1", #9
  957. "NV2g020623000.1", "NV2g025991000.1", "NV2g021230000.1", "NV2g018400000.1", "NV2g023233000.1", #10
  958. "NV2g000160000.1", "NV2g005192000.1", "NV2g004530000.1", "NV2g014941000.1", "NV2g024159000.1", #11
  959. "NV2g011608000.1", "NV2g016903000.1", "NV2g008266000.1", "NV2g017107000.1", "NV2g001999000.1", #12
  960. "NV2g002620000.1", "NV2g002627000.1", "NV2g015930000.1", "NV2g002623000.1", "NV2g009018000.1", #13
  961. "NV2g015362000.1", "NV2g018588000.1", "NV2g003875000.1", "NV2g024629000.1", "NV2g006966000.1", #14
  962. "NV2g000450000.1", "NV2g005143000.1", "NV2g005142000.1", "NV2g006552000.1", "NV2g005046000.1", #15
  963. "NV2g001204000.1", "NV2g001205000.1", "NV2g007735000.1", "NV2g018362000.1", "NV2g022444000.1", #16
  964. "NV2g010279000.1", "NV2g009708000.1", "NV2g014290000.1", "NV2g023195000.1", "NV2g015846000.1", #17
  965. "NV2g021058000.1", "NV2g009262000.1", "NV2g013128000.1", "NV2g014663000.1", "NV2g016346000.1" #18
  966. )
  968. plt4 <- DotPlot(ashA_neural_sub3_15dims, features = list_of_genes, scale = FALSE) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  969. plt4 <- plt4 + scale_x_discrete(labels = c("NV2g021840000.1", "NV2g021302000.1", "NV2g024129000.1", "NV2g001446000.1", "NV2g012348000.1", #0
  970. "NV2g007811000.1", "NV2g012793000.1", "NV2g018662000.1", "NV2g014260000.1", "NV2g010809000.1", #1
  971. "NV2g014348000.1", "NV2g014774000.1", "NV2g001200000.1", "NV2g016299000.1", "NV2g019361000.1", #2
  972. "NV2g018827000.1", "NV2g018825000.1", "NV2g002299000.1", "NV2g002596000.1", "NV2g020327000.1", #3
  973. "NV2g011343000.1", "NV2g023614000.1", "NV2g012875000.1", "NV2g024863000.1", "NV2g003856000.1", #4
  974. "NV2g011441000.1", "NV2g008165000.1", "NV2g011806000.1", "NV2g011442000.1", "NV2g000333000.1", #5
  975. "NV2g011772000.1", "NV2g014153000.1", "NV2g025073000.1", "NV2g004915000.1", "NV2g021248000.1", #6
  976. "NV2g006813000.1", "NV2g006810000.1", "NV2g008789000.1", "NV2g017754000.1", "NV2g006809000.1", #7
  977. "NV2g007753000.1", "NV2g001757000.1", "NV2g004704000.1", "NV2g020120000.1", "NV2g004096000.1", #8
  978. "NV2g013863000.1", "NV2g002898000.1", "NV2g014792000.1", "NV2g006333000.1", "NV2g004360000.1", #9
  979. "NV2g020623000.1", "NV2g025991000.1", "NV2g021230000.1", "NV2g018400000.1", "NV2g023233000.1", #10
  980. "NV2g000160000.1", "NV2g005192000.1", "NV2g004530000.1", "NV2g014941000.1", "NV2g024159000.1", #11
  981. "NV2g011608000.1", "NV2g016903000.1", "NV2g008266000.1", "NV2g017107000.1", "NV2g001999000.1", #12
  982. "NV2g002620000.1", "NV2g002627000.1", "NV2g015930000.1", "NV2g002623000.1", "NV2g009018000.1", #13
  983. "NV2g015362000.1", "NV2g018588000.1", "NV2g003875000.1", "NV2g024629000.1", "NV2g006966000.1", #14
  984. "NV2g000450000.1", "NV2g005143000.1", "NV2g005142000.1", "NV2g006552000.1", "NV2g005046000.1", #15
  985. "NV2g001204000.1", "NV2g001205000.1", "NV2g007735000.1", "NV2g018362000.1", "NV2g022444000.1", #16
  986. "NV2g010279000.1", "NV2g009708000.1", "NV2g014290000.1", "NV2g023195000.1", "NV2g015846000.1", #17
  987. "NV2g021058000.1", "NV2g009262000.1", "NV2g013128000.1", "NV2g014663000.1", "NV2g016346000.1" #18
  988. ))
  989. plt4 <- plt4 + ggtitle('top 5 subset markers on NvashA neural subset')
  990. plt4
  991. ggsave(plt4, file = "/top_markersNeural_dotplot.jpg", width = 20, height = 10)
  992. ```
  994. # Nvasha+ neural cluster markers on Cole et al. 2024's neuroglandular dataset
  995. ```{r, fig.width = 15, fig.height = 10}
  997. genes <- c("NV2.21840", "PRGa-R.199", "SHB-like-6", "NV2.1446", "PR1-like-1", #0
  998. "NV2.7811", "OtxC", "ADA1A-like-16", "NV2.14260", "ACH10-like-9", #1
  999. "PAPSS-like-1", "NV2.14774", "THIO-like-1", "PRGa", #"HLSa/SNIa/NTVa", #2
  1000. "myc1", "Myc3", "NV2.2299", "YPF17-like-2", "ADRB2-like-31", #3
  1001. "AA2BR-like-2", "NV2.23614", "MRP4-like-1", "MUC1-like-2", "NV2.3856", #4
  1002. "FoxA", "VKT52-like-1", "PPR27-like-12", "NV2.11442", "RSPH1-like-4", #5
  1003. "BLML2-like-1", "NV2.14153", "AA3R-like-7", "TAAR6-like-2", "NV2.21248", #6
  1004. "NV2.6813", "MLRP2-like-30", "NV2.8789", "NV2.17754", "GR101-like-11", #7
  1005. "VIAAT-like-40", "fez", "CNTP1-like-14", "TMPS6-like-2", "SC6A1-like-9", #8
  1006. "NV2.13863", "NV2.2898", "RL12-like-1", "RL7-like-1", "RL13-like-1", #9
  1007. "RPAC1-like-1", "FOXD1-like-1", "ATS6-like-7", "DD3-like-2", "NV2.23233", #10
  1008. "NV2.160", "NV2.5192", "SEGN-like-4", "BIR7B-like-1", "MCTP1-like-1", #11
  1009. "BHA15-like-1", "GRIA3-like-1", "NV2.8266", "NR1BA-like-1", "VIAAT-like-31", #12
  1010. "NV2.2620", "NV2.2627", "NV2.15930", "NV2.2623", "NV2.9018", #13
  1011. "HCE2-like-1", "ACH10-like-11", "QRFPR-like-31", "ADRB1-like-1", "TLR1-like-4", #14
  1012. "NV2.450", "NV2.5143", "NV2.5142", "DMBT1-like-29", "NV2.5046", #15
  1013. "NV2.1204", "NV2.1205", "OPN4B-like-6", "ADRB2-like-40", "AGRL3-like-1", #16
  1014. "NV2.10279", "GBRB1-like-5", "NV2.14290", "LOX5-like-3", "CNTP5-like-7", #17
  1015. "NV2.21058", "NV2.9262", "AA2AR-like-3", "NV2.14663", "OPN4-like-13" #18
  1016. )
  1018. plt9 <- DotPlot(neurogland.all, features = genes) + RotatedAxis() + theme(axis.text.x = element_text(angle = 90)) + FontSize(10)
  1019. plt9 <- plt9 + ggtitle("ashA neural cluster markers on neurogland")
  1020. #plt9 <- plt9 + scale_x_discrete(labels = expression(
  1021. # italic("")))
  1022. plt9
  1023. ggsave(plt9, file = "/ashAneuralmarkersNGdata.jpg", width = 15, height = 10)
  1025. ```
  1027. ```{r}
  1028. dim(ashA_neural_sub3_15dims)
  1029. ncol(ashA_neural_sub3_15dims)
  1030. ```
  1033. #trajectory and pseudotime
  1034. trajectory analysis with node selection
  1035. ```{r}
  1036. library(Seurat)
  1037. library(monocle3)
  1038. library(SeuratWrappers)
  1039. library(ggplot2)
  1040. library(dplyr)
  1041. library(viridis)
  1043. #----------------------------------------------------------
  1044. # 1) Start from your Seurat object
  1045. #----------------------------------------------------------
  1046. seu <- ashA_neural_sub3_15dims
  1048. #----------------------------------------------------------
  1049. # 2) Convert Seurat -> Monocle3 CDS
  1050. #----------------------------------------------------------
  1051. cds <- SeuratWrappers::as.cell_data_set(seu)
  1053. #----------------------------------------------------------
  1054. # 3) Cluster cells & define partitions (Monocle3)
  1055. #----------------------------------------------------------
  1056. cds <- monocle3::cluster_cells(
  1057. cds = cds,
  1058. reduction_method = "UMAP",
  1059. resolution = 0.001 # ~Seurat res 0.1
  1060. )
  1062. # (optional sanity checks)
  1063. monocle3::plot_cells(cds) # clusters
  1064. monocle3::plot_cells(cds, color_cells_by = "partition") # partitions
  1066. #----------------------------------------------------------
  1067. # 4) Learn the trajectory graph using those partitions
  1068. #----------------------------------------------------------
  1069. cds <- monocle3::learn_graph(
  1070. cds,
  1071. use_partition = TRUE,
  1072. close_loop = TRUE
  1073. )
  1075. #----------------------------------------------------------
  1076. # 5) Choose roots interactively & compute pseudotime
  1077. # - A plot appears: click the node(s) where you want
  1078. #----------------------------------------------------------
  1079. cds <- monocle3::order_cells(cds)
  1081. #----------------------------------------------------------
  1082. # 6) Extract pseudotime and CLEAN it (no Inf)
  1083. #----------------------------------------------------------
  1084. pt <- monocle3::pseudotime(cds)
  1086. # replace Inf / -Inf / NaN with NA so Seurat can plot it
  1087. pt_clean <- ifelse(is.finite(pt), pt, NA_real_)
  1089. # put cleaned pseudotime into Seurat
  1090. seu$cds_pseudotime <- pt_clean
  1092. # quick check
  1093. summary(seu$cds_pseudotime)
  1095. #----------------------------------------------------------
  1096. # 7) Summarize which clusters are early/late (optional)
  1097. #----------------------------------------------------------
  1098. meta <- as.data.frame(colData(cds))
  1099. meta$partition_m3 <- monocle3::partitions(cds)
  1100. meta$cluster_m3 <- monocle3::clusters(cds)
  1101. meta$pt <- pt_clean
  1103. meta_summary <- meta %>%
  1104. group_by(partition = partition_m3, cluster = cluster_m3) %>%
  1105. summarise(
  1106. n_cells = n(),
  1107. min_pt = min(pt, na.rm = TRUE),
  1108. mean_pt = mean(pt, na.rm = TRUE),
  1109. .groups = "drop"
  1110. ) %>%
  1111. arrange(partition, min_pt)
  1113. meta_summary # <-- table tells you what’s earliest in each partition
  1115. #----------------------------------------------------------
  1116. # 8) Pseudotime-only monocle plot (blue-yellow viridis) + mirror
  1117. #----------------------------------------------------------
  1118. # helper to flip left-right
  1119. mirror_x <- function(p) p + scale_x_reverse()
  1121. p_pseudo <- monocle3::plot_cells(
  1122. cds,
  1123. color_cells_by = "pseudotime",
  1124. show_trajectory_graph = FALSE, # no lines, just cells
  1125. label_groups_by_cluster = FALSE,
  1126. graph_label_size = 0,
  1127. cell_size = 0.5
  1128. )
  1130. p_pseudo_colored <- p_pseudo +
  1131. scale_color_viridis_c(option = "C") + # blue->yellow palette
  1132. theme_classic()
  1134. p_pseudo_mirrored <- mirror_x(p_pseudo_colored)
  1136. p_pseudo_mirrored # <-- this should be “pseudotime only” mirrored panel
  1138. #----------------------------------------------------------
  1139. # 9) FeaturePlot pseudotime on Seurat UMAP
  1140. #----------------------------------------------------------
  1141. FeaturePlot(
  1142. seu,
  1143. features = "cds_pseudotime",
  1144. reduction = "umap",
  1145. pt.size = 0.1
  1146. )
  1148. #----------------------------------------------------------
  1149. # 10) Save rooted CDS
  1150. #----------------------------------------------------------
  1151. saveRDS(cds, file = "ashA_neural_sub3_monocle_rooted.rds")
  1153. ```
  1155. ```{r}
  1156. # Trajectory + partitions, not mirrored
  1157. p_part_traj <- monocle3::plot_cells(
  1158. cds,
  1159. color_cells_by = "partition",
  1160. show_trajectory_graph = TRUE,
  1161. label_groups_by_cluster = FALSE, # no cluster labels
  1162. label_roots = TRUE,
  1163. label_leaves = TRUE,
  1164. label_branch_points = TRUE
  1165. )
  1167. p_part_traj
  1169. # Mirrored version (to match other panels)
  1170. p_part_traj_mirrored <- mirror_x(p_part_traj)
  1171. p_part_traj_mirrored
  1173. ```
  1176. ```{r}
  1177. p_part_traj_clean <- monocle3::plot_cells(
  1178. cds,
  1179. color_cells_by = "partition",
  1180. show_trajectory_graph = TRUE,
  1181. label_groups_by_cluster = FALSE,
  1182. label_roots = FALSE,
  1183. label_leaves = FALSE,
  1184. label_branch_points = FALSE
  1185. )
  1187. p_part_traj_clean_mirrored <- mirror_x(p_part_traj_clean)
  1188. p_part_traj_clean_mirrored
  1190. ```
  1192. ```{r}
  1193. # pseudotime-only (new rooting)
  1194. pC_new <- monocle3::plot_cells(
  1195. cds,
  1196. color_cells_by = "pseudotime",
  1197. show_trajectory_graph = FALSE,
  1198. label_groups_by_cluster = FALSE,
  1199. graph_label_size = 0,
  1200. cell_size = 0.5
  1201. ) +
  1202. scale_x_reverse() # mirrored to match the others
  1204. pC_new
  1206. # pseudotime + trajectory (new rooting)
  1207. pD_new <- monocle3::plot_cells(
  1208. cds,
  1209. color_cells_by = "pseudotime",
  1210. show_trajectory_graph = TRUE,
  1211. label_roots = TRUE,
  1212. label_leaves = TRUE,
  1213. label_branch_points = TRUE,
  1214. label_groups_by_cluster = FALSE,
  1215. cell_size = 0.5
  1216. ) +
  1217. scale_x_reverse()
  1219. pD_new
  1221. ```
  1223. ```{r}
  1224. p_traj_clean <- monocle3::plot_cells(
  1225. cds,
  1226. color_cells_by = "pseudotime",
  1227. show_trajectory_graph = TRUE,
  1228. label_roots = FALSE,
  1229. label_leaves = FALSE,
  1230. label_branch_points = FALSE,
  1231. label_groups_by_cluster = FALSE,
  1232. cell_size = 0.5
  1233. )
  1235. p_traj_clean
  1237. p_traj_clean_mirrored <- p_traj_clean + scale_x_reverse()
  1238. p_traj_clean_mirrored
  1240. ```

05_NvashA_subsetting_trajectory.Rmd at commit 07b4710, no license · at the source

Overview

Authors: Jamie A Havrilak1,2, MingHe Cheng1,2, Layla Al-Shaer1,2, Whitney B Leach1, Mia Yagodich1, Dylan Faltine-Gonzalez1,3, Michael J Layden1,2
  1. Department of Biological Sciences, Lehigh University, Bethlehem, PA USA
  2. Lehigh Oceans Research Center, Lehigh University, Bethlehem, PA USA
  3. Present Address: Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD USA
Institutions: Lehigh University (United States); Johns Hopkins University (United States)
Journal: Scientific reports, volume 16, issue 1, article 12151
Dates: received 8 August 2025; accepted 20 February 2026; published online 5 March 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41598-026-41460-z · PMID 41781584 · PMCID PMC13076895 · OpenAlex W7133754179
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: other (organism), developmental (subfield)
Keywords: Neurogenesis, Nematostella, NvashA, Neuronal patterning, Nerve-net, Neural development, Developmental biology, Evolution, Neuroscience
MeSH: Cnidaria*, Neurogenesis*, Neurons*, Sea Anemones*, Animals, Body Patterning, Gene Expression Regulation, Developmental, Larva, Nervous System (* major topic)
Topic: Marine Invertebrate Physiology and Ecology (Paleontology, Earth and Planetary Sciences), according to OpenAlex
Funding: National Science Foundation (1942777); NIGMS NIH HHS (R01 GM127615); Department of Health and Human Services National Institute of General Medical Sciences (R01GM127615)
Citations: cited by 2 papers (Europe PMC); 59 references in the paper
Research resources: RRID:SCR_024675

Abstract

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

Repository

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

LaydenLab/NvashA_scRNAseq

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 07b47105178678bfa348e9c3474025a3246234ec, 5 December 2025
Languages: R (5)
Size: 6 files, 5 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 5 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (5 files), patchwork (5 files), Seurat (5 files), tidyverse (5 files), igraph (1 file), Monocle 3 (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
6 files
At the source: github.com/LaydenLab/NvashA_scRNAseq

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;
  • 5 scripts, each with its path and the digest of its content;
  • 5 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

Code and data availability statement

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

Read it in the paper: doi.org/10.1038/s41598-026-41460-z.

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, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 9 keywords, 9 MeSH terms, 3 funders, 49 references, 1 RRID.

Cite

This paper

Havrilak, J. A., Cheng, M., Al-Shaer, L., Leach, W. B., Yagodich, M., Faltine-Gonzalez, D., & Layden, M. J. (2026). NvashA function reveals temporal differences in neural subtype generation in cnidarians. Scientific reports, 16(1), 12151. https://doi.org/10.1038/s41598-026-41460-z

BibTeX

@article{havrilak2026nvasha,
author = {Havrilak, Jamie A and Cheng, MingHe and Al-Shaer, Layla and Leach, Whitney B and Yagodich, Mia and Faltine-Gonzalez, Dylan and Layden, Michael J},
title = {{NvashA function reveals temporal differences in neural subtype generation in cnidarians}},
journal = {Scientific reports},
year = {2026},
month = mar,
volume = {16},
number = {1},
pages = {12151},
publisher = {Nature Publishing Group},
issn = {2045-2322},
doi = {10.1038/s41598-026-41460-z},
url = {https://doi.org/10.1038/s41598-026-41460-z},
pmid = {41781584},
pmcid = {PMC13076895}
}

RIS

TY - JOUR
AU - Havrilak, Jamie A
AU - Cheng, MingHe
AU - Al-Shaer, Layla
AU - Leach, Whitney B
AU - Yagodich, Mia
AU - Faltine-Gonzalez, Dylan
AU - Layden, Michael J
TI - NvashA function reveals temporal differences in neural subtype generation in cnidarians
T2 - Scientific reports
J2 - Sci Rep
PY - 2026
DA - 2026/03/05
VL - 16
IS - 1
SP - 12151
SN - 2045-2322
PB - Nature Publishing Group
DO - 10.1038/s41598-026-41460-z
UR - https://doi.org/10.1038/s41598-026-41460-z
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41598-026-41460-z",
"type": "article-journal",
"title": "NvashA function reveals temporal differences in neural subtype generation in cnidarians",
"container-title": "Scientific reports",
"author": [
{
"family": "Havrilak",
"given": "Jamie A"
},
{
"family": "Cheng",
"given": "MingHe"
},
{
"family": "Al-Shaer",
"given": "Layla"
},
{
"family": "Leach",
"given": "Whitney B"
},
{
"family": "Yagodich",
"given": "Mia"
},
{
"family": "Faltine-Gonzalez",
"given": "Dylan"
},
{
"family": "Layden",
"given": "Michael J"
}
],
"container-title-short": "Sci Rep",
"volume": "16",
"issue": "1",
"page": "12151",
"DOI": "10.1038/s41598-026-41460-z",
"PMID": "41781584",
"PMCID": "PMC13076895",
"ISSN": "2045-2322",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41598-026-41460-z",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
5
]
]
}
}

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.1371/journal.pbio.3003803 [code]
The anti-neural role of BMP signaling is a consequence of its ancestral function in dorsoventral patterning.
Journal: PLoS biology
In common: Seurat, ggplot2, tidyverse, other, 12 references
[2] doi:10.1371/journal.pbio.3003931 [code]
A single-cell transcriptomic atlas reveals the emergence of medusa-specific cell states in the scyphozoan Aurelia coerulea.
Journal: PLoS biology
In common: Seurat, patchwork, ggplot2, 1 other tool, 7 references
[3] 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: Monocle 3, igraph, Seurat, 3 other tools, 2 references
[4] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Monocle 3, igraph, Seurat, 3 other tools, 2 references
[5] doi:10.1016/j.stem.2026.05.005 [code]
Generation of human appetite-regulating neurons and tanycytes from pluripotent stem cells.
Journal: Cell stem cell
In common: Monocle 3, igraph, Seurat, 3 other tools, 1 reference
[6] doi:10.3390/ijms27104466 [code]
Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: Monocle 3, igraph, Seurat, 3 other tools, 1 reference
[7] doi:10.1038/s41467-026-69944-6 [code]
Multi-modal dissection of cell-type specific TDP-43 pathology in the motor cortex.
Journal: Nature communications
In common: Monocle 3, igraph, Seurat, 3 other tools, 1 reference
[8] doi:10.1038/s41420-026-02971-w [code]
Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells.
Journal: Cell death discovery
In common: Monocle 3, igraph, Seurat, 3 other tools, 1 reference
[9] doi:10.1038/s41467-026-75273-5 [code]
Spinal cord regeneration deploys cell-type specific developmental and non-developmental strategies to restore neuron diversity.
Journal: Nature communications
In common: Monocle 3, Seurat, patchwork, 2 other tools, other, developmental, 1 reference
[10] doi:10.1002/jsp2.70200 [code]
A Porcine Model of Intervertebral Disc Injury Recapitulates Human Discogenic Pain Via Notochordal Cell Loss and Pain-Inducing Nucleus Pulposus Cell Emergence.
Journal: JOR spine
In common: Monocle 3, igraph, Seurat, 3 other tools, other

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.