A single-nucleus RNA-seq dataset of the colon in Pink1-deficient and wild-type mice.
The 7 matches
- [1] § Technical Validation › Anchor-based label transfer and validation ↔ Code.R, lines 361–448 · score 0.89 · FindTransferAnchors, MapQuery, reference.reduction, wt_pink1ko, LogNormalize, pcaproject
- [2] § Methods › Single-nucleus RNA-seq data processing and analysis ↔ Code.R, lines 1–63 · score 0.84 · FindNeighbors, IntegrateLayers, variable features, Seurat, QC, resolution
- [3] § Background & Summary ↔ Code.R, lines 201–258 · score 0.67 · intestinal epithelial cells, goblet cells, stem cells, colonocytes, clustering, genes
- [4] § Technical Validation › Anchor-based label transfer and validation ↔ Code.R, lines 260–310 · score 0.59 · Col1a1, Acta2, Lyz2, Pdgfra, Ptprc, fibroblasts
- [5] § Technical Validation › Anchor-based label transfer and validation ↔ Code.R, lines 739–816 · score 0.57 · decrease Gini, decrease accuracy, MDA, MDG, class, predicted
- [6] § Technical Validation › Anchor-based label transfer and validation ↔ Code.R, lines 739–816 · score 0.56 · Decrease Gini, Decrease Accuracy, Random forest, MDA, MDG, predicted cell
- [7] § Background & Summary ↔ Code.R, lines 201–258 · score 0.55 · intestinal epithelium, stem cells, DEGs, colonocytes, immune, gene
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 · 816 lines · 28 KB · CC-BY-4.0 · 7 matches
- library(Seurat)
- library(SeuratObject)
- library(dplyr)
- library(Matrix)
- library(ggplot2)
- library(patchwork)
- # Loading datasets
- WT_data <- Read10X(data.dir = "WT")
- WT <- CreateSeuratObject(counts = WT_data, project = "WT")
- PINK1_data <- Read10X(data.dir = "PINK1")
- PINK1 <- CreateSeuratObject(counts = PINK1_data, project = "PINK1")
- # Assigning the dataset metadata
- WT$dataset <- "WT"
- PINK1$dataset <- "PINK1"
- WT <- JoinLayers(WT)
- PINK1 <- JoinLayers(PINK1)
- # Merging datasets
- merged <- merge(WT, y = PINK1, add.cell.ids = c("WT", "PINK1"), project = "Gut_tissue")
- Idents(object = merged) <- "Gut_tissue"
- merged[["percent.mt"]] <- PercentageFeatureSet(merged, pattern = "^mt-")
- # QC metrics
- VlnPlot(merged, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
- merged <- subset(merged, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 &
- nCount_RNA > 500 & nCount_RNA < 10000 &
- percent.mt < 10)
- FeatureScatter(merged, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = FALSE)
- dim(WT)
- dim(merged)
- head(merged)
- # preprocessing
- merged <- NormalizeData(merged)
- merged <- FindVariableFeatures(merged)
- merged <- ScaleData(merged)
- merged <- RunPCA(merged)
- DimPlot(merged, reduction = "pca", group.by = "orig.ident") + ggtitle("PCA After Integration")
- ElbowPlot(merged)
- merged <- FindNeighbors(merged, dims = 1:20, reduction = "pca")
- merged <- FindClusters(merged, resolution = 2, cluster.name = "unintegrated_clusters")
- merged <- RunUMAP(merged, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
- merged <- RunTSNE(merged, dims = 1:30, reduction = "pca", reduction.name = "tsne.unintegrated")
- DimPlot(merged, reduction = "umap.unintegrated", group.by = "orig.ident") + ggtitle("umap Before Integration")
- DimPlot(merged, reduction = "tsne.unintegrated", group.by = "orig.ident") +
- ggtitle("t-SNE Before Integration")
- merged <- IntegrateLayers(
- object = merged, method = CCAIntegration,
- orig.reduction = "pca", new.reduction = "integrated.cca",
- verbose = FALSE
- )
- merged <- FindNeighbors(merged, reduction = "integrated.cca", dims = 1:20)
- merged <- FindClusters(merged, resolution = 0.2, cluster.name = "cca_clusters")
- merged <- RunUMAP(merged, reduction = "integrated.cca", dims = 1:20, reduction.name = "umap.cca")
- DimPlot(merged, reduction = "umap.cca", group.by = "orig.ident") + ggtitle("umap After Integration")
- merged <- RunTSNE(merged, reduction = "integrated.cca", dims = 1:20, reduction.name = "tsne.cca")
- DimPlot(merged, reduction = "tsne.cca", group.by = "orig.ident") + ggtitle("tsne After Integration")
- # Create a vector of different resolutions
- #resolutions <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.5)
- # Loop through resolutions and find clusters at each level
- #for (res in resolutions) {
- #merged <- FindClusters(merged, resolution = res, cluster.name = paste0("res_", res))
- #}
- #library(patchwork) # For combining plots
- # Generate UMAP plots for different resolutions
- #plots <- lapply(resolutions, function(res) {
- #DimPlot(merged, reduction = "umap.cca", group.by = paste0("res_", res)) + ggtitle(paste("Resolution:", res))
- #})
- # Combine all plots
- #wrap_plots(plots, ncol = 2)
- #table(merged$seurat_clusters) # Default clustering (last applied resolution)
- #for (res in resolutions) {
- #print(table([email hidden][[paste0("res_", res)]]))
- #}
- DimPlot(merged, reduction = "umap.cca", label = TRUE)
- DimPlot(merged, reduction = "tsne.cca", label = TRUE)
- merged<-JoinLayers(merged)
- [email hidden]
- Idents(merged) <- "cca_clusters"
- # colors selection
- cluster_colors <- c(
- "#1f78b4", "#33a02c", "#e31a1c", "#ff7f00", "#6a3d9a",
- "#b15928", "#b2df8a", "#fb9a99", "#fdbf6f", "#cab2d6",
- "#a6cee3", "#ffff99", "#8c510a", "#d73027", "#4575b4",
- "#f46d43", "#74add1", "#d53e4f"
- )
- # UMAP for WT (Keeping Cluster Colors)
- p1 <- DimPlot(
- merged,
- reduction = "tsne.cca",
- group.by = "cca_clusters",
- cells = WhichCells(merged, expression = orig.ident == "WT"),
- cols = cluster_colors,
- pt.size = 0.1
- ) +
- ggtitle("WT") +
- theme_minimal() +
- theme(
- plot.title = element_text(size = 18, face = "bold"),
- legend.position = "right",
- axis.text = element_text(size = 14),
- axis.title = element_text(size = 16)
- )
- # UMAP for PINK1 KO (Keeping Cluster Colors)
- p2 <- DimPlot(
- merged,
- reduction = "tsne.cca",
- group.by = "cca_clusters",
- cells = WhichCells(merged, expression = orig.ident == "PINK1"),
- cols = cluster_colors,
- pt.size = 0.1
- ) +
- ggtitle("PINK1 KO") +
- theme_minimal() +
- theme(
- plot.title = element_text(size = 18, face = "bold"),
- legend.position = "right",
- axis.text = element_text(size = 14),
- axis.title = element_text(size = 16)
- )
- # Combining both plots side by side
- final_plot <- p1 + p2
- table(merged$cca_clusters, merged$orig.ident)
- DEGs <- FindAllMarkers(merged, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
- View(DEGs)
- # Count total cells per condition (WT & PINK1 KO)
- total_cells <- table(merged$orig.ident)
- table(merged$orig.ident)
- # Count the number of cells in each cluster per condition
- cluster_counts <- as.data.frame(table(merged$cca_clusters, merged$orig.ident))
- # Rename columns
- colnames(cluster_counts) <- c("cca_clusters", "dataset", "cluster_counts")
- # Calculate relative abundance (% of total cells in each condition)
- cluster_counts <- cluster_counts %>%
- group_by(dataset) %>%
- mutate(Relative_Abundance = (cluster_counts / sum(cluster_counts)) * 100)
- head(cluster_counts)
- # Plotting the relative abundance of each cluster in WT vs. PINK1 KO
- ggplot(cluster_counts, aes(x = cca_clusters, y = Relative_Abundance, fill = dataset)) +
- geom_bar(stat = "identity", position = "dodge", width = 0.7) +
- theme_minimal() +
- labs(title = "",
- x = "Cluster",
- y = "Relative Abundance (%)") +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1, size = 14, face = "bold"),
- axis.text.y = element_text(size = 12, face = "bold"),
- axis.title = element_text(size = 14, face = "bold"),
- plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
- legend.title = element_text(size = 14, face = "bold"),
- legend.text = element_text(size = 12)
- ) +
- scale_fill_manual(values = c("WT" = "#1f78b4", "PINK1" = "#e31a1c"))
- VlnPlot(merged, features = c("Pink1", "Map2", "Syt1"), raster = FALSE, pt.size = 0, group.by = "cca_clusters")
- top_DEGs <- DEGs %>%
- group_by(cluster) %>%
- top_n(n = 3, wt = avg_log2FC) %>%
- pull(gene)
- print(top_DEGs)
- View(DEGs)
- # color palette
- custom_colors <- c("#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854",
- "#ffd92f", "#e5c494", "#b3b3b3", "#1f78b4", "#33a02c")
- # dot plot
- dot_plot <- DotPlot(merged, features = top_DEGs, cols = custom_colors, dot.scale = 6) +
- theme_minimal() +
- labs(title = "Top DEGs Across Clusters", x = "Clusters", y = "Genes") +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1, size = 9, face = "bold"),
- axis.text.y = element_text(size = 9, face = "bold"),
- plot.title = element_text(size = 14, face = "bold"),
- legend.position = "right"
- )
- print(dot_plot)
- ############ cell type assingnment
- merged$cell_type <- as.character(Idents(merged))
- merged$cell_type[merged$seurat_clusters == 0] <- "Stem cell"
- merged$cell_type[merged$seurat_clusters == 1] <- "Stem cell"
- merged$cell_type[merged$seurat_clusters == 2] <- "Goblet cells"
- merged$cell_type[merged$seurat_clusters == 3] <- "Goblet cells"
- merged$cell_type[merged$seurat_clusters == 4] <- "Colonocytes"
- merged$cell_type[merged$seurat_clusters == 5] <- "Myofibroblasts"
- merged$cell_type[merged$seurat_clusters == 6] <- "Pericytes"
- merged$cell_type[merged$seurat_clusters == 7] <- "Stem cell"
- merged$cell_type[merged$seurat_clusters == 8] <- "Intestinal epithelial cells"
- merged$cell_type[merged$seurat_clusters == 9] <- "Immune cells"
- merged$cell_type[merged$seurat_clusters == 10] <- "Immune cells"
- merged$cell_type[merged$seurat_clusters == 11] <- "Colonocytes"
- merged$cell_type[merged$seurat_clusters == 12] <- "Enteroendocrine cells"
- merged$cell_type[merged$seurat_clusters == 13] <- "Endothelial cells"
- merged$cell_type[merged$seurat_clusters == 14] <- "Endothelial cells"
- merged$cell_type[merged$seurat_clusters == 15] <- "Schwann cells"
- merged$cell_type[merged$seurat_clusters == 16] <- "Stem cell"
- merged$cell_type[merged$seurat_clusters == 17] <- "Schwann cells"
- Idents(merged) <- merged$cell_type
- my_colors <- c(
- "#E64B35", "#4DBBD5", "#00A087", "#3C5488", "#F39B7F", "#8491B4", "#91D1C2",
- "#DC0000", "#7E6148", "#B09C85", "#FFDC91", "#00A6D6", "#FF61CC", "#1B9E77",
- "#D95F02", "#7570B3", "#E7298A", "#66A61E", "#E6AB02", "#A6761D", "#A6CEE3",
- "#FB9A99", "#FDBF6F", "#CAB2D6", "#B3DE69", "#BC80BD", "#9E9AC8",
- "#8DD3C7", "#FCCDE5", "#D9D9D9", "#FFB3AB", "#A1D99B", "#B3CDE3", "#C49C94",
- "#F2B447", "#C7CEEA", "#FF6F91", "#88CCEE", "#DDCC77", "#A594F9" # NEW FINAL COLOR
- )
- umap_plot <- DimPlot(merged, reduction = "tsne.cca", label = FALSE, cols = cluster_colors) +
- theme_classic(base_size = 14) +
- theme(
- axis.text = element_blank(),
- axis.ticks = element_blank(),
- axis.title = element_blank(),
- legend.title = element_text(size = 12, face = "bold"),
- legend.text = element_text(size = 10)
- ) +
- ggtitle("")
- ggsave("wt_vs_pink1.png", plot = umap_plot, width = 5.1, height = 3.1, units = "in", dpi = 300, bg = "transparent")
- # Define Ordered Marker Genes
- features_ordered <- c(
- "Fcgbp", "Zg16", "Muc2", # Goblet Cells
- "Ptprc", "Lyz2", # Myeloid Immune Cells
- "Cd19", "Ms4a1", # Lymphoid Immune Cells
- "Plp1", "Ngfr", "Sox10", # Schwann Cells
- "Pecam1", "Lyve1", "Flt4", # Endothelial Cells
- "Chga", "Chgb", "Tph1", "Neurod1", # Enteroendocrine Cells
- "Col1a1", "Pdgfra", "Acta2", "Vim", # Fibroblasts
- "Mki67", "Pcna", # Proliferating Cells
- "Prom1", # Adult Stem Cell
- "Tfrc", "Clock", "Ccnd2", "Fut9", # Stem Cells
- "mt-Nd4", "Cox6a1", "Rplp2", # Mitochondria-Rich Cells
- "Tex9", "Cluh", # Migrating Cells
- "Slc26a3", "Slc15a1", "Epcam" # Colonocytes/Enterocytes
- )
- library(Seurat)
- library(ggplot2)
- # Define ordered marker genes
- features_ordered <- c(
- "Chga", "Neurod1", # Enteroendocrine
- "Lgr5", "Axin2", # Stem Cells
- "Prom1", "Car1", # Colonocytes
- "Slc15a1", "Epcam", # Enterocytes
- "Col1a1", "Pdgfra", # Fibroblasts
- "Zg16", "Muc2", # Goblet Cells
- "Pecam1", "Lyve1", # Endothelial
- "Ptprc", "Lyz2", # Immune Cells
- "Acta2", "Vim", # Myofibroblasts
- "Plp1", "Ngfr" # Schwann Cells
- )
- # Create DotPlot
- dot_plot <- DotPlot(merged, features = features_ordered) +
- RotatedAxis() +
- ggtitle("Marker Genes Across Gut Cell Types") +
- scale_color_gradientn(
- colors = c("#e0ecf4", "#9ecae1", "#fdae61", "#f46d43", "#a50026"), # Blue to red custom palette
- name = "Expression"
- ) +
- scale_size(range = c(1, 6), name = "% Cells") + # Adjust dot sizes
- theme_classic(base_size = 14) +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1, size = 10, face = "bold"),
- axis.text.y = element_text(size = 10, face = "bold"),
- plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
- panel.border = element_rect(color = "black", fill = NA, size = 1) # Add box line
- )
- # UMAP for WT
- p1 <- DimPlot(
- merged,
- reduction = "tsne.cca",
- cells = WhichCells(merged, expression = orig.ident == "WT"),
- cols = cluster_colors,
- pt.size = 0.1
- ) +
- ggtitle("WT") +
- theme_classic(base_size = 14) + # Clean base theme
- theme(
- plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
- axis.text = element_text(size = 12),
- axis.title = element_text(size = 14),
- legend.position = "none", # Hide legend
- panel.border = element_rect(color = "black", fill = NA, size = 1), # Add box line
- panel.grid = element_blank() # Remove grid lines
- )
- # UMAP for PINK1 KO
- p2 <- DimPlot(
- merged,
- reduction = "tsne.cca",
- cells = WhichCells(merged, expression = orig.ident == "PINK1"),
- cols = cluster_colors,
- pt.size = 0.1
- ) +
- ggtitle("PINK1 KO") +
- theme_classic(base_size = 14) +
- theme(
- plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
- axis.text = element_text(size = 12),
- axis.title = element_text(size = 14),
- legend.position = "none",
- panel.border = element_rect(color = "black", fill = NA, size = 1),
- panel.grid = element_blank()
- )
- # Combine plots
- library(patchwork)
- final_plot <- p1 + p2
- final_plot
- ggsave("wt_vs_pink1_diff.png", plot = final_plot, width = 5.1, height = 2.8, units = "in", dpi = 300, bg = "transparent")
- library(dplyr)
- library(ggplot2)
- # Creating a data frame of counts
- metadata_df <- as.data.frame([email hidden])
- summary_df <- metadata_df %>%
- group_by(orig.ident, cell_type) %>%
- summarize(count = n(), .groups = "drop")
- # Calculating proportions
- summary_df <- summary_df %>%
- group_by(orig.ident) %>%
- mutate(proportion = count / sum(count))
- # colors for cell types
- my_colors <- cluster_colors
- # Stacked bar plot
- ggplot(summary_df, aes(x = orig.ident, y = proportion, fill = cell_type)) +
- geom_col(color = "black", size = 0.3) + # Thin black outline
- scale_y_continuous(labels = scales::percent_format()) +
- scale_fill_manual(values = my_colors) +
- labs(
- title = "",
- x = "Condition",
- y = "Cell Type Proportion",
- fill = "Cell Type"
- ) +
- theme_classic(base_size = 14) +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1, size = 12, face = "bold"),
- axis.text.y = element_text(size = 12),
- plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
- panel.border = element_rect(color = "black", fill = NA, size = 1) # Add box line
- )
- #####################################label transfer validation######################
- reference <- mouse_colon_male
- query <- wt_pink1ko
- reference <- NormalizeData(reference)
- reference <- FindVariableFeatures(reference)
- query <- NormalizeData(query)
- query <- FindVariableFeatures(query)
- anchors <- FindTransferAnchors(
- reference = reference,
- query = query,
- normalization.method = "LogNormalize",
- reduction = "pcaproject",
- reference.reduction = "pca"
- )
- reference <- RunUMAP(reference, reduction = "pca", dims = 1:20, return.model = T)
- reference[["umap.new"]] <- CreateDimReducObject(
- embeddings = reference[["umap"]]@cell.embeddings,
- key = "UMAPnew_", # Set the key with an underscore
- assay = "RNA"
- )
- Validation <- MapQuery(
- anchorset = anchors,
- query = query,
- reference = reference,
- refdata = list(celltype = reference$Annotation),
- reference.reduction = "pca",
- reduction.model = "umap"
- )
- colors <- c(
- "#e6194B", "#3cb44b", "#ffe119", "#4363d8", "#f58231", "#911eb4",
- "#46f0f0", "#f032e6", "#bcf60c", "#fabebe", "#008080", "#e6beff",
- "#9A6324", "#fffac8", "#800000", "#aaffc3", "#808000", "#ffd8b1",
- "#0000FF", "#FF1493", "#00FF7F", "#FF4500", "#1E90FF", "#FF00FF",
- "#ADFF2F", "#FF6347", "#20B2AA"
- )
- table(Validation$predicted.celltype) # See label distribution
- p1 <- DimPlot(Validation, group.by = "cell_type", reduction = "tsne.cca", label = FALSE, cols = colors) +
- ggtitle("Actual cell type")
- p2 <- DimPlot(Validation, group.by = "predicted.celltype", reduction = "tsne.cca", label = FALSE, cols = colors) +
- ggtitle("Predicted cell type")
- p1 + p2
- ########### validation dataet to show predictic score distrtibution in density plot
- library(ggplot2)
- library(viridis)
- # Create a data frame of all predicted cell type scores
- scores <- [email hidden]$predicted.celltype.score
- cell_types <- [email hidden]$predicted.celltype
- data <- data.frame(cell_type = cell_types, score = scores)
- # Set the y-axis scale limits to be between 0 and the maximum density value
- max_density <- max(density(data$score)$y)
- ylim_max <- max_density + max_density*0.2 # Set a little extra space at the top
- ggplot(data, aes(x = score, fill = cell_type)) +
- geom_density(alpha = 0.5) +
- xlab("Predicted cell type score") +
- ylab("Density") +
- ggtitle("Distribution of predicted cell type scores") +
- ylim(0, ylim_max)
- # Generate a slightly darker palette
- num_types <- length(unique(data$cell_type))
- palette_colors <- viridis(num_types, option = "C")
- # Create a histogram of all predicted cell type scores
- validation_score <- ggplot(data, aes(x = score, fill = cell_type)) +
- geom_histogram(alpha = 0.5, position = "identity", bins = 30) +
- xlab("Predicted cell type score") +
- ylab("Count") +
- scale_fill_manual(values = palette_colors) +
- labs(
- title = "Distribution of Predicted Cell Type Scores",
- x = "Predicted Cell Type Score",
- y = "Cell Count",
- fill = "Predicted Cell Type"
- ) +
- theme_classic(base_size = 14) + # No background, no grid
- theme(
- plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
- axis.title = element_text(face = "bold"),
- axis.text = element_text(color = "black"),
- legend.title = element_text(face = "bold"),
- legend.key.size = unit(0.4, "cm"),
- legend.text = element_text(size = 10),
- legend.position = "right",
- panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.7) # Add outer box
- )
- ggsave("validation_score.png", plot = validation_score, width = 6.19, height = 3.42, units = "in", dpi = 300, bg = "transparent")
- ggplot(data, aes(x = score, fill = cell_type)) +
- geom_histogram(bins = 30, color = "black", fill = "#FFB3BA") +
- facet_wrap(~ cell_type, scales = "free_y", ncol = 5) +
- theme_minimal() +
- labs(
- title = "Predicted Cell Type Score Distribution per Cell Type",
- x = "Predicted Score", y = "Cell Count"
- ) +
- theme(strip.text = element_text(face = "bold"))
- table(Validation$predicted.celltype)
- table(Validation_high$predicted.celltype)
- Validation_high <- subset(Validation, subset = predicted.celltype.score > 0.5)
- table(Validation$predicted.celltype.score > 0.5)
- DimPlot(Validation_high, group.by = "predicted.celltype", reduction = "tsne.cca", label = FALSE, cols = colors) +
- ggtitle("Predicted cell type")
- p1 <- DimPlot(Validation_high, group.by = "cell_type", reduction = "tsne.cca", label = FALSE, cols = colors) +
- ggtitle("Actual cell type")
- p2 <- DimPlot(Validation_high, group.by = "predicted.celltype", reduction = "tsne.cca", label = TRUE, cols = colors) +
- ggtitle("Predicted cell type")
- p1 + p2
- table(Validation_high$predicted.celltype)
- table(Validation$predicted.celltype)
- library(randomForest)
- library(caret)
- library(pheatmap)
- library(ggplot2)
- library(multiROC)
- #1. data
- expr_mat <- GetAssayData(Validation_high, slot = "data") %>% as.matrix() %>% t()
- labels <- as.factor([email hidden]$predicted.celltype) # Use predicted label
- num_features <- 2000
- var_genes <- FindVariableFeatures(Validation_high, nfeatures = num_features) %>% VariableFeatures()
- expr_mat_reduced <- expr_mat[, var_genes]
- # 2. Build ML dataframe
- ml_data <- as.data.frame(expr_mat_reduced)
- ml_data$cell_type <- labels
- colnames(ml_data) <- make.names(colnames(ml_data))
- # Train/Test Split
- set.seed(123)
- trainIndex <- createDataPartition(ml_data$cell_type, p = 0.7, list = FALSE)
- train_data <- ml_data[trainIndex, ]
- test_data <- ml_data[-trainIndex, ]
- # Train RF Classifier
- set.seed(123)
- rf_model <- randomForest(cell_type ~ ., data = train_data, ntree = 500, importance = TRUE)
- # Predict on test set
- test_predictions <- predict(rf_model, newdata = test_data)
- test_probs <- predict(rf_model, newdata = test_data, type = "prob")
- # ----------------------------
- # Confusion Matrix & Enhanced Visualization
- # ----------------------------
- conf_mat <- confusionMatrix(test_predictions, test_data$cell_type)
- print(conf_mat)
- # Convert to matrix for heatmap
- cm <- conf_mat$table
- cm_df <- as.data.frame(cm)
- colnames(cm_df) <- c("Prediction", "Reference", "Freq")
- cm_df$Freq_log <- log1p(cm_df$Freq)
- # ggplot2 heatmap (log scale for low counts)
- ggplot(cm_df, aes(x = Reference, y = Prediction, fill = Freq_log)) +
- geom_tile(color = "grey90", linewidth = 0.3) +
- geom_text(aes(label = Freq), size = 3.5, color = "black") +
- scale_fill_gradientn(colors = c("white", "lightyellow", "orange", "red", "darkred"),
- name = "log(Freq+1)",
- guide = guide_colorbar(barwidth = 10, barheight = 1)) +
- labs(title = "Random Forest Confusion Matrix",
- x = "Actual Cell Type", y = "Predicted Cell Type") +
- theme_minimal(base_size = 13) +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1),
- plot.title = element_text(hjust = 0.5, face = "bold"),
- legend.position = "top"
- )
- # ----------------------------
- # Filter Confusion Matrix: remove all-zero rows/columns
- # ----------------------------
- # Extract raw confusion matrix
- cm <- conf_mat$table
- # Keep only classes that appear in both predicted and actual
- keep_classes <- intersect(rownames(cm)[rowSums(cm) > 0],
- colnames(cm)[colSums(cm) > 0])
- cm_filtered <- cm[keep_classes, keep_classes]
- # Convert to dataframe for ggplot
- cm_df <- as.data.frame(cm_filtered)
- colnames(cm_df) <- c("Prediction", "Reference", "Freq")
- cm_df$Freq_log <- log1p(cm_df$Freq)
- ggplot(cm_df, aes(x = Reference, y = Prediction, fill = Freq_log)) +
- geom_tile(color = "grey90", linewidth = 0.3) +
- geom_text(aes(label = Freq), size = 3.5, color = "black") +
- scale_fill_gradientn(colors = c("white", "grey90", "lightyellow","red", "darkred"),
- name = "log(Freq+1)",
- guide = guide_colorbar(barwidth = 10, barheight = 1)) +
- labs(title = "Random Forest Confusion Matrix",
- x = "Actual Cell Type", y = "Predicted Cell Type") +
- theme_minimal(base_size = 13) +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1),
- plot.title = element_text(hjust = 0.5, face = "bold"),
- legend.position = "top"
- )
- # Filtered cm_df is already prepared
- cm_df$Freq_log <- log1p(cm_df$Freq)
- random <- ggplot(cm_df, aes(x = Reference, y = Prediction, fill = Freq_log)) +
- # Draw tiles only where Freq > 0
- geom_tile(data = subset(cm_df, Freq > 0), color = "grey90", linewidth = 0.3) +
- # Add text labels only for nonzero cells
- geom_text(data = subset(cm_df, Freq > 0),
- aes(label = Freq), size = 3.5, color = "black") +
- scale_fill_gradientn(colors = c("white", "grey90", "lightyellow","red", "darkred"),
- name = "log(Freq+1)",
- guide = guide_colorbar(barwidth = 10, barheight = 1),
- na.value = "white") + # Make zero cells blank
- labs(title = "Random Forest Confusion Matrix",
- x = "Actual Cell Type", y = "Predicted Cell Type") +
- theme_minimal(base_size = 13) +
- theme(
- axis.text.x = element_text(angle = 45, hjust = 1),
- plot.title = element_text(hjust = 0.5, face = "bold"),
- legend.position = "top"
- )
- # assumes: cm_df with columns Reference, Prediction, Freq
- cm_df$Freq_log <- log1p(cm_df$Freq)
- random <- ggplot(subset(cm_df, Freq > 0),
- aes(x = Reference, y = Prediction, fill = Freq_log)) +
- geom_tile(color = "grey90", linewidth = 0.3) +
- scale_fill_gradientn(
- colors = c("white", "grey90", "lightyellow","red", "darkred"),
- name = "log(Freq+1)",
- guide = guide_colorbar(barwidth = 10, barheight = 1),
- na.value = "white"
- ) +
- coord_fixed(expand = FALSE) +
- labs(
- title = "Random Forest Confusion Matrix",
- x = "Actual Cell Type", y = "Predicted Cell Type"
- ) +
- theme_classic(base_size = 13) +
- theme(
- legend.position = "top",
- plot.title = element_text(hjust = 0.5, face = "bold"),
- axis.text.x = element_text(angle = 45, hjust = 1),
- panel.border = element_rect(color = "black", fill = NA, linewidth = 0.7) # boxed frame
- # no panel.grid in theme_classic(), so gridlines are gone
- )
- random
- ggsave("rf.png", plot = random, width = 6.00, height = 6.13, units = "in", dpi = 300, bg = "transparent")
- getwd()
- # Extract confusion matrix
- cm <- conf_mat$table
- # Calculate per-class accuracy (Recall)
- class_accuracy <- diag(cm) / rowSums(cm)
- # Remove NA (where row sum = 0, i.e., no true cells of that type)
- class_accuracy <- class_accuracy[!is.na(class_accuracy)]
- # Make dataframe
- acc_df <- data.frame(
- CellType = names(class_accuracy),
- Accuracy = class_accuracy
- )
- # Plot as horizontal barplot
- library(ggplot2)
- ggplot(acc_df, aes(x = reorder(CellType, Accuracy), y = Accuracy)) +
- geom_col(fill = "#9ABB50") +
- geom_text(aes(label = sprintf("%.2f", Accuracy)),
- hjust = -0.2, size = 3.5, color = "black") +
- scale_y_continuous(limits = c(0, 1.05), expand = c(0,0)) +
- coord_flip() +
- labs(title = "Per-Class prediction accuracy (random forest)",
- x = "Cell Type", y = "Accuracy (recall)") +
- theme_minimal(base_size = 14) +
- theme(
- axis.text.x = element_text(size = 12),
- axis.text.y = element_text(size = 12),
- plot.title = element_text(hjust = 0.5, face = "bold")
- )
- library(ggplot2)
- library(ggplot2)
- ggplot(acc_df, aes(x = reorder(CellType, Accuracy), y = Accuracy)) +
- geom_col(fill = "#9ABB50", width = 0.6) + # Sleek bar width
- scale_y_continuous(limits = c(0, 1.05), expand = c(0,0)) +
- coord_flip() +
- labs(
- title = "Per-Class Prediction Accuracy",
- x = NULL,
- y = "Accuracy (Recall)"
- ) +
- theme_classic(base_size = 10) +
- theme(
- axis.text.x = element_text(size = 10, color = "black"),
- axis.text.y = element_text(size = 10, color = "black"),
- axis.ticks.y = element_blank(),
- plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
- panel.grid.major.y = element_blank(),
- panel.grid.minor = element_blank()
- )
- library(dplyr)
- library(ggplot2)
- # Extract importance from trained RF
- imp_mat <- randomForest::importance(rf_model)
- imp_df <- data.frame(Gene = rownames(imp_mat), imp_mat, row.names = NULL)
- ## ---- Top 20 by MeanDecreaseAccuracy (MDA) ----
- top20_mda <- imp_df %>%
- filter(!is.na(MeanDecreaseAccuracy)) %>%
- arrange(desc(MeanDecreaseAccuracy)) %>%
- slice_head(n = 20)
- p_mda <- ggplot(top20_mda,
- aes(x = MeanDecreaseAccuracy,
- y = reorder(Gene, MeanDecreaseAccuracy))) +
- geom_point(size = 4, color = "#E69F00") + # orange dots
- labs(title = "Top 20 Features by Mean Decrease Accuracy",
- x = "Mean Decrease Accuracy", y = NULL) +
- theme_classic(base_size = 12) +
- theme(
- axis.text.y = element_text(size = 9),
- plot.title = element_text(hjust = 0.5, face = "bold"),
- panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8)
- )
- ggsave("rf_top20_MDA_dotplot.png", p_mda, width = 7, height = 6, dpi = 300, bg = "white")
- ## ---- Top 20 by MeanDecreaseGini (MDG) ----
- top20_mdg <- imp_df %>%
- filter(!is.na(MeanDecreaseGini)) %>%
- arrange(desc(MeanDecreaseGini)) %>%
- slice_head(n = 20)
- p_mdg <- ggplot(top20_mdg,
- aes(x = MeanDecreaseGini,
- y = reorder(Gene, MeanDecreaseGini))) +
- geom_point(size = 4, color = "#56B4E9") + # blue dots
- labs(title = "Top 20 Features by Mean Decrease Gini",
- x = "Mean Decrease Gini", y = NULL) +
- theme_classic(base_size = 12) +
- theme(
- axis.text.y = element_text(size = 9),
- plot.title = element_text(hjust = 0.5, face = "bold"),
- panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8)
- )
- ggsave("rf_top20_MDG_dotplot.png", p_mdg, width = 7, height = 6, dpi = 300, bg = "white")
- ########################
- # build a tiny data frame of diagonal positions (all class levels)
- lev <- levels(factor(cm_df$Reference, levels = unique(cm_df$Reference)))
- diag_df <- data.frame(
- Reference = factor(lev, levels = lev),
- Prediction = factor(lev, levels = lev)
- )
- diag_color <- "#1f1f1f" # <- change this to any accent (e.g., "#9C27B0", "#FF5722")
- random <- ggplot(subset(cm_df, Freq > 0),
- aes(x = Reference, y = Prediction, fill = log1p(Freq))) +
- geom_tile(color = "grey90", linewidth = 0.3) +
- # make the diagonal pop: thin white halo + colored stroke
- geom_tile(data = diag_df, fill = NA, color = "white", linewidth = 2.0) +
- geom_tile(data = diag_df, fill = NA, color = diag_color, linewidth = 1.1) +
- scale_fill_gradientn(
- colors = c("#FFFFFF", "#EEEEEE", "#F4F1D0", "#E57373", "#B71C1C"),
- name = "log(Freq+1)",
- guide = guide_colorbar(barwidth = 10, barheight = 1),
- na.value = "white"
- ) +
- coord_fixed(expand = FALSE) +
- labs(title = "Random Forest Confusion Matrix",
- x = "Actual Cell Type", y = "Predicted Cell Type") +
- theme_classic(base_size = 13) +
- theme(
- legend.position = "top",
- plot.title = element_text(hjust = 0.5, face = "bold"),
- axis.text.x = element_text(angle = 45, hjust = 1),
- panel.border = element_rect(color = "black", fill = NA, linewidth = 0.7)
- )
Code.R, under CC-BY-4.0 · at the source
Overview
- Department of Biochemistry & Molecular Biology, Ajou University School of Medicine,Suwon, 16499 South Korea
- Department of Biomedical Science, Graduate School of Ajou University,Suwon, 16499 South Korea
- BK21 R&E initiative for Advanced Precision Medicine, Ajou University School of Medicine,Suwon, 16499 South Korea
- Department of Brain Science, Ajou University School of Medicine,Suwon, 16499 South Korea
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 7 matches between paragraphs and lines of code.
figshare 30702425
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
- 29 September 2026: the link answers (HTTP 200)
1 file
- Code.R, R, 816 lines, 7 matches
Code availability statement
The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:
- it points to the authors' code: figshare 30702425
Read it in the paper: doi.org/10.1038/s41597-026-07193-4.
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;
- 1 script, each with its path and the digest of its content;
- 7 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
No dataset and no data link were found in the paper.
Data availability statement
The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:
- no repository, dataset or request procedure was recognized in it
Read it in the paper: doi.org/10.1038/s41597-026-07193-4.
Versions
The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.
Version 1, 29 September 2026: the first record
Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 9 MeSH terms, 2 funders, 21 references.
Cite
This paper
Junaid, M., Park, S. J., Bae, Y., Lee, E. J., & Bin Lim, S. (2026). A single-nucleus RNA-seq dataset of the colon in Pink1-deficient and wild-type mice. Scientific data, 13(1), 819. https://
BibTeX
@article{junaid2026singl
author = {Junaid, Muhammad and Park, Soo Jung and Bae, Yiseul and Lee, Eun Jeong and Bin Lim, Su},
title = {{A single-nucleus RNA-seq dataset of the colon in Pink1-deficient and wild-type mice}},
journal = {Scientific data},
year = {2026},
month = apr,
volume = {13},
number = {1},
pages = {819},
publisher = {Nature Publishing Group},
issn = {2052-4463},
doi = {10.1038/
url = {https://
pmid = {41942478},
pmcid = {PMC13230531}
}
RIS
TY - JOUR
AU - Junaid, Muhammad
AU - Park, Soo Jung
AU - Bae, Yiseul
AU - Lee, Eun Jeong
AU - Bin Lim, Su
TI - A single-nucleus RNA-seq dataset of the colon in Pink1-deficient and wild-type mice
T2 - Scientific data
J2 - Sci Data
PY - 2026
DA - 2026/
VL - 13
IS - 1
SP - 819
SN - 2052-4463
PB - Nature Publishing Group
DO - 10.1038/
UR - https://
LA - en
ER -
CSL-JSON
{
"id": "10.1038/
"type": "article-journal",
"title": "A single-nucleus RNA-seq dataset of the colon in Pink1-deficient and wild-type mice",
"container-title": "Scientific data",
"author": [
{
"family": "Junaid",
"given": "Muhammad"
},
{
"family": "Park",
"given": "Soo Jung"
},
{
"family": "Bae",
"given": "Yiseul"
},
{
"family": "Lee",
"given": "Eun Jeong"
},
{
"family": "Bin Lim",
"given": "Su"
}
],
"container-title-short":
"volume": "13",
"issue": "1",
"page": "819",
"DOI": "10.1038/
"PMID": "41942478",
"PMCID": "PMC13230531",
"ISSN": "2052-4463",
"publisher": "Nature Publishing Group",
"URL": "https://
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
6
]
]
}
}
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.7717/peerj.21426 [code]
- Integrated transcriptomic identification and validation reveal key autophagy-associated biomarkers in sleep deprivation.Journal: PeerJIn common: randomForest, caret, pheatmap, 4 other tools, genetics / omics, mouse
- [2] doi:10.3390/ijms27156925 [code]
- XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis.Journal: International journal of molecular sciencesIn common: randomForest, caret, pheatmap, 4 other tools, genetics / omics
- [3] doi:10.1038/s42003-026-10957-8 [code]
- Brain defence by the extracellular matrix protein Cochlin.Journal: Communications biologyIn common: randomForest, caret, pheatmap, 4 other tools, mouse
- [4] doi:10.1038/s41467-026-75723-0 [code]
- Spatial transcriptomics reveals distinct cell type dynamics following opioid dependence in female mice with the common human μ-opioid receptor variant Oprm1 A118G.Journal: Nature communicationsIn common: randomForest, pheatmap, Seurat, 3 other tools, genetics / omics, mouse
- [5] doi:10.1002/advs.76205 [code]
- CHCHD10 Mitigates Alzheimer's Disease-Related Phenotypes in Association With Epigenetic Remodeling in Directly Reprogrammed Neurons.Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)In common: randomForest, pheatmap, Seurat, 3 other tools, genetics / omics
- [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 sciencesIn common: caret, pheatmap, Seurat, 3 other tools, genetics / omics - [7] doi:10.1038/s41467-026-77170-3 [code]
- DNA methylation profiling identifies long-range epigenetic silencing of clustered protocadherins as a key determinant of meningioma progression.Journal: Nature communicationsIn common: randomForest, caret, pheatmap, 2 other tools, genetics / omics
- [8] doi:10.1016/j.celrep.2026.117298 [code]
- Midbrain endocannabinoids actuate dopamine-based action selection.Journal: Cell reportsIn common: randomForest, caret, patchwork, 2 other tools, mouse
- [9] doi:10.1038/s41593-026-02387-w [code]
- Developing mouse inhibitory neuron single-cell transcriptomes reveal distinct modes of cell-type diversification.Journal: Nature neuroscienceIn common: randomForest, Seurat, patchwork, 2 other tools, genetics / omics, mouse
- [10] doi:10.1016/j.nbd.2026.107379 [code]
- DYRK1A and Parkinson's disease, facts and hypotheses.Journal: Neurobiology of diseaseIn common: ggplot2, tidyverse, Parkinson's, 3 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.
Claim this paper
Correct its record
Say what each link of this record is, remove the ones that are not the paper's, add the ones that are missing. The correction becomes a new version of the record, in its Versions section.
Validate its tracing map
You validate the map as this page shows it: 1 repository of the authors' code, each at its verified commit and with its license, 1 script, and 7 matches between paragraphs and code (see the Code and Map sections). It then receives a DOI on Zenodo, with you (your ORCID iD) and OSCR as its creators; the code itself is not deposited.
The map's fingerprint: sha256:15f78efe6d59780d…
Add the badge to its README
The badge links the code to this page. Copy one of these into the README of the paper's code: only you decide where it goes, and nothing is changed for you.
Markdown
[.
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.
