FC Scanner files for Analysis-part2

1. load libraries

#Differential Expression Analysis

2. load seurat object and save metadata file for samples

#Load Seurat Object L7
load("/home/bioinfo/0-imp_Robj/Harmony_integrated_All_samples_Merged_with_PBMC10x_with_harmony_clustering.Robj")
# sct_data <- GetAssayData(All_samples_Merged, assay = "SCT", layer = "data")

# memory.limit(size = 64000)
# 
# # Transpose the data so that cells are rows and genes are columns
# transposed_data <- t(as.data.frame(sct_data))
# 
# # Specify the file name and save as CSV
# write.csv(transposed_data, file = "table/SCT_data_All_samples_Merged_transposed.csv", row.names = TRUE)
# 
# 
# 
# 
# 
# # Extract metadata from Seurat object
# metadata <- All_samples_Merged@meta.data
# 
# # Write metadata to CSV
# write.csv(metadata, file = "Extra/Metadata_All_samples_Merged.csv", row.names = TRUE)

2. create barplot of CD4 cells in all samples

# Load necessary libraries
library(Seurat)
library(dplyr)
library(ggplot2)

# Load your CSV file
data <- read.csv("Extra/Step_by_step_DE_Analysis/1-Metadata_All_samples_Merged/All_Samples_metadata/Metadata_All_samples_Merged.csv")

# Filter rows where all three predicted columns contain "CD4 T"
filtered_data <- data %>%
  filter(grepl("CD4", predicted.celltype.l1) & 
         grepl("CD4", predicted.celltype.l2) & 
         grepl("CD4", predicted.celltype.l3))

# Count the number of CD4 T cells in each orig.ident group
cd4_t_counts <- filtered_data %>%
  group_by(orig.ident) %>%
  summarise(count = n())

# Print the counts
print(cd4_t_counts)

# Define colors for the samples (optional)
sample_colors <- c("L1" = "#1f77b4", "L2" = "#ff7f0e", "L3" = "#2ca02c", "L4" = "#d62728", 
                   "L5" = "#9467bd", "L6" = "#8c564b", "L7" = "#e377c2", "PBMC" = "#7f7f7f", 
                   "PBMC10x" = "#bcbd22")

# Create the bar plot
celltype_plot <- ggplot(cd4_t_counts, aes(x = orig.ident, y = count, fill = orig.ident)) +
  geom_col() +
  theme_classic() +
  geom_text(aes(label = count), 
            position = position_dodge(width = 0.9), 
            vjust = -0.25) +
  scale_fill_manual(values = sample_colors) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5)) +
  ggtitle("Count of CD4  Cells by Sample") +
  xlab("Sample") +
  ylab("Number of CD4  Cells")

# Print the plot
print(celltype_plot)

2. create PBMC CD4 cells (Control) file

# Load necessary libraries
library(Seurat)
library(dplyr)
library(ggplot2)

# Load your CSV file
data <- read.csv("Extra/NewFiles/PBMC/pbmc_METADATA.csv")

# Filter rows where all three predicted columns contain "CD4 T"
filtered_data <- data %>%
  filter(grepl("CD4", predicted.celltype.l1) & 
         grepl("CD4", predicted.celltype.l2) & 
         grepl("CD4", predicted.celltype.l3))

write.csv(filtered_data, "Extra/Step_by_step_DE_Analysis/1-Metadata_All_samples_Merged/Control_CD4Tcells/CD4_cells_PBMC_Control.csv", row.names = FALSE)

# Save the first column (PBMC cells) to a txt file without header
write.table(filtered_data[, 1], "Extra/Step_by_step_DE_Analysis/1-Metadata_All_samples_Merged/Control_CD4Tcells/CD4_Control.txt", row.names = FALSE, col.names = FALSE, quote = FALSE)

#Differential Expression Analysis

3. FC SCanner for DE

library(foreach)
library(doParallel)

setwd("/home/bioinfo/17-SingleCellFCscanner/")

source("scFCscanner_028_in_list.r")

Script to calculate logFC, Filter logFC >2 & logfc <-1.5, exclude lines when mean_L1_L7 <0.2 & control_group <0.2

