From cf01a8c0727251b4323d417b9e74cc9ee44506e4 Mon Sep 17 00:00:00 2001 From: Charmie Vang <147016287+charmvang@users.noreply.github.com> Date: Thu, 12 Feb 2026 07:41:14 -0700 Subject: [PATCH 1/2] Add unsupervised clustering script for genomes This script performs unsupervised clustering of genomes to investigate the separation of light ecotypes using various clustering techniques and visualizations. --- scripts/unsup_clust_genomes.R | 114 ++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 scripts/unsup_clust_genomes.R diff --git a/scripts/unsup_clust_genomes.R b/scripts/unsup_clust_genomes.R new file mode 100644 index 0000000..1f39579 --- /dev/null +++ b/scripts/unsup_clust_genomes.R @@ -0,0 +1,114 @@ +# Not sure what number script yet, still finalizing and will update later + +# Purpose: Does unsupervised clustering of the genomes show separation of light ecotypes? +# Based on github issue: https://github.com/JRaviLab/cyano_adaptation/issues/2 + +# Use the PanTA pangenome files from +# https://ondemand-rmacc.rc.colorado.edu/pun/sys/dashboard/files/fs//pl/active/jravilab/BRaVE/manuscript/data/PanTA_All_NoAlignments + +library(data.table) +library(pheatmap) +library(vegan) # For calculating jacard +library(ape) +library(ggplot2) +library(uwot) +library(dbscan) + +pa_matrix <- fread("/pl/active/jravilab/BRaVE/manuscript/data/PanTA_All_NoAlignments/gene_presence_absence.Rtab") + +# Extract gene names +gene_ids <- pa_matrix[[1]] + +# Extract numeric matrix (genes × genomes) +mat_gene_genome <- as.matrix(pa_matrix[, -1, with=FALSE]) +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] +) + +ggplot(pcoa_df, aes(PC1, PC2)) + + geom_point(size=2) + + theme_classic() + + ggtitle("PCoA of PanTA gene-content (Jaccard)") + +set.seed(1) +umap_res <- 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] +) +ggplot(umap_df, aes(UMAP1, UMAP2)) + + geom_point(size=2) + + theme_classic() + + ggtitle("UMAP of PanTA gene-content") + +hc <- hclust(dist_jaccard, method = "ward.D2") +plot(hc, labels = FALSE) +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) +ggplot(umap_df, aes(UMAP1, UMAP2, color = cluster)) + + geom_point(size=2) + + theme_classic() + + ggtitle("UMAP + HDBSCAN clusters") + +# 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) +pheatmap( + mat_filt[, top_genes], + cluster_rows = TRUE, + cluster_cols = TRUE, + show_rownames = FALSE, + show_colnames = FALSE, + annotation_row = annotation_row +) + +cluster_labels <- clusters_hc +genes <- colnames(mat_filt) +enrich <- lapply(unique(cluster_labels), function(cl){ + sapply(genes, function(g){ + a <- sum(mat_filt[cluster_labels == cl, g] == 1) + b <- sum(mat_filt[cluster_labels == cl, g] == 0) + c <- sum(mat_filt[cluster_labels != cl, g] == 1) + d <- sum(mat_filt[cluster_labels != cl, g] == 0) + fisher.test(matrix(c(a,b,c,d), nrow=2))$p.value + }) +}) + +# Adjust p-values +enrich_adj <- lapply(enrich, p.adjust, method="BH") From 1c99bb44cd89e866fd995471fc91589056cc94e3 Mon Sep 17 00:00:00 2001 From: Janani Ravi Date: Wed, 15 Apr 2026 20:29:55 -0600 Subject: [PATCH 2/2] minor edits - sapply to purrr map function - alpine path - file saves - pkg::function calls in some cases - switched to here::here and readr for consistency - c is for concatenate. d is inbuilt function too. risky to set a variable name to `c`!! check new names nad runs. --- scripts/unsup_clust_genomes.R | 107 +++++++++++++++++++++++----------- 1 file changed, 74 insertions(+), 33 deletions(-) diff --git a/scripts/unsup_clust_genomes.R b/scripts/unsup_clust_genomes.R index 1f39579..235dea6 100644 --- a/scripts/unsup_clust_genomes.R +++ b/scripts/unsup_clust_genomes.R @@ -1,26 +1,35 @@ -# Not sure what number script yet, still finalizing and will update later +### 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 +# Based on github issue: https://github.com/JRaviLab/cyano_adaptation/issues/2 -# Use the PanTA pangenome files from -# https://ondemand-rmacc.rc.colorado.edu/pun/sys/dashboard/files/fs//pl/active/jravilab/BRaVE/manuscript/data/PanTA_All_NoAlignments - -library(data.table) +library(here) +library(purrr) +library(readr) library(pheatmap) -library(vegan) # For calculating jacard +library(vegan) # For calculating Jaccard distance library(ape) library(ggplot2) library(uwot) library(dbscan) -pa_matrix <- fread("/pl/active/jravilab/BRaVE/manuscript/data/PanTA_All_NoAlignments/gene_presence_absence.Rtab") +# 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, with=FALSE]) +mat_gene_genome <- as.matrix(pa_matrix[, -1]) rownames(mat_gene_genome) <- gene_ids # Transpose → genomes × genes @@ -40,36 +49,45 @@ keep_genes <- which( mat_filt <- mat_bin[, keep_genes] dim(mat_filt) -# Compute Jaccard distance +# 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] + PC1 = pcoa_res$vectors[, 1], + PC2 = pcoa_res$vectors[, 2] ) -ggplot(pcoa_df, aes(PC1, PC2)) + - geom_point(size=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) -set.seed(1) -umap_res <- umap(mat_filt, n_neighbors = 15, min_dist = 0.1) +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] + UMAP1 = umap_res[, 1], + UMAP2 = umap_res[, 2] ) -ggplot(umap_df, aes(UMAP1, UMAP2)) + - geom_point(size=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) @@ -77,38 +95,61 @@ 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) -ggplot(umap_df, aes(UMAP1, UMAP2, color = cluster)) + - geom_point(size=2) + +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] +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, + 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 <- lapply(unique(cluster_labels), function(cl){ - sapply(genes, function(g){ - a <- sum(mat_filt[cluster_labels == cl, g] == 1) - b <- sum(mat_filt[cluster_labels == cl, g] == 0) - c <- sum(mat_filt[cluster_labels != cl, g] == 1) - d <- sum(mat_filt[cluster_labels != cl, g] == 0) - fisher.test(matrix(c(a,b,c,d), nrow=2))$p.value +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 <- lapply(enrich, p.adjust, method="BH") +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") +)