ind4_replication.Rmd

---
title: "paper replication"
author: "Grigoris Ntoulaveris"
date: "`r Sys.Date()`"
output: html_document
---

The following is an attempt at replicating the results of the study https://europepmc.org/article/MED/29795293#abstract, specifically for the scRNA seq data of the "Individual 4" dataset. Differences with the paper's methods can be found only in the clustering section. Besides density based clustering (DBSCAN), which corresponded to the paper's methods, two more clustering approaches are presented (GMMs, SNN-Louvain). A comparison with the paper's results concering the marker genes of interest ensues the two new algorithms.


```{r}
source("scRNA_R_pipeline.R")

```


# Data exploration


```{r}

file_path_id4 <- "data/GSM3099846_Ind4_Expression_Matrix.txt.gz"
data <- fread(file_path_id4, sep = "\t", header = TRUE, data.table = TRUE)
data <- transpose(data, keep.names = "cells", make.names = "V1")

dataset_name <- "GSM3099846_Ind4_Expression_Matrix"


```


```{r}
#melted_data <- melt_dataset(data)

```

```{r}

#library_size <- get_library_size(melted_data)

```

```{r}

#plot_GeneExpression_heatmap(melted_data)


```


# Data preprocessing


```{r}
data_copy <- data.frame(data)

```

```{r}
data_seurat <- build_seurat_object(data_copy, dataset_name)

```


# Quality Control


```{r}
data_seurat <- find_mitochondrial_genes(data_seurat)
print(data_seurat@meta.data)

```

In the following plots the nCount_RNA corresponds to the total number of RNA molecules detected in a cell, which can be considered as a measure of the cell's RNA content or overall expression level. nFeature_RNA represents the count of unique genes or features with non-zero expression in a cell and is an indicator of the diversity of expressed genes in a cell. percent_mito corresponds to the percentage of mitochondrial genes that are expressed in a cell and a high number is indicative of cells of low-quality.


```{r}
plot_Features_violin(data_seurat)

```

In the following plot the total number of RNA molecules and the total number of unique genes for each cell are plotted together, to ensure that the previous visualization is not misleading in terms of cell quality. Good quality cells usually have a high value in both metrics. If many cells are clustered towards the bottom right of the plot it would mean that only a few number of unique genes were detected and those are sequenced continuously providing a misleading high number of transcripts. If many cells are clustered towards the top left of the plot it would mean that the sequencer had discovered many unique genes but they wouldn't be deeply sequenced enough to provide meaningful results.


```{r}
plot_Features_scatter(data_seurat)

```



# Filtering

In the filtering step low-quality cells are removed from the dataset, based on the results of the quality control.

```{r}
data_seurat <- filter_cells(data_seurat, 
                             nFeature_RNA_min = 500, nFeature_RNA_max = 6000, percent_mito_max = 10)

data_seurat
```



# Normalization
In order to be able to compare the levels of expression across multiple cells the data need to be normalized. The gene expression measurements for each cell is divided by the total expression of all cells, is then multiplied by a scaling factor. The result is also transformed into log space.

```{r}
data_seurat <- NormalizeData(data_seurat, normalization.method = "LogNormalize")

```



# Regress Out Normalization

```{r}
# Access the RNA assay
assay_data <- data_seurat@assays$RNA

# Retrieve the counts matrix
counts_matrix <- assay_data@counts

# Find the UMI values
umi_values <- colSums(counts_matrix)

```



```{r}

data_seurat <- perform_regress_out(data_seurat)

```



# Highly variable genes identification


```{r}

# find all high variable genes
data_seurat <- FindPlot_variable_genes(data_seurat)

```



# PCA


```{r}
data_seurat <- reduce_dimensions(data_seurat, "pca")

```



# Save - Load seurat

```{r}
saveRDS(data_seurat, file="ind4_seurat")

```


```{r}
data_seurat <- readRDS(file = "ind4_seurat")

```



# Clustering

## GMM clustering

```{r}
clustering_results <- GMM_clustering(data_seurat, "pca")

data_seurat <- clustering_results[[1]]
gmm_model <- clustering_results[[2]]
cell_embeddings <- clustering_results[[3]]

```


```{r}
data_seurat@meta.data

```

```{r}
# with tSNE as projection for plots
clustering_plots <- plot_clustering(dataset_name, data_seurat, gmm_model)

```


```{r}
# Find marker genes for all clusters
marker_genes <- FindAllMarkers(data_seurat, only.pos = TRUE, random.seed = 42)

# Print the marker genes for each cluster
print(marker_genes)


```


```{r}

# Group marker_genes by cluster and select the top 10 genes based on avg_log2FC
top_genes <- marker_genes %>%
  group_by(cluster) %>%
  top_n(10, avg_log2FC) %>%
  ungroup()


```


