Mice_Data.R

# Load Required Libraries ---------------------------------------------------
library(tidyverse)      # Data manipulation and visualization
library(dplyr)          # Data manipulation
library(Seurat)         # Single-cell RNA-seq analysis
library(SeuratObject)   # Data structures for Seurat
library(patchwork)      # Combining plots
library(here)           # File path management
library(Matrix)         # Sparse matrix handling
library(biomaRt)        # Accessing BioMart for gene annotations

# Load and Preprocess Datasets ---------------------------------------------


# Load Paper 24 Dataset (10X Genomics HDF5 format)
Data <- Read10X_h5(filename = "GSM6997574_Control_Hep_raw_feature_bc_matrix.h5") 
Data24_seu <- CreateSeuratObject(counts = Data, project = 'Paper_24', min.cells = 3, min.features = 200) # create Seurat object
Data24_seu[["percent_mt"]] <- PercentageFeatureSet(Data_seu, pattern = "^mt-") # Calculate and add mitochondrial content (mouse gene symbols start with 'mt-')

# Save as Seurat object for further analysis
saveRDS(Data24_seu, file = "Data_seu_Paper24.rds")
Data_seu <- readRDS(file = "Data_seu_Paper24.rds") # During anlysis


# Load Paper 27 Dataset (Matrix Market format with separate features and barcodes)
barcodes <- read.table(gzfile("C:/Users/natal/Desktop/Disseration/Dissertation_Project_ScRNAseq/Paper_27/GSM4321524_Liver_Aged_barcodes.tsv.gz"), # Raw counts availble in gz format
                       header = FALSE, sep = "\t")
features <- read.table("C:/Users/natal/Desktop/Disseration/Dissertation_Project_ScRNAseq/Paper_27/features.tsv",
                       header = FALSE, sep = "\t")
matrix <- readMM("C:/Users/natal/Desktop/Disseration/Dissertation_Project_ScRNAseq/Paper_27/matrix.mtx")

# Assign row names (genes) to the matrix from the features file
rownames(matrix) <- features$V1
# Assign column names (barcodes) to the matrix
colnames(matrix) <- barcodes$V1

Data27_seu <- CreateSeuratObject(counts = matrix) # Create Seurat object

# Annotate gene names with biomaRt to correctly identify mitochondrial genes
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
gene_ids <- rownames(Data27_seu)
gene_annotations <- getBM(attributes = c("ensembl_gene_id", "external_gene_name"),
                          filters = "ensembl_gene_id",
                          values = gene_ids,
                          mart = ensembl)

mt_genes <- gene_annotations$ensembl_gene_id[grepl("^mt-", gene_annotations$external_gene_name, ignore.case = TRUE)]

Data27_seu[["percent_mt"]] <- PercentageFeatureSet(Data27_seu, features = mt_genes)


saveRDS(Data27_seu, file = "Data_seu_Paper27.rds") #Save as seurat object for downstream analysis
Data_seu <- readRDS(file = "Data_seu_Paper27.rds")



# Load Paper 56 Dataset (10X Genomics HDF5 format)
Data56 <- Read10X_h5(here("Paper_56", "GSM5354106_RS025_filtered_gene_bc_matrices_h5.h5"))

Data56_seu <- CreateSeuratObject(counts = Data56, project = 'Paper_56', min.cells = 3, min.features = 200) #features are genes

Data56_seu[["percent_mt"]] <- PercentageFeatureSet(Data56_seu, pattern = "^mt-")

saveRDS(Data56_seu, file = "Data_seu_Paper56.rds") # Save as Seurat object
Data_seu <- readRDS(file = "Data_seu_Paper56.rds")



# VIOLIN PLOTS

# Reload saved datasets
Data24_seu <- readRDS(file = "Data_seu_Paper24.rds")
Data27_seu <- readRDS(file = "Data_seu_Paper27.rds")
Data56_seu <- readRDS(file = "Data_seu_Paper56.rds")

