OSCR

MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos.

Code ↔ Paper

20 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 20 matches
  1. [1] § Results › MIC-Drop-seq enables multiplexed embryo phenotyping ↔ Figure1_MIC-Drop-seq.Rmd, lines 907–958 · score 0.81 · hoxa9a, hoxb9a, hoxc6b, spinal cord, identified cdx4, clusters
  2. [2] § Methods › Dimensionality reduction, clustering and cell cluster labeling ↔ Figure1_MIC-Drop-seq.Rmd, lines 551–607 · score 0.81 · transfer anchors, 21–26, Predicted cell, reference atlas, Daniocell, hpf
  3. [3] § Methods › Determination of CRISPR editing frequency ↔ Figure1_MIC-Drop-seq.Rmd, lines 442–484 · score 0.80 · rhAmpSeq, cut site, editing frequency, Genomic DNA, indels, sequencing
  4. [4] § Results › MIC-Drop-seq reveals transcriptional and cellular phenotypes ↔ Figure3_MIC-Drop-seq.Rmd, lines 473–498 · score 0.78 · optic tectum, nkx2.7, gata2a, sox11b, dlx1a expression, hypothalamus
  5. [5] § Methods › Differential cell abundance analysis ↔ Figure2_MIC-Drop-seq.Rmd, lines 951–1019 · score 0.76 · Poisson Lognormal models, biological replicate, cell abundance, Monocle3, Hooke, genotype
  6. [6] § Methods › Statistics and reproducibility ↔ Figure4_MIC-Drop-seq.Rmd, lines 912–959 · score 0.75 · pairwise comparisons, post hoc, mRNA, Tukey, rescue, ANOVA
  7. [7] § Results › MIC-Drop-seq reveals transcriptional and cellular phenotypes ↔ Figure3_MIC-Drop-seq.Rmd, lines 473–498 · score 0.74 · optic tectum, nkx2.7, gata2a, sox11b, dlx1a expression, hypothalamus
  8. [8] § Methods › In situ imaging and quantification ↔ Figure3_MIC-Drop-seq.Rmd, lines 285–312 · score 0.71 · pixel area, dlx1a, gata2a, Stained, trunk, head
  9. [9] § Results › MIC-Drop-seq enables multiplexed embryo phenotyping ↔ Figure1_MIC-Drop-seq.Rmd, lines 442–484 · score 0.71 · cut sites, genomic DNA, editing frequency, genotypic composition, bulk, indel
  10. [10] § Methods › Dimensionality reduction, clustering and cell cluster labeling ↔ Figure1_MIC-Drop-seq.Rmd, lines 529–549 · score 0.71 · variable features, remove cells, dimensionality reduction, Seurat, doublets, gRNA
  11. [11] § Methods › Differential gene expression analysis ↔ Figure2_MIC-Drop-seq.Rmd, lines 664–749 · score 0.70 · edgeR, glmQLFit, aggregate, pseudobulk, FDR, gene expression
  12. [12] § Results › MIC-Drop-seq enables multiplexed embryo phenotyping ↔ Figure1_MIC-Drop-seq.Rmd, lines 698–727 · score 0.65 · presomitic mesoderm, tbx16 mutant, tbxta mutant, tyr, accumulation, depletion
  13. [13] § Results › MIC-Drop-seq detects indirect developmental effects ↔ Figure4_MIC-Drop-seq.Rmd, lines 912–959 · score 0.64 · post hoc, mRNA, way ANOVA, CVP, Tukey, outliers
  14. [14] § Methods › Gene expression matrix generation and genotype assignment ↔ Figure2_MIC-Drop-seq.Rmd, lines 159–241 · score 0.62 · gRNAs, multiple gene, mutant genotype, protospacer, doublet, detection
  15. [15] § Results › MIC-Drop-seq reveals transcriptional and cellular phenotypes ↔ Figure4_MIC-Drop-seq.Rmd, lines 169–292 · score 0.61 · presomitic mesoderm, optic tectum, myotome, sclerotome, dorsal, secreted
  16. [16] § Results › MIC-Drop-seq reveals transcriptional and cellular phenotypes ↔ Figure4_MIC-Drop-seq.Rmd, lines 169–292 · score 0.60 · optic tectum, hindbrain, neurons, ventral, MIC Drop seq, precursors
  17. [17] § Results › MIC-Drop-seq detects indirect developmental effects ↔ Figure4_MIC-Drop-seq.Rmd, lines 839–888 · score 0.59 · phenotype class, DEG thresholds, cell extrinsic, lineage, classification, seq
  18. [18] § Results › MIC-Drop-seq detects indirect developmental effects ↔ Figure4_MIC-Drop-seq.Rmd, lines 839–888 · score 0.58 · lineage intrinsic, cell intrinsic, cell extrinsic, threshold, DEGs, classified
  19. [19] § Methods › Determination of cell extrinsic and cell intrinsic effects ↔ Figure4_MIC-Drop-seq.Rmd, lines 635–666 · score 0.52 · cell intrinsic, cell extrinsic, threshold, lineage, classified
  20. [20] § Results › MIC-Drop-seq enables multiplexed embryo phenotyping ↔ Figure1_MIC-Drop-seq.Rmd, lines 698–727 · score 0.52 · hoxb1b, foxa2, hand2, tbx16, tyr, tbxta

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,137 lines · 37 KB · MIT · 7 matches

  1. ---
  2. title: "Figure 1: MIC-Drop-seq Proof-of-Concept Validation"
  3. author: "Clay Carey"
  4. date: "2024-05-08"
  5. output: html_document
  6. ---
  9. First- set your base directory to the folder with the downloaded 'input_data" folder. This contains all source data to run the code to reprodce figures.
  12. ```{r setup}
  13. # Set working directory for analysis- you MUST download the input dat aand set this directory manually for the code to work
  14. knitr::opts_knit$set(root.dir = "")
  15. ```
  18. # =============================================================================
  19. # SECTION 1: LOAD LIBRARIES
  20. # =============================================================================
  22. ```{r load-libraries, message=FALSE, warning=FALSE}
  23. # Core data manipulation
  24. library(Seurat)
  25. library(dplyr)
  26. library(tidyverse)
  27. library(stringr)
  29. # Single-cell analysis
  30. library(scCustomize)
  31. library(dittoSeq)
  32. library(scDblFinder)
  33. library(SingleCellExperiment)
  35. # Visualization
  36. library(ggplot2)
  37. library(ggridges)
  38. library(viridis)
  39. library(EnhancedVolcano)
  41. # Differential expression
  42. library(edgeR)
  43. library(presto)
  45. # Table formatting
  46. library(gtExtras)
  47. ```
  49. # =============================================================================
  50. # SECTION 2: LOAD AND PROCESS gRNA DATA
  51. # =============================================================================
  53. ## Load Cell Ranger CRISPR output
  55. Cell Ranger generates a "protospacer_calls_per_cell.csv" file that contains:
  56. - cell_barcode: 10X cell barcode
  57. - feature_call: pipe-delimited list of detected gRNA IDs
  58. - num_umis: pipe-delimited list of UMI counts per gRNA
  60. ```{r load-grna-data}
  61. # Load gRNA detection data from Cell Ranger
  62. pilot_CRISPR_data <- read.csv(file = "input_data/protospacer_calls_per_cell.csv")
  64. # Load experimental configuration (maps gRNA IDs to target genes)
  65. exp_config <- read.csv(file = "input_data/feature_reference.csv")
  67. head(pilot_CRISPR_data)
  68. head(exp_config)
  69. ```
  71. ## Parse gRNA detection data
  73. Cell Ranger stores multiple gRNAs per cell as pipe-delimited strings.
  74. We split these into separate columns for downstream classification.
  76. ```{r parse-grna-function}
  77. #' Process Cell Ranger CRISPR output
  78. #'
  79. #' Splits pipe-delimited feature_call and num_umis columns into
  80. #' separate columns (up to 20 features per cell)
  81. #'
  82. #' @param data Dataframe from Cell Ranger protospacer_calls_per_cell.csv
  83. #' @return Dataframe with feature_call_1:20 and num_umis_1:20 columns
  84. process_data <- function(data) {
  86. max_features <- 20
  87. feature_cols <- paste0("feature_call_", 1:max_features)
  89. # Split feature_call column (gRNA IDs)
  90. split_features <- strsplit(as.character(data$feature_call), "\\|")
  92. for (i in 1:max_features) {
  93. data[, feature_cols[i]] <- sapply(split_features, function(x) {
  94. if (length(x) >= i) return(x[i]) else return(0)
  95. })
  96. }
  98. # Split num_umis column (UMI counts)
  99. max_umis <- 20
  100. umis_cols <- paste0("num_umis_", 1:max_umis)
  101. split_umis <- strsplit(as.character(data$num_umis), "\\|")
  103. for (i in 1:max_umis) {
  104. data[, umis_cols[i]] <- sapply(split_umis, function(x) {
  105. if (length(x) >= i) return(x[i]) else return(0)
  106. })
  107. }
  109. return(data)
  110. }
  112. # Apply parsing function
  113. CRdata <- process_data(pilot_CRISPR_data)
  114. ```
  116. ## Classify cells by genotype
  118. For each cell, determine if gRNAs target a single gene or multiple genes.
  119. - Single gene → Assign that gene as genotype
  120. - Multiple genes → Classify as "Multiple" (doublet)
  121. - No gRNA → Will be classified as "ND" (Not Detected) later
  123. ```{r classify-genotypes}
  124. #' Classify cells by gRNA detection pattern
  125. #'
  126. #' Checks each cell across all feature_call columns and determines
  127. #' whether detected gRNAs target a single gene or multiple genes
  128. #'
  129. #' @param df Parsed CRISPR dataframe with feature_call_1:20 columns
  130. #' @param proto_ids Vector of gRNA IDs from experimental config
  131. #' @return Dataframe with 'class' column (gene name or "Multiple")
  132. classifier <- function(df, proto_ids) {
  134. # Create binary detection columns for each gRNA
  135. for (proto in proto_ids) {
  136. df <- df %>%
  137. mutate(!!paste0(proto, "_count") := case_when(
  138. feature_call_1 == proto ~ 1,
  139. feature_call_2 == proto ~ 1,
  140. feature_call_3 == proto ~ 1,
  141. feature_call_4 == proto ~ 1,
  142. feature_call_5 == proto ~ 1,
  143. feature_call_6 == proto ~ 1,
  144. feature_call_7 == proto ~ 1,
  145. feature_call_8 == proto ~ 1,
  146. feature_call_9 == proto ~ 1,
  147. feature_call_10 == proto ~ 1,
  148. feature_call_11 == proto ~ 1,
  149. feature_call_12 == proto ~ 1,
  150. feature_call_13 == proto ~ 1,
  151. feature_call_14 == proto ~ 1,
  152. feature_call_15 == proto ~ 1,
  153. feature_call_16 == proto ~ 1,
  154. feature_call_17 == proto ~ 1,
  155. feature_call_18 == proto ~ 1,
  156. feature_call_19 == proto ~ 1,
  157. feature_call_20 == proto ~ 1,
  158. TRUE ~ 0
  159. ))
  160. }
  162. # Extract detection counts and identify detected gRNAs per cell
  163. temp_df <- select(df, cell_barcode, ends_with("_count"))
  164. temp_df <- temp_df %>%
  165. pivot_longer(-cell_barcode) %>%
  166. filter(value == 1)
  168. # Summarize detected gRNAs per cell
  169. temp_df <- temp_df %>%
  170. group_by(cell_barcode) %>%
  171. summarise(class = paste(unique(name), collapse = ","))
  173. # Clean up gRNA names to gene names
  174. temp_df$all_protos <- temp_df$class
  175. temp_df$all_protos <- gsub("_count", "", temp_df$all_protos)
  176. temp_df$class <- gsub("-\\d+_count", "", temp_df$class)
  178. # Separate into columns to check for multiple targets
  179. temp_df <- separate(temp_df,
  180. class,
  181. into = c("count1", "count2", "count3", "count4", "count5"),
  182. sep = ",",
  183. extra = "merge")
  185. # Classify based on target gene consistency
  186. temp_df <- mutate(temp_df, class = case_when(
  187. is.na(count2) ~ count1,
  188. count1 == count2 & is.na(count3) ~ count1,
  189. count1 == count2 & count2 == count3 & is.na(count4) ~ count1,
  190. count1 == count2 & count2 == count3 & count3 == count4 & is.na(count5) ~ count1,
  191. TRUE ~ "Multiple" # Different genes detected = doublet
  192. ))
  194. # Merge back with original data
  195. temp_df <- select(temp_df, cell_barcode, all_protos, class)
  196. df <- left_join(df, temp_df, by = "cell_barcode")
  197. df <- select(df, cell_barcode, num_features, feature_call, num_umis, all_protos, class)
  199. return(df)
  200. }
  202. # Apply classification
  203. CRdata_classified <- classifier(CRdata, unique(exp_config$id))
  204. head(CRdata_classified)
  205. ```
  207. # =============================================================================
  208. # SECTION 3: CREATE SEURAT OBJECT AND ADD METADATA
  209. # =============================================================================
  211. ## Load 10X data
  213. Data contains two assays:
  214. 1. Gene Expression (RNA)
  215. 2. CRISPR Guide Capture (gRNA counts)
  217. ```{r load-10x-data}
  218. # Read both assays from Cell Ranger output
  219. tenx_data <- Read10X(data.dir = "input_data/micdrop_overload_matrix")
  220. tenx_data.gex <- tenx_data$`Gene Expression`
  221. tenx_data.crispr <- tenx_data$`CRISPR Guide Capture`
  223. # Create Seurat object with RNA assay
  224. mdp <- CreateSeuratObject(counts = tenx_data.gex, min.cells = 3)
  226. # Add CRISPR assay
  227. mdp[['CRISPR']] <- CreateAssayObject(counts = tenx_data.crispr)
  228. ```
  230. ## Merge gRNA classifications into Seurat metadata
  232. ```{r add-grna-metadata}
  233. # Extract metadata and add cell barcodes
  234. mdp_meta <- [email hidden]
  235. mdp_meta$cell_barcode <- rownames(mdp_meta)
  237. # Merge with gRNA classifications
  238. mdp_meta <- left_join(mdp_meta, CRdata_classified, by = "cell_barcode")
  239. rownames(mdp_meta) <- mdp_meta$cell_barcode
  241. # Add back to Seurat object
  242. [email hidden] <- mdp_meta
  244. # Clean up classifications
  245. [email hidden] <- [email hidden] %>%
  246. mutate(
  247. class = ifelse(is.na(class), "ND", class), # Cells without gRNA
  248. class = ifelse(class == "Non-Targeting", "tyr", class) # Rename control
  249. )
  251. table(mdp$class)
  252. ```
  254. # =============================================================================
  255. # SECTION 4: QUALITY CONTROL
  256. # =============================================================================
  258. ## Calculate QC metrics
  260. ```{r calculate-qc}
  261. # Calculate mitochondrial percentage
  262. mdp[["percent.mt"]] <- PercentageFeatureSet(mdp, pattern = "^mt-")
  264. # Visualize QC metrics
  265. DefaultAssay(mdp) <- "RNA"
  266. VlnPlot(mdp, features = c("percent.mt", "nCount_RNA", "nFeature_RNA"))
  267. ```
  269. ## Filter low-quality cells
  271. ```{r filter-cells}
  272. # Apply quality thresholds
  273. mdp <- subset(mdp, subset = nFeature_RNA > 200 &
  274. nFeature_RNA < 9000 &
  275. percent.mt < 15)
  277. cat("Cells after filtering:", ncol(mdp), "\n")
  278. ```
  280. # =============================================================================
  281. # SECTION 5: FIGURE 1B - gRNA UMI DISTRIBUTION
  282. # =============================================================================
  284. Distribution of gRNA UMI counts in cells with detected gRNAs
  286. ```{r figure1B, fig.width=2, fig.height=2}
  287. mdp_meta <- [email hidden]
  289. mdp_meta %>%
  290. filter(class != "ND") %>%
  291. ggplot(aes(x = nCount_CRISPR)) +
  292. geom_histogram(binwidth = 1, fill = "black") +
  293. scale_x_continuous(limits = c(0, 75)) +
  294. theme_classic() +
  295. theme(
  296. text = element_text(family = "Helvetica", size = 12, color = "black"),
  297. axis.text = element_text(size = 12, color = "black", family = "Helvetica"),
  298. axis.title = element_text(size = 12)
  299. ) +
  300. xlab("gRNA UMIs") +
  301. ylab("Cell Count")
  303. ggsave(filename = "outputs_fig1/Figure_1B.eps",
  304. device = "ps",
  305. width = 2, height = 2, units = 'in',
  306. bg = "transparent")
  307. ```
  309. # =============================================================================
  310. # SECTION 6: FIGURE 1C - gRNA TARGET DISTRIBUTION
  311. # =============================================================================
  313. Cells classified by number of gene targets detected
  314. - 0: No gRNA detected
  315. - 1: Single gene target (expected for MIC-Drop)
  316. - 2+: Multiple gene targets (doublets)
  318. ```{r figure1C, fig.width=1.5, fig.height=2}
  319. # Categorize cells by number of targets
  320. mdp_meta <- [email hidden] %>%
  321. mutate(type = case_when(
  322. class == "ND" ~ "0",
  323. class == "Multiple" ~ "2+",
  324. TRUE ~ "1"
  325. ))
  327. # Plot distribution
  328. mdp_meta %>%
  329. dplyr::count(type) %>%
  330. ggplot() +
  331. geom_bar(aes(x = type, y = n / 1000), stat = "identity", fill = "black") +
  332. theme_classic() +
  333. xlab("gRNA targets") +
  334. ylab("Cells (thousands)") +
  335. scale_y_continuous(labels = scales::comma) +
  336. theme(
  337. axis.text = element_text(size = 12, family = "Helvetica", color = "black"),
  338. axis.title = element_text(size = 12, family = "Helvetica")
  339. )
  341. ggsave(filename = "outputs_fig1/Figure_1C.eps",
  342. device = "ps",
  343. width = 1.5, height = 2, units = 'in',
  344. bg = "transparent")
  345. ```
  347. # =============================================================================
  348. # SECTION 7: DOUBLET DETECTION WITH scDblFinder
  349. # =============================================================================
  351. Use computational doublet detection to validate that cells with
  352. multiple gRNA targets are enriched for doublets
  354. ```{r doublet-detection}
  355. # Prepare for scDblFinder (requires matching order)
  356. correctorder <- dimnames(mdp@assays$RNA$counts)[[2]]
  357. mdp <- mdp[, correctorder]
  358. [email hidden] <- [email hidden][match(correctorder, rownames([email hidden])), ]
  360. # Normalize and scale for doublet detection
  361. mdp <- NormalizeData(mdp)
  362. mdp <- ScaleData(mdp)
  364. # Convert to SingleCellExperiment and run scDblFinder
  365. mdp_sce <- as.SingleCellExperiment(mdp)
  366. mdp_sce <- scDblFinder(mdp_sce)
  367. mdp_sce <- as.Seurat(mdp_sce)
  369. # Extract doublet information
  370. dub_info <- data.frame([email hidden]) %>%
  371. select(scDblFinder.class, scDblFinder.score,
  372. scDblFinder.weighted, scDblFinder.cxds_score, class) %>%
  373. mutate(type = case_when(
  374. class == "Multiple" ~ "2+",
  375. TRUE ~ "0-1"
  376. ))
  377. ```
  379. ## Figure S2A: Doublet score distribution
  381. Cells with multiple gRNA targets show higher doublet scores
  383. ```{r figureS2A, fig.width=4, fig.height=3}
  384. ggplot(dub_info, aes(x = scDblFinder.weighted, fill = type)) +
  385. geom_density(alpha = 0.8, position = "identity", color = NA) +
  386. scale_fill_manual(values = c("0-1" = "black", "2+" = "grey")) +
  387. xlab("scDblFinder Doublet Score") +
  388. ylab("Density") +
  389. theme_classic() +
  390. guides(fill = guide_legend(title = "gRNA targets")) +
  391. theme(
  392. axis.text = element_text(size = 12, family = "Helvetica", color = "black"),
  393. axis.title = element_text(size = 12, family = "Helvetica"),
  394. legend.text = element_text(size = 12, family = "Helvetica")
  395. )
  397. ggsave(filename = "outputs_fig1/Figure_S2A.eps",
  398. device = "ps",
  399. width = 4, height = 3, units = 'in',
  400. bg = "transparent")
  401. ```
  403. ## Figure S2B: Doublet classification by gRNA pattern
  405. ```{r figureS2B, fig.width=3.5, fig.height=3}
  406. # Calculate doublet proportions
  407. class_freq <- dub_info %>%
  408. select(type, scDblFinder.class) %>%
  409. group_by(type) %>%
  410. dplyr::count(scDblFinder.class) %>%
  411. left_join(
  412. dub_info %>%
  413. select(type, scDblFinder.class) %>%
  414. group_by(type) %>%
  415. dplyr::count(scDblFinder.class) %>%
  416. summarize(tot = sum(n)),
  417. by = "type"
  418. ) %>%
  419. mutate(pct = n / tot)
  421. # Plot proportions
  422. ggplot(class_freq, aes(x = type, fill = scDblFinder.class, y = pct)) +
  423. scale_fill_manual(values = c("singlet" = "black", "doublet" = "grey")) +
  424. geom_bar(position = "stack", stat = "identity") +
  425. theme_classic() +
  426. theme(
  427. axis.text = element_text(size = 12, family = "Helvetica", color = "black"),
  428. axis.title = element_text(size = 12, family = "Helvetica"),
  429. legend.text = element_text(size = 12, family = "Helvetica")
  430. ) +
  431. guides(fill = guide_legend(title = "scDblFinder Class")) +
  432. ylab("Proportion of Cells") +
  433. xlab("gRNA targets recovered") +
  434. labs(fill = "")
  436. ggsave(filename = "outputs_fig1/Figure_S2B.eps",
  437. device = "ps",
  438. width = 3.5, height = 3, units = 'in',
  439. bg = "transparent")
  440. ```
  442. # =============================================================================
  443. # SECTION 8: FIGURE 1E - MUTAGENESIS EFFICIENCY
  444. # =============================================================================
  446. Compare CRISPR editing frequency (from genomic DNA) to
  447. genotype composition (from gRNA detection)
  449. ```{r load-rhampseq-data}
  450. # Load rhAmpSeq results (bulk DNA sequencing of cut sites)
  451. rhamp <- read.csv(file = "input_data/rhamp_results.csv")
  453. # Extract gene names from target column
  454. rhamp <- mutate(rhamp, class = case_when(
  455. str_detect(target, "^cdx4") ~ "cdx4",
  456. str_detect(target, "^rx3") ~ "rx3",
  457. str_detect(target, "^tbx16") ~ "tbx16",
  458. str_detect(target, "^tbxta") ~ "tbxta",
  459. str_detect(target, "^hand2") ~ "hand2",
  460. str_detect(target, "^hoxb1b") ~ "hoxb1b",
  461. str_detect(target, "^foxa2") ~ "foxa2",
  462. str_detect(target, "^tyr") ~ "tyr",
  463. TRUE ~ "Other"
  464. ))
  466. # Calculate percentage of cells with each genotype
  467. cellfreq <- data.frame(table(mdp$class)) %>%
  468. filter(Var1 != "ND", Var1 != "Multiple") %>%
  469. mutate(pct = (Freq / sum(Freq)) * 100) %>%
  470. select(Var1, pct)
  472. names(cellfreq) <- c("class", "cell_pct")
  473. cellfreq$class <- gsub("Non-Targeting", "tyr", cellfreq$class)
  475. # Merge with rhAmpSeq data, keep highest edited site per gene
  476. rhamp <- rhamp %>%
  477. left_join(cellfreq, by = "class") %>%
  478. select(class, target, pct_indel, wt_indel, cell_pct) %>%
  479. group_by(class) %>%
  480. filter(pct_indel == max(pct_indel)) %>%
  481. ungroup() %>%
  482. mutate(class = factor(class, levels = c("hoxb1b", "rx3", "cdx4", "hand2",
  483. "tbx16", "tbxta", "foxa2", "tyr")))
  484. ```
  486. ```{r figure1E, fig.width=2.5, fig.height=3}
  487. # Normalize indel rate by subtracting background
  488. rhamp_normalized <- rhamp %>%
  489. mutate(indel_normalized = pct_indel - wt_indel)
  491. # Dumbbell plot comparing DNA editing to gRNA detection
  492. ggplot(rhamp_normalized) +
  493. geom_segment(aes(x = class, xend = class,
  494. y = indel_normalized, yend = cell_pct),
  495. color = "grey") +
  496. geom_point(aes(x = class, y = indel_normalized,
  497. color = factor("% DNA with indel",
  498. levels = c("% DNA with indel", "% Cells with gRNA"))),
  499. size = 3) +
  500. geom_point(aes(x = class, y = cell_pct,
  501. color = factor("% Cells with gRNA",
  502. levels = c("% DNA with indel", "% Cells with gRNA"))),
  503. size = 3) +
  504. scale_color_manual(values = c("black", "#888888")) +
  505. coord_flip() +
  506. scale_y_continuous(limits = c(0, 17), breaks = c(0, 5, 10, 15)) +
  507. theme_classic() +
  508. theme(
  509. text = element_text(family = "Helvetica"),
  510. axis.text.x = element_text(size = 12),
  511. axis.text.y = element_text(size = 12, face = "italic", color = "black"),
  512. axis.title = element_text(size = 12),
  513. legend.text = element_text(size = 12),
  514. legend.title = element_blank(),
  515. legend.position = c(0.4, -0.25),
  516. plot.margin = margin(b = 70),
  517. legend.background = element_rect(fill = alpha("white", 0))
  518. ) +
  519. xlab("") +
  520. ylab(element_blank()) +
  521. guides(color = guide_legend(title = "", ncol = 1))
  523. ggsave(filename = "outputs_fig1/Figure_1E.eps",
  524. device = "ps",
  525. width = 2.5, height = 3, units = 'in',
  526. bg = "transparent")
  527. ```
  529. # =============================================================================
  530. # SECTION 9: DIMENSIONALITY REDUCTION AND CLUSTERING
  531. # =============================================================================
  533. Remove doublets and perform standard Seurat workflow
  535. ```{r clustering}
  536. # Remove cells with multiple gRNA targets (doublets)
  537. mdp_filt <- subset(mdp, subset = class != "Multiple")
  539. # Standard Seurat workflow
  540. mdp_filt <- NormalizeData(mdp_filt)
  541. mdp_filt <- FindVariableFeatures(mdp_filt, selection.method = "vst", nfeatures = 2000)
  542. mdp_filt <- ScaleData(mdp_filt)
  543. mdp_filt <- RunPCA(mdp_filt, features = VariableFeatures(object = mdp_filt))
  544. mdp_filt <- FindNeighbors(mdp_filt, dims = 1:30)
  545. mdp_filt <- FindClusters(mdp_filt, resolution = 0.5)
  546. mdp_filt <- RunUMAP(mdp_filt, dims = 1:30, min.dist = 0.4)
  548. DimPlot(mdp_filt)
  549. ```
  551. # =============================================================================
  552. # SECTION 10: CELL TYPE ANNOTATION VIA LABEL TRANSFER
  553. # =============================================================================
  555. Use Daniocell reference atlas to annotate cell types
  557. ```{r label-transfer}
  558. # Load Daniocell reference (21-26 hpf timepoints)
  559. dcell_ref <- readRDS(file = "input_data/dcell_21_26.rds")
  561. # Prepare reference
  562. dcell_ref <- RunPCA(dcell_ref, dims = 1:30)
  563. dcell_ref <- RunUMAP(dcell_ref, dims = 1:30, return.model = TRUE)
  565. # Separate cell type hierarchy into columns
  566. [email hidden] <- separate([email hidden],
  567. full_ident,
  568. into = c("tissue", "cell_type", "sub_type"),
  569. sep = "[|>]",
  570. remove = FALSE)
  572. # Find transfer anchors and transfer labels
  573. xfer_anchors <- FindTransferAnchors(
  574. reference = dcell_ref,
  575. query = mdp_filt,
  576. dims = 1:30,
  577. reference.assay = "RNA",
  578. query.assay = "RNA",
  579. reference.reduction = "pca"
  580. )
  582. predictions <- TransferData(
  583. anchorset = xfer_anchors,
  584. refdata = dcell_ref$full_ident,
  585. dims = 1:30
  586. )
  588. # Add predictions to query
  589. mdp_filt <- AddMetaData(mdp_filt, metadata = predictions)
  591. # Map query to reference UMAP
  592. mdp_filt <- MapQuery(
  593. anchorset = xfer_anchors,
  594. reference = dcell_ref,
  595. query = mdp_filt,
  596. refdata = "full_ident",
  597. reference.reduction = "pca",
  598. reduction.model = "umap"
  599. )
  601. # Separate predicted cell type hierarchy
  602. [email hidden] <- separate([email hidden],
  603. predicted.id,
  604. into = c("tissue", "cell_type", "sub_type"),
  605. sep = "[|>]",
  606. remove = FALSE)
  607. ```
  609. ## Figure S4: Label transfer validation
  611. ```{r figureS4A, fig.width=4, fig.height=4}
  612. # S4A: Reference atlas tissue types
  613. p1 <- DimPlot(dcell_ref, group.by = "tissue", label = FALSE, reduction = 'umap')
  614. LabelClusters(p1, id = "tissue", repel = TRUE) +
  615. NoLegend() + NoAxes() +
  616. ggtitle("Daniocell reference tissue types") +
  617. theme(plot.title = element_text(size = 12, family = "Helvetica"))
  619. ggsave(filename = "outputs_fig1/Figure_S4A.eps",
  620. device = "ps",
  621. width = 4, height = 4, units = 'in',
  622. bg = "transparent")
  623. ```
  625. ```{r figureS4B, fig.width=4, fig.height=4}
  626. # S4B: Query cells projected to reference
  627. DimPlot(mdp_filt, group.by = "tissue", label = FALSE, reduction = 'ref.umap') +
  628. NoLegend() + NoAxes() +
  629. ggtitle("MIC-Drop-seq cells projected to reference") +
  630. theme(plot.title = element_text(size = 12, family = "Helvetica"))
  632. ggsave(filename = "outputs_fig1/Figure_S4B.eps",
  633. device = "ps",
  634. width = 4, height = 4, units = 'in',
  635. bg = "transparent")
  636. ```
  638. ## Assign consensus cell types to clusters
  640. Take the most frequent predicted cell type for each cluster
  642. ```{r consensus-labels}
  643. consensus_assignments <- [email hidden] %>%
  644. group_by(seurat_clusters) %>%
  645. dplyr::count(predicted.id) %>%
  646. top_n(n = 1, wt = n) %>%
  647. select(seurat_clusters, predicted.id) %>%
  648. separate(predicted.id,
  649. into = c("consensus_tissue", "consensus_cell_type", "consensus_sub_type"),
  650. sep = "[|>]",
  651. remove = FALSE) %>%
  652. dplyr::rename("consensus_predicted.id" = "predicted.id")
  654. # Merge back with metadata
  655. [email hidden] <- [email hidden] %>%
  656. left_join(consensus_assignments, by = "seurat_clusters")
  658. rownames([email hidden]) <- [email hidden]$cell_barcode
  659. ```
  661. ```{r figureS4C, fig.width=4, fig.height=4}
  662. # S4C: Consensus tissue labels
  663. p1 <- DimPlot(mdp_filt, group.by = "consensus_tissue", label = FALSE)
  664. LabelClusters(p1, id = "consensus_tissue", repel = TRUE, size = 4) +
  665. NoLegend() + NoAxes() +
  666. ggtitle("MIC-Drop-seq consensus tissue labels") +
  667. theme(plot.title = element_text(size = 12, family = "Helvetica"))
  669. ggsave(filename = "outputs_fig1/Figure_S4C.eps",
  670. device = "ps",
  671. width = 4, height = 4, units = 'in',
  672. bg = "transparent")
  673. ```
  675. ## Figure 1D: Cluster UMAP
  677. ```{r figure1D, fig.width=5, fig.height=5}
  678. # Custom color palette for 35 clusters
  679. palette <- c('#AADAEEFF', '#FCDACAFF', '#DCD4E6FF', '#92C8E0FF', '#FFBDF2FF',
  680. '#7BB6D3FF', '#64A4C5FF', '#C4DFB6FF', '#FCE4A6FF', '#EC85D8FF',
  681. '#ACCF9FFF', '#EBC576FF', '#F0ABAEFF', '#458DB3FF', '#D69E3BFF',
  682. '#AF7623FF', '#95C089FF', '#2675A1FF', '#1A6088FF', '#D35C61FF',
  683. '#7EB173FF', '#619F57FF', '#9A70ABFF', '#545454FF', '#8B1713FF',
  684. '#104B6FFF', '#428B39FF', '#D747BBFF', '#B90497FF', '#44236EFF',
  685. '#277822FF', '#1D6623FF', '#7C000CFF', '#770262FF', '#135524FF')
  687. p1 <- DimPlot(mdp_filt, cols = palette, group.by = 'seurat_clusters') +
  688. NoLegend() + NoAxes() + ggtitle("")
  690. LabelClusters(p1, id = "seurat_clusters", bg.color = "white")
  692. ggsave(filename = "outputs_fig1/Figure_1D.eps",
  693. device = "ps",
  694. width = 5, height = 5, units = 'in',
  695. bg = "transparent")
  696. ```
  698. # =============================================================================
  699. # SECTION 11: PHENOTYPE VALIDATION
  700. # =============================================================================
  702. ## Presomitic Mesoderm (tbx16, tbxta)
  704. Known phenotypes:
  705. - tbx16 mutants: Accumulation of presomitic mesoderm
  706. - tbxta mutants: Depletion of presomitic mesoderm
  708. ```{r presomitic-mesoderm}
  709. # Add presomitic mesoderm annotation
  710. [email hidden] <- mutate([email hidden], psm = case_when(
  711. seurat_clusters == 15 ~ "Presomitic Mesoderm",
  712. TRUE ~ "Other"
  713. ))
  715. # Calculate percentage of cells per genotype in PSM
  716. psmstats <- data.frame(table(mdp_filt$psm, mdp_filt$class))
  717. totals <- psmstats %>% group_by(Var2) %>% summarise(total = sum(Freq))
  719. psmstats <- left_join(psmstats, totals, by = "Var2") %>%
  720. mutate(pct = (Freq / total) * 100) %>%
  721. filter(Var1 == "Presomitic Mesoderm", Var2 != "ND", Var2 != "Multiple")
  723. psmstats$Var2 <- gsub("Non-Targeting", "tyr", psmstats$Var2)
  724. psmstats$Var2 <- factor(psmstats$Var2,
  725. levels = c("cdx4", "foxa2", "hand2", "hoxb1b",
  726. "rx3", "tbx16", "tbxta", "tyr"))
  727. ```
  729. ### Figure 1F: PSM UMAP split by genotype
  731. ```{r figure1F, fig.width=8, fig.height=1.5}
  732. Idents(mdp_filt) <- "seurat_clusters"
  733. PSM <- subset(mdp_filt, idents = "15")
  735. [email hidden] <- mutate([email hidden], psmplot = case_when(
  736. class == "tbx16" ~ "tbx16",
  737. class == "tyr" ~ "tyr",
  738. class == "tbxta" ~ "tbxta",
  739. TRUE ~ "Other"
  740. ))
  742. dittoDimPlot(PSM, "psmplot", split.by = "psmplot",
  743. color.panel = c("#000000", "black", "black", "black"),
  744. size = 0.5, split.nrow = 1) +
  745. NoAxes() + NoLegend() +
  746. coord_cartesian(xlim = c(-14, -9), ylim = c(-1, 2.5)) +
  747. theme(strip.text = element_text(size = 12, face = "italic", family = "Helvetica"),
  748. strip.background = element_rect(fill = "white")) +
  749. ggtitle("")
  751. ggsave(filename = "outputs_fig1/Figure_1F.eps",
  752. device = "ps",
  753. width = 8, height = 1.5, units = 'in',
  754. bg = "transparent")
  755. ```
  757. ### Figure 1G: PSM percentage barplot
  759. ```{r figure1G, fig.width=3, fig.height=2}
  760. ggplot(data = psmstats, aes(x = Var2, y = pct)) +
  761. geom_bar(stat = "identity", fill = "black") +
  762. theme_classic() +
  763. ylab("% cells in cluster") +
  764. xlab("Mutation") +
  765. theme(
  766. axis.text.x = element_text(size = 12, family = "Helvetica",
  767. color = "black", face = "italic",
  768. angle = 45, hjust = 1),
  769. axis.text.y = element_text(size = 12, family = "Helvetica", color = "black"),
  770. axis.title = element_text(size = 12, family = "Helvetica")
  771. )
  773. ggsave(filename = "outputs_fig1/Figure_1G.eps",
  774. device = "ps",
  775. width = 3, height = 2, units = 'in',
  776. bg = "transparent")
  777. ```
  779. ## Optic Primordia (rx3)
  781. Known phenotype: rx3 mutants lack eyes (depletion of retinal cells)
  783. ```{r optic-primordia}
  784. # Add optic primordia annotation
  785. [email hidden] <- mutate([email hidden], eye = case_when(
  786. seurat_clusters %in% c(2, 22) ~ "Optic Primordia",
  787. TRUE ~ "Other"
  788. ))
  790. # Calculate percentage of cells per genotype in optic primordia
  791. eyestats <- data.frame(table(mdp_filt$eye, mdp_filt$class))
  792. totals <- eyestats %>% group_by(Var2) %>% summarise(total = sum(Freq))
  794. eyestats <- left_join(eyestats, totals, by = "Var2") %>%
  795. mutate(pct = (Freq / total) * 100) %>%
  796. filter(Var1 == "Optic Primordia", Var2 != "ND", Var2 != "Multiple")
  798. eyestats$Var2 <- gsub("Non-Targeting", "tyr", eyestats$Var2)
  799. eyestats$Var2 <- factor(eyestats$Var2,
  800. levels = c("cdx4", "foxa2", "hand2", "hoxb1b",
  801. "rx3", "tbx16", "tbxta", "tyr"))
  802. ```
  804. ### Figure S5A: Optic primordia UMAP
  806. ```{r figureS5A, fig.width=6, fig.height=1.5}
  807. optic <- subset(mdp_filt, idents = c(2, 22))
  809. [email hidden] <- mutate([email hidden], opticplot = case_when(
  810. class == "rx3" ~ "rx3",
  811. class == "tyr" ~ "tyr",
  812. TRUE ~ "Other"
  813. ))
  815. dittoDimPlot(optic, "opticplot", split.by = "opticplot",
  816. color.panel = c("#000000", "black", "black", "black"),
  817. size = 0.5, split.nrow = 1) +
  818. NoAxes() + NoLegend() +
  819. coord_cartesian(xlim = c(2, 11), ylim = c(-7, -2)) +
  820. theme(strip.text = element_text(size = 12, face = "italic", family = "Helvetica"),
  821. strip.background = element_rect(fill = "white")) +
  822. ggtitle("")
  824. ggsave(filename = "outputs_fig1/Figure_S5A.eps",
  825. device = "ps",
  826. width = 6, height = 1.5, units = 'in',
  827. bg = "transparent")
  828. ```
  830. ### Figure S5B: Optic primordia percentage
  832. ```{r figureS5B, fig.width=3, fig.height=2}
  833. ggplot(data = eyestats, aes(x = Var2, y = pct)) +
  834. geom_bar(stat = "identity", fill = 'black') +
  835. theme_classic() +
  836. ylab("% cells in cluster") +
  837. xlab("Mutation") +
  838. theme(
  839. axis.text.x = element_text(size = 12, family = "Helvetica",
  840. color = "black", face = "italic",
  841. angle = 45, hjust = 1),
  842. axis.text.y = element_text(size = 12, family = "Helvetica", color = "black"),
  843. axis.title = element_text(size = 12, family = "Helvetica")
  844. )
  846. ggsave(filename = "outputs_fig1/Figure_S5B.eps",
  847. device = "ps",
  848. width = 3, height = 2, units = 'in',
  849. bg = "transparent")
  850. ```
  852. ## Spinal Cord (cdx4)
  854. Known phenotype: cdx4 regulates hox gene expression in spinal cord
  856. ```{r spinal-cord-deg}
  857. # Add cdx4 classification
  858. [email hidden] <- [email hidden] %>%
  859. mutate(cdx = ifelse(class == "cdx4", "cdx4", "other"))
  860. ```
  862. ### Figure 1H: Spinal cord UMAP
  864. ```{r figure1H, fig.width=3, fig.height=1.5}
  865. Idents(mdp_filt) <- "consensus_cell_type"
  866. sc <- subset(mdp_filt, idents = "spinal cord")
  868. [email hidden] <- mutate([email hidden], cdx = ifelse(class == "cdx4", "cdx4", "other"))
  870. dittoDimPlot(sc, "cdx", split.by = "cdx",
  871. color.panel = c("#000000", "black"),
  872. size = 0.5, split.nrow = 1) +
  873. NoAxes() + NoLegend() +
  874. coord_cartesian(xlim = c(-5, 0), ylim = c(-5, 2)) +
  875. theme(strip.text = element_text(size = 12, face = c("italic", "plain"), family = "Helvetica"),
  876. strip.background = element_rect(fill = "white")) +
  877. ggtitle("")
  879. ggsave(filename = "outputs_fig1/Figure_1H.eps",
  880. device = "ps",
  881. width = 3, height = 1.5, units = 'in',
  882. bg = "transparent")
  883. ```
  885. ### Figure 1I: Hox gene expression in spinal cord
  887. ```{r figure1I, fig.width=3, fig.height=2.5}
  888. Idents(mdp_filt) <- "consensus_tissue"
  890. VlnPlot(mdp_filt,
  891. group.by = "cdx",
  892. features = c('elavl3', 'hoxb9a', 'hoxc6b', 'hoxa9a'),
  893. idents = "spinal cord",
  894. stack = TRUE,
  895. flip = TRUE,
  896. cols = rep("darkgrey", 5)) +
  897. NoLegend() +
  898. scale_x_discrete(labels = c(expression(italic(cdx4), "other"))) +
  899. xlab("Mutation")
  901. ggsave(filename = "outputs_fig1/Figure_1I.eps",
  902. device = "ps",
  903. width = 3, height = 2.5, units = 'in',
  904. bg = "transparent")
  905. ```
  907. ### Figure S5C: cdx4 volcano plot
  909. ```{r figureS5C, fig.width=3, fig.height=3}
  910. # Identify cdx4 and control cells in spinal cord clusters
  911. cdx_spinal_cells <- [email hidden] %>%
  912. filter(seurat_clusters %in% c(4, 28), class == 'cdx4') %>%
  913. pull(cell_barcode)
  915. control_cells <- [email hidden] %>%
  916. filter(seurat_clusters %in% c(4, 28), class != "cdx4") %>%
  917. pull(cell_barcode)
  919. # Find differentially expressed genes
  920. cdx_deg <- FindMarkers(mdp_filt,
  921. ident.1 = cdx_spinal_cells,
  922. ident.2 = control_cells)
  924. # Prepare for plotting
  925. genes_of_interest <- c("hoxb9a", "hoxa9a", "hoxc6b")
  927. cdx_deg <- cdx_deg %>%
  928. rownames_to_column("gene") %>%
  929. mutate(
  930. is_highlight = gene %in% genes_of_interest,
  931. log_p = -log10(p_val_adj)
  932. )
  934. highlight_subset <- cdx_deg %>% filter(is_highlight)
  936. # Volcano plot
  937. ggplot(cdx_deg, aes(x = avg_log2FC, y = log_p)) +
  938. geom_point(color = "grey70", alpha = 0.6, size = 1.5) +
  939. geom_point(data = highlight_subset, color = "black", size = 3) +
  940. geom_label_repel(data = highlight_subset,
  941. aes(label = gene),
  942. fontface = "italic",
  943. force = 2) +
  944. theme_classic() +
  945. labs(
  946. title = expression(italic("cdx4") ~ "- Spinal Cord DEGs"),
  947. x = "Average Log2 Fold Change",
  948. y = "-Log10 Adjusted P-value"
  949. ) +
  950. geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  951. geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  952. theme(plot.title = element_text(hjust = 0.5))
  954. ggsave(filename = "outputs_fig1/Figure_S5C.eps",
  955. device = "ps",
  956. width = 3, height = 3, units = 'in',
  957. bg = "transparent")
  958. ```
  960. # =============================================================================
  961. # SECTION 12: SUPPLEMENTARY FIGURES
  962. # =============================================================================
  964. ## Figure S3A: Cells per genotype
  966. ```{r figureS3A, fig.width=4, fig.height=3}
  967. [email hidden] <- [email hidden] %>%
  968. mutate(class = factor(class,
  969. levels = c("hoxb1b", "rx3", "cdx4", "hand2",
  970. "tbx16", "tbxta", "foxa2", "tyr", "ND")))
  972. [email hidden] %>%
  973. dplyr::count(class) %>%
  974. ggplot(aes(x = class, y = n)) +
  975. geom_bar(stat = "identity", fill = "black") +
  976. theme_classic() +
  977. labs(y = "Assigned Cells", x = "Genotype") +
  978. scale_y_continuous(breaks = seq(0, 9000, by = 1000)) +
  979. theme(
  980. axis.title = element_text(size = 12, family = "Helvetica"),
  981. axis.text.y = element_text(size = 12, family = "Helvetica"),
  982. axis.text.x = element_text(angle = 45, hjust = 1, size = 12,
  983. face = c(rep("italic", 8), 'plain'),
  984. family = "Helvetica")
  985. )
  987. ggsave(filename = "outputs_fig1/Figure_S3A.eps",
  988. device = "ps",
  989. width = 4, height = 3, units = 'in',
  990. bg = "transparent")
  991. ```
  993. ## Figure S3B: UMAP split by genotype
  995. ```{r figureS3B, fig.width=6, fig.height=6}
  996. dittoDimPlot(mdp_filt, "class", split.by = "class") +
  997. NoLegend() + NoAxes() +
  998. theme(strip.text = element_text(size = 12,
  999. face = c(rep("italic", 8), "plain"),
  1000. family = "Helvetica")) +
  1001. ggtitle("")
  1003. ggsave(filename = "outputs_fig1/Figure_S3B.eps",
  1004. device = "ps",
  1005. width = 6, height = 6, units = 'in',
  1006. bg = "transparent")
  1007. ```
  1011. Figure S3D
  1012. ```{r}
  1013. #read amplican results
  1014. amp_res <- read.csv(file = "input_data/amplican_results_summary.csv")
  1017. amp_res_normalized <- amp_res %>%
  1018. #calculate the percentage of reads frameshifted
  1019. mutate(frameshift_pct = Reads_Frameshifted / Reads_Filtered * 100) %>%
  1020. #calculate the proportion of frameshifted reads relative to the percentage of cells assigned the genotype
  1021. mutate(non_frameshift_pr_normalized = 1 - (frameshift_pct / Cell_pct),
  1022. frameshift_pr_normalized = (frameshift_pct / Cell_pct))
  1024. cumul_frameshift <- amp_res_normalized %>%
  1025. #Now we can group the rows by mutation, and then we summarize to a new column to show the cumulative frameshift efficiency
  1026. group_by(Group) %>%
  1027. summarise(cumulative_frameshift_efficiency = 1 - prod(non_frameshift_pr_normalized) )
  1029. amp_results_final <- amp_res_normalized %>%
  1030. left_join(cumul_frameshift, by = "Group") %>%
  1031. mutate(Group = gsub("Non-Targeting","tyr",.$Group))
  1033. amp_results_plot_data <- amp_results_final %>%
  1034. #select only relevant columns
  1035. dplyr::select(ID, Group, frameshift_pr_normalized, cumulative_frameshift_efficiency ) %>%
  1036. pivot_longer(names_to = "name", cols = -c(Group, ID))
  1039. p1 <- ggplot(amp_results_plot_data, aes(x = value, y = Group, color = name)) +
  1040. geom_point(size = 3) +
  1041. scale_x_continuous(limits = c(0,1)) +
  1042. scale_color_manual(
  1043. values = c("dodgerblue", "darkgrey"),
  1044. labels = c("Cumulative Frameshift Probability", "Individual Guide Frameshift Frequency")
  1045. ) +
  1046. labs(
  1047. x = "Normalized Frameshift Frequency",
  1048. y = "Target Gene",
  1049. color = "" # Blank legend title
  1050. ) +
  1051. theme_minimal() +
  1052. theme(
  1053. text = element_text(size = 12),
  1054. axis.text = element_text(size = 12), # Axis numbers/labels
  1055. axis.title = element_text(size = 12), # Axis titles
  1056. legend.text = element_text(size = 12), # Legend text
  1057. plot.title = element_text(size = 12), # Plot title (if you add one)
  1058. panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
  1059. legend.position = "bottom",
  1060. legend.direction = "vertical"
  1061. )
  1063. p1
  1065. ggsave(filename = "outputs_fig1/Figure_S3D.eps",
  1066. device = "ps",
  1067. width = 3, height = 3, units = 'in',
  1068. bg = "transparent")
  1071. ```
  1073. Figure S1C plot
  1074. ```{r}
  1075. # Load libraries
  1076. library(ggplot2)
  1077. library(dplyr)
  1079. # 1. Data Setup (based on your table)
  1080. data <- data.frame(
  1081. gRNA = rep(c("Unmodified gRNA", "Modified gRNA"), each = 6),
  1082. phenotype = rep(c("Full", "Full", "Partial", "Partial", "WT", "WT"), 2),
  1083. replicate = rep(c(1, 2), 6),
  1084. percentage = c(88, 94.7, 12, 5.3, 0, 0, # Unmodified
  1085. 90.4, 94.1, 9.6, 5.9, 0, 0) # Modified
  1086. )
  1088. # 2. Data Processing
  1089. # Set 'Full' as the base layer
  1090. data$phenotype <- factor(data$phenotype, levels = c("Full", "Partial", "WT"))
  1092. # Calculate means for the bars
  1093. bar_data <- data %>%
  1094. group_by(gRNA, phenotype) %>%
  1095. summarize(mean_perc = mean(percentage), .groups = 'drop')
  1097. # Filter dots to show ONLY the 'Full' phenotype variance
  1098. dot_data <- data %>%
  1099. filter(phenotype == "Full")
  1101. # 3. Create the Plot
  1102. p1 <- ggplot() +
  1103. # Bars: Show the average composition
  1104. geom_col(data = bar_data, aes(x = gRNA, y = mean_perc, fill = phenotype),
  1105. color = "black", width = 0.6, position = position_stack(reverse = TRUE)) +
  1107. # Dots: Show individual variance for the 'Full' phenotype only
  1108. geom_point(data = dot_data,
  1109. aes(x = gRNA, y = percentage),
  1110. shape = 21,
  1111. fill = "black", # White fill makes dots pop against black bars
  1112. color = "white",
  1113. size = 3.5,
  1114. stroke = 1,
  1115. position = position_dodge(width = 0.4),
  1116. show.legend = FALSE) +
  1118. # Aesthetics and Styling
  1119. scale_fill_manual(values = c("Full" = "black",
  1120. "Partial" = "grey70",
  1121. "WT" = "white")) +
  1122. scale_y_continuous(limits = c(0, 105), expand = c(0, 0), breaks = seq(0, 100, 20)) +
  1123. labs(y = "% embryos with phenotype", x = NULL, fill = NULL) +
  1124. theme_classic() +
  1125. theme(
  1126. axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 12, color = "black"),
  1127. axis.text.y = element_text(size = 12, color = "black"),
  1128. legend.position = "right"
  1129. )
  1130. p1
  1131. ggsave(filename = "outputs_fig1/Figure_S1C.eps",
  1132. device = "ps",
  1133. width = 3, height = 4, units = 'in',
  1134. bg = "transparent")
  1136. ```