```{r}

# Specify the genes to include in the new data frame
genes_to_include <- c("LTF", "SERPINB4", "SERPINB3", "WFDC2", "LCN2", "BTG1", "CLDN4", "ANXA1", "HMGA1", "S100A2", "TIMP1", "MMP3", "TAGLN", "ACTA2", "ACTG2", "EIF5A", "CAV1", "VIM", "STC2", "AGR2", "AREG", "TNFSF10", "SERPINA1", "PIP", "APOD", "SFRP4", "IGFBP7", "GNG11", "ANGPT2", "SERPINE1", "HSPA6", "ZNF90", "CORO1A")

# Filter the top_genes data frame based on the specified genes
top_genes_small <- top_genes[top_genes$gene %in% genes_to_include, ]

# Convert the gene column to a factor with the desired order
top_genes_small$gene <- factor(top_genes_small$gene, levels = genes_to_include)


ggplot(top_genes_small, aes(x = cluster, y = gene, fill = avg_log2FC)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Cluster", y = "Gene")


```




## Density based clustering


```{r}
data_seurat <- perform_dbscan_clustering(data_seurat, "pca")

```

```{r}
DimPlot(data_seurat, group.by = "cluster_label")

```

```{r}
# only with PCA
clustering_plots <- plot_clustering(dataset_name, data_seurat, gmm_model)

```




## KNN graphs clustering

```{r}

data_seurat <- FindNeighbors(data_seurat, dims = 1:10)
data_seurat <- FindClusters(data_seurat, resolution = 0.5, algorithm = 2, random.seed = 42)

```

```{r}
data_seurat <- RunTSNE(data_seurat, dims = 1:10)


```

```{r}
DimPlot(data_seurat, reduction = "tsne")

```

Next, we find the marker genes that were specified in the paper inside the clusters that were created with Louvain. After that, based again on the paper's results, according to their marker genes, some clusters will be combined into one to reveal the three main cell types.


```{r}
# Find marker genes for all clusters
marker_genes <- FindAllMarkers(data_seurat, only.pos = TRUE, random.seed = 42)

# Print the marker genes for each cluster
print(marker_genes)


```

```{r}

# Group marker_genes by cluster and select the top 10 genes based on avg_log2FC
top_genes <- marker_genes %>%
  group_by(cluster) %>%
  top_n(10, avg_log2FC) %>%
  ungroup()


```

```{r}
# Convert top_genes to a data.table
top_genes_dt <- as.data.table(top_genes)

# Reshape the data table for heatmap plotting
heatmap_data <- dcast(top_genes_dt, gene ~ cluster, value.var = "avg_log2FC")

# Reorder the rows based on gene names
heatmap_data <- heatmap_data[order(gene)]

# Generate the heatmap
ggplot(top_genes_dt, aes(x = cluster, y = gene, fill = avg_log2FC)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Cluster", y = "Gene")

```


```{r}

# Specify the genes to include in the new data frame
genes_to_include <- c("LTF", "SERPINB4", "SERPINB3", "WFDC2", "LCN2", "BTG1", "CLDN4", "ANXA1", "HMGA1", "S100A2", "TIMP1", "MMP3", "TAGLN", "ACTA2", "ACTG2", "EIF5A", "CAV1", "VIM", "STC2", "AGR2", "AREG", "TNFSF10", "SERPINA1", "PIP", "APOD", "SFRP4", "IGFBP7", "GNG11", "ANGPT2", "SERPINE1", "HSPA6", "ZNF90", "CORO1A")

# Filter the top_genes data frame based on the specified genes
top_genes_small <- top_genes[top_genes$gene %in% genes_to_include, ]

# Convert the gene column to a factor with the desired order
top_genes_small$gene <- factor(top_genes_small$gene, levels = genes_to_include)


ggplot(top_genes_small, aes(x = cluster, y = gene, fill = avg_log2FC)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Cluster", y = "Gene") +
  ggtitle("Marker genes of interest")


```

```{r}

# Find genes not in top_genes
genes_not_found <- setdiff(genes_to_include, top_genes_small$gene)

# Find the cluster for genes not found in top_genes from marker_genes
missing_genes_clusters <- marker_genes[marker_genes$gene %in% genes_not_found, c("gene", "cluster")]

# Print the genes and their corresponding clusters
cat(genes_not_found)
print(missing_genes_clusters)


```



### Comparison with paper

Finally, the classified genes for Individual 4 listed in the paper's Supplementary material are compared with our own classified genes to access similarity of results.

```{r}


excel_file <- "genes_labels_ind4_7.xlsx"
paper_data <- read_excel(excel_file)

filtered_data <- paper_data[paper_data$individual == "I4", ]

# Create a new data frame with updated labels
updated_marker_genes <- marker_genes

# Update labels in the "cluster" column
updated_marker_genes$cluster <- ifelse(updated_marker_genes$cluster %in% c(0, 1), "L1",
                                       ifelse(updated_marker_genes$cluster %in% c(2, 3, 5), "B",
                                              ifelse(updated_marker_genes$cluster %in% c(4, 6), "L2", "X")))

updated_paper_data <- paper_data

updated_paper_data$cluster <- ifelse(updated_paper_data$cluster %in% c("Basal", "Basal_Myoepithelial"), "B",
                                       ifelse(updated_paper_data$cluster %in% c("Luminal_1_1", "Luminal_1_2"), "L1",
                                              ifelse(updated_paper_data$cluster %in% c("Luminal_2"), "L2", "X")))


```