# # Load necessary libraries
# library(dplyr)
# library(readr)
# 
# # Define the file processing function
# process_file <- function(file_path, output_dir) {
#   # Read the TSV file
#   Exp_Allsample <- read_tsv(file_path)
#   
#   # Check if the first column name is missing, and if so, set it to "gene"
#   if (colnames(Exp_Allsample)[1] == "...1") {
#     colnames(Exp_Allsample)[1] <- "gene"
#   }
#   
#   # Sanitize column names to handle spaces
#   colnames(Exp_Allsample) <- make.names(colnames(Exp_Allsample), unique = TRUE)
#   
#   # Identify the FC and mean columns dynamically
#   fc_column <- grep("^FC\\.", colnames(Exp_Allsample), value = TRUE)  # FC column with spaces replaced by dots
#   mean_columns <- grep("^mean\\.", colnames(Exp_Allsample), value = TRUE)  # Mean columns with spaces replaced by dots
#   
#   # Ensure the required columns are present
#   if (length(fc_column) != 1 || length(mean_columns) != 2) {
#     stop(paste("Required columns not found in file:", file_path))
#   }
#   
#   # Dynamically assign mean column names for clarity
#   mean_P2 <- mean_columns[1]
#   mean_P3 <- mean_columns[2]
#   
#   # Calculate log2 fold-change directly within the pipe
#   Exp_Allsample <- Exp_Allsample %>%
#     mutate(log2FC = log2(!!sym(fc_column)))  # Add log2FC as a new column
#   
#   # Define thresholds
#   threshold_positive <- 1.5    # Positive log2FC threshold
#   threshold_negative <- -1 # Negative log2FC threshold
#   
#   # Filter rows based on log2FC and mean thresholds
#   filtered_data_final <- Exp_Allsample %>%
#     filter(log2FC > threshold_positive | log2FC < threshold_negative) %>%
#     filter(!(!!sym(mean_P2) < 0.2 & !!sym(mean_P3) < 0.2))
#   
#   # Generate output file path
#   output_file <- file.path(output_dir, paste0("2-filtered_", basename(file_path)))
#   
#   # Write the filtered data to a CSV file
#   write.csv(filtered_data_final, output_file, row.names = FALSE)
#   
#   # Return the output file path for logging
#   return(output_file)
# }
# 
# # Set the input directory, output directory, and file list
# input_dir <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/"
# output_dir <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/Filtered/"
# file_list <- list.files(input_dir, pattern = "\\.csv$", full.names = TRUE)
# 
# # Create the output directory if it doesn't exist
# if (!dir.exists(output_dir)) {
#   dir.create(output_dir, recursive = TRUE)
# }
# 
# # Apply the process_file function to each file
# output_files <- lapply(file_list, process_file, output_dir = output_dir)
# 
# # Print completion message
# cat("Processing complete. Filtered files saved in the directory:", output_dir, "\n")

Script to calculate logFC, Filter logFC >2 & logfc <-1.5, exclude lines when mean_L1_L7 <0.2 & control_group <0.2

# # Load necessary libraries
# library(dplyr)
# library(readr)
# 
# # Enhanced file processing function
# process_file <- function(file_path, output_dir) {
#   tryCatch({
#     # Read the TSV file with explicit column types
#     Exp_Allsample <- read_tsv(file_path, col_types = cols(.default = col_guess()))
#     
#     # Check if the first column name is missing, and if so, set it to "gene"
#     if (colnames(Exp_Allsample)[1] == "...1") {
#       colnames(Exp_Allsample)[1] <- "gene"
#     }
#     
#     # Sanitize column names to handle spaces and special characters
#     colnames(Exp_Allsample) <- make.names(colnames(Exp_Allsample), unique = TRUE)
#     
#     # Identify the FC column dynamically
#     fc_column <- grep("^FC\\.", colnames(Exp_Allsample), value = TRUE)  # FC column with spaces replaced by dots
#     
#     # Dynamically identify mean columns by matching specific patterns (e.g., columns starting with "mean")
#     mean_columns <- grep("^mean\\.", colnames(Exp_Allsample), value = TRUE)  # Mean columns with spaces replaced by dots
#     
#     # If the mean column names are inconsistent, use alternative logic
#     if (length(mean_columns) < 2) {
#       # Attempt to identify mean columns by finding numeric columns with distinct naming patterns
#       numeric_columns <- colnames(Exp_Allsample)[sapply(Exp_Allsample, is.numeric)]
#       mean_columns <- numeric_columns[grep("Mean|mean|Average|avg", numeric_columns, ignore.case = TRUE)]
#     }
#     
#     # Ensure the required columns are present
#     if (length(fc_column) != 1 || length(mean_columns) < 2) {
#       stop(paste("Required columns not found in file:", file_path))
#     }
#     
#     # Select the first two mean columns dynamically
#     mean_P2 <- mean_columns[1]
#     mean_P3 <- mean_columns[2]
#     
#     # Calculate log2 fold-change directly within the pipe
#     Exp_Allsample <- Exp_Allsample %>%
#       mutate(log2FC = log2(!!sym(fc_column)))  # Add log2FC as a new column
#     
#     # Define thresholds
#     threshold_positive <- 1.5    # Positive log2FC threshold
#     threshold_negative <- -1 # Negative log2FC threshold
#     
#     # Filter rows based on log2FC and mean thresholds
#     filtered_data_final <- Exp_Allsample %>%
#       filter(log2FC > threshold_positive | log2FC < threshold_negative) %>%
#       filter(!(!!sym(mean_P2) < 0.2 & !!sym(mean_P3) < 0.2))
#     
#     # Extract the portion of the filename after the prefix "tab_" (if present)
#     base_filename <- basename(file_path)
#     new_filename <- if (grepl("tab_", base_filename)) {
#       sub("^.*tab_", "", base_filename)
#     } else {
#       base_filename
#     }
#     
#     # Generate output file path
#     output_file <- file.path(output_dir, paste0("2-filtered_", new_filename))
#     
#     # Write the filtered data to a CSV file
#     write.csv(filtered_data_final, output_file, row.names = FALSE)
#     
#     # Return the output file path for logging
#     return(output_file)
#   }, error = function(e) {
#     # Log error message
#     message(paste("Error processing file:", file_path, "\n", e))
#     return(NULL)
#   })
# }
# 
# # Set the input directory, output directory, and file list
# input_dir <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/"
# output_dir <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/Filtered2/"
# file_list <- list.files(input_dir, pattern = "\\.csv$", full.names = TRUE)
# 
# # Create the output directory if it doesn't exist
# if (!dir.exists(output_dir)) {
#   dir.create(output_dir, recursive = TRUE)
# }
# 
# # Apply the process_file function to each file and log any errors
# output_files <- lapply(file_list, function(file) {
#   message("Processing file:", file)
#   process_file(file, output_dir)
# })
# 
# # Print completion message
# cat("Processing complete. Filtered files saved in the directory:", output_dir, "\n")
# 
# 

