adding Mean Expression_P1_vs_P2 - Filtering and Visualization_Patients

1. Load Libraries

2. load seurat object

#Load Seurat Object L7
load("../../0-robj/5-Harmony_Integrated_All_samples_Merged_CD4Tcells_final_Resolution_Selected_0.8_ADT_Normalized_cleaned_mt.robj")



All_samples_Merged

An object of class Seurat 
62900 features across 49305 samples within 6 assays 
Active assay: SCT (26176 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 5 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap, harmony

3. Load Data

3.1 Load CSV Files with Mean Expression

# Load the CSV files with mean expression data
P1_vs_P2 <- read.csv("comparison_P1_vs_P2.csv", row.names = 1)

4. Compute Mean Expression for Groups

library(Seurat)
library(dplyr)

# Define cell groups
group1_cells <- WhichCells(All_samples_Merged, idents = c("5", "1", "9"))  # Malignant CD4+ T cells
group2_cells <- WhichCells(All_samples_Merged, idents = c("2", "6", "8"))  # Normal CD4+ T cells

# Extract normalized expression values for both assays
expression_data_SCT <- GetAssayData(All_samples_Merged, assay = "SCT", layer = "data")  # SCT-normalized counts (not Pearson residuals)

# Function to calculate mean expression for each group
calculate_mean_expression <- function(markers, group1_cells, group2_cells, expression_data) {
  group1_mean <- rowMeans(expression_data[, group1_cells, drop = FALSE], na.rm = TRUE)
  group2_mean <- rowMeans(expression_data[, group2_cells, drop = FALSE], na.rm = TRUE)
  
  markers <- markers %>%
    rownames_to_column("gene") %>%
    mutate(mean_expr_group1 = group1_mean[gene],
           mean_expr_group2 = group2_mean[gene])
  
  return(markers)
}

# Compute mean expression P1_vs_P2
P1_vs_P2 <- calculate_mean_expression(P1_vs_P2, group1_cells, group2_cells, expression_data_SCT)

# Save updated results with mean expression

write.csv(P1_vs_P2, file = "comparison_P1_vs_P2_with_mean_expression.csv", row.names = TRUE)

5. Summarize Results

5.1 Summary Function

summarize_markers <- function(markers) {
  num_pval0 <- sum(markers$p_val_adj == 0)
  num_pval1 <- sum(markers$p_val_adj == 1)
  num_significant <- sum(markers$p_val_adj < 0.05)
  num_upregulated <- sum(markers$avg_log2FC > 1)
  num_downregulated <- sum(markers$avg_log2FC < -1)
  
  cat("Number of genes with p_val_adj = 0:", num_pval0, "\n")
  cat("Number of genes with p_val_adj = 1:", num_pval1, "\n")
  cat("Number of significant genes (p_val_adj < 0.05):", num_significant, "\n")
  cat("Number of upregulated genes (avg_log2FC > 1):", num_upregulated, "\n")
  cat("Number of downregulated genes (avg_log2FC < -1):", num_downregulated, "\n")
}

5.2 Summarize Markers1 (P1_vs_P2)

cat("Markers1 Summary (P1_vs_P2):\n")

Markers1 Summary (P1_vs_P2):

summarize_markers(P1_vs_P2)

Number of genes with p_val_adj = 0: 0 
Number of genes with p_val_adj = 1: 5134 
Number of significant genes (p_val_adj < 0.05): 7299 
Number of upregulated genes (avg_log2FC > 1): 2486 
Number of downregulated genes (avg_log2FC < -1): 2772 

6. Filter and Summarize Results

6.1 Apply Expression Filter First

# Apply the expression filter first
P1_vs_P2 <- P1_vs_P2 %>%
  filter(!(mean_expr_group1 < 0.2 & mean_expr_group2 < 0.2))

6.2 Summary Function

summarize_markers <- function(markers) {
  num_pval0 <- sum(markers$p_val_adj == 0)
  num_pval1 <- sum(markers$p_val_adj == 1)
  num_significant <- sum(markers$p_val_adj < 0.05)
  num_upregulated <- sum(markers$avg_log2FC > 1)
  num_downregulated <- sum(markers$avg_log2FC < -1)
  
  cat("Number of genes with p_val_adj = 0:", num_pval0, "\n")
  cat("Number of genes with p_val_adj = 1:", num_pval1, "\n")
  cat("Number of significant genes (p_val_adj < 0.05):", num_significant, "\n")
  cat("Number of upregulated genes (avg_log2FC > 1):", num_upregulated, "\n")
  cat("Number of downregulated genes (avg_log2FC < -1):", num_downregulated, "\n")
}

5.3 Summarize Markers1 (P1_vs_P2)

cat("Markers1 Summary (P1_vs_P2):\n")

Markers1 Summary (P1_vs_P2):

summarize_markers(P1_vs_P2)

Number of genes with p_val_adj = 0: 0 
Number of genes with p_val_adj = 1: 1347 
Number of significant genes (p_val_adj < 0.05): 4015 
Number of upregulated genes (avg_log2FC > 1): 722 
Number of downregulated genes (avg_log2FC < -1): 756