```{r}

# Create a scatter plot to compare cluster labels
ggplot(comparison_df, aes(x = cluster_paper, y = cluster_marker)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(x = "Cluster Label (Paper Data)", y = "Cluster Label (Marker Genes)") +
  ggtitle("Comparison of Cluster Labels")


```



```{r}

# Find common genes between updated_paper_data and updated_marker_genes
common_genes <- intersect(updated_paper_data$gene, updated_marker_genes$gene)

# Filter updated_paper_data and updated_marker_genes for common genes
common_genes_paper_data <- updated_paper_data[updated_paper_data$gene %in% common_genes, ]
common_genes_marker_genes <- updated_marker_genes[updated_marker_genes$gene %in% common_genes, ]

# Merge the two data frames to compare cluster labels
comparison_df <- merge(common_genes_paper_data, common_genes_marker_genes, by = "gene", suffixes = c("_paper", "_marker"))

# Calculate the number of differences in cluster labels
comparison_df$label_diff <- ifelse(comparison_df$cluster_paper != comparison_df$cluster_marker, 1, 0)

# Calculate the number of differences in cluster labels
num_differences <- sum(comparison_df$cluster_paper != comparison_df$cluster_marker)
total_gene_entries <- nrow(comparison_df)
percentage_of_differences <- round((num_differences / total_gene_entries)*100, 2)

# Print the number of differences
#print(paste("Number of Differences:", num_differences))

ggplot(comparison_df, aes(x = gene, fill = factor(label_diff))) +
  geom_bar() +
  scale_fill_manual(values = c("0" = "green", "1" = "red"), labels = c("No Difference", "Difference")) +
  theme_minimal() +
  labs(x = "Genes", y = "Number of Differences", fill= "") +
  ggtitle(sprintf("Misclassified genes according to the paper's data\n(%.2f%% of difference)", percentage_of_differences)) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

```





# Save the labeled results in a csv

Save the expression matrix with the cluster labels of only the top 2000 highly variable genes in an attempt to evade the curse of dimensionality (and the excessive overfitting) in the classification process.

```{r}
var_genes <- VariableFeatures(data_seurat)
top_var_genes <- head(var_genes, 2000)


```

```{r}
expression_matrix <- t(GetAssayData(data_seurat, assay = "RNA")[top_var_genes, ])

expression_df <- as.data.frame(expression_matrix)
expression_dt <- data.table(expression_df)

cluster_labels <- Idents(data_seurat)
cluster_labels_dt <- data.table(cluster_labels)

joined_dt <- data.table(expression_dt, cluster = cluster_labels_dt)
setnames(joined_dt, old = "cluster.cluster_labels", new = "cluster_label")

write.csv(joined_dt, file=gzfile("ind4_var_genes_with_clusters.csv.gz"), row.names = FALSE)

```


Save the initial expression matrix but with labeled cells
```{r}
#install.packages("zip")
library(zip)

expression_matrix <- t(GetAssayData(data_seurat, assay = "RNA"))
expression_df <- as.data.frame(expression_matrix)
expression_dt <- data.table(expression_df)

cluster_labels <- Idents(data_seurat)
cluster_labels_dt <- data.table(cluster_labels)

joined_dt <- data.table(expression_dt, cluster = cluster_labels_dt)
setnames(joined_dt, old = "cluster.cluster_labels", new = "cluster_label")

# Save the data frame as a CSV file
#write.csv(joined_dt, "ind4_exp_matrix_with_clusters.csv", row.names = FALSE)
write.csv(joined_dt, file=gzfile("ind4_exp_matrix_with_clusters.csv.gz"), row.names = FALSE)

```






Save all PCs

```{r}
# Extract the PCA embeddings
pca_embeddings <- data_seurat@reductions$pca@cell.embeddings

# Get the cluster labels for each cell
cluster_labels <- Idents(data_seurat)
cluster_labels_dt <- data.table(cluster_labels)

# Create a data frame with PCA embeddings and cluster labels
pca_data <- data.table(pca_embeddings)

# Add the cluster labels as a new column
joined_dt <- data.table(pca_data, cluster = cluster_labels_dt)
setnames(joined_dt, old = "cluster.cluster_labels", new = "cluster_label")


# Save the data frame as a CSV file
write.csv(joined_dt, "ind4_pca_embeddings_with_clusters.csv", row.names = FALSE)


```

Save only the first 10 PCs
```{r}
# Extract the PCA embeddings
pca_embeddings <- data_seurat@reductions$pca@cell.embeddings[, 1:10]

# Get the cluster labels for each cell
cluster_labels <- Idents(data_seurat)
cluster_labels_dt <- data.table(cluster_labels)

# Create a data frame with PCA embeddings and cluster labels
pca_data <- data.table(pca_embeddings)

# Add the cluster labels as a new column
joined_dt <- data.table(pca_data, cluster = cluster_labels_dt)
setnames(joined_dt, old = "cluster.cluster_labels", new = "cluster_label")


# Save the data frame as a CSV file
write.csv(joined_dt, "ind4_pca_embeddings_only_10PCs.csv", row.names = FALSE)


```