# Visualise basic QC metrics for each dataset
VlnPlot(Data_seu, features=c("nCount_RNA", "nFeature_RNA", "percent_mt"))
VlnPlot(Data27_seu, features=c("nCount_RNA", "nFeature_RNA", "percent_mt"))
VlnPlot(Data56_seu, features=c("nCount_RNA", "nFeature_RNA", "percent_mt"))

# Add dataset labels for merging
Data24_seu$dataset <- "Paper 24"
Data27_seu$dataset <- "Paper 27"
Data56_seu$dataset <- "Paper 56"


# Merge all datasets for joint visualisation
merged_seu <- merge(Data24_seu, y = c(Data27_seu, Data56_seu), add.cell.ids = c("Data24_seu", "Data27_seu", "Data56_seu"))

# Generate grouped violin plots for key metrics
vln_genes <- VlnPlot(merged_seu, features = "nFeature_RNA", group.by = "dataset") + ggtitle("Genes per Cell")
vln_umi <- VlnPlot(merged_seu, features = "nCount_RNA", group.by = "dataset") + ggtitle("UMIs per Cell")
vln_mt <- VlnPlot(merged_seu, features = "percent_mt", group.by = "dataset") + ggtitle("Mitochondrial Percentage")



# APPLYING QC FILTERING, CLUSTERING AND DEG

# Function to filter cells based on QC metrics and run clustering + DEG detection
QC_filter <- function(seurat_object, 
                       min_counts = 0, max_counts = Inf, 
                       min_genes = 0, max_genes = Inf,
                       max_percent_mt = Inf) {
  try(
    seurat_object_filtered <- subset(seurat_object, subset = 
                                       nFeature_RNA > min_genes & nFeature_RNA < max_genes & 
                                       nCount_RNA > min_counts & nCount_RNA < max_counts &
                                       percent_mt < max_percent_mt))
  if(exists("seurat_object_filtered")) {
    ncells <- ncol(seurat_object_filtered)
    clusters_DEG <- cluster_cells_and_DEG(seurat_object_filtered)
    
    return(list(ncells = ncells,
                clusters = clusters_DEG$num_clusters,
                DEG = clusters_DEG$num_DEG))  
  } else {
    return(list(ncells = 0,
                clusters = 0,
                DEG = 0))  
  }
}


# Function to normalise, cluster, and compute DEGs
cluster_cells_and_DEG <- function(seurat_object, resolution = 0.4) {
  # Normalise using SCT
  seurat_object <- seurat_object %>% 
    NormalizeData(normalization.method = "LogNormalize", scale.factor = 10000) %>% 
    ScaleData() %>% 
    # Find the 2000 top variable features, used to calculate PCA
    FindVariableFeatures(nfeatures = 2000) %>% 
    # Run Principal Component Analysis, returning the first 50 principal components
    RunPCA(npcs=50) %>% 
    FindNeighbors() %>% 
    FindClusters(resolution = resolution) 
  markers <- FindAllMarkers(seurat_object, 
                            logfc.threshold = 0.25, min.pct = 0.1,
                            only.pos = TRUE)
  markers <- markers %>% 
    filter(p_val_adj < 0.05)
  
  #table(markers$cluster)
  
  return(list(
    num_clusters = length(unique(Idents(seurat_object))),
    num_DEG = nrow(markers)))
}


#######################################################################################################################
# QC strategy 1: Changing Genes/cell
G_cells <- NULL #make sure this is empty everytime i run the code

i <- 1
for(g in seq(0,3000,500)){
  print(g)
  cells<- QC_filter2(Data_seu, 
                     min_counts = 750, max_counts = 30000,
                     min_genes = g,
                     max_percent_mt = 10)
  
  G_cells[[i]] <- cells
  i <- i + 1
}   

# Collect results
n_cells <- sapply(G_cells, function(x){return(x$ncells)})
n_clusters <- sapply(G_cells, function(x){return(x$clusters)})
DEG <- sapply(G_cells, function(x){return(x$DEG)})

