Add unsupervised clustering of genomes

scripts/unsup_clust_genomes.R

@@ -0,0 +1,155 @@
+### Run after pangenome generation (script 02) to perform unsupervised clustering
+
+# Purpose: Does unsupervised clustering of the genomes show separation of light ecotypes?
+# Based on github issue: https://github.com/JRaviLab/cyano_adaptation/issues/2
+
+library(here)
+library(purrr)
+library(readr)
+library(pheatmap)
+library(vegan)   # For calculating Jaccard distance
+library(ape)
+library(ggplot2)
+library(uwot)
+library(dbscan)
+
+# Paths
+data_dir    <- here::here("data")
+results_dir <- here::here("results")
+plots_dir   <- here::here("plots")
+dir.create(results_dir, recursive = TRUE, showWarnings = FALSE)
+dir.create(plots_dir,   recursive = TRUE, showWarnings = FALSE)
+
+pa_matrix <- readr::read_tsv(
+  file.path(data_dir, "PanTA_All_NoAlignments", "gene_presence_absence.Rtab"),
+  show_col_types = FALSE
+)
+
+# Extract gene names
+gene_ids <- pa_matrix[[1]]
+
+# Extract numeric matrix (genes × genomes)
+mat_gene_genome <- as.matrix(pa_matrix[, -1])
+rownames(mat_gene_genome) <- gene_ids
+
+# Transpose → genomes × genes
+mat_genome_gene <- t(mat_gene_genome)
+
+# Convert to binary presence/absence
+mat_bin <- (mat_genome_gene > 0) * 1
+dim(mat_bin)   # genomes × genes
+
+min_present <- 3                 # present in at least 3 genomes
+max_frac    <- 0.98              # present in ≤ 98% of genomes
+gene_counts <- colSums(mat_bin)
+keep_genes <- which(
+  gene_counts >= min_present &
+  gene_counts <= nrow(mat_bin) * max_frac
+)
+mat_filt <- mat_bin[, keep_genes]
+dim(mat_filt)
+
+# Compute Jaccard distance
+dist_jaccard <- vegdist(mat_filt, method = "jaccard", binary = TRUE)
+
+# Ordination for visualization
+pcoa_res <- pcoa(dist_jaccard)
+pcoa_df <- data.frame(
+  Genome = rownames(mat_filt),
+  PC1 = pcoa_res$vectors[, 1],
+  PC2 = pcoa_res$vectors[, 2]
+)
+
+p_pcoa <- ggplot(pcoa_df, aes(PC1, PC2)) +
+  geom_point(size = 2) +
+  theme_classic() +
+  ggtitle("PCoA of PanTA gene-content (Jaccard)")
+ggsave(file.path(plots_dir, "unsup_pcoa_jaccard.png"), p_pcoa, width = 7, height = 5.5, dpi = 300)
+ggsave(file.path(plots_dir, "unsup_pcoa_jaccard.svg"), p_pcoa, width = 7, height = 5.5)
+
+umap_seed <- 42L
+set.seed(umap_seed)
+umap_res <- uwot::umap(mat_filt, n_neighbors = 15, min_dist = 0.1)
+umap_df <- data.frame(
+  Genome = rownames(mat_filt),
+  UMAP1 = umap_res[, 1],
+  UMAP2 = umap_res[, 2]
+)
+
+p_umap <- ggplot(umap_df, aes(UMAP1, UMAP2)) +
+  geom_point(size = 2) +
+  theme_classic() +
+  ggtitle("UMAP of PanTA gene-content")
+ggsave(file.path(plots_dir, "unsup_umap.png"), p_umap, width = 7, height = 5.5, dpi = 300)
+ggsave(file.path(plots_dir, "unsup_umap.svg"), p_umap, width = 7, height = 5.5)
+
+hc <- hclust(dist_jaccard, method = "ward.D2")
+png(file.path(plots_dir, "unsup_dendrogram.png"), width = 2400, height = 1800, res = 300)
+plot(hc, labels = FALSE)
+dev.off()
+
+k <- 6    # choose by dendrogram or silhouette
+clusters_hc <- cutree(hc, k = k)
+table(clusters_hc)
+
+hdb <- hdbscan(umap_res, minPts = 8)
+clusters_hdb <- hdb$cluster   # 0 = noise
+table(clusters_hdb)
+
+# Add to UMAP plot
+umap_df$cluster <- factor(clusters_hdb)
+p_umap_clust <- ggplot(umap_df, aes(UMAP1, UMAP2, color = cluster)) +
+  geom_point(size = 2) +
+  theme_classic() +
+  ggtitle("UMAP + HDBSCAN clusters")
+ggsave(file.path(plots_dir, "unsup_umap_hdbscan.png"), p_umap_clust, width = 7, height = 5.5, dpi = 300)
+ggsave(file.path(plots_dir, "unsup_umap_hdbscan.svg"), p_umap_clust, width = 7, height = 5.5)
+
+# Save cluster assignments
+cluster_df <- data.frame(
+  Genome          = rownames(mat_filt),
+  cluster_hclust  = clusters_hc,
+  cluster_hdbscan = clusters_hdb
+)
+readr::write_tsv(cluster_df, file.path(results_dir, "unsup_cluster_assignments.tsv"))
+
+# Take top 200 most variable genes
+gene_var <- apply(mat_filt, 2, var)
+top_genes <- names(sort(gene_var, decreasing = TRUE))[1:200]
+annotation_row <- data.frame(Cluster = factor(clusters_hc))
+rownames(annotation_row) <- rownames(mat_filt)
+
+png(file.path(plots_dir, "unsup_heatmap_top200.png"), width = 3600, height = 3600, res = 300)
+pheatmap(
+  mat_filt[, top_genes],
+  cluster_rows  = TRUE,
+  cluster_cols  = TRUE,
+  show_rownames = FALSE,
+  show_colnames = FALSE,
+  annotation_row = annotation_row
+)
+dev.off()
+
+cluster_labels <- clusters_hc
+genes <- colnames(mat_filt)
+enrich <- purrr::map(sort(unique(cluster_labels)), function(cl) {
+  purrr::map_dbl(purrr::set_names(genes), function(g) {
+    present_in  <- sum(mat_filt[cluster_labels == cl, g] == 1)
+    absent_in   <- sum(mat_filt[cluster_labels == cl, g] == 0)
+    present_out <- sum(mat_filt[cluster_labels != cl, g] == 1)
+    absent_out  <- sum(mat_filt[cluster_labels != cl, g] == 0)
+    fisher.test(matrix(c(present_in, absent_in, present_out, absent_out), nrow = 2))$p.value
+  })
+})
+names(enrich) <- paste0("cluster_", sort(unique(cluster_labels)))
+
+# Adjust p-values
+enrich_adj <- purrr::map(enrich, p.adjust, method = "BH")
+
+# Save enrichment results
+enrich_df <- as.data.frame(do.call(cbind, enrich_adj))
+enrich_df$Gene <- rownames(enrich_df)
+readr::write_tsv(
+  enrich_df[, c("Gene", names(enrich_adj))],
+  file.path(results_dir, "unsup_enrichment_adj_pvalues.tsv")
+)