Barplot for up and down regulated genes

# Load necessary libraries
library(readr)
library(dplyr)
library(stringr)

# Define file paths
input_file <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/Filtered2/2-filtered_Cluster19_cells_vs_CD4_Control.csv"
output_folder <- "Comparisons_Postprocessing/Step_by_step_DE_Analysis/3-FC-Comparison_to_Control/6-All_Files_Together_for_Filtering/TEST/Filtered2/Enrichment_Analysis_Files_for_Autocompare/C19_cells_vs_CD4_Control/"

# Ensure output folder exists
if (!dir.exists(output_folder)) {
  dir.create(output_folder, recursive = TRUE)
}

data <- tryCatch({
  read.csv(input_file, stringsAsFactors = FALSE)
}, error = function(e) {
  stop("Error reading input file. Check the file path or format.")
})



# Check if the necessary columns exist
if (!all(c("log2FC", "gene") %in% colnames(data))) {
  stop("Required columns ('log2FC' and 'gene') are missing in the data.")
}

# Choose your own log2FC threshold for both positive and negative values
threshold_positive <- 3.5 # Set your chosen threshold for positive log2FC
threshold_negative <- -1  # Set your chosen threshold for negative log2FC

# Filter the data based on log2FC criteria
filtered_data <- data %>%
  filter(log2FC > threshold_positive | log2FC < threshold_negative)

# Filter for upregulated and downregulated genes based on the filtered data
upregulated_genes <- filtered_data %>% filter(log2FC > threshold_positive)
downregulated_genes <- filtered_data %>% filter(log2FC < threshold_negative)

# Extract only gene names as vectors
upregulated_gene_names <- upregulated_genes %>% pull(gene)
downregulated_gene_names <- downregulated_genes %>% pull(gene)

# Extract base name for the input file (without extension) and remove the "2-filtered_" prefix
base_name <- tools::file_path_sans_ext(basename(input_file))  # Extract base name without extension
base_name <- gsub("^2-filtered_", "", base_name)  # Remove "2-filtered_" prefix

# Create output file paths with dynamic names
upregulated_file <- file.path(output_folder, paste0(base_name, "_upregulated_genes.txt"))
downregulated_file <- file.path(output_folder, paste0(base_name, "_downregulated_genes.txt"))

# Save the gene names to text files
write_lines(upregulated_gene_names, upregulated_file)
write_lines(downregulated_gene_names, downregulated_file)

# Print the number of upregulated and downregulated genes
cat("Number of upregulated genes:", length(upregulated_gene_names), "\n")
cat("Number of downregulated genes:", length(downregulated_gene_names), "\n")
cat("Gene names saved to:\n")
cat("  Upregulated genes: ", upregulated_file, "\n")
cat("  Downregulated genes: ", downregulated_file, "\n")