Figure1_MIC-Drop-seq.Rmd at commit 9de3c19, under MIT · at the source

Overview

Authors: Clayton M. Carey1, Saba Parvez2,3, Zachary J. Brandt2, Brent W. Bisgrove1,2, Christopher J. Yates2, Randall T. Peterson2, James A. Gagnon1,4
ORCID iDs: Brent W. Bisgrove, Christopher J. Yates, Randall T. Peterson, James A. Gagnon
  1. School of Biological Sciences, University of Utah,Salt Lake City, UT USA
  2. Department of Pharmacology and Toxicology, College of Pharmacy, University of Utah,Salt Lake City, UT USA
  3. Department of Cell & Developmental Biology, Feinberg School of Medicine, Northwestern University,Chicago, IL USA
  4. Henry Eyring Center for Cell & Genome Science, University of Utah,Salt Lake City, UT USA
Institutions: University of Utah (United States); Northwestern University (United States)
Journal: Nature communications, volume 17, issue 1, article 4738
Dates: received 9 June 2025; accepted 10 March 2026; published online 1 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-70989-w · PMID 41922342 · PMCID PMC13216532 · OpenAlex W7147444157
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: zebrafish (organism)
Methods: Statistics, Smoothing, state filtering, decompositions, Evoked potentials, fMRI & imaging
Keywords: Cell lineage, Development, Genomic engineering, CRISPR-Cas systems
MeSH: Embryo, Nonmammalian*, Single-Cell Analysis*, Zebrafish*, Animals, CRISPR-Cas Systems, Gene Expression Regulation, Developmental, Gene Regulatory Networks, Mesoderm, Mutation, Phenotype, Sequence Analysis, RNA, Single-Cell Gene Expression Analysis, Transcription Factors, Zebrafish Proteins (* major topic)
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: NIH (F32HL156644, K01HG013682, R00HG012593, R01GM134069, R35GM142950, R24OD035409)
Citations: not cited yet (Europe PMC); 65 references in the paper