QC_results_df <- data.frame(
  Genecount = seq(0,3000,500),  # The varying mitochondrial percentage
  Cells_Filtered = n_cells,      # Number of cells retained after filtering
  Clusters = n_clusters,         # Number of clusters
  DEG_Count = DEG,                 # Number of differentially expressed genes (or a placeholder)
  DEG_Per= DEG/max(DEG, na.rm = TRUE),
  Cells_Per = n_cells/max(n_cells, na.rm = TRUE),
  Clusters_Per = n_clusters/max(n_clusters, na.rm = TRUE)
)

write.csv(QC_results_df, file = "QC_resultsXX.csv", row.names = FALSE) # Correct Paper number


#####################################################################################################
#QC STRAT 2: Changing UMI

#Reload csvs corresponding to UMI

U_cells <- NULL #make sure this is empty everytime i run the code

i <- 1
for(u in seq(0,6000,1000)){
  print(u)
  cells<- QC_filter2(Data_seu, 
                     min_counts = u,
                     min_genes = 200, max_genes = 6000,
                     max_percent_mt = 10)
  
  U_cells[[i]] <- cells
  i <- i + 1
}   

n_cells <- sapply(U_cells, function(x){return(x$ncells)})
n_clusters <- sapply(U_cells, function(x){return(x$clusters)})
DEG <- sapply(U_cells, function(x){return(x$DEG)})

QC_results_df <- data.frame(
  UMI = seq(0,6000,1000),  # The varying mitochondrial percentage
  Cells_Filtered = n_cells,      # Number of cells retained after filtering
  Clusters = n_clusters,         # Number of clusters
  DEG_Count = DEG,
  DEG_Per= DEG/max(DEG, na.rm = TRUE),# Number of differentially expressed genes (or a placeholder)
  Clusters_Per = n_clusters/max(n_clusters, na.rm = TRUE),
  Cells_Per = n_cells/max(n_cells, na.rm = TRUE)
)


write.csv(QC_results_df, file = "QC_resultsXX.csv", row.names = FALSE)



###################################################################################################
# QC Strategy 3: Changing Mitochondrial %

# Reload csvs corresponding to Mt%

M_cells <- NULL #make sure this is empty everytime i run the code

i <- 1
for(m in seq(0,100,25)){
  print(m)
  cells<- QC_filter2(Data_seu, 
                     min_counts = 750, max_counts = 30000,
                     min_genes = 200, max_genes = 6000,
                     max_percent_mt = m)
  
  M_cells[[i]] <- cells
  i <- i + 1
}   

n_cells <- sapply(M_cells, function(x){return(x$ncells)})
n_clusters <- sapply(M_cells, function(x){return(x$clusters)})
DEG <- sapply(M_cells, function(x){return(x$DEG)})

QC_results_df <- data.frame(
  MITO = seq(0,100,25),  # The varying mitochondrial percentage
  Cells_Filtered = n_cells,      # Number of cells retained after filtering
  Clusters = n_clusters,         # Number of clusters
  DEG_Count = DEG,
  DEG_Per= DEG/max(DEG, na.rm = TRUE),# Number of differentially expressed genes (or a placeholder)
  Clusters_Per = n_clusters/max(n_clusters, na.rm = TRUE),
  Cells_Per = n_cells/max(n_cells, na.rm = TRUE)
)


write.csv(QC_results_df, file = "QC_resultsXX.csv", row.names = FALSE)


# By the end of this process, there will be 9 csv files



################################################################################################
# Plots (3 QC STRATS(Gene/cell, UMI/cell, Mt%) X 3 Outputs(Filteredcells, Clusters, DEGs) = 9 Plots)

# 3 PLOTS FOR GENE/CELL

#Read in required csv files
QC_results24 <- read_csv("QC_results24GENE.csv")
QC_results27 <- read_csv("QC_results27GENE.csv")
QC_results56 <- read_csv("QC_results56GENE.csv")

