# Load necessary libraries
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.3.3
## Loading required package: SeuratObject
## Warning: package 'SeuratObject' was built under R version 4.3.3
## Loading required package: sp
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
# Define a function to perform analysis on Seurat objects
analyze_seurat <- function(seurat_obj_path, experiment_name, metadata_path, display_name) {
# Print experiment name as a markdown heading
cat(sprintf("\n# %s (%s)\n", display_name, experiment_name))
cat("\n---\nAnalyzing experiment:", experiment_name, "\n")
# Load metadata
cat("Loading metadata...\n")
if (!file.exists(metadata_path)) {
stop("Metadata file not found: ", metadata_path)
}
metadata <- tryCatch({
read.csv(metadata_path, header = TRUE, stringsAsFactors = FALSE)
}, error = function(e) {
stop("Error reading metadata file: ", e$message)
})
colnames(metadata)[1] <- "Barcode"
# Subset metadata for the specific experiment
metadata_subset <- subset(metadata, Experiment == experiment_name)
cat("Total metadata entries for", experiment_name, ":", nrow(metadata_subset), "\n")
if (nrow(metadata_subset) == 0) {
stop("No metadata found for experiment:", experiment_name)
}
# Load the Seurat object
cat("Loading Seurat object from:", seurat_obj_path, "\n")
if (!file.exists(seurat_obj_path)) {
stop("Seurat object file not found: ", seurat_obj_path)
}
seurat_obj <- tryCatch({
readRDS(seurat_obj_path)
}, error = function(e) {
stop("Error reading Seurat object file: ", e$message)
})
# Extract the counts matrix
counts_matrix <- seurat_obj@assays$RNA@counts
counts_df <- as.data.frame(as.matrix(counts_matrix))
# Find common barcodes between the counts matrix and the metadata
common_barcodes <- intersect(colnames(counts_df), metadata_subset$Barcode)
cat("Number of common barcodes:", length(common_barcodes), "\n")
if (length(common_barcodes) == 0) {
stop("No common barcodes found between counts matrix and metadata for experiment:", experiment_name)
}
# Add common metadata to the Seurat object
metadata_common <- metadata_subset %>% filter(Barcode %in% common_barcodes)
metadata_common <- metadata_common[match(common_barcodes, metadata_common$Barcode), ]
for (col in colnames(metadata_common)[-1]) {
seurat_obj@meta.data[common_barcodes, col] <- metadata_common[, col]
}
cat("Barcodes in Seurat object after adding metadata:", nrow(seurat_obj@meta.data), "\n")
# Filter Seurat object to remove NAs in Subclass
if ("Subclass" %in% colnames(seurat_obj@meta.data)) {
cells_to_keep <- rownames(seurat_obj@meta.data[!is.na(seurat_obj@meta.data$Subclass), ])
cat("Cells to keep after filtering non-NA Subclass:", length(cells_to_keep), "\n")
if (length(cells_to_keep) == 0) {
stop("No cells with non-NA Subclass for experiment:", experiment_name)
}
seurat_obj <- subset(seurat_obj, cells = cells_to_keep)
} else {
stop("Subclass column not found in metadata for experiment:", experiment_name)
}
# Filter and normalize the Seurat object
cat("Performing log-normalization\n")
seurat_obj <- seurat_obj %>%
PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mito") %>%
subset(nFeature_RNA > 500 & nFeature_RNA < 8000 & percent.mito < 0.5) %>%
NormalizeData(normalization.method = "LogNormalize", scale.factor = 10000) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData(features = rownames(.)) %>%
RunPCA(features = VariableFeatures(object = .))
# Perform UMAP and clustering
cat("Running UMAP and clustering\n")
seurat_obj <- seurat_obj %>%
FindNeighbors(dims = 1:18) %>%
FindClusters(resolution = 0.5) %>%
RunUMAP(dims = 1:18)
# Plot UMAP
umap_plot <- DimPlot(seurat_obj, reduction = "umap", group.by = "Subclass", label = TRUE, label.size = 5, repel = TRUE, pt.size = 0.5) + NoLegend()
# Define genes of interest
genes_of_interest <- c("WPRE", "TC66T", "bGHpolyA", "oG", "Cre")
# Initialize lists to store plots
gene_plots_linear <- list()
gene_plots_log <- list()
# Perform analysis for each gene of interest
for (gene in genes_of_interest) {
if (gene %in% rownames(counts_matrix)) {
cat(sprintf("Performing %s analysis...\n", gene))
gene_counts <- counts_matrix[gene, ]
gene_df <- data.frame(Barcode = colnames(counts_matrix), Gene_reads = as.numeric(gene_counts))
seurat_obj@meta.data$Barcode <- colnames(seurat_obj)
metadata_with_gene <- merge(seurat_obj@meta.data, gene_df, by = "Barcode", all.x = TRUE)
subclasses_of_interest <- unique(metadata_subset$Subclass)
filtered_metadata <- metadata_with_gene %>% filter(Subclass %in% subclasses_of_interest)
# Plot gene distribution by subclass (linear scale)
gene_plot_linear <- ggplot(filtered_metadata, aes(x = Gene_reads, fill = Subclass)) +
geom_histogram(binwidth = 1, color = "black") +
facet_wrap(~ Subclass, scales = "free_y", nrow = 2, strip.position = "top") +
labs(title = paste(experiment_name, sprintf("Distribution of %s reads per barcode by subclass (linear scale)", gene)),
x = sprintf("Number of %s reads", gene), y = "Frequency") +
theme_minimal(base_size = 15)
# Plot gene distribution by subclass (log scale)
gene_plot_log <- ggplot(filtered_metadata, aes(x = Gene_reads, fill = Subclass)) +
geom_histogram(binwidth = 1, color = "black") +
scale_x_log10() +
facet_wrap(~ Subclass, scales = "free_y", nrow = 2, strip.position = "top") +
labs(title = paste(experiment_name, sprintf("Distribution of %s reads per barcode by subclass (log scale)", gene)),
x = sprintf("Number of %s reads (log scale)", gene), y = "Frequency") +
theme_minimal(base_size = 15)
gene_plots_linear[[gene]] <- gene_plot_linear
gene_plots_log[[gene]] <- gene_plot_log
} else {
cat(sprintf("%s gene not found in the dataset.\n", gene))
gene_plots_linear[[gene]] <- ggplot() + labs(title = sprintf("%s gene not found in the dataset.", gene))
gene_plots_log[[gene]] <- ggplot() + labs(title = sprintf("%s gene not found in the dataset.", gene))
}
}
list(umap_plot = umap_plot, gene_plots_linear = gene_plots_linear, gene_plots_log = gene_plots_log)
}
# Metadata path
metadata_path <- "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/patino2024_metadata/patino_metadata.csv"
# Paths to Seurat objects
seurat_paths <- list(
ST.S11 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_S11.rds",
ST.SC5 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_SC5.rds",
ST.T12 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T12.rds",
ST.T13 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T13.rds",
ST.NT6 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_NT6.rds"
)
# Experiment display names
experiment_display_names <- list(
ST.S11 = "Sepw1",
ST.SC5 = "Scnn1a",
ST.T12 = "Tlx3",
ST.T13 = "Tlx3",
ST.NT6 = "Ntsr1"
)
# Perform analysis for each Seurat object
for (experiment_name in names(seurat_paths)) {
seurat_obj_path <- seurat_paths[[experiment_name]]
display_name <- experiment_display_names[[experiment_name]]
results <- analyze_seurat(seurat_obj_path, experiment_name, metadata_path, display_name)
# Display UMAP plot
print(results$umap_plot)
for (gene in c("WPRE", "TC66T", "bGHpolyA", "oG", "Cre")) {
cat(sprintf("\n## %s (%s) - %s\n", display_name, experiment_name, gene))
# Display gene plots
print(results$gene_plots_linear[[gene]])
print(results$gene_plots_log[[gene]])
}
}
##
## # Sepw1 (ST.S11)
##
## ---
## Analyzing experiment: ST.S11
## Loading metadata...
## Total metadata entries for ST.S11 : 265
## Loading Seurat object from: /Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_S11.rds
## Number of common barcodes: 265
## Barcodes in Seurat object after adding metadata: 288
## Cells to keep after filtering non-NA Subclass: 265
## Performing log-normalization
## Centering and scaling data matrix
## PC_ 1
## Positive: Grm5, Robo2, Setbp1, Ppp2r2b, Clstn2, Ptprg, Cntn3, Nrg1, Fgf12, Ctnnd2
## Epha6, Xkr4, Lsamp, Cntn4, Cadm1, Rgs7, Grik2, Cntn1, Dync1i1, Zfp385b
## Fat3, Nol4, 5730522E02Rik, Dpp10, Dync1i2, Dab1, Rbfox1, Cntnap5a, Pdzrn3, Asic2
## Negative: 4930562C15Rik, Jhy, Cecr2, Or2h2c, Stk33, Adamtsl1, Matn2, Or10ac1, Pabpc1l, Pycard
## Col23a1, Ush2a, Stab2, Creb5, Tcaf3, Xrra1, Rnls, Tcaf2, Rprd2, Ccdc7a
## Fchsd2, Dnah12, Csf1r, Slc5a1, Jun, Or2h1, Prdm9, Ccdc162, Farp2, Il12rb2
## PC_ 2
## Positive: Hs6st3, Lingo2, Tafa1, Slc1a2, Nell2, Cacna1e, Kcnq5, Nrgn, Pde4b, Sv2b
## Oprm1, Khdrbs2, Slit3, Rorb, Nrn1, Arc, Mapk4, Ptn, Nrg1, Cracdl
## Ptprk, Cadps2, Nwd2, Pdzd2, Khdrbs3, Rgs6, Scube1, Kcnt2, Mlip, Brinp3
## Negative: Erbb4, Nxph1, Abtb3, Sox6, Zfp536, Kcnip1, Luzp2, Gad1, 6330411D24Rik, Cacna2d2
## Cntnap4, Rbms3, Eya4, Slc6a1, Grip1, Nek7, Cntnap3, Coro6, Tox3, Alk
## Gad2, Vwc2, Cntnap5c, Robo1, Kcns3, Gabrg3, Maf, Kcnmb2, Cemip, Ankrd55
## PC_ 3
## Positive: Fyb, Cx3cr1, Hk2, Fli1, Slfn8, Arhgap45, Mrgpra1, Ly86, Lpcat2, Insyn2b
## Maml2, Fcer1g, C1qc, Tgfbr1, Pik3r5, Blnk, Pik3ap1, Slfn5, Ptprc, Mertk
## Apbb1ip, Fam180a, Adap2, Dock8, Inpp5d, Ifi209, Spi1, Ifi204, Or2c1, Ifi214
## Negative: Apaf1, Mpp7, Polk, Abca6, Slc25a3, Ifrd1, Plod2, Usp29, Mdm2, Arhgap12
## Ddit3, Sh3tc2, Slc9b1, Faxdc2, Tpm1, Nhs, Rps15a, Rhoq, Tmc1, Abcb1b
## Grik1, Ccdc162, Tubb4b, Map3k19, Mt1, Vps36, Eya4, Odad2, Scg2, Abtb3
## PC_ 4
## Positive: Mro, Cfap61, Ear1, Stac, 4930579F01Rik, Pde7b, Veph1, Pebp4, Obscn, Ntsr1
## Cfap95, Pam, Zscan10, Tnfrsf9, Plekhg1, Aspa, Bmp3, Thegl, Or6x1, Adgrl3
## Adamts12, Dnah3, Tnr, Inhca, Awat1, Ttc21a, Catspere1, Greb1l, Lrcol1, Rbks
## Negative: Hk2, Ptprc, Insyn2b, Fcer1g, Dock8, Il10ra, Prkcd, Ifi209, Inpp5d, Fyb
## Ifi214, Ly86, Slfn8, Dock2, Plcg2, Blnk, Oasl1, Pik3r5, Pik3ap1, Skap2
## F11r, Apbb1ip, Cxcl16, Adap2, Ifi207, Ifi204, Bcl2a1b, Ccl3, Lpcat2, Tlr7
## PC_ 5
## Positive: Selplg, Stac, Ankrd36, Iqca, 4930579F01Rik, Cfap95, Trcg1, Pkp1, Lmntd1, Parm1
## Itgad, Esrrb, Pole, Batf, Il11, Il23r, Dnah14, Tgfb1, Il6ra, Obscn
## Ank1, Acod1, Myo7b, Crtam, Mamdc2, Ccdc136, Syt2, Ctf2, Endou, Abcc2
## Negative: Ptgds, Mag, Plp1, St18, Aspa, Ermn, Gatm, Cldn11, Fa2h, Gjc3
## Mobp, Prr5l, Car2, Tspan2, Ugt8a, Pstpip2, Megf10, Gpc5, Myrf, Sec14l5
## Dscaml1, Chrm3, Grb14, Mbp, Apod, Sox10, Bcas1, Fign, Fbxw15, Csrp1
## Running UMAP and clustering
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 256
## Number of edges: 7037
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7468
## Number of communities: 5
## Elapsed time: 0 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 22:32:57 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:32:57 Read 256 rows and found 18 numeric columns
## 22:32:57 Using Annoy for neighbor search, n_neighbors = 30
## 22:32:57 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:32:57 Writing NN index file to temp file /var/folders/m2/v2yg6yh92818c42pf0tfmbwr0000gn/T//Rtmp9fhWbf/file9ca874774486
## 22:32:57 Searching Annoy index using 1 thread, search_k = 3000
## 22:32:57 Annoy recall = 100%
## 22:32:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:32:57 Initializing from normalized Laplacian + noise (using RSpectra)
## 22:32:57 Commencing optimization for 500 epochs, with 8928 positive edges
## 22:32:58 Optimization finished
## Performing WPRE analysis...
## Performing TC66T analysis...
## Performing bGHpolyA analysis...
## Performing oG analysis...
## Performing Cre analysis...
##
## ## Sepw1 (ST.S11) - WPRE
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 106 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Sepw1 (ST.S11) - TC66T
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 250 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Sepw1 (ST.S11) - bGHpolyA
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 180 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Sepw1 (ST.S11) - oG
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 228 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Sepw1 (ST.S11) - Cre
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 254 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## # Scnn1a (ST.SC5)
##
## ---
## Analyzing experiment: ST.SC5
## Loading metadata...
## Total metadata entries for ST.SC5 : 3025
## Loading Seurat object from: /Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_SC5.rds
## Number of common barcodes: 3016
## Barcodes in Seurat object after adding metadata: 3171
## Cells to keep after filtering non-NA Subclass: 3016
## Performing log-normalization
## Centering and scaling data matrix
## PC_ 1
## Positive: Stk33, Ush2a, 4930562C15Rik, Jhy, Creb5, Stab2, Cecr2, Dnah12, Ccdc7a, Pebp4
## Adamtsl1, Pabpc1l, Il12rb2, Tcaf3, Hormad2, Col23a1, Or10ac1, Pld4, Jun, Pycard
## Veph1, Slc5a1, Myo3b, Dnah3, Nup210l, Crtam, Abtb2, Csf1r, Tmem132c, 1110017D15Rik
## Negative: Dpp10, Lrfn5, Nrg1, Grm5, 5730522E02Rik, Kcnq5, Brinp3, Sv2b, Khdrbs2, Dync1i1
## Nyap2, Pdzrn3, Fat3, Kcnb2, Rorb, Car10, Sgcz, Prr16, Grid2, Gria4
## Pde4b, Pbx1, Kcnip4, Grm3, Chsy3, Ncam2, Fut9, Cdh8, Slit3, Egfem1
## PC_ 2
## Positive: Rorb, Lingo2, Kcnh5, Nrg1, Tafa1, Prr16, Car10, Cadps2, Sorcs1, Mlip
## Pdzrn3, Egr1, Cdh6, Xirp2, Pde7b, Brinp3, Slc1a2, Pde4b, Prkg1, Kcnq5
## Trpm3, Mapk4, Marchf1, Arc, Junb, Cdh8, Ccdc7a, Khdrbs2, Dynll1, Rasl11b
## Negative: Erbb4, Nxph1, Sox6, Kcnip1, Gabrg3, Robo1, Zfp536, Luzp2, Rbms3, Abtb3
## Grip1, Slc6a1, Gad1, Kcnmb2, Vwc2, Cacna2d2, Eya4, Mgat4c, 6330411D24Rik, Zfp804a
## Cntnap4, Cntnap3, Kcns3, Cemip, Sox5, Afap1, Cntn6, Tenm1, Gad2, Ubash3b
## PC_ 3
## Positive: Scn9a, Abtb3, 6330411D24Rik, Vwc2, Gad1, Kcnip1, Kcnmb2, Erbb4, Kcnc1, Cacna2d2
## Pde5a, Sox6, Zfp536, Cemip, Eya4, Kcns3, Pvalb, Slc6a1, Cntnap3, Aldh1l2
## WPRE, Thsd7a, Ankrd55, Mybpc1, Grik1, Tafa2, Grip1, Rnf144b, Nxph1, Cntnap4
## Negative: Rgs6, Pde1a, Nfib, Nfia, Dgkb, Mertk, Grm8, Dpp10, Ldb2, Nkain2
## Khdrbs2, Meis2, Csgalnact1, Wdr17, 5730522E02Rik, Pex5l, Chrm3, Pde4d, Prex2, Slc1a3
## Lrrtm4, Slit3, Sgcz, Gpc5, Sv2b, Nrp1, Marchf1, Pdzrn4, Gpc6, Pde4b
## PC_ 4
## Positive: Rasgrf2, Pam, Ntrk3, Lingo2, Kctd1, Pde7b, Kcnip4, Slc4a4, Ccbe1, Pdzrn3
## Nectin3, Cdh7, Ifit1, Mgat5, Col19a1, Epb41l2, Calb1, Alcam, Pparg, Robo3
## Gsg1l, Gucy1a1, Stard8, Kctd8, Lrrtm4, Enpp2, Chrm3, Kcnc2, Mlip, Slc9a9
## Negative: Hs3st4, Hs3st2, Slc35f1, Il1rapl2, Kcnk2, Pex5l, Sdk2, Fras1, Celf4, Hs3st5
## Deptor, Grik3, Agbl4, Galnt17, Nalf1, Plcxd3, Entrep2, Etv1, Sulf2, Olfm3
## Foxp2, Gng12, Thsd7b, Cdh18, Rab3c, Esrrg, Prkcg, Pde1a, Galnt9, Htr1f
## PC_ 5
## Positive: Chrm3, Gria3, Rgs6, Lrrtm4, Ntrk3, Pde4d, Dgkb, Lrrtm3, Rbfox1, Pde1a
## Cntn4, Lhfpl3, Sorcs3, Tnr, Gria1, Grm1, Ccbe1, Nos1ap, Pebp4, Grm8
## Dscaml1, Enpp2, Edil3, Csgalnact1, Atp2b4, Me3, Gsg1l, Pam, Kcnip4, Ldb2
## Negative: Hk2, Ly86, Dock8, Inpp5d, Mertk, Fyb, Ikzf1, Dock2, Fgd2, Unc93b1
## Pik3ap1, Adap2, Ifi213, Ifi204, Arhgap45, Cx3cr1, Ctss, Tgfbr1, Epsti1, Sp100
## Cyth4, Cd37, Ifi207, Rbm47, Siglech, Ifi206, Plcg2, C1qb, Cd84, Apbb1ip
## Running UMAP and clustering
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 3012
## Number of edges: 94130
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8953
## Number of communities: 16
## Elapsed time: 0 seconds
## 22:33:09 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:33:09 Read 3012 rows and found 18 numeric columns
## 22:33:09 Using Annoy for neighbor search, n_neighbors = 30
## 22:33:09 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:33:09 Writing NN index file to temp file /var/folders/m2/v2yg6yh92818c42pf0tfmbwr0000gn/T//Rtmp9fhWbf/file9ca84233f526
## 22:33:09 Searching Annoy index using 1 thread, search_k = 3000
## 22:33:10 Annoy recall = 100%
## 22:33:10 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:33:10 Initializing from normalized Laplacian + noise (using RSpectra)
## 22:33:10 Commencing optimization for 500 epochs, with 120072 positive edges
## 22:33:12 Optimization finished
## Performing WPRE analysis...
## Performing TC66T analysis...
## Performing bGHpolyA analysis...
## Performing oG analysis...
## Performing Cre analysis...
##
## ## Scnn1a (ST.SC5) - WPRE
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 851 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Scnn1a (ST.SC5) - TC66T
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2965 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Scnn1a (ST.SC5) - bGHpolyA
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2167 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Scnn1a (ST.SC5) - oG
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2899 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Scnn1a (ST.SC5) - Cre
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2992 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## # Tlx3 (ST.T12)
##
## ---
## Analyzing experiment: ST.T12
## Loading metadata...
## Total metadata entries for ST.T12 : 2956
## Loading Seurat object from: /Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T12.rds
## Number of common barcodes: 2940
## Barcodes in Seurat object after adding metadata: 3035
## Cells to keep after filtering non-NA Subclass: 2940
## Performing log-normalization
## Centering and scaling data matrix
## PC_ 1
## Positive: 4930562C15Rik, Jhy, Stk33, Ush2a, Cecr2, Creb5, Or10ac1, Col23a1, Adamtsl1, Dnah12
## Matn2, Stab2, Or2h2c, Tcaf3, Slc5a1, Il12rb2, Pebp4, Jun, Pycard, Myo3b
## Tcaf2, Nup210l, Veph1, Rftn2, Prdm9, Xrra1, Pabpc1l, Crtam, Bicdl2, Tmem132c
## Negative: Lrfn5, Ctnna2, Dpp10, Cntn3, Nrg1, Lsamp, Kcnq5, Fat3, Brinp3, Fstl4
## Car10, Xkr4, Robo2, Grid2, Pbx1, Gria4, Pdzrn3, Sv2b, Dync1i1, Frmd4a
## 5730522E02Rik, Slc1a2, Lingo2, Fut9, Slc2a13, Kcnip4, Rorb, Prr16, Chsy3, Khdrbs2
## PC_ 2
## Positive: Erbb4, Nxph1, Sox6, Robo1, Gabrg3, Zfp536, Kcnip1, Abtb3, Luzp2, Rbms3
## 6330411D24Rik, Gad1, Grip1, Cacna2d2, Vwc2, Mgat4c, Sox5, Ano4, Eya4, Cntnap4
## Cntn4, Kcnmb2, Rbfox1, Slc6a1, Alk, Cntnap3, Zfp804a, Tenm1, Kcns3, Ubash3b
## Negative: Rorb, Nrg1, Kcnh5, Lingo2, Fstl4, Tafa1, Prr16, Arc, Cadps2, Mapk4
## Slc1a2, Pdzrn3, Nell1, Dlg1, Pdp1, Egr1, Car10, Dynll1, Kcnq5, Pde7b
## Chrm2, Ier2, Sorcs1, Scube1, Egr3, Plk2, Junb, Brinp3, Adam2, Pde4b
## PC_ 3
## Positive: Abtb3, Sox6, Gad1, Scn9a, Vwc2, Cacna2d2, Erbb4, 6330411D24Rik, Eya4, Kcnmb2
## Zfp536, Nxph1, Slc6a1, Cntnap5c, Afap1, Kcnip1, Nhs, Cemip, Grip1, Cntnap4
## Ankrd55, Lhx6, Frmpd1, Cntnap3, Gad2, Kcns3, Npas3, Nek7, Kcnc1, Slit2
## Negative: Dpp10, 5730522E02Rik, Khdrbs2, Rgs6, Pde1a, Sv2b, Pex5l, Slit3, Chsy3, Nfib
## Prkg1, Oprm1, Marchf1, Ldb2, Sgcz, Nkain2, Car10, Pde4b, Grm8, Thsd7b
## Mctp1, Slc1a2, Foxp2, Fyb, Nos1ap, Meis2, Egfem1, Kcnn2, Pde4d, Garnl3
## PC_ 4
## Positive: Kcnb2, Pdzrn3, Pde7b, Kcnc2, Prr16, Fstl4, Tmem108, Tafa1, Myo16, Prkg1
## Lingo2, Tafa2, Rorb, Sphkap, Frmd5, Pparg, Gria4, Lrfn5, Rcan2, Epha6
## Thsd7a, Kcnh5, Thrb, Brinp3, Dock4, Nrg1, Sorbs2, Slc2a13, Ptprk, Chrm2
## Negative: Foxp2, Grm8, Cdh18, Zfpm2, Grik3, Nfia, Gng12, Sdk2, Nrp1, Nfib
## Sema3e, Syt6, Hs3st2, Thsd7b, Bcl11b, Slc35f1, Hs3st4, Tle4, Nos1ap, Rai14
## Dlc1, Peak1, Htr1f, Hs3st5, Arhgap25, Mgat4c, Prkcg, A830018L16Rik, Pde1a, Pcsk5
## PC_ 5
## Positive: Gm6040, Obscn, Parm1, Ankrd36, Garnl3, Abcc2, Nalf1, Slc35f1, Chat, Tenm1
## Ttc21a, Hs3st4, Pde1a, Cdh18, Veph1, Lgr5, Pou6f2, Lmo3, Sdk2, Grik3
## Pex5l, Crtam, Mrgpra2b, Ms4a4b, Spink13, Gm5724, Foxp2, Oosp3, Vwc2l, Esrrg
## Negative: Hk2, Fyb, Dock8, Ly86, Fgd2, Mertk, Ikzf1, Arhgap45, Epsti1, Pik3ap1
## Inpp5d, Insyn2b, Bin2, Ifi204, Adap2, Ifi213, Unc93b1, Dock2, Sp100, Apobec1
## Ctss, C1qa, C1qc, Apobec3, Ifi207, Apbb1ip, Ptprc, Cd84, Lpcat2, Tgfbr2
## Running UMAP and clustering
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2927
## Number of edges: 93856
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8841
## Number of communities: 13
## Elapsed time: 0 seconds
## 22:33:26 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:33:26 Read 2927 rows and found 18 numeric columns
## 22:33:26 Using Annoy for neighbor search, n_neighbors = 30
## 22:33:26 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:33:26 Writing NN index file to temp file /var/folders/m2/v2yg6yh92818c42pf0tfmbwr0000gn/T//Rtmp9fhWbf/file9ca86123b3be
## 22:33:26 Searching Annoy index using 1 thread, search_k = 3000
## 22:33:27 Annoy recall = 100%
## 22:33:27 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:33:27 Initializing from normalized Laplacian + noise (using RSpectra)
## 22:33:27 Commencing optimization for 500 epochs, with 118518 positive edges
## 22:33:30 Optimization finished
## Performing WPRE analysis...
## Performing TC66T analysis...
## Performing bGHpolyA analysis...
## Performing oG analysis...
## Performing Cre analysis...
##
## ## Tlx3 (ST.T12) - WPRE
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 1118 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T12) - TC66T
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2882 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T12) - bGHpolyA
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2366 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T12) - oG
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2793 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T12) - Cre
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2922 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## # Tlx3 (ST.T13)
##
## ---
## Analyzing experiment: ST.T13
## Loading metadata...
## Total metadata entries for ST.T13 : 3272
## Loading Seurat object from: /Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T13.rds
## Number of common barcodes: 3244
## Barcodes in Seurat object after adding metadata: 3428
## Cells to keep after filtering non-NA Subclass: 3244
## Performing log-normalization
## Centering and scaling data matrix
## PC_ 1
## Positive: 4930562C15Rik, Jhy, Stab2, Stk33, Dnah12, Cecr2, Rftn2, Ush2a, Creb5, Matn2
## Jun, Il12rb2, Veph1, Crtam, Adamtsl1, Pebp4, Nup210l, Pycard, Col23a1, Xrra1
## Ccdc7a, Pabpc1l, Dnah3, Or10ac1, Ankrd36, Or2h2c, Abtb2, Myo3b, Ms4a4d, Ms4a6c
## Negative: Dpp10, Ctnna2, Cntn3, Nrg1, Lsamp, Brinp3, 5730522E02Rik, Xkr4, Sv2b, Grid2
## Fat3, Frmd4a, Fstl4, Slc2a13, Slc1a2, Gria4, Kcnj3, Robo2, Pdzrn3, Oprm1
## Rorb, Chsy3, Grik2, Kcnip4, Fut9, Cntn5, Nalf1, Cadm1, Khdrbs2, Pde4b
## PC_ 2
## Positive: Rorb, Kcnh5, Lingo2, Nrg1, Xirp2, Tafa1, Fstl4, Homer1, Mapk4, Slc1a2
## Prr16, Nell1, Pde7b, Cadps2, Pdzrn3, Brinp3, Pde4b, Chrm2, Nwd2, Tacc1
## Cdh8, Pamr1, Scube1, Egr1, Cpne8, Trpm3, Arc, Mlip, Tanc1, Plk2
## Negative: Erbb4, Sox6, Nxph1, Gabrg3, Abtb3, Robo1, Zfp536, Kcnip1, Slc6a1, 6330411D24Rik
## Rbms3, Luzp2, Vwc2, Gad1, Grip1, Cntnap4, Sox5, Cntnap3, Cacna2d2, Kcnmb2
## Kcns3, Mgat4c, Ano4, Gad2, Eya4, Cemip, Nek7, Alk, Zfp804a, Kazn
## PC_ 3
## Positive: Abtb3, Aldh1l2, Scn9a, Gad1, Zfp536, Kcnip1, Cntnap5c, Nxph1, Sox6, Cntnap3
## Vwc2, Cacna2d2, Nhs, Grip1, Kcnc1, Pde5a, Ankrd55, Apaf1, Slit2, Lhx6
## Slc6a1, Gad2, Slc32a1, Maf, Kcns3, Erbb4, Mpp7, Dynll2, Eya4, Nek7
## Negative: Rgs6, Dpp10, Pde1a, 5730522E02Rik, Slit3, Khdrbs2, Sv2b, Prkg1, Sgcz, Asic2
## Wdr17, Marchf1, Kcnip4, Csgalnact1, Oprm1, Ldb2, Nfib, Pde7b, Kcnn2, Dgkb
## Pex5l, Nkain2, Veph1, Chsy3, Lhfpl3, Pde4d, Fli1, Obscn, Dab1, Nos1ap
## PC_ 4
## Positive: Grm8, Foxp2, Nrp1, Grik3, Cdh18, Sema3e, Sdk2, Gng12, Zfpm2, Nfib
## Nfia, Hs3st2, Thsd7b, Syt6, Nos1ap, Pde1a, Rai14, Hs3st5, Pcsk5, Peak1
## Tle4, Bcl11b, Slc35f1, Prkcg, Grik4, Dlc1, Meis2, Pcp4, Igsf21, Sox5
## Negative: Pdzrn3, Pde7b, Kcnc2, Kcnb2, Tafa2, Prr16, Tmem108, Rorb, Fstl4, Thsd7a
## Prkg1, Frmd5, Grik1, Pparg, Gria4, Rcan2, Kcnh5, Nrg1, Chrm2, Erbb4
## 6330411D24Rik, Col19a1, Thrb, Lingo2, Slc2a13, Limch1, Dpp6, Brinp3, Rasgrf2, Sorbs2
## PC_ 5
## Positive: Mertk, Fyb, Hk2, Ly86, Dock8, Inpp5d, C1qc, Ifi213, Adap2, Pik3ap1
## Dock2, Tcf7l2, Epsti1, Fgd2, Insyn2b, Unc93b1, Ptprc, Arhgap45, Gpc5, Ikzf1
## Csf3r, Epb41l2, Runx1, Sp100, Cx3cr1, Ifi204, C1qb, C1qa, Ifi207, Lair1
## Negative: Veph1, Hs3st4, Obscn, Ankrd36, Garnl3, Gm6040, Slc35f1, Parm1, Tenm1, Ttc21a
## Cntn5, Galntl6, Il23r, 1700022I11Rik, Abcc2, Fras1, Mrgpra2b, Pou6f2, Srrm3, Cidec
## Iqca, Plcxd3, Oosp3, Cdh18, Tnfrsf10b, Acot12, Crtam, Galnt17, Lmo3, Kazn
## Running UMAP and clustering
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 3231
## Number of edges: 99634
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8889
## Number of communities: 13
## Elapsed time: 0 seconds
## 22:33:42 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:33:42 Read 3231 rows and found 18 numeric columns
## 22:33:42 Using Annoy for neighbor search, n_neighbors = 30
## 22:33:42 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:33:43 Writing NN index file to temp file /var/folders/m2/v2yg6yh92818c42pf0tfmbwr0000gn/T//Rtmp9fhWbf/file9ca86eedb0c4
## 22:33:43 Searching Annoy index using 1 thread, search_k = 3000
## 22:33:43 Annoy recall = 100%
## 22:33:43 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:33:44 Initializing from normalized Laplacian + noise (using RSpectra)
## 22:33:44 Commencing optimization for 500 epochs, with 128664 positive edges
## 22:33:46 Optimization finished
## Performing WPRE analysis...
## Performing TC66T analysis...
## Performing bGHpolyA analysis...
## Performing oG analysis...
## Performing Cre analysis...
##
## ## Tlx3 (ST.T13) - WPRE
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 1515 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T13) - TC66T
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 3192 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T13) - bGHpolyA
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 2370 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T13) - oG
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 3140 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Tlx3 (ST.T13) - Cre
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 3223 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## # Ntsr1 (ST.NT6)
##
## ---
## Analyzing experiment: ST.NT6
## Loading metadata...
## Total metadata entries for ST.NT6 : 1786
## Loading Seurat object from: /Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_NT6.rds
## Number of common barcodes: 1768
## Barcodes in Seurat object after adding metadata: 1854
## Cells to keep after filtering non-NA Subclass: 1768
## Performing log-normalization
## Centering and scaling data matrix
## PC_ 1
## Positive: Grm5, Ank2, Fut9, Lsamp, Pbx1, Xkr4, Lrfn5, Etl4, Grik2, Clstn2
## Robo2, A830018L16Rik, Cntn3, Cntn4, Sv2b, Lrrtm3, Khdrbs2, 5730522E02Rik, Chsy3, Mmp16
## Cntn5, Tmeff2, Frmd4a, Ptprg, Nrg1, Ifit1, Rbfox1, Flrt2, Slc1a2, Pex5l
## Negative: Jhy, 4930562C15Rik, Creb5, Stab2, Stk33, Col23a1, Cecr2, Ush2a, Nup210l, Adamtsl1
## Pycard, Pabpc1l, Il12rb2, Ccdc7a, Slc5a1, Abtb2, Or10ac1, Nlrp4c, Or2h2c, Jun
## Myo3b, Matn2, Pebp4, Xrra1, Dnah3, Dnah12, Crtam, Tcaf3, Rftn2, Hormad2
## PC_ 2
## Positive: Mgat4c, Robo1, Sox5, Erbb4, Cdh18, Tenm1, Ano4, Gabrg3, Dlc1, Cntn4
## Sox6, A830018L16Rik, Me3, Zfpm2, Nxph1, Zfp536, Bcl11b, Thsd7b, Rbfox1, Sh3rf3
## Kcnmb2, Rbms3, Vwc2, Slc35f1, Cntnap4, Foxp2, Gad1, Fmn1, Kcnip1, Slc35f3
## Negative: Cacna2d3, Hs6st3, Lingo2, Kcnh5, Rorb, Tafa1, Fstl4, Pdzrn3, Prr16, Mapk4
## Xirp2, Nell1, Cux2, Ptn, Tafa2, Nwd2, Pamr1, Kcnq5, Kcnb2, Sorcs1
## Sphkap, Cux1, Adam2, Cdh6, Serpine2, Ptprt, Nrg1, Cdh8, Scube1, Bdnf
## PC_ 3
## Positive: Kcnc2, Erbb4, Zfp536, Galntl6, Epha6, Sox6, Limch1, Dpp6, Vwc2, Hs6st3
## Pdzrn3, Abtb3, Maf, Ptprm, Frmd5, Kcnb2, Gad1, Cntnap5a, Myo16, Slc6a1
## Gad2, Kcns3, Gria4, Scn9a, Prr16, Luzp2, Eya4, Thrb, Cacna2d3, Grip1
## Negative: Foxp2, Zfpm2, Hs3st4, Thsd7b, Nfia, Cdh18, Grm8, Nrp1, Syt6, Nfib
## Grik3, Tle4, Rai14, Gng12, Slc35f1, Pde1a, Sdk2, Sema3e, Ephb1, Hs3st5
## Rgs6, Vxn, Grik4, Ifi203, Garnl3, Mctp1, Slc35f3, Igsf21, Nos1ap, Dlc1
## PC_ 4
## Positive: Pde7b, Prkg1, Kcnb2, Gm6040, Nrg1, Veph1, Lhfpl3, Fstl4, Tafa1, Rgs7
## Slit3, Cacna2d3, Sgpp2, Pebp4, Gabra2, Ptn, Gm8369, Hs6st3, Ms4a4b, Ms4a6c
## Sv2b, Mrgpra2b, Dnah11, Ptprk, Nyap2, Mlip, Grm5, Nkain2, Fat3, St6galnac3
## Negative: Apaf1, Mdm2, Mpp7, Polk, Hcn1, Ckap2, Cacna2d2, Map3k20, Nxph1, Abtb3
## Col18a1, Cntnap5c, Npas3, Zfp536, Sox6, bGHpolyA, Gad1, Nek7, Vwc2, Hspa5
## Rpl37, Grip1, Hsp90b1, Hspe1, Ccnb1ip1, Sltm, Afap1, Plxdc2, Pde5a, Eya4
## PC_ 5
## Positive: Ly86, Dock8, Mertk, Apbb1ip, Csf3r, Pik3ap1, Hk2, Inpp5d, Fgd2, Adap2
## Fyb, Unc93b1, C1qc, Ifi204, Siglech, Ctss, Tns3, Ikzf1, Bin2, Arhgap45
## Cx3cr1, Epsti1, C1qb, Sp100, Dock2, Cd84, Runx1, Apobec1, Ifi213, C1qa
## Negative: 6330411D24Rik, Kcnmb2, Abtb3, Nxph1, Gad1, Erbb4, Tmem132c, Cemip, Pde5a, Kcns3
## Cacna2d2, Pvalb, Sox6, Vwc2, Cntnap5c, Alk, Lhx6, Gad2, Cntnap3, Slc6a1
## Kcnip1, Cntnap4, Frmpd1, Grip1, Ankrd55, Afap1, Kcnc1, Gabrg3, Syt2, Adra1a
## Running UMAP and clustering
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1764
## Number of edges: 51408
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8821
## Number of communities: 12
## Elapsed time: 0 seconds
## 22:33:56 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:33:56 Read 1764 rows and found 18 numeric columns
## 22:33:56 Using Annoy for neighbor search, n_neighbors = 30
## 22:33:56 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:33:56 Writing NN index file to temp file /var/folders/m2/v2yg6yh92818c42pf0tfmbwr0000gn/T//Rtmp9fhWbf/file9ca836f04cc1
## 22:33:56 Searching Annoy index using 1 thread, search_k = 3000
## 22:33:57 Annoy recall = 100%
## 22:33:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:33:58 Initializing from normalized Laplacian + noise (using RSpectra)
## 22:33:58 Commencing optimization for 500 epochs, with 68814 positive edges
## 22:33:59 Optimization finished
## Performing WPRE analysis...
## Performing TC66T analysis...
## Performing bGHpolyA analysis...
## Performing oG analysis...
## Performing Cre analysis...
##
## ## Ntsr1 (ST.NT6) - WPRE
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 327 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Ntsr1 (ST.NT6) - TC66T
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 1655 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Ntsr1 (ST.NT6) - bGHpolyA
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 535 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Ntsr1 (ST.NT6) - oG
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 1521 rows containing non-finite outside the scale range
## (`stat_bin()`).
##
## ## Ntsr1 (ST.NT6) - Cre
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 1756 rows containing non-finite outside the scale range
## (`stat_bin()`).
# Function to add metadata from Maribel's project to Seurat object
add_metadata_to_seurat <- function(seurat_obj, metadata_path, experiment_name) {
# Load metadata
metadata <- read.csv(metadata_path, header = TRUE, stringsAsFactors = FALSE)
# Ensure the Barcode column is named correctly
colnames(metadata)[1] <- "Barcode"
# Subset metadata for the specific experiment
metadata_subset <- subset(metadata, Experiment == experiment_name)
# Ensure barcodes are the same format
rownames(metadata_subset) <- metadata_subset$Barcode
# Add metadata to Seurat object
common_barcodes <- intersect(colnames(seurat_obj), metadata_subset$Barcode)
seurat_obj <- AddMetaData(seurat_obj, metadata = metadata_subset[common_barcodes, ])
return(seurat_obj)
}
# Function to add gene reads to metadata
add_gene_reads_to_metadata <- function(seurat_obj, gene) {
counts_matrix <- seurat_obj@assays$RNA@counts
if (gene %in% rownames(counts_matrix)) {
gene_counts <- counts_matrix[gene, ]
gene_df <- data.frame(Barcode = colnames(counts_matrix), Gene_reads = as.numeric(gene_counts))
# Ensure metadata has Barcode column
seurat_obj@meta.data$Barcode <- rownames(seurat_obj@meta.data)
# Merge gene reads with metadata
metadata_with_gene <- merge(seurat_obj@meta.data, gene_df, by = "Barcode", all.x = TRUE)
rownames(metadata_with_gene) <- metadata_with_gene$Barcode
seurat_obj@meta.data <- metadata_with_gene
} else {
stop(sprintf("%s gene not found in the dataset.", gene))
}
return(seurat_obj)
}
# Function to calculate statistics for gene reads
calculate_gene_statistics <- function(seurat_obj, experiment_name, gene) {
# Ensure gene reads are added to metadata
seurat_obj <- add_gene_reads_to_metadata(seurat_obj, gene)
# Extract the metadata from the Seurat object
metadata <- seurat_obj@meta.data
# Check if Subclass column exists
if (!"Subclass" %in% colnames(metadata)) {
stop("Subclass column not found in metadata for experiment:", experiment_name)
}
# Group by Subclass and calculate statistics
gene_stats <- metadata %>%
group_by(Subclass) %>%
summarise(
mean_Gene = mean(Gene_reads, na.rm = TRUE),
median_Gene = median(Gene_reads, na.rm = TRUE),
sd_Gene = sd(Gene_reads, na.rm = TRUE),
var_Gene = var(Gene_reads, na.rm = TRUE),
min_Gene = min(Gene_reads, na.rm = TRUE),
max_Gene = max(Gene_reads, na.rm = TRUE),
count = n()
) %>%
mutate(Gene = gene)
return(gene_stats)
}
# Perform the analysis for each Seurat object and gene
genes_of_interest <- c("WPRE", "TC66T", "bGHpolyA", "oG", "Cre")
for (experiment_name in names(seurat_paths)) {
seurat_obj_path <- seurat_paths[[experiment_name]]
display_name <- experiment_display_names[[experiment_name]]
# Load the Seurat object
seurat_obj <- readRDS(seurat_obj_path)
# Add metadata from Maribel's project
seurat_obj <- add_metadata_to_seurat(seurat_obj, metadata_path, experiment_name)
# Ensure Subclass column exists in metadata
if (!"Subclass" %in% colnames(seurat_obj@meta.data)) {
cat(sprintf("Subclass column not found in metadata for experiment: %s\n", experiment_name))
next
}
for (gene in genes_of_interest) {
# Calculate gene statistics
gene_stats <- calculate_gene_statistics(seurat_obj, experiment_name, gene)
# Print the statistics
cat(sprintf("\n# %s (%s) - %s\n", display_name, experiment_name, gene))
print(gene_stats)
}
}
##
## # Sepw1 (ST.S11) - WPRE
## # A tibble: 14 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 18.5 5 42.6 1815. 0 268 68 WPRE
## 2 "L4 IT" 3.47 0 7.33 53.7 0 52 90 WPRE
## 3 "L5 IT" 8.02 1 14.6 214. 0 62 52 WPRE
## 4 "L5 NP" 0 0 0 0 0 0 2 WPRE
## 5 "L5 PT" 0 0 0 0 0 0 3 WPRE
## 6 "L6 CT" 0.2 0 0.447 0.2 0 1 5 WPRE
## 7 "L6 IT" 2.75 0 5.5 30.2 0 11 4 WPRE
## 8 "L6b" 0 0 NA NA 0 0 1 WPRE
## 9 "Lamp5" 21.4 15 17.4 301. 4 50 5 WPRE
## 10 "Pvalb" 36.6 28.5 33.2 1105. 0 120 20 WPRE
## 11 "Sncg" 0 0 NA NA 0 0 1 WPRE
## 12 "Sst" 13.4 3.5 20.5 421. 0 64 10 WPRE
## 13 "Vip " 3.25 1.5 4.57 20.9 0 10 4 WPRE
## 14 <NA> 58.3 1 138. 19036. 0 577 23 WPRE
##
## # Sepw1 (ST.S11) - TC66T
## # A tibble: 14 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0 0 0 0 0 0 68 TC66T
## 2 "L4 IT" 0.0111 0 0.105 0.0111 0 1 90 TC66T
## 3 "L5 IT" 0.0385 0 0.194 0.0377 0 1 52 TC66T
## 4 "L5 NP" 0 0 0 0 0 0 2 TC66T
## 5 "L5 PT" 0 0 0 0 0 0 3 TC66T
## 6 "L6 CT" 0 0 0 0 0 0 5 TC66T
## 7 "L6 IT" 0 0 0 0 0 0 4 TC66T
## 8 "L6b" 0 0 NA NA 0 0 1 TC66T
## 9 "Lamp5" 0 0 0 0 0 0 5 TC66T
## 10 "Pvalb" 0.15 0 0.489 0.239 0 2 20 TC66T
## 11 "Sncg" 0 0 NA NA 0 0 1 TC66T
## 12 "Sst" 0.2 0 0.422 0.178 0 1 10 TC66T
## 13 "Vip " 0 0 0 0 0 0 4 TC66T
## 14 <NA> 0.391 0 1.67 2.79 0 8 23 TC66T
##
## # Sepw1 (ST.S11) - bGHpolyA
## # A tibble: 14 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 1.84 0 4.63 21.5 0 29 68 bGHp…
## 2 "L4 IT" 0.411 0 1.07 1.14 0 7 90 bGHp…
## 3 "L5 IT" 0.558 0 1.43 2.06 0 6 52 bGHp…
## 4 "L5 NP" 0.5 0.5 0.707 0.5 0 1 2 bGHp…
## 5 "L5 PT" 0 0 0 0 0 0 3 bGHp…
## 6 "L6 CT" 0 0 0 0 0 0 5 bGHp…
## 7 "L6 IT" 0.5 0 1 1 0 2 4 bGHp…
## 8 "L6b" 0 0 NA NA 0 0 1 bGHp…
## 9 "Lamp5" 2.4 2 2.07 4.3 1 6 5 bGHp…
## 10 "Pvalb" 3.9 2 3.81 14.5 0 11 20 bGHp…
## 11 "Sncg" 0 0 NA NA 0 0 1 bGHp…
## 12 "Sst" 1.1 0.5 1.45 2.1 0 4 10 bGHp…
## 13 "Vip " 0.25 0 0.5 0.25 0 1 4 bGHp…
## 14 <NA> 7.09 0 16.5 274. 0 59 23 bGHp…
##
## # Sepw1 (ST.S11) - oG
## # A tibble: 14 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.412 0 1.43 2.04 0 10 68 oG
## 2 "L4 IT" 0.111 0 0.529 0.280 0 4 90 oG
## 3 "L5 IT" 0.0769 0 0.269 0.0724 0 1 52 oG
## 4 "L5 NP" 0 0 0 0 0 0 2 oG
## 5 "L5 PT" 0 0 0 0 0 0 3 oG
## 6 "L6 CT" 0 0 0 0 0 0 5 oG
## 7 "L6 IT" 0 0 0 0 0 0 4 oG
## 8 "L6b" 0 0 NA NA 0 0 1 oG
## 9 "Lamp5" 0.2 0 0.447 0.2 0 1 5 oG
## 10 "Pvalb" 0.7 0 1.49 2.22 0 5 20 oG
## 11 "Sncg" 0 0 NA NA 0 0 1 oG
## 12 "Sst" 0.4 0 0.699 0.489 0 2 10 oG
## 13 "Vip " 0.25 0 0.5 0.25 0 1 4 oG
## 14 <NA> 1.43 0 5.65 31.9 0 27 23 oG
##
## # Sepw1 (ST.S11) - Cre
## # A tibble: 14 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0294 0 0.170 0.0290 0 1 68 Cre
## 2 "L4 IT" 0 0 0 0 0 0 90 Cre
## 3 "L5 IT" 0 0 0 0 0 0 52 Cre
## 4 "L5 NP" 0 0 0 0 0 0 2 Cre
## 5 "L5 PT" 0 0 0 0 0 0 3 Cre
## 6 "L6 CT" 0 0 0 0 0 0 5 Cre
## 7 "L6 IT" 0 0 0 0 0 0 4 Cre
## 8 "L6b" 0 0 NA NA 0 0 1 Cre
## 9 "Lamp5" 0 0 0 0 0 0 5 Cre
## 10 "Pvalb" 0 0 0 0 0 0 20 Cre
## 11 "Sncg" 0 0 NA NA 0 0 1 Cre
## 12 "Sst" 0 0 0 0 0 0 10 Cre
## 13 "Vip " 0 0 0 0 0 0 4 Cre
## 14 <NA> 0.0435 0 0.209 0.0435 0 1 23 Cre
##
## # Scnn1a (ST.SC5) - WPRE
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 5.25 2 11.9 141. 0 177 774 WPRE
## 2 "L4 IT" 3.88 2 9.19 84.5 0 181 1005 WPRE
## 3 "L5 IT" 6.63 2 14.6 214. 0 236 639 WPRE
## 4 "L5 NP" 1.76 1 4.40 19.4 0 28 42 WPRE
## 5 "L5 PT" 9.46 3 16.4 268. 0 80 63 WPRE
## 6 "L6 CT" 3.89 1 9.25 85.6 0 69 133 WPRE
## 7 "L6 IT" 5.32 1 14.0 197. 0 97 72 WPRE
## 8 "L6b" 2.62 1 4.44 19.7 0 13 8 WPRE
## 9 "Lamp5" 13.0 8 21.0 440. 0 101 23 WPRE
## 10 "Pvalb" 23.7 9 54.7 2994. 0 449 164 WPRE
## 11 "Serpin… 2 2 NA NA 2 2 1 WPRE
## 12 "Sncg" 3.5 3.5 4.95 24.5 0 7 2 WPRE
## 13 "Sst" 109. 23.5 372. 138553. 0 2706 60 WPRE
## 14 "Vip " 2.1 1 3.97 15.7 0 19 30 WPRE
## 15 <NA> 40.3 1 336. 113221. 0 3994 155 WPRE
##
## # Scnn1a (ST.SC5) - TC66T
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0103 0 0.101 0.0102 0 1 774 TC66T
## 2 "L4 IT" 0.0109 0 0.113 0.0128 0 2 1005 TC66T
## 3 "L5 IT" 0.0188 0 0.157 0.0247 0 2 639 TC66T
## 4 "L5 NP" 0 0 0 0 0 0 42 TC66T
## 5 "L5 PT" 0.0635 0 0.246 0.0604 0 1 63 TC66T
## 6 "L6 CT" 0.00752 0 0.0867 0.00752 0 1 133 TC66T
## 7 "L6 IT" 0.0139 0 0.118 0.0139 0 1 72 TC66T
## 8 "L6b" 0 0 0 0 0 0 8 TC66T
## 9 "Lamp5" 0 0 0 0 0 0 23 TC66T
## 10 "Pvalb" 0.0854 0 0.406 0.164 0 4 164 TC66T
## 11 "Serpin… 0 0 NA NA 0 0 1 TC66T
## 12 "Sncg" 0 0 0 0 0 0 2 TC66T
## 13 "Sst" 0.183 0 0.911 0.830 0 5 60 TC66T
## 14 "Vip " 0 0 0 0 0 0 30 TC66T
## 15 <NA> 0.0323 0 0.265 0.0704 0 3 155 TC66T
##
## # Scnn1a (ST.SC5) - bGHpolyA
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.494 0 1.19 1.42 0 19 774 bGHp…
## 2 "L4 IT" 0.337 0 0.937 0.877 0 15 1005 bGHp…
## 3 "L5 IT" 0.510 0 1.36 1.84 0 22 639 bGHp…
## 4 "L5 NP" 0.190 0 0.671 0.451 0 4 42 bGHp…
## 5 "L5 PT" 0.746 0 1.45 2.10 0 7 63 bGHp…
## 6 "L6 CT" 0.346 0 0.708 0.501 0 4 133 bGHp…
## 7 "L6 IT" 0.472 0 1.15 1.32 0 7 72 bGHp…
## 8 "L6b" 0.25 0 0.707 0.5 0 2 8 bGHp…
## 9 "Lamp5" 1.09 0 2.02 4.08 0 9 23 bGHp…
## 10 "Pvalb" 1.96 1 4.72 22.3 0 37 164 bGHp…
## 11 "Serpin… 0 0 NA NA 0 0 1 bGHp…
## 12 "Sncg" 0.5 0.5 0.707 0.5 0 1 2 bGHp…
## 13 "Sst" 8.07 1.5 26.8 720. 0 193 60 bGHp…
## 14 "Vip " 0.133 0 0.346 0.120 0 1 30 bGHp…
## 15 <NA> 2.93 0 24.1 583. 0 288 155 bGHp…
##
## # Scnn1a (ST.SC5) - oG
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0271 0 0.178 0.0316 0 2 774 oG
## 2 "L4 IT" 0.0199 0 0.140 0.0195 0 1 1005 oG
## 3 "L5 IT" 0.0360 0 0.210 0.0442 0 2 639 oG
## 4 "L5 NP" 0 0 0 0 0 0 42 oG
## 5 "L5 PT" 0.0317 0 0.177 0.0312 0 1 63 oG
## 6 "L6 CT" 0.0301 0 0.171 0.0294 0 1 133 oG
## 7 "L6 IT" 0.0278 0 0.165 0.0274 0 1 72 oG
## 8 "L6b" 0 0 0 0 0 0 8 oG
## 9 "Lamp5" 0.130 0 0.344 0.119 0 1 23 oG
## 10 "Pvalb" 0.262 0 0.717 0.514 0 5 164 oG
## 11 "Serpin… 0 0 NA NA 0 0 1 oG
## 12 "Sncg" 0 0 0 0 0 0 2 oG
## 13 "Sst" 1 0 2.92 8.51 0 17 60 oG
## 14 "Vip " 0 0 0 0 0 0 30 oG
## 15 <NA> 0.342 0 2.96 8.73 0 36 155 oG
##
## # Scnn1a (ST.SC5) - Cre
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0 0 0 0 0 0 774 Cre
## 2 "L4 IT" 0.0149 0 0.121 0.0147 0 1 1005 Cre
## 3 "L5 IT" 0.00626 0 0.0789 0.00623 0 1 639 Cre
## 4 "L5 NP" 0 0 0 0 0 0 42 Cre
## 5 "L5 PT" 0 0 0 0 0 0 63 Cre
## 6 "L6 CT" 0 0 0 0 0 0 133 Cre
## 7 "L6 IT" 0 0 0 0 0 0 72 Cre
## 8 "L6b" 0 0 0 0 0 0 8 Cre
## 9 "Lamp5" 0 0 0 0 0 0 23 Cre
## 10 "Pvalb" 0 0 0 0 0 0 164 Cre
## 11 "Serpin… 0 0 NA NA 0 0 1 Cre
## 12 "Sncg" 0 0 0 0 0 0 2 Cre
## 13 "Sst" 0.0167 0 0.129 0.0167 0 1 60 Cre
## 14 "Vip " 0 0 0 0 0 0 30 Cre
## 15 <NA> 0.00645 0 0.0803 0.00645 0 1 155 Cre
##
## # Tlx3 (ST.T12) - WPRE
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 4.82 1 32.0 1023. 0 696 718 WPRE
## 2 "L4 IT" 2.29 1 5.94 35.3 0 92 921 WPRE
## 3 "L5 IT" 7.31 1 25.2 634. 0 522 652 WPRE
## 4 "L5 NP" 1.05 0.5 1.77 3.14 0 10 56 WPRE
## 5 "L5 PT" 10.8 3 20.2 408. 0 96 55 WPRE
## 6 "L6 CT" 3.36 1 8.21 67.5 0 74 225 WPRE
## 7 "L6 IT" 3.14 1 7.03 49.4 0 40 74 WPRE
## 8 "L6b" 6.33 0.5 12.7 162. 0 32 6 WPRE
## 9 "Lamp5" 32.4 3 66.4 4414. 0 223 13 WPRE
## 10 "Pvalb" 30.7 5 113. 12681. 0 1007 141 WPRE
## 11 "Sst" 18.8 4.5 40.7 1660. 0 226 46 WPRE
## 12 "Vip " 4.67 1 8.45 71.5 0 32 33 WPRE
## 13 <NA> 2.74 0 9.33 87.1 0 71 95 WPRE
##
## # Tlx3 (ST.T12) - TC66T
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0209 0 0.170 0.0289 0 3 718 TC66T
## 2 "L4 IT" 0 0 0 0 0 0 921 TC66T
## 3 "L5 IT" 0.00920 0 0.0956 0.00913 0 1 652 TC66T
## 4 "L5 NP" 0 0 0 0 0 0 56 TC66T
## 5 "L5 PT" 0.0727 0 0.325 0.106 0 2 55 TC66T
## 6 "L6 CT" 0.0133 0 0.115 0.0132 0 1 225 TC66T
## 7 "L6 IT" 0.0270 0 0.163 0.0267 0 1 74 TC66T
## 8 "L6b" 0 0 0 0 0 0 6 TC66T
## 9 "Lamp5" 0.692 0 1.55 2.40 0 5 13 TC66T
## 10 "Pvalb" 0.0922 0 0.314 0.0986 0 2 141 TC66T
## 11 "Sst" 0.0870 0 0.354 0.126 0 2 46 TC66T
## 12 "Vip " 0 0 0 0 0 0 33 TC66T
## 13 <NA> 0.0316 0 0.176 0.0309 0 1 95 TC66T
##
## # Tlx3 (ST.T12) - bGHpolyA
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.323 0 1.98 3.92 0 43 718 bGHp…
## 2 "L4 IT" 0.187 0 0.513 0.263 0 4 921 bGHp…
## 3 "L5 IT" 0.468 0 1.51 2.29 0 26 652 bGHp…
## 4 "L5 NP" 0.0893 0 0.345 0.119 0 2 56 bGHp…
## 5 "L5 PT" 0.764 0 1.14 1.29 0 6 55 bGHp…
## 6 "L6 CT" 0.204 0 0.592 0.351 0 5 225 bGHp…
## 7 "L6 IT" 0.243 0 0.592 0.351 0 3 74 bGHp…
## 8 "L6b" 0.167 0 0.408 0.167 0 1 6 bGHp…
## 9 "Lamp5" 2.31 0 4.40 19.4 0 15 13 bGHp…
## 10 "Pvalb" 1.74 0 5.88 34.5 0 48 141 bGHp…
## 11 "Sst" 0.957 0 2.49 6.22 0 15 46 bGHp…
## 12 "Vip " 0.333 0 0.990 0.979 0 4 33 bGHp…
## 13 <NA> 0.232 0 0.643 0.414 0 4 95 bGHp…
##
## # Tlx3 (ST.T12) - oG
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0613 0 0.447 0.200 0 8 718 oG
## 2 "L4 IT" 0.0195 0 0.153 0.0235 0 2 921 oG
## 3 "L5 IT" 0.0966 0 0.416 0.173 0 4 652 oG
## 4 "L5 NP" 0.0179 0 0.134 0.0179 0 1 56 oG
## 5 "L5 PT" 0.164 0 0.501 0.251 0 3 55 oG
## 6 "L6 CT" 0.0222 0 0.175 0.0308 0 2 225 oG
## 7 "L6 IT" 0.0541 0 0.281 0.0792 0 2 74 oG
## 8 "L6b" 0.167 0 0.408 0.167 0 1 6 oG
## 9 "Lamp5" 0.308 0 0.855 0.731 0 3 13 oG
## 10 "Pvalb" 0.340 0 0.869 0.755 0 6 141 oG
## 11 "Sst" 0.391 0 1.18 1.40 0 6 46 oG
## 12 "Vip " 0 0 0 0 0 0 33 oG
## 13 <NA> 0.0947 0 0.485 0.236 0 4 95 oG
##
## # Tlx3 (ST.T12) - Cre
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0 0 0 0 0 0 718 Cre
## 2 "L4 IT" 0 0 0 0 0 0 921 Cre
## 3 "L5 IT" 0.00767 0 0.0873 0.00762 0 1 652 Cre
## 4 "L5 NP" 0 0 0 0 0 0 56 Cre
## 5 "L5 PT" 0 0 0 0 0 0 55 Cre
## 6 "L6 CT" 0 0 0 0 0 0 225 Cre
## 7 "L6 IT" 0 0 0 0 0 0 74 Cre
## 8 "L6b" 0 0 0 0 0 0 6 Cre
## 9 "Lamp5" 0 0 0 0 0 0 13 Cre
## 10 "Pvalb" 0 0 0 0 0 0 141 Cre
## 11 "Sst" 0 0 0 0 0 0 46 Cre
## 12 "Vip " 0 0 0 0 0 0 33 Cre
## 13 <NA> 0.0105 0 0.103 0.0105 0 1 95 Cre
##
## # Tlx3 (ST.T13) - WPRE
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 2.31 0 10.4 108. 0 194 757 WPRE
## 2 "L4 IT" 1.53 0 6.23 38.8 0 135 1078 WPRE
## 3 "L5 IT" 5.31 1 15.6 243. 0 172 755 WPRE
## 4 "L5 NP" 0.804 0 1.11 1.24 0 4 51 WPRE
## 5 "L5 PT" 5.48 1 11.7 137. 0 84 81 WPRE
## 6 "L6 CT" 5.55 2 10.8 116. 0 91 198 WPRE
## 7 "L6 IT" 7.67 0.5 24.7 608. 0 122 70 WPRE
## 8 "Lamp5" 12.4 1 29.8 887. 0 114 17 WPRE
## 9 "Pvalb" 29.9 5 186. 34740. 0 2320 159 WPRE
## 10 "Sncg" 145 145 NA NA 145 145 1 WPRE
## 11 "Sst" 11.2 2 33.0 1088. 0 218 44 WPRE
## 12 "Vip " 1.94 1 3.67 13.4 0 19 33 WPRE
## 13 <NA> 6.35 0 54.3 2952. 0 731 184 WPRE
##
## # Tlx3 (ST.T13) - TC66T
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.00661 0 0.0811 0.00657 0 1 757 TC66T
## 2 "L4 IT" 0.00278 0 0.0527 0.00278 0 1 1078 TC66T
## 3 "L5 IT" 0.0119 0 0.109 0.0118 0 1 755 TC66T
## 4 "L5 NP" 0 0 0 0 0 0 51 TC66T
## 5 "L5 PT" 0.0123 0 0.111 0.0123 0 1 81 TC66T
## 6 "L6 CT" 0.0505 0 0.242 0.0583 0 2 198 TC66T
## 7 "L6 IT" 0.0143 0 0.120 0.0143 0 1 70 TC66T
## 8 "Lamp5" 0.0588 0 0.243 0.0588 0 1 17 TC66T
## 9 "Pvalb" 0.0818 0 0.405 0.164 0 3 159 TC66T
## 10 "Sncg" 0 0 NA NA 0 0 1 TC66T
## 11 "Sst" 0.0682 0 0.255 0.0650 0 1 44 TC66T
## 12 "Vip " 0 0 0 0 0 0 33 TC66T
## 13 <NA> 0.0272 0 0.194 0.0375 0 2 184 TC66T
##
## # Tlx3 (ST.T13) - bGHpolyA
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.419 0 1.45 2.11 0 28 757 bGHp…
## 2 "L4 IT" 0.278 0 0.916 0.840 0 19 1078 bGHp…
## 3 "L5 IT" 0.774 0 2.10 4.39 0 23 755 bGHp…
## 4 "L5 NP" 0.235 0 0.473 0.224 0 2 51 bGHp…
## 5 "L5 PT" 0.852 0 2.10 4.43 0 15 81 bGHp…
## 6 "L6 CT" 0.823 0 1.57 2.45 0 12 198 bGHp…
## 7 "L6 IT" 1.07 0 3.22 10.4 0 17 70 bGHp…
## 8 "Lamp5" 1.65 0 3.95 15.6 0 16 17 bGHp…
## 9 "Pvalb" 3.96 1 24.7 608. 0 308 159 bGHp…
## 10 "Sncg" 22 22 NA NA 22 22 1 bGHp…
## 11 "Sst" 1.39 0 3.98 15.8 0 26 44 bGHp…
## 12 "Vip " 0.394 0 0.704 0.496 0 2 33 bGHp…
## 13 <NA> 0.967 0 7.39 54.6 0 99 184 bGHp…
##
## # Tlx3 (ST.T13) - oG
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0198 0 0.181 0.0327 0 3 757 oG
## 2 "L4 IT" 0.0102 0 0.146 0.0213 0 3 1078 oG
## 3 "L5 IT" 0.105 0 0.652 0.425 0 12 755 oG
## 4 "L5 NP" 0.0196 0 0.140 0.0196 0 1 51 oG
## 5 "L5 PT" 0.0370 0 0.190 0.0361 0 1 81 oG
## 6 "L6 CT" 0.0152 0 0.122 0.0150 0 1 198 oG
## 7 "L6 IT" 0.0286 0 0.239 0.0571 0 2 70 oG
## 8 "Lamp5" 0.176 0 0.393 0.154 0 1 17 oG
## 9 "Pvalb" 0.321 0 1.77 3.13 0 21 159 oG
## 10 "Sncg" 0 0 NA NA 0 0 1 oG
## 11 "Sst" 0.114 0 0.493 0.243 0 3 44 oG
## 12 "Vip " 0 0 0 0 0 0 33 oG
## 13 <NA> 0.0924 0 0.570 0.325 0 7 184 oG
##
## # Tlx3 (ST.T13) - Cre
## # A tibble: 13 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.00132 0 0.0363 0.00132 0 1 757 Cre
## 2 "L4 IT" 0 0 0 0 0 0 1078 Cre
## 3 "L5 IT" 0.00795 0 0.0888 0.00789 0 1 755 Cre
## 4 "L5 NP" 0 0 0 0 0 0 51 Cre
## 5 "L5 PT" 0 0 0 0 0 0 81 Cre
## 6 "L6 CT" 0.00505 0 0.0711 0.00505 0 1 198 Cre
## 7 "L6 IT" 0 0 0 0 0 0 70 Cre
## 8 "Lamp5" 0 0 0 0 0 0 17 Cre
## 9 "Pvalb" 0 0 0 0 0 0 159 Cre
## 10 "Sncg" 0 0 NA NA 0 0 1 Cre
## 11 "Sst" 0 0 0 0 0 0 44 Cre
## 12 "Vip " 0 0 0 0 0 0 33 Cre
## 13 <NA> 0.00543 0 0.0737 0.00543 0 1 184 Cre
##
## # Ntsr1 (ST.NT6) - WPRE
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 13.9 3 27.7 770. 0 229 327 WPRE
## 2 "L4 IT" 6.92 3 14.0 197. 0 136 375 WPRE
## 3 "L5 IT" 22.6 6 43.4 1882. 0 445 364 WPRE
## 4 "L5 NP" 1.83 2 1.94 3.76 0 8 63 WPRE
## 5 "L5 PT" 38.9 13 71.9 5163. 0 430 46 WPRE
## 6 "L6 CT" 7.97 2 24.8 614. 0 340 289 WPRE
## 7 "L6 IT" 5.85 2 11.0 121. 0 79 170 WPRE
## 8 "L6b" 7.6 2 19.9 394. 0 79 15 WPRE
## 9 "Lamp5" 8 4 8.72 76 2 18 3 WPRE
## 10 "Pvalb" 63.5 15 185. 34289. 0 1540 77 WPRE
## 11 "Serpin… 0 0 NA NA 0 0 1 WPRE
## 12 "Sncg" 1 1 NA NA 1 1 1 WPRE
## 13 "Sst" 89.6 7 194. 37599. 0 841 27 WPRE
## 14 "Vip " 25.4 2 43.1 1857. 0 128 10 WPRE
## 15 <NA> 7.81 1 24.7 609. 0 201 86 WPRE
##
## # Ntsr1 (ST.NT6) - TC66T
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.0765 0 0.346 0.120 0 3 327 TC66T
## 2 "L4 IT" 0.056 0 0.263 0.0690 0 2 375 TC66T
## 3 "L5 IT" 0.0962 0 0.330 0.109 0 2 364 TC66T
## 4 "L5 NP" 0 0 0 0 0 0 63 TC66T
## 5 "L5 PT" 0.217 0 0.513 0.263 0 2 46 TC66T
## 6 "L6 CT" 0.0415 0 0.247 0.0608 0 3 289 TC66T
## 7 "L6 IT" 0.0235 0 0.187 0.0349 0 2 170 TC66T
## 8 "L6b" 0.133 0 0.516 0.267 0 2 15 TC66T
## 9 "Lamp5" 0.667 0 1.15 1.33 0 2 3 TC66T
## 10 "Pvalb" 0.221 0 0.681 0.464 0 5 77 TC66T
## 11 "Serpin… 0 0 NA NA 0 0 1 TC66T
## 12 "Sncg" 0 0 NA NA 0 0 1 TC66T
## 13 "Sst" 0.407 0 0.797 0.635 0 3 27 TC66T
## 14 "Vip " 0 0 0 0 0 0 10 TC66T
## 15 <NA> 0.0465 0 0.212 0.0449 0 1 86 TC66T
##
## # Ntsr1 (ST.NT6) - bGHpolyA
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 4.95 1 9.61 92.3 0 68 327 bGHp…
## 2 "L4 IT" 2.62 1 4.92 24.2 0 48 375 bGHp…
## 3 "L5 IT" 8.54 3 16.5 274. 0 180 364 bGHp…
## 4 "L5 NP" 0.889 1 0.969 0.939 0 4 63 bGHp…
## 5 "L5 PT" 15.1 4.5 23.7 563. 0 130 46 bGHp…
## 6 "L6 CT" 3.00 1 8.60 73.9 0 110 289 bGHp…
## 7 "L6 IT" 2.17 1 4.12 17.0 0 36 170 bGHp…
## 8 "L6b" 2.53 1 7.09 50.3 0 28 15 bGHp…
## 9 "Lamp5" 2.33 1 3.21 10.3 0 6 3 bGHp…
## 10 "Pvalb" 23.7 6 74.4 5534. 0 631 77 bGHp…
## 11 "Serpin… 0 0 NA NA 0 0 1 bGHp…
## 12 "Sncg" 3 3 NA NA 3 3 1 bGHp…
## 13 "Sst" 33.1 3 73.4 5389. 0 329 27 bGHp…
## 14 "Vip " 8.2 2 14.5 211. 0 47 10 bGHp…
## 15 <NA> 2.97 0 9.44 89.1 0 78 86 bGHp…
##
## # Ntsr1 (ST.NT6) - oG
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0.211 0 0.904 0.817 0 10 327 oG
## 2 "L4 IT" 0.173 0 0.633 0.400 0 6 375 oG
## 3 "L5 IT" 0.321 0 1.55 2.40 0 25 364 oG
## 4 "L5 NP" 0.0159 0 0.126 0.0159 0 1 63 oG
## 5 "L5 PT" 0.652 0 1.46 2.14 0 8 46 oG
## 6 "L6 CT" 0.218 0 0.789 0.622 0 9 289 oG
## 7 "L6 IT" 0.100 0 0.387 0.150 0 3 170 oG
## 8 "L6b" 0 0 0 0 0 0 15 oG
## 9 "Lamp5" 2.67 0 4.62 21.3 0 8 3 oG
## 10 "Pvalb" 0.870 0 2.00 3.98 0 10 77 oG
## 11 "Serpin… 0 0 NA NA 0 0 1 oG
## 12 "Sncg" 1 1 NA NA 1 1 1 oG
## 13 "Sst" 1.11 0 1.97 3.87 0 7 27 oG
## 14 "Vip " 0.3 0 0.483 0.233 0 1 10 oG
## 15 <NA> 0.209 0 0.616 0.379 0 4 86 oG
##
## # Ntsr1 (ST.NT6) - Cre
## # A tibble: 15 × 9
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 "L2/3 I… 0 0 0 0 0 0 327 Cre
## 2 "L4 IT" 0 0 0 0 0 0 375 Cre
## 3 "L5 IT" 0 0 0 0 0 0 364 Cre
## 4 "L5 NP" 0 0 0 0 0 0 63 Cre
## 5 "L5 PT" 0 0 0 0 0 0 46 Cre
## 6 "L6 CT" 0.0208 0 0.143 0.0204 0 1 289 Cre
## 7 "L6 IT" 0 0 0 0 0 0 170 Cre
## 8 "L6b" 0.0667 0 0.258 0.0667 0 1 15 Cre
## 9 "Lamp5" 0 0 0 0 0 0 3 Cre
## 10 "Pvalb" 0.0130 0 0.114 0.0130 0 1 77 Cre
## 11 "Serpin… 0 0 NA NA 0 0 1 Cre
## 12 "Sncg" 0 0 NA NA 0 0 1 Cre
## 13 "Sst" 0 0 0 0 0 0 27 Cre
## 14 "Vip " 0 0 0 0 0 0 10 Cre
## 15 <NA> 0 0 0 0 0 0 86 Cre
# Function to calculate statistics for a gene and add to the combined data frame
calculate_and_add_stats <- function(seurat_obj, gene, experiment_name, combined_stats) {
# Ensure gene reads are added to metadata
seurat_obj <- add_gene_reads_to_metadata(seurat_obj, gene)
# Extract the metadata from the Seurat object
metadata <- seurat_obj@meta.data
# Check if Subclass column exists
if (!"Subclass" %in% colnames(metadata)) {
stop("Subclass column not found in metadata for experiment:", experiment_name)
}
# Group by Subclass and calculate statistics
gene_stats <- metadata %>%
group_by(Subclass) %>%
summarise(
mean_Gene = mean(Gene_reads, na.rm = TRUE),
median_Gene = median(Gene_reads, na.rm = TRUE),
sd_Gene = sd(Gene_reads, na.rm = TRUE),
var_Gene = var(Gene_reads, na.rm = TRUE),
min_Gene = min(Gene_reads, na.rm = TRUE),
max_Gene = max(Gene_reads, na.rm = TRUE),
count = n()
) %>%
mutate(Gene = gene, Experiment = experiment_name)
# Append to the combined statistics data frame
combined_stats <- rbind(combined_stats, gene_stats)
return(combined_stats)
}
# Initialize an empty data frame to store the combined statistics
combined_stats <- data.frame()
# Define the genes of interest
genes_of_interest <- c("WPRE", "TC66T", "bGHpolyA", "oG", "Cre")
# Paths to Seurat objects
seurat_paths <- list(
ST.S11 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_S11.rds",
ST.SC5 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_SC5.rds",
ST.T12 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T12.rds",
ST.T13 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_T13.rds",
ST.NT6 = "/Users/saikrishnabhamidipati/Library/Mobile Documents/com~apple~CloudDocs/grad school/research/data/processed/patino2024_wpre_analyses/patino2024_wpre_analysis_ver3/rds files/ST_NT6.rds"
)
# Experiment display names
experiment_display_names <- list(
ST.S11 = "Sepw1",
ST.SC5 = "Scnn1a",
ST.T12 = "Tlx3",
ST.T13 = "Tlx3",
ST.NT6 = "Ntsr1"
)
# Perform analysis for each Seurat object and each gene
for (experiment_name in names(seurat_paths)) {
seurat_obj_path <- seurat_paths[[experiment_name]]
display_name <- experiment_display_names[[experiment_name]]
# Load the Seurat object
seurat_obj <- readRDS(seurat_obj_path)
# Add metadata from Maribel's project
seurat_obj <- add_metadata_to_seurat(seurat_obj, metadata_path, experiment_name)
# Ensure Subclass column exists in metadata
if (!"Subclass" %in% colnames(seurat_obj@meta.data)) {
cat(sprintf("Subclass column not found in metadata for experiment: %s\n", experiment_name))
next
}
for (gene in genes_of_interest) {
# Calculate gene statistics and add to the combined data frame
combined_stats <- calculate_and_add_stats(seurat_obj, gene, experiment_name, combined_stats)
}
}
# Print the combined statistics table
print(combined_stats)
## # A tibble: 350 × 10
## Subclass mean_Gene median_Gene sd_Gene var_Gene min_Gene max_Gene count Gene
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 L2/3 IT 18.5 5 42.6 1815. 0 268 68 WPRE
## 2 L4 IT 3.47 0 7.33 53.7 0 52 90 WPRE
## 3 L5 IT 8.02 1 14.6 214. 0 62 52 WPRE
## 4 L5 NP 0 0 0 0 0 0 2 WPRE
## 5 L5 PT 0 0 0 0 0 0 3 WPRE
## 6 L6 CT 0.2 0 0.447 0.2 0 1 5 WPRE
## 7 L6 IT 2.75 0 5.5 30.2 0 11 4 WPRE
## 8 L6b 0 0 NA NA 0 0 1 WPRE
## 9 Lamp5 21.4 15 17.4 301. 4 50 5 WPRE
## 10 Pvalb 36.6 28.5 33.2 1105. 0 120 20 WPRE
## # ℹ 340 more rows
## # ℹ 1 more variable: Experiment <chr>