Abstract

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

Repositories

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

clay-carey/MIC-Drop-seq

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 9de3c197ea1f6a13758e52b6d2b419ceb3cb0582, 11 February 2026
Languages: R (4)
Size: 8 files, 4 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 4 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (4 files), Seurat (4 files), tidyverse (4 files), ComplexHeatmap (3 files), edgeR (2 files), ggpubr (2 files), patchwork (2 files), rstatix (2 files), circlize (1 file), cowplot (1 file), igraph (1 file), Monocle 3 (1 file), SingleCellExperiment (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
6 files
At the source: github.com/clay-carey/MIC-Drop-seq

Zenodo 18615355

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (4 files), Seurat (4 files), tidyverse (4 files), ComplexHeatmap (3 files), edgeR (2 files), ggpubr (2 files), patchwork (2 files), rstatix (2 files), circlize (1 file), cowplot (1 file), igraph (1 file), Monocle 3 (1 file), SingleCellExperiment (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
6 files
At the source:

Code availability statement

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

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

Tracing map

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

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 8 scripts, each with its path and the digest of its content;
  • 20 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/s41467-026-70989-w.

Versions

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

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 4 keywords, 14 MeSH terms, 1 funder, 64 references.

Cite

This paper

Carey, C. M., Parvez, S., Brandt, Z. J., Bisgrove, B. W., Yates, C. J., Peterson, R. T., & Gagnon, J. A. (2026). MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos. Nature communications, 17(1), 4738. https://doi.org/10.1038/s41467-026-70989-w

BibTeX

@article{carey2026mic,
author = {Carey, Clayton M. and Parvez, Saba and Brandt, Zachary J. and Bisgrove, Brent W. and Yates, Christopher J. and Peterson, Randall T. and Gagnon, James A.},
title = {{MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {4738},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-70989-w},
url = {https://doi.org/10.1038/s41467-026-70989-w},
pmid = {41922342},
pmcid = {PMC13216532}
}

RIS

TY - JOUR
AU - Carey, Clayton M.
AU - Parvez, Saba
AU - Brandt, Zachary J.
AU - Bisgrove, Brent W.
AU - Yates, Christopher J.
AU - Peterson, Randall T.
AU - Gagnon, James A.
TI - MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/01
VL - 17
IS - 1
SP - 4738
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-70989-w
UR - https://doi.org/10.1038/s41467-026-70989-w
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-70989-w",
"type": "article-journal",
"title": "MIC-Drop-seq: scalable single-cell phenotyping of mutant vertebrate embryos",
"container-title": "Nature communications",
"author": [
{
"family": "Carey",
"given": "Clayton M."
},
{
"family": "Parvez",
"given": "Saba"
},
{
"family": "Brandt",
"given": "Zachary J."
},
{
"family": "Bisgrove",
"given": "Brent W."
},
{
"family": "Yates",
"given": "Christopher J."
},
{
"family": "Peterson",
"given": "Randall T."
},
{
"family": "Gagnon",
"given": "James A."
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "4738",
"DOI": "10.1038/s41467-026-70989-w",
"PMID": "41922342",
"PMCID": "PMC13216532",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-70989-w",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
1
]
]
}
}

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: Monocle 3, SingleCellExperiment, edgeR, 9 other tools, 2 references
[2] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Monocle 3, SingleCellExperiment, edgeR, 9 other tools, 2 references
[3] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Monocle 3, SingleCellExperiment, edgeR, 10 other tools
[4] doi:10.1186/s12967-026-08266-z [code]
Single-cell multi-omic integration analysis prioritizes druggable genes and reveals cell-type-specific causal effects in glioblastomagenesis.
Journal: Journal of translational medicine
In common: Monocle 3, edgeR, igraph, 8 other tools, 2 references
[5] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Monocle 3, SingleCellExperiment, edgeR, 9 other tools
[6] doi:10.1038/s41467-026-76341-6 [code]
Neonatal inflammation disrupts a temporally restricted postnatal Numb-enriched microglial state in mice.
Journal: Nature communications
In common: Monocle 3, SingleCellExperiment, igraph, 8 other tools, 1 reference
[7] doi:10.1038/s41398-026-04200-5 [code]
Postmortem brain single-nucleus and bulk gene expression analyses identify shared and distinct abnormalities in bipolar disorder and major depressive disorder.
Journal: Translational psychiatry
In common: SingleCellExperiment, edgeR, rstatix, 8 other tools, 1 reference
[8] doi:10.1038/s44318-026-00806-z [code]
Interspecific diversity in the neuronal composition of the mammalian cortex arises from heterochrony in neurogenesis.
Journal: The EMBO journal
In common: Monocle 3, SingleCellExperiment, igraph, 8 other tools, 1 reference
[9] doi:10.1016/j.celrep.2026.117500 [code]
Spatio-molecular gene expression reflects dorsal anterior cingulate cortex structure and function in the human brain.
Journal: Cell reports
In common: SingleCellExperiment, edgeR, igraph, 8 other tools, 1 reference
[10] doi:10.1038/s41467-026-71595-6 [code]
A single-cell and spatial atlas of early human olfactory development.
Journal: Nature communications
In common: SingleCellExperiment, igraph, circlize, 7 other tools, 2 references

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.