QC_results24$Dataset <- "Paper 24"
QC_results27$Dataset <- "Paper 27"
QC_results56$Dataset <- "Paper 56"

QC_all <- bind_rows(QC_results24, QC_results27, QC_results56)

max<- max(QC_all$Clusters) # change depending on what is being plot

plot <- ggplot(QC_all, aes(x = Genecount, y = Clusters, color = Dataset)) +  #Change Y according to what is being plot
  geom_line(linewidth = 1) + 
  geom_point(size = 2.5) +
  theme_classic() +
  labs(
    title = "Clusters After Filtering with changing minimum gene thresholds",
    x = "Gene/cell threshold",
    y = "Number of Clusters") +
  scale_x_continuous(breaks = seq(0,3000,500))+
  scale_y_continuous(breaks = seq(0,max(QC_all$Clusters),2))+
  # scale_y_continuous(labels = scales::label_percent(accuracy = 1))+  # Required if plotting a percentage
  theme(
    axis.title.x = element_text(size = 14, face = "bold"),  # X-axis label size
    axis.title.y = element_text(size = 14, face = "bold"),  # Y-axis label size
    axis.text.x = element_text(size = 12),  # X-axis tick label size
    axis.text.y = element_text(size = 12),  # Y-axis tick label size
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # Title size and center alignment
    legend.text = element_text(size = 10),  # Legend text size
    legend.title = element_text(size = 10, face = "bold")  # Legend title size
  )

# save with species and QC strategy
ggsave("combined_plot_CLustersGENEMice.png", plot = plot, width = 7, height = 7.5, units = "in", dpi = 300) #save image



# 3 PLOTS FOR UMI/CELL

plot <- ggplot(QC_all, aes(x = UMI, y = Cells_Per, color = Dataset)) +  #Change Y according to what is being plot
  geom_line(linewidth = 1) + 
  geom_point(size = 2.5) +
  theme_classic() +
  labs(
    title = "Percentage of Filtered Cells with Increasing Minimum UMI Thresholds",
    x = "UMI/cell threshold",
    y = "Percentage of Filtered Cells") +
  scale_x_continuous(breaks = seq(0,6000,1000))+
  # scale_y_continuous(labels = scales::label_percent(accuracy = 1))+  # Required if plotting a percentage
  theme(
    axis.title.x = element_text(size = 14, face = "bold"),  # X-axis label size
    axis.title.y = element_text(size = 14, face = "bold"),  # Y-axis label size
    axis.text.x = element_text(size = 12),  # X-axis tick label size
    axis.text.y = element_text(size = 12),  # Y-axis tick label size
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # Title size and center alignment
    legend.text = element_text(size = 10),  # Legend text size
    legend.title = element_text(size = 10, face = "bold")  # Legend title size
  )

ggsave("combined_plot_CLustersUMIMice.png", plot = plot, width = 7, height = 7.5, units = "in", dpi = 300) #save image



# 3 PLOTS FOR Mt%

plot <- ggplot(QC_all, aes(x = MITO, y = Clusters, color = Dataset)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2.5) +
  theme_classic() +
  labs(
    title = "Number of DEGs After Filtering by Mt_% Cutoff",
    x = "mt%",
    y = "Percentage (Relative to Max)") +
  scale_x_continuous(breaks = seq(0,100,25))+
  # scale_y_continuous(labels = scales::label_percent(accuracy = 1))+
  theme(
    axis.title.x = element_text(size = 14, face = "bold"),  # X-axis label size
    axis.title.y = element_text(size = 14, face = "bold"),  # Y-axis label size
    axis.text.x = element_text(size = 12),  # X-axis tick label size
    axis.text.y = element_text(size = 12),  # Y-axis tick label size
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # Title size and center alignment
    legend.text = element_text(size = 10),  # Legend text size
    legend.title = element_text(size = 10, face = "bold")  # Legend title size
  )

ggsave("combined_plot_CLustersMITOMice.png", plot = plot, width = 7, height = 7.5, units = "in", dpi = 300)