Phylogeny_based on RNA Expression

1. load libraries

2. Load Seurat Object

load("/home/bioinfo/Cluster_to_Computer_Transfer_files_folder/All_Normal-PBMC_Abnormal-cellLines_T_cells_Merged_Annotated_UMAP_on_Clusters_to_USE.Robj")

3. Phylogeney of RNA based Expression-eucleadian

DefaultAssay(All_samples_Merged) <- 'SCT'

# 1. Calculate average expression for each gene in each cell line
avg_expression <- AggregateExpression(All_samples_Merged, 
                                      group.by = "cell_line", 
                                      assays = "SCT", 
                                      slot = "data",
                                      return.seurat = FALSE)$SCT

# 2. Select top variable genes
var_genes <- VariableFeatures(All_samples_Merged)
avg_expression_var <- avg_expression[var_genes,]

# Convert to matrix and handle potential issues
avg_expression_mat <- as.matrix(avg_expression_var)
print(paste("Dimensions of expression matrix:", paste(dim(avg_expression_mat), collapse = "x")))

[1] "Dimensions of expression matrix: 3000x8"

# Remove any genes with zero variance
gene_var <- apply(avg_expression_mat, 1, var)
avg_expression_mat <- avg_expression_mat[gene_var > 0, ]
print(paste("Dimensions after removing zero variance genes:", paste(dim(avg_expression_mat), collapse = "x")))

[1] "Dimensions after removing zero variance genes: 3000x8"

# 3. Calculate distance matrix
dist_matrix <- dist(t(avg_expression_mat))

Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'

# Check distance matrix
print(paste("Any NA in distance matrix:", any(is.na(dist_matrix))))

[1] "Any NA in distance matrix: FALSE"

print(paste("Any infinite values in distance matrix:", any(is.infinite(dist_matrix))))

[1] "Any infinite values in distance matrix: FALSE"

print("Summary of distance values:")

[1] "Summary of distance values:"

print(summary(as.vector(dist_matrix)))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 592041 1136261 1515296 1702792 2290386 2801467 

# 4. Create phylogenetic tree
if (any(is.na(dist_matrix))) {
  print("Distance matrix contains NA values, cannot create phylogenetic tree")
  # Create a simple dendrogram instead
  hc <- hclust(dist_matrix)
  plot(hc, main="Dendrogram of Cell Lines Based on Gene Expression")
} else {
  tree <- nj(dist_matrix)
  plot(tree, main="Phylogenetic Tree of Cell Lines Based on Gene Expression")
}

# 5. Create a heatmap
pheatmap(avg_expression_mat, 
         cluster_rows = TRUE, 
         cluster_cols = TRUE, 
         show_rownames = FALSE,
         main = "Heatmap of Top Variable Genes Across Cell Lines")

NA
NA

4. Phylogeney of RNA based Expression-Pearson

# Check if SCT assay exists and perform SCTransform if needed
DefaultAssay(All_samples_Merged) <- 'SCT'


# 1. Verify that cell_line metadata exists
if (!"cell_line" %in% colnames(All_samples_Merged@meta.data)) {
  stop("cell_line column not found in metadata. Ensure each cell has a cell line label.")
}

# 2. Calculate average expression for each gene in each cell line
avg_expression <- AggregateExpression(All_samples_Merged, 
                                      group.by = "cell_line", 
                                      assays = "SCT", 
                                      slot = "data",
                                      return.seurat = FALSE)$SCT

# Ensure that the aggregation is correct
print(paste("Number of genes:", nrow(avg_expression)))

[1] "Number of genes: 25901"

print(paste("Number of cell lines:", ncol(avg_expression)))

[1] "Number of cell lines: 8"

# 3. Select top variable genes for phylogenetic analysis
var_genes <- VariableFeatures(All_samples_Merged)
avg_expression_var <- avg_expression[var_genes, ]

# Convert to matrix: genes as rows, cell lines as columns
avg_expression_mat <- as.matrix(avg_expression_var)
print(paste("Dimensions of expression matrix:", paste(dim(avg_expression_mat), collapse = "x")))

[1] "Dimensions of expression matrix: 3000x8"

# Remove genes with zero variance
gene_var <- apply(avg_expression_mat, 1, var)
avg_expression_mat <- avg_expression_mat[gene_var > 0, ]
print(paste("Dimensions after removing zero variance genes:", paste(dim(avg_expression_mat), collapse = "x")))

[1] "Dimensions after removing zero variance genes: 3000x8"

# 4. Calculate Pearson correlation distance between **cell lines**
# Transpose the matrix so that cell lines are compared
cor_matrix <- cor(avg_expression_mat, method = "pearson")
cor_dist <- as.dist(1 - cor_matrix)

Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by 'BiocGenerics'

# Check distance matrix
print(paste("Any NA in distance matrix:", any(is.na(cor_dist))))

[1] "Any NA in distance matrix: FALSE"

print(paste("Any infinite values in distance matrix:", any(is.infinite(cor_dist))))

[1] "Any infinite values in distance matrix: FALSE"

print("Summary of distance values:")

[1] "Summary of distance values:"

print(summary(as.vector(cor_dist)))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02325 0.08369 0.13251 0.13627 0.19588 0.24904 

# 5. Create phylogenetic tree using NJ or dendrogram as fallback
library(ape)  # Ensure ape is loaded for nj()
if (any(is.na(cor_dist))) {
  print("Distance matrix contains NA values, cannot create phylogenetic tree")
  # Create a dendrogram as fallback
  hc <- hclust(cor_dist)
  plot(hc, main="Dendrogram of Cell Lines Based on Correlation Distance")
} else {
  tree <- nj(cor_dist)
  plot(tree, main="Phylogenetic Tree of Cell Lines Based on Correlation Distance")
}

# 6. Create heatmap for visualization (ensure it clusters cell lines)
library(pheatmap)
pheatmap(avg_expression_mat, 
         cluster_rows = FALSE,    # Do not cluster genes
         cluster_cols = TRUE,     # Cluster cell lines (columns)
         show_rownames = FALSE,
         main = "Heatmap of Top Variable Genes Across Cell Lines")

NA
NA

6. Save the Seurat object as an Robj file

#save(All_samples_Merged, file = "../5-SS_ScRNA_Data_Analysis/4-ScSS_MyAnalysis_on_SS/0-Important_R_OBJ/All_samples_Merged_WNN_correct_on_HPC.Robj")