---
title: "13.00-multiomics-barnacle"
author: "Sam White"
date: "2025-10-03"
output: 
  github_document:
    toc: true
    number_sections: true
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
  html_document:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
bibliography: references.bib
citeproc: true
---

# BACKGROUND

This analysis prepares three-species ortholog expression matrices, normalizes the data, constructs a multi‑omics tensor (genes × combined species-samples × timepoints), and runs a sparse tensor decomposition using the [barnacle](https://github.com/blasks/barnacle) [@blaskowski2024] workflow to discover shared gene/time/species patterns.

## Workflow overview

1. **Data Preparation** - Extract and normalize ortholog expression data
2. **Tensor Construction** - Build 3D tensor combining species, samples, and timepoints
3. **Rank Comparison** - Test multiple ranks systematically to find optimal components
4. **Interpretation** - Choose best rank and examine factors, visualizations, and biological meaning

## Rank selection approach

**This analysis uses a systematic rank comparison approach.**

The workflow:
1. Configure `ranks_to_test` in the comparison chunk (e.g., `[5, 8, 10, 12, 15, 20]`)
2. Run decompositions for all specified ranks (10-30 minutes)
3. Review comparison metrics and elbow plot
4. Choose optimal rank based on:
   - Variance explained (aim for >60%)
   - Diminishing returns in elbow plot
   - Sparsity patterns
   - Component interpretability
5. Use chosen rank for final biological interpretation

**Key metrics for comparison:**
- **Variance explained**: How much of data variation is captured
- **Relative error**: Reconstruction accuracy
- **Sparsity**: How many genes/samples load strongly per component
- **Convergence**: Whether optimization completed successfully
- **Component weights**: Distribution of importance across components

## Key steps performed by this script:

- Load annotated ortholog groups and per-species count matrices.
- Filter ortholog groups to retain complete three‑way matches with available expression data.
- Extract per‑species expression matrices mapped to ortholog group IDs.
- Normalize counts (preferred: `sctransform`; fallback: log2(CPM + 1)).
- Build a 3D tensor combining species/sample and timepoint dimensions.
- Run [barnacle](https://github.com/blasks/barnacle) [@blaskowski2024] sparse CP decomposition and save factor matrices, metadata and figures.
- Assess rank appropriateness through multiple metrics and visualizations.

## Input files

- `../output/12-ortho-annot/ortholog_groups_annotated.csv` : annotated ortholog groups with columns for `group_id`, `apul`, `peve`, `ptua`, `type`, etc.
- `../../D-Apul/output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv` : Apul gene count matrix (contains `gene_id` plus sample columns).
- `../../E-Peve/output/02.20-E-Peve-RNAseq-alignment-HiSat2/peve-gene_count_matrix.csv` : Peve gene count matrix.
- `../../F-Ptua/output/02.20-F-Ptua-RNAseq-alignment-HiSat2/ptua-gene_count_matrix.csv` : Ptua gene count matrix.
- (Also referenced) transcript-level matrices: `apul-transcript_count_matrix.csv`, `peve-transcript_count_matrix.csv`, `ptua-transcript_count_matrix.csv` when needed.

## Output files

- `apul_ortholog_expression.csv`, `peve_ortholog_expression.csv`, `ptua_ortholog_expression.csv` : per‑species expression matrices aligned to ortholog `group_id`.
- `apul_normalized_expression.csv`, `peve_normalized_expression.csv`, `ptua_normalized_expression.csv` : normalized expression matrices (sctransform or log2(CPM+1)).
- `multiomics_tensor.npy` : saved NumPy array of the 3D tensor used for decomposition (genes × combined_samples × timepoints).
- `barnacle_factors/` directory containing:
    - `gene_factors.csv` : gene loadings per component (genes × components).
    - `sample_factors.csv` : combined sample (species_sample) loadings per component with `Species` and `Sample_ID` metadata.
    - `time_factors.csv` : timepoint loadings per component.
    - `component_weights.csv` : component weights / importance.
    - `sample_mapping.csv` : mapping of combined sample indices to species and sample IDs.
    - `metadata.json` : analysis parameters and tensor metadata (shape, rank, lambdas, convergence, etc.).
    - `figures/` : generated visualizations (component weights, time loadings, sample heatmap, PCA, top ortholog plots).

## Notes / assumptions:

- Sample column names are parsed expecting a dot-separated format with a `TP#` timepoint token (e.g., `ACR.139.TP1`).
- Apul ortholog IDs in the ortholog table include transcript suffixes (e.g., `-T1`) which are removed for matching to the gene count matrix.
- Missing values in the tensor are handled by substitution (current workflow fills NaNs with zeros before decomposition).
- Samples must have all four timepoints (TP1, TP2, TP3, TP4) to be included in the analysis.


# SETUP

## Libraries

```{r setup, include=FALSE}
library(knitr)
library(tidyverse)
library(reticulate)
library(matrixStats)  # Required for sctransform internal functions
library(sctransform)
library(glmGamPoi)

knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
    eval = FALSE,        # Evaluate code chunks
    results = "hold",   # Hold outputs and show them after the full code chunk
  warning = FALSE,     # Hide warnings
        collapse = FALSE,    # Keep code and output in separate blocks
        warning = FALSE,     # Hide warnings
        message = FALSE,     # Hide messages
        comment = "##"      # Prefix output lines with '##' so output is visually distinct
)
```

## Set R variables

```{r R-variables, eval=TRUE}
# OUTPUT DIRECTORY
output_dir <- "../output/13.00-multiomics-barnacle"

#INPUT FILE(S)
ortholog_groups_file <- "../output/12-ortho-annot/ortholog_groups_annotated.csv"

# Transcript count matrices for each species
apul_transcript_matrix_file <- "../../D-Apul/output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-transcript_count_matrix.csv"
peve_transcript_matrix_file <- "../../E-Peve/output/02.20-E-Peve-RNAseq-alignment-HiSat2/peve-transcript_count_matrix.csv"
ptua_transcript_matrix_file <- "../../F-Ptua/output/02.20-F-Ptua-RNAseq-alignment-HiSat2/ptua-transcript_count_matrix.csv"

# CONDA
conda_env_name <- c("/home/sam/programs/mambaforge/envs/barnacle_py311_env")
conda_path <- c("/opt/anaconda/anaconda3/bin/conda")
```

## Load [barnacle](https://github.com/blasks/barnacle) conda environment

If this is successful, the first line of output should show that the Python being used is the one in your [barnacle](https://github.com/blasks/barnacle) [@blaskowski2024] conda environment path.

E.g.

`python:         /home/sam/programs/mambaforge/envs/barnacle_py311_env/bin/python`

```{r load-barnacle-conda-env, eval=TRUE}
use_condaenv(condaenv = conda_env_name, conda = conda_path)
py_config()
```

# DATA PREP

## Load ortholog groups data

```{r load-ortholog-data, eval=TRUE}
# Read in the ortholog groups data
ortholog_groups <- read.csv(ortholog_groups_file)



# Display basic info about the data
cat("Dimensions of ortholog groups data:", dim(ortholog_groups), "\n\n")
cat("Column names:", colnames(ortholog_groups), "\n\n")
head(ortholog_groups)
str(ortholog_groups)
```

## Extract ortholog expression data

Now let's extract expression data for genes that are present in the ortholog groups. We'll use the gene count matrices with gene IDs to properly map the data.

### Load gene count matrices

```{r load-gene-matrices, eval=TRUE}
# Define file paths for gene count matrices
apul_gene_matrix_file <- "../../D-Apul/output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv"
peve_gene_matrix_file <- "../../E-Peve/output/02.20-E-Peve-RNAseq-alignment-HiSat2/peve-gene_count_matrix.csv"
ptua_gene_matrix_file <- "../../F-Ptua/output/02.20-F-Ptua-RNAseq-alignment-HiSat2/ptua-gene_count_matrix.csv"

# Load gene count matrices for each species
cat("Loading gene count matrices...\n\n")

apul_gene_matrix <- read.csv(apul_gene_matrix_file)
cat("Apul gene matrix dimensions:", dim(apul_gene_matrix), "\n")

peve_gene_matrix <- read.csv(peve_gene_matrix_file)
cat("Peve gene matrix dimensions:", dim(peve_gene_matrix), "\n")

ptua_gene_matrix <- read.csv(ptua_gene_matrix_file)
cat("Ptua gene matrix dimensions:", dim(ptua_gene_matrix), "\n\n")
```

### Filter ortholog groups for complete three-way matches

```{r filter-ortholog-groups, eval=TRUE}
cat("Filtering for complete three-way ortholog groups...\n")

# Keep only rows where all three species have entries (no NA values or empty strings)
complete_ortholog_groups <- ortholog_groups[nzchar(ortholog_groups$apul) & 
                                          nzchar(ortholog_groups$peve) & 
                                          nzchar(ortholog_groups$ptua), ]

cat("Total ortholog groups:", nrow(ortholog_groups), "\n")
cat("Complete three-way ortholog groups:", nrow(complete_ortholog_groups), "\n")

```

### Filter for expression data availability

```{r filter-expression-availability, eval=TRUE}
cat("Filtering ortholog groups to ensure all genes have expression data...\n")

# Clean gene IDs to check against expression data
# For Apul: remove -T[n] suffix from ortholog groups to match gene matrix format
apul_ortholog_genes_check <- gsub("-T[0-9]+$", "", complete_ortholog_groups$apul)

# For Peve and Ptua: use as-is (will clean gene matrix IDs later)
peve_ortholog_genes_check <- complete_ortholog_groups$peve
ptua_ortholog_genes_check <- complete_ortholog_groups$ptua

# Check which genes are present in expression data
# (Note: We need to clean gene matrix IDs to match)
apul_gene_matrix_ids <- gsub("^gene-", "", apul_gene_matrix$gene_id)  # Remove gene- prefix if present
peve_gene_matrix_ids <- gsub("^gene-", "", peve_gene_matrix$gene_id)  # Remove gene- prefix
ptua_gene_matrix_ids <- gsub("^gene-", "", ptua_gene_matrix$gene_id)  # Remove gene- prefix

# Find which ortholog genes are present in each species' expression data
apul_present <- apul_ortholog_genes_check %in% apul_gene_matrix_ids
peve_present <- peve_ortholog_genes_check %in% peve_gene_matrix_ids
ptua_present <- ptua_ortholog_genes_check %in% ptua_gene_matrix_ids

# Keep only ortholog groups where all three species have expression data
expression_complete_mask <- apul_present & peve_present & ptua_present
complete_ortholog_groups <- complete_ortholog_groups[expression_complete_mask, ]

cat("Ortholog groups after filtering for expression data availability:", nrow(complete_ortholog_groups), "\n")

```

### Gene ID cleaning examples

```{r gene-id-examples-testing, eval=TRUE}
cat("\n=== GENE ID CLEANING EXAMPLES ===\n")
cat("Apul (clean ortholog groups to match gene matrix):\n")
cat("Ortholog groups original:", head(complete_ortholog_groups$apul, 3), "\n")
cat("Ortholog groups cleaned:", head(gsub("-T[0-9]+$", "", complete_ortholog_groups$apul), 3), "\n")
cat("Gene matrix (target format):", head(apul_gene_matrix$gene_id, 3), "\n\n")

cat("Peve (clean gene matrix to match ortholog groups):\n")
cat("Ortholog groups (target format):", head(complete_ortholog_groups$peve, 3), "\n") 
cat("Gene matrix original:", head(peve_gene_matrix$gene_id, 3), "\n")
cat("Gene matrix cleaned:", head(peve_gene_matrix$gene_id_clean, 3), "\n\n")

cat("Ptua (clean gene matrix to match ortholog groups):\n")
cat("Ortholog groups (target format):", head(complete_ortholog_groups$ptua, 3), "\n")
cat("Gene matrix original:", head(ptua_gene_matrix$gene_id, 3), "\n")
cat("Gene matrix cleaned:", head(ptua_gene_matrix$gene_id_clean, 3), "\n\n")

```

### Clean gene IDs for matching

```{r clean-gene-ids, eval=TRUE}
cat("Cleaning gene matrix IDs to match ortholog group format...\n")

# For Apul: ortholog groups have "FUN_000185-T1", gene matrix has "FUN_002326"
# We need to remove "-T1" from ortholog groups to match gene matrix
apul_ortholog_genes <- unique(gsub("-T[0-9]+$", "", complete_ortholog_groups$apul))

# For Peve and Ptua: keep ortholog groups as-is and clean gene matrix
peve_ortholog_genes <- unique(complete_ortholog_groups$peve)
ptua_ortholog_genes <- unique(complete_ortholog_groups$ptua)

# Clean gene matrix IDs accordingly
# Apul: gene matrix already in correct format (no cleaning needed)
apul_gene_matrix$gene_id_clean <- apul_gene_matrix$gene_id

# Peve: gene matrix has "gene-Peve_00000032", ortholog groups have "Peve_00037402"  
# So we need to remove "gene-" prefix from gene matrix
peve_gene_matrix$gene_id_clean <- gsub("^gene-", "", peve_gene_matrix$gene_id)

# Ptua: gene matrix has "gene-Pocillopora_meandrina_HIv1___RNAseq.g20905.t1", 
# ortholog groups have "Pocillopora_meandrina_HIv1___RNAseq.g28886.t1"
# So we just need to remove "gene-" prefix from gene matrix
ptua_gene_matrix$gene_id_clean <- gsub("^gene-", "", ptua_gene_matrix$gene_id)

cat("Apul ortholog genes (complete groups only):", length(apul_ortholog_genes), "\n")
cat("Peve ortholog genes (complete groups only):", length(peve_ortholog_genes), "\n")
cat("Ptua ortholog genes (complete groups only):", length(ptua_ortholog_genes), "\n\n")
```

### Define sample filtering function

```{r define-sample-filter-function, eval=TRUE}
# Function to filter samples that have all four timepoints (TP1, TP2, TP3, TP4)
filter_complete_samples <- function(gene_matrix, species_name) {
  cat("  Filtering samples with complete timepoints for", species_name, "...\n")
  
  # Get expression columns (exclude gene_id column, gene_id_clean may not exist yet)
  all_cols <- colnames(gene_matrix)
  id_cols <- c("gene_id", "gene_id_clean")
  expr_cols <- setdiff(all_cols, id_cols)
  
  if(length(expr_cols) == 0) {
    cat("  ERROR: No expression columns found! Returning original matrix.\n")
    return(gene_matrix)
  }
  
  # Parse sample column names to identify sample groups and timepoints
  sample_timepoint_map <- list()
  
  for(col in expr_cols) {
    # Expected formats: ACR.139.TP1, POR.216.TP1, POC.201.TP1, etc.
    # Note: Using dots (.) as separators, not hyphens (-)
    parts <- strsplit(col, "\\.")[[1]]  # Split on dots, not hyphens
    
    if(length(parts) >= 3) {
      # Extract timepoint from last part (e.g., "TP1", "TP2", etc.)
      tp_part <- parts[length(parts)]
      
      if(grepl("^TP[0-9]+$", tp_part)) {
        timepoint <- as.numeric(gsub("TP", "", tp_part))
        
        # Sample ID is everything except the timepoint part, joined with dots
        sample_id <- paste(parts[1:(length(parts)-1)], collapse = ".")
        
        if(!sample_id %in% names(sample_timepoint_map)) {
          sample_timepoint_map[[sample_id]] <- list()
        }
        sample_timepoint_map[[sample_id]][[as.character(timepoint)]] <- col
      }
    }
  }
  
  # Identify samples with all four timepoints (1, 2, 3, 4)
  complete_samples <- c()
  incomplete_samples <- c()
  required_timepoints <- c("1", "2", "3", "4")
  
  for(sample_id in names(sample_timepoint_map)) {
    available_timepoints <- names(sample_timepoint_map[[sample_id]])
    
    if(all(required_timepoints %in% available_timepoints)) {
      complete_samples <- c(complete_samples, sample_id)
    } else {
      incomplete_samples <- c(incomplete_samples, sample_id)
      missing_tps <- setdiff(required_timepoints, available_timepoints)
      cat("    Sample", sample_id, "missing timepoints:", paste(paste0("TP", missing_tps), collapse = ", "), "\n")
    }
  }
  
  cat("  Complete samples (all 4 timepoints):", length(complete_samples), "\n")
  cat("  Incomplete samples:", length(incomplete_samples), "\n")
  
  if(length(incomplete_samples) > 0) {
    cat("  Removing incomplete samples:", paste(incomplete_samples, collapse = ", "), "\n")
  }
  
  # Build list of columns to keep (include gene_id columns plus ONLY complete sample columns)
  if(length(complete_samples) == 0) {
    cat("  WARNING: No complete samples found! Keeping all expression columns.\n")
    # Fallback: keep all expression columns
    keep_cols <- c("gene_id")
    if("gene_id_clean" %in% colnames(gene_matrix)) {
      keep_cols <- c(keep_cols, "gene_id_clean")
    }
    keep_cols <- c(keep_cols, expr_cols)
  } else {
    keep_cols <- c("gene_id")
    if("gene_id_clean" %in% colnames(gene_matrix)) {
      keep_cols <- c(keep_cols, "gene_id_clean")
    }
    
    # ONLY add columns from complete samples (this excludes incomplete sample columns)
    for(sample_id in complete_samples) {
      for(tp in required_timepoints) {
        if(tp %in% names(sample_timepoint_map[[sample_id]])) {
          keep_cols <- c(keep_cols, sample_timepoint_map[[sample_id]][[tp]])
        }
      }
    }
    
    # Double-check: ensure we're not accidentally including incomplete sample columns
    incomplete_cols <- c()
    for(sample_id in incomplete_samples) {
      for(tp in names(sample_timepoint_map[[sample_id]])) {
        incomplete_cols <- c(incomplete_cols, sample_timepoint_map[[sample_id]][[tp]])
      }
    }
    
    if(length(incomplete_cols) > 0) {
      # Remove any incomplete sample columns that might have been included
      keep_cols <- setdiff(keep_cols, incomplete_cols)
    }
  }
  
  # Filter the gene matrix to keep only complete samples
  keep_cols <- intersect(keep_cols, colnames(gene_matrix))  # Ensure columns exist
  filtered_matrix <- gene_matrix[, keep_cols, drop = FALSE]
  
  # Count ID vs expression columns in final result
  final_expr_cols <- setdiff(colnames(filtered_matrix), c("gene_id", "gene_id_clean"))
  cat("  Final expression columns:", length(final_expr_cols), "\n")
  
  return(filtered_matrix)
}
```

### Define expression extraction function

```{r define-expression-extraction-function, eval=TRUE}
# Function to extract expression data for a species using ortholog group mapping
extract_species_expression <- function(ortholog_groups_df, gene_matrix, species_col, species_name) {
  cat("Processing", species_name, "...\n")
  
  # First, filter the gene matrix to keep only samples with complete timepoints
  filtered_gene_matrix <- filter_complete_samples(gene_matrix, species_name)
  
  # Remove duplicate group_ids, keeping first occurrence of each
  unique_groups <- ortholog_groups_df[!duplicated(ortholog_groups_df$group_id), ]
  cat("  Unique ortholog groups:", nrow(unique_groups), "\n")
  
  # Create results data frame starting with group_id
  result_df <- data.frame(group_id = unique_groups$group_id, stringsAsFactors = FALSE)
  
  # Get expression columns (exclude gene_id and gene_id_clean columns)
  expr_cols <- setdiff(colnames(filtered_gene_matrix), c("gene_id", "gene_id_clean"))
  
  # Initialize expression columns with NA
  for(col in expr_cols) {
    result_df[[col]] <- NA
  }
  
  # For each ortholog group, find the corresponding gene and extract expression
  for(i in seq_len(nrow(unique_groups))) {
    target_gene <- unique_groups[[species_col]][i]
    
    # Clean target gene for matching
    if(species_name == "Apul") {
      # Remove transcript suffix for Apul
      target_gene_clean <- gsub("-T[0-9]+$", "", target_gene)
      matching_rows <- which(filtered_gene_matrix$gene_id_clean == target_gene_clean)
    } else {
      # For Peve and Ptua, match cleaned gene_id
      matching_rows <- which(filtered_gene_matrix$gene_id_clean == target_gene)
    }
    
    if(length(matching_rows) == 1) {
      # Single match - copy expression data
      for(col in expr_cols) {
        result_df[i, col] <- filtered_gene_matrix[matching_rows, col]
      }
    } else if(length(matching_rows) > 1) {
      # Multiple matches - take first gene (no averaging)
      first_match <- matching_rows[1]
      for(col in expr_cols) {
        result_df[i, col] <- filtered_gene_matrix[first_match, col]
      }
    } else {
      # No match - leave as NA (will be removed later)
    }
  }
  
  # Remove rows with all NA expression values
  expr_na_mask <- apply(result_df[, expr_cols, drop = FALSE], 1, function(x) all(is.na(x)))
  result_df <- result_df[!expr_na_mask, ]
  
  cat("  Final dimensions:", nrow(result_df), "ortholog groups x", ncol(result_df)-1, "samples\n")
  cat("  Removed", sum(expr_na_mask), "groups with no expression data\n\n")
  
  return(result_df)
}
```

### Extract ortholog expression data

```{r extract-ortholog-expression-data, eval=TRUE}
cat("Creating ortholog expression data with proper group_id mapping...\n")

# Extract expression data for each species using the ortholog group mapping
apul_ortholog_expression <- extract_species_expression(complete_ortholog_groups, apul_gene_matrix, "apul", "Apul")
peve_ortholog_expression <- extract_species_expression(complete_ortholog_groups, peve_gene_matrix, "peve", "Peve")
ptua_ortholog_expression <- extract_species_expression(complete_ortholog_groups, ptua_gene_matrix, "ptua", "Ptua")

cat("=== FINAL ORTHOLOG EXPRESSION DATA DIMENSIONS ===\n")
cat("Apul:", nrow(apul_ortholog_expression), "ortholog groups x", ncol(apul_ortholog_expression)-1, "samples\n")
cat("Peve:", nrow(peve_ortholog_expression), "ortholog groups x", ncol(peve_ortholog_expression)-1, "samples\n")
cat("Ptua:", nrow(ptua_ortholog_expression), "ortholog groups x", ncol(ptua_ortholog_expression)-1, "samples\n\n")
```

### Write ortholog expression data

```{r write-expression-data, eval=TRUE}
cat("Exporting ortholog expression data to CSV files...\n")

# Define output file paths
apul_output_file <- file.path(output_dir, "apul_ortholog_expression.csv")
peve_output_file <- file.path(output_dir, "peve_ortholog_expression.csv")
ptua_output_file <- file.path(output_dir, "ptua_ortholog_expression.csv")

# Write CSV files without quotes
write.csv(apul_ortholog_expression, file = apul_output_file, quote = FALSE, row.names = FALSE)
cat("Exported Apul ortholog expression data to:", apul_output_file, "\n")

write.csv(peve_ortholog_expression, file = peve_output_file, quote = FALSE, row.names = FALSE)
cat("Exported Peve ortholog expression data to:", peve_output_file, "\n")

write.csv(ptua_ortholog_expression, file = ptua_output_file, quote = FALSE, row.names = FALSE)
cat("Exported Ptua ortholog expression data to:", ptua_output_file, "\n")

cat("\nAll ortholog expression data exported successfully!\n\n")
```

### Column structure analysis

```{r column-structure-analysis, eval=TRUE}
cat("=== COLUMN STRUCTURE ANALYSIS ===\n")
cat("Apul columns:", ncol(apul_ortholog_expression), "\n")
cat("Apul column names (first 10):", paste(head(colnames(apul_ortholog_expression), 10), collapse = ", "), "\n")
cat("Apul column names (last 10):", paste(tail(colnames(apul_ortholog_expression), 10), collapse = ", "), "\n\n")

cat("Peve columns:", ncol(peve_ortholog_expression), "\n")
cat("Peve column names (first 10):", paste(head(colnames(peve_ortholog_expression), 10), collapse = ", "), "\n")
cat("Peve column names (last 10):", paste(tail(colnames(peve_ortholog_expression), 10), collapse = ", "), "\n\n")

cat("Ptua columns:", ncol(ptua_ortholog_expression), "\n")
cat("Ptua column names (first 10):", paste(head(colnames(ptua_ortholog_expression), 10), collapse = ", "), "\n")
cat("Ptua column names (last 10):", paste(tail(colnames(ptua_ortholog_expression), 10), collapse = ", "), "\n\n")

```

### Summary statistics

```{r summary-statistics, eval=TRUE}
cat("=== LINE COUNTS FOR ORTHOLOG EXPRESSION DATA ===\n")
cat("Apul ortholog expression with info: ", nrow(apul_ortholog_expression), " rows\n")
cat("Peve ortholog expression with info: ", nrow(peve_ortholog_expression), " rows\n")
cat("Ptua ortholog expression with info: ", nrow(ptua_ortholog_expression), " rows\n\n")
```

### Sample filtering verification

```{r sample-filtering-verification, eval=TRUE}
cat("=== SAMPLE FILTERING VERIFICATION ===\n")
cat("This analysis filters samples to include ONLY those with all four timepoints (TP1, TP2, TP3, TP4)\n\n")

# Function to analyze sample completeness in the filtered data
analyze_sample_completeness <- function(data_df, species_name) {
  cat("Analyzing", species_name, "sample completeness:\n")
  
  # Get expression columns (exclude group_id)
  expr_cols <- setdiff(colnames(data_df), "group_id")
  
  # Parse sample column names to identify sample groups and timepoints
  sample_timepoint_map <- list()
  
  for(col in expr_cols) {
    parts <- strsplit(col, "\\.")[[1]]  # Split on dots to match actual format
    
    if(length(parts) >= 3) {
      tp_part <- parts[length(parts)]
      
      if(grepl("^TP[0-9]+$", tp_part)) {
        timepoint <- as.numeric(gsub("TP", "", tp_part))
        sample_id <- paste(parts[1:(length(parts)-1)], collapse = ".")
        
        if(!sample_id %in% names(sample_timepoint_map)) {
          sample_timepoint_map[[sample_id]] <- c()
        }
        sample_timepoint_map[[sample_id]] <- c(sample_timepoint_map[[sample_id]], timepoint)
      }
    }
  }
  
  # Check completeness
  complete_samples <- 0
  incomplete_samples <- 0
  required_timepoints <- c(1, 2, 3, 4)
  
  for(sample_id in names(sample_timepoint_map)) {
    available_timepoints <- sort(sample_timepoint_map[[sample_id]])
    
    if(all(required_timepoints %in% available_timepoints)) {
      complete_samples <- complete_samples + 1
    } else {
      incomplete_samples <- incomplete_samples + 1
      missing_tps <- setdiff(required_timepoints, available_timepoints)
      cat("  WARNING: Sample", sample_id, "missing timepoints:", paste(paste0("TP", missing_tps), collapse = ", "), "\n")
    }
  }
  
  cat("  Total samples:", length(sample_timepoint_map), "\n")
  cat("  Complete samples (all 4 timepoints):", complete_samples, "\n")
  cat("  Incomplete samples:", incomplete_samples, "\n")
  
  if(incomplete_samples > 0) {
    cat("  ERROR: Filtering did not work correctly!\n")
  } else {
    cat("  SUCCESS: All samples have complete timepoints\n")
  }
  
  cat("\n")
  
  return(list(
    total = length(sample_timepoint_map),
    complete = complete_samples,
    incomplete = incomplete_samples
  ))
}

# Verify each species
apul_stats <- analyze_sample_completeness(apul_ortholog_expression, "Apul")
peve_stats <- analyze_sample_completeness(peve_ortholog_expression, "Peve")
ptua_stats <- analyze_sample_completeness(ptua_ortholog_expression, "Ptua")

cat("=== FILTERING SUMMARY ===\n")
cat("All samples in the filtered dataset should have exactly 4 timepoints (TP1, TP2, TP3, TP4)\n")
cat("Apul: ", apul_stats$complete, "/", apul_stats$total, " complete samples\n")
cat("Peve: ", peve_stats$complete, "/", peve_stats$total, " complete samples\n")
cat("Ptua: ", ptua_stats$complete, "/", ptua_stats$total, " complete samples\n")

total_incomplete <- apul_stats$incomplete + peve_stats$incomplete + ptua_stats$incomplete
if(total_incomplete == 0) {
  cat("SUCCESS: All samples across all species have complete timepoints!\n")
} else {
  cat("ERROR: Found", total_incomplete, "incomplete samples - filtering needs to be fixed!\n")
}
cat("\n")
```

### Diagnostic analysis

```{r diagnostic-analysis, eval=TRUE}
cat("\n=== VERIFICATION: THREE-WAY ORTHOLOG COUNTS ===\n")
cat("After filtering, all species should have identical three_way ortholog counts.\n\n")

# Check how many genes we have for each species
cat("Genes found in expression data by species:\n")
cat("Apul ortholog expression (before adding info):", nrow(apul_ortholog_expression), "\n")
cat("Peve ortholog expression (before adding info):", nrow(peve_ortholog_expression), "\n")
cat("Ptua ortholog expression (before adding info):", nrow(ptua_ortholog_expression), "\n\n")

# Verify all ortholog groups are three-way
cat("Complete three-way ortholog groups available:", nrow(complete_ortholog_groups), "\n")
cat("All should be type 'three_way'? Check:", table(complete_ortholog_groups$type), "\n\n")

# Verify perfect gene coverage (should be 100% for all species now)
apul_coverage <- length(intersect(apul_ortholog_genes, apul_gene_matrix$gene_id_clean))
peve_coverage <- length(intersect(peve_ortholog_genes, peve_gene_matrix$gene_id_clean))
ptua_coverage <- length(intersect(ptua_ortholog_genes, ptua_gene_matrix$gene_id_clean))

cat("Gene coverage in expression data (should be 100% for all):\n")
cat("Apul: ", apul_coverage, "/", length(apul_ortholog_genes), " (", round(apul_coverage/length(apul_ortholog_genes)*100, 1), "%)\n")
cat("Peve: ", peve_coverage, "/", length(peve_ortholog_genes), " (", round(peve_coverage/length(peve_ortholog_genes)*100, 1), "%)\n")
cat("Ptua: ", ptua_coverage, "/", length(ptua_ortholog_genes), " (", round(ptua_coverage/length(ptua_ortholog_genes)*100, 1), "%)\n\n")
```

### Preview expression data

```{r preview-expression-data, eval=TRUE}
cat("\n=== PREVIEW OF ORTHOLOG EXPRESSION DATA ===\n")

cat("Apul ortholog expression:\n\n")
str(apul_ortholog_expression)
cat("\n\n")

cat("\nPeve ortholog expression:\n")
str(peve_ortholog_expression)
cat("\n\n")

cat("\nPtua ortholog expression:\n")
str(ptua_ortholog_expression)
cat("\n\n")
```

# BARNACLE ANALYSIS

This workflow focuses on comparing multiple ranks to find the optimal number of components.

Based on the barnacle workflow, we: 1. Normalize count data with `sctransform` 2. Create tensors for multiomics analysis 3. Compare multiple ranks systematically 4. Choose optimal rank based on metrics

## Load expression data

```{r load-normalized-data, eval=TRUE}
# Read the exported ortholog expression data
apul_expr <- read.csv(file.path(output_dir, "apul_ortholog_expression.csv"))
peve_expr <- read.csv(file.path(output_dir, "peve_ortholog_expression.csv"))  
ptua_expr <- read.csv(file.path(output_dir, "ptua_ortholog_expression.csv"))

cat("Loaded expression data:\n")
cat("Apul:", nrow(apul_expr), "genes x", ncol(apul_expr)-1, "samples\n")
cat("Peve:", nrow(peve_expr), "genes x", ncol(peve_expr)-1, "samples\n") 
cat("Ptua:", nrow(ptua_expr), "genes x", ncol(ptua_expr)-1, "samples\n")
```

## Normalize data with sctransform

Following the barnacle manuscript approach, we'll use `sctransform` to normalize each species' data. We'll use a bulk RNA-seq appropriate approach.

```{r normalize-sctransform, eval=TRUE}
# Function to normalize count data with sctransform for bulk RNA-seq
normalize_with_sctransform <- function(count_data, species_name) {
  cat("Normalizing", species_name, "data with sctransform...\n")
  
  # Check if we have group_id or gene_id column
  id_col <- if("group_id" %in% colnames(count_data)) "group_id" else "gene_id"
  cat("Using", id_col, "as identifier column\n")
  
  # Check for and handle duplicate group_ids/gene_ids
  duplicate_ids <- count_data[[id_col]][duplicated(count_data[[id_col]])]
  if(length(duplicate_ids) > 0) {
    cat("Found", length(unique(duplicate_ids)), "duplicate", id_col, "values:\n")
    cat("  Examples:", head(unique(duplicate_ids), 10), "\n")
    cat("  Note: sctransform may fail due to duplicate row names, will fall back to log2(CPM + 1)\n")
  } else {
    cat("No duplicate", id_col, "values found\n")
  }
  
  # Use original data without aggregation - let sctransform handle duplicates or fail
  agg_data <- count_data
  
  # Extract count matrix (genes as rows, samples as columns)
  count_matrix <- as.matrix(agg_data[, -1])  # Remove id column
  rownames(count_matrix) <- agg_data[[id_col]]
  
  # Check for and handle problematic values
  cat("Checking data quality...\n")
  cat("  - Zero values:", sum(count_matrix == 0), "/", length(count_matrix), "\n")
  cat("  - NA values:", sum(is.na(count_matrix)), "\n")
  cat("  - Infinite values:", sum(is.infinite(count_matrix)), "\n")
  cat("  - Min value:", min(count_matrix, na.rm = TRUE), "\n")
  cat("  - Max value:", max(count_matrix, na.rm = TRUE), "\n")
  
  # Remove genes with all zeros or very low expression
  gene_sums <- rowSums(count_matrix)
  keep_genes <- gene_sums > 10  # Keep genes with total counts > 10
  count_matrix_filtered <- count_matrix[keep_genes, , drop = FALSE]
  
  cat("  - Filtered to", nrow(count_matrix_filtered), "ortholog groups (from", nrow(count_matrix), ")\n")
  
  # Ensure data is numeric and handle any potential issues
  count_matrix_filtered <- apply(count_matrix_filtered, c(1,2), function(x) {
    if(is.na(x) || !is.finite(x)) return(0)
    return(as.numeric(x))
  })
  
  # Ensure row and column names are character strings (not factors)
  rownames(count_matrix_filtered) <- as.character(rownames(count_matrix_filtered))
  colnames(count_matrix_filtered) <- as.character(colnames(count_matrix_filtered))
  
  # sctransform expects genes as rows and cells (samples) as columns —
  # our matrix is already genes x samples (rows=genes, cols=samples), so do NOT transpose.
  count_matrix_for_vst <- count_matrix_filtered
  
  # Apply sctransform normalization with bulk RNA-seq appropriate parameters
  normalized_df <- tryCatch({
    # First try with minimal parameters to isolate the issue
    cat("  Attempting sctransform with minimal parameters...\n")
    # Pass genes x samples matrix directly
    normalized <- sctransform::vst(
      count_matrix_for_vst,
      verbosity = 1
    )

    # Extract normalized data (should be genes x samples)
    normalized_data <- normalized$y
    
    # Create output data frame with original gene set (fill missing with zeros)
    full_normalized_data <- matrix(0, nrow = nrow(count_matrix), ncol = ncol(count_matrix))
    rownames(full_normalized_data) <- rownames(count_matrix)
    colnames(full_normalized_data) <- colnames(count_matrix)
    
    # Fill in normalized values for kept genes
    full_normalized_data[rownames(normalized_data), ] <- normalized_data
    
    result_df <- data.frame(
      group_id = rownames(full_normalized_data),
      full_normalized_data,
      stringsAsFactors = FALSE
    )
    
    cat("sctransform normalization successful for", species_name, "\n")
    return(result_df)
    
  }, error = function(e) {
    cat("sctransform with minimal parameters failed for", species_name, ":", e$message, "\n")
    cat("Trying sctransform with glmGamPoi method...\n")
    
    # Try again with glmGamPoi method
    tryCatch({
      normalized <- sctransform::vst(
        count_matrix_for_vst, 
        method = "glmGamPoi",
        n_genes = min(2000, nrow(count_matrix_for_vst)),
        return_cell_attr = TRUE,
        verbosity = 1
      )

      # Extract normalized data (genes x samples)
      normalized_data <- normalized$y
      
      # Create output data frame with original gene set (fill missing with zeros)
      full_normalized_data <- matrix(0, nrow = nrow(count_matrix), ncol = ncol(count_matrix))
      rownames(full_normalized_data) <- rownames(count_matrix)
      colnames(full_normalized_data) <- colnames(count_matrix)
      
      # Fill in normalized values for kept genes
      full_normalized_data[rownames(normalized_data), ] <- normalized_data
      
      result_df <- data.frame(
        group_id = rownames(full_normalized_data),
        full_normalized_data,
        stringsAsFactors = FALSE
      )
      
      cat("sctransform normalization with glmGamPoi successful for", species_name, "\n")
      return(result_df)
      
    }, error = function(e2) {
      cat("sctransform failed for", species_name, ":", e2$message, "\n")
      cat("Falling back to log2(CPM + 1) normalization...\n")
      
  # Fallback: log2(CPM + 1) normalization
  # Calculate CPM (Counts Per Million) using the filtered count matrix
  lib_sizes <- colSums(count_matrix_filtered)
  cpm_matrix <- sweep(count_matrix_filtered, 2, lib_sizes/1e6, FUN = "/")
      
      # Log2 transform with pseudocount
      normalized_data <- log2(cpm_matrix + 1)
      
      result_df <- data.frame(
        group_id = rownames(normalized_data),
        normalized_data,
        stringsAsFactors = FALSE
      )
      
      cat("Log2(CPM + 1) normalization complete for", species_name, "\n")
      return(result_df)
    })
  })
  
  cat("Input dimensions:", nrow(count_data), "rows x", ncol(count_data)-1, "samples\n")
  cat("Output dimensions:", nrow(normalized_df), "ortholog groups x", ncol(normalized_df)-1, "samples\n\n")
  
  return(normalized_df)
}

# Normalize each species
cat("=== STARTING NORMALIZATION ===\n\n")
apul_normalized <- normalize_with_sctransform(apul_expr, "Apul")
peve_normalized <- normalize_with_sctransform(peve_expr, "Peve") 
ptua_normalized <- normalize_with_sctransform(ptua_expr, "Ptua")
cat("=== NORMALIZATION COMPLETE ===\n\n")
```

## Export normalized data for Python analysis

```{r export-normalized-data, eval=TRUE}
# Export normalized data for Python processing
apul_norm_file <- file.path(output_dir, "apul_normalized_expression.csv")
peve_norm_file <- file.path(output_dir, "peve_normalized_expression.csv")
ptua_norm_file <- file.path(output_dir, "ptua_normalized_expression.csv")

write.csv(apul_normalized, apul_norm_file, row.names = FALSE, quote = FALSE)
write.csv(peve_normalized, peve_norm_file, row.names = FALSE, quote = FALSE)
write.csv(ptua_normalized, ptua_norm_file, row.names = FALSE, quote = FALSE)

cat("Exported normalized data:\n")
cat("Apul:", apul_norm_file, "\n")
cat("Peve:", peve_norm_file, "\n")
cat("Ptua:", ptua_norm_file, "\n\n")
```

## Create tensor dataset in Python

Now we'll switch to Python to create the multiomics tensor and run barnacle analysis.

```{python setup-logging, eval=TRUE}
import sys
import os
from datetime import datetime

# Set up logging to both console and file
class Logger:
    def __init__(self, log_file):
        self.terminal = sys.stdout
        self.log = open(log_file, 'w')
        
    def write(self, message):
        self.terminal.write(message)
        self.log.write(message)
        self.log.flush()
        
    def flush(self):
        self.terminal.flush()
        self.log.flush()

# Create log file path
output_dir = r.output_dir
log_file = os.path.join(output_dir, f'barnacle_analysis_log_{datetime.now().strftime("%Y%m%d_%H%M%S")}.txt')

# Redirect stdout to both console and log file
sys.stdout = Logger(log_file)
sys.stderr = sys.stdout  # Also capture error messages

print("="*60)
print(f"BARNACLE ANALYSIS LOG")
print(f"Started: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
print(f"Log file: {log_file}")
print("="*60)
print()
```

```{python create-multiomics-tensor, eval=TRUE}
import pandas as pd
import numpy as np
import os
from pathlib import Path

# Set up paths
output_dir = r.output_dir
print(f"Working in output directory: {output_dir}")

# Load normalized data
apul_norm = pd.read_csv(os.path.join(output_dir, "apul_normalized_expression.csv"))
peve_norm = pd.read_csv(os.path.join(output_dir, "peve_normalized_expression.csv"))
ptua_norm = pd.read_csv(os.path.join(output_dir, "ptua_normalized_expression.csv"))

print("Loaded normalized data:")
print(f"Apul: {apul_norm.shape}")
print(f"Peve: {peve_norm.shape}")  
print(f"Ptua: {ptua_norm.shape}")

# Check which genes are common across all species
apul_genes = set(apul_norm['group_id'])
peve_genes = set(peve_norm['group_id'])
ptua_genes = set(ptua_norm['group_id'])

common_genes = apul_genes & peve_genes & ptua_genes
print(f"\nCommon genes across all species: {len(common_genes)}")

# Filter to common genes and align gene order
common_genes_list = sorted(list(common_genes))

apul_common = apul_norm[apul_norm['group_id'].isin(common_genes_list)].set_index('group_id').reindex(common_genes_list)
peve_common = peve_norm[peve_norm['group_id'].isin(common_genes_list)].set_index('group_id').reindex(common_genes_list)
ptua_common = ptua_norm[ptua_norm['group_id'].isin(common_genes_list)].set_index('group_id').reindex(common_genes_list)

print(f"\nFiltered to common genes:")
print(f"Apul: {apul_common.shape}")
print(f"Peve: {peve_common.shape}")
print(f"Ptua: {ptua_common.shape}")
```

## Parse sample information

```{python parse-sample-info, eval=TRUE}
# Parse sample names to extract sample information for each species independently
def parse_species_samples(columns, species_name):
    """Parse sample column names for a specific species"""
    sample_map = {}
    sample_ids = []
    timepoints = set()
    
    for col in columns:
        # Expected format: PREFIX-NUMBER-TP# (e.g., ACR-139-TP1, POR-216-TP1, POC-201-TP1)
        # Also support legacy format: PREFIX.NUMBER.TP# for backwards compatibility
        parts = col.split('-') if '-' in col else col.split('.')
        
        if len(parts) >= 3 and parts[-1].startswith('TP'):
            # Extract timepoint from last part (e.g., "TP1", "TP2", etc.)
            tp_part = parts[-1]
            timepoint = int(tp_part[2:])  # e.g., 1 from TP1
            
            # Sample ID is everything except the timepoint part
            sample_id_parts = parts[:-1]
            sample_id = '-'.join(sample_id_parts) if '-' in col else '.'.join(sample_id_parts)
            
            sample_map[(sample_id, timepoint)] = col
            if sample_id not in sample_ids:
                sample_ids.append(sample_id)
            timepoints.add(timepoint)
        else:
            print(f"Warning: Could not parse column name: {col}")
    
    return sample_map, sample_ids, sorted(timepoints)

# Parse sample information for each species independently
print("Parsing sample information for each species...")

species_data = {
    'apul': apul_common,
    'peve': peve_common, 
    'ptua': ptua_common
}

species_info = {}
all_timepoints = set()

for species, data in species_data.items():
    sample_map, sample_ids, timepoints = parse_species_samples(data.columns, species)
    species_info[species] = {
        'sample_map': sample_map,
        'sample_ids': sample_ids,
        'timepoints': timepoints,
        'n_samples': len(sample_ids)
    }
    all_timepoints.update(timepoints)
    
    print(f"{species}:")
    print(f"  Samples: {len(sample_ids)} ({sample_ids[:3]}...)")
    print(f"  Timepoints: {timepoints}")

common_timepoints = sorted(list(all_timepoints))
print(f"\nTimepoints found across all species: {common_timepoints}")

# Find the maximum number of samples to determine tensor dimensions
max_samples = max(info['n_samples'] for info in species_info.values())
print(f"Maximum samples in any species: {max_samples}")

# Print detailed sample structure
print(f"\nDetailed sample structure:")
for species, info in species_info.items():
    print(f"{species}: {info['n_samples']} samples × {len(info['timepoints'])} timepoints")
```

## Create 3D tensor (genes × species_samples × timepoints)

```{python create-3d-tensor, eval=TRUE}
# Create a 3D tensor: genes × (species_samples) × timepoints
# This flattens species and samples into a single dimension that Barnacle can handle

print("Creating 3D tensor by combining species and samples...")
print("Note: R-level filtering has already removed samples without all 4 timepoints")

# First, collect all actual sample-timepoint combinations that have data
all_sample_columns = []
sample_labels = []  # Track which sample belongs to which species
species_sample_map = {}  # Map from combined index to (species, sample_idx, sample_id)

sample_idx = 0
for species in ['apul', 'peve', 'ptua']:
    data = species_data[species]
    info = species_info[species]
    
    print(f"\nProcessing {species}:")
    for sample_id in info['sample_ids']:
        # Check if this sample has data for ALL required timepoints
        available_timepoints = []
        sample_timepoint_cols = []
        
        for timepoint in common_timepoints:
            if (sample_id, timepoint) in info['sample_map']:
                col_name = info['sample_map'][(sample_id, timepoint)]
                sample_timepoint_cols.append(col_name)
                available_timepoints.append(timepoint)
        
        # Only include samples that have ALL timepoints (should be all samples due to R filtering)
        if len(available_timepoints) == len(common_timepoints):
            all_sample_columns.extend(sample_timepoint_cols)
            sample_labels.append(f"{species}_{sample_id}")
            species_sample_map[sample_idx] = {
                'species': species,
                'sample_id': sample_id,
                'sample_idx_in_species': info['sample_ids'].index(sample_id)
            }
            sample_idx += 1
            print(f"  Added {sample_id} with {len(sample_timepoint_cols)} timepoints: {sorted(available_timepoints)}")
        else:
            print(f"  Skipped {sample_id} - incomplete timepoints: {sorted(available_timepoints)} (expected: {sorted(common_timepoints)})")

n_genes = len(common_genes_list)
n_combined_samples = len(sample_labels)
n_timepoints = len(common_timepoints)

print(f"\nCreating 3D tensor with shape: ({n_genes}, {n_combined_samples}, {n_timepoints})")
print(f"Combined samples from all species: {n_combined_samples}")

# Initialize tensor
tensor_3d = np.full((n_genes, n_combined_samples, n_timepoints), np.nan)

# Fill tensor
filled_count = 0
missing_count = 0

for combined_idx, sample_label in enumerate(sample_labels):
    species_info_map = species_sample_map[combined_idx]
    species = species_info_map['species']
    sample_id = species_info_map['sample_id']
    
    data = species_data[species]
    info = species_info[species]
    
    for time_idx, timepoint in enumerate(common_timepoints):
        if (sample_id, timepoint) in info['sample_map']:
            col_name = info['sample_map'][(sample_id, timepoint)]
            tensor_3d[:, combined_idx, time_idx] = data[col_name].values
            filled_count += 1
        else:
            missing_count += 1

# Check tensor statistics
n_missing = np.sum(np.isnan(tensor_3d))
n_total = tensor_3d.size
n_finite = np.sum(np.isfinite(tensor_3d))

print(f"\n=== TENSOR STATISTICS ===")
print(f"Tensor shape: {tensor_3d.shape}")
print(f"Total elements: {n_total}")
print(f"Finite values: {n_finite}")
print(f"Missing/NaN values: {n_missing}")
print(f"Missing percentage: {n_missing / n_total * 100:.2f}%")
print(f"Filled {filled_count} sample-timepoint combinations")
print(f"Missing {missing_count} sample-timepoint combinations")

# Check non-zero values among finite values
finite_values = tensor_3d[np.isfinite(tensor_3d)]
n_nonzero = np.sum(finite_values != 0)
print(f"Non-zero finite values: {n_nonzero}")
print(f"Zero finite values: {len(finite_values) - n_nonzero}")
print(f"Sparsity among finite values: {(len(finite_values) - n_nonzero) / len(finite_values) * 100:.2f}%")

# Save sample mapping for later interpretation
sample_mapping = pd.DataFrame([
    {
        'combined_index': i,
        'sample_label': label,
        'species': species_sample_map[i]['species'],
        'sample_id': species_sample_map[i]['sample_id']
    }
    for i, label in enumerate(sample_labels)
])
print(f"\nSample mapping:")
print(sample_mapping.head(10))
```



## Rank Comparison

**This workflow systematically compares multiple ranks to find the optimal number of components.**

The comparison runs decompositions for each specified rank and generates:
1. Metrics comparison table
2. Elbow plot (variance explained vs rank)
3. Multi-panel summary plots
4. Automated recommendation

**Results:** Each rank is saved to `rank_XX/` directory with complete factor matrices.

### Define helper functions for rank comparison

```{python rank-comparison-functions, eval=TRUE}
def run_single_rank_decomposition(tensor, rank, lambdas, n_iter_max, random_state=42):
    """
    Run sparse CP decomposition for a given rank.
    Returns dict with results and metrics.
    """
    from barnacle.decomposition import SparseCP
    
    print(f"\n{'='*60}")
    print(f"Running decomposition with rank = {rank}")
    print(f"{'='*60}")
    
    # Handle missing values
    tensor_filled = np.nan_to_num(tensor, nan=0.0)
    
    # Create model
    model = SparseCP(
        rank=rank,
        lambdas=lambdas,
        nonneg_modes=[0],
        norm_constraint=True,
        init='random',
        tol=1e-5,
        n_iter_max=n_iter_max,
        random_state=random_state,
        n_initializations=1
    )
    
    # Fit model
    try:
        decomposition = model.fit_transform(tensor_filled, verbose=0)
        
        # Calculate metrics
        metrics = calculate_decomposition_metrics(tensor_filled, decomposition, model)
        
        results = {
            'rank': rank,
            'decomposition': decomposition,
            'model': model,
            'metrics': metrics,
            'success': True
        }
        
        print(f"✓ Rank {rank} complete")
        print(f"  Final loss: {metrics['final_loss']:.6f}")
        print(f"  Variance explained: {metrics['variance_explained']*100:.2f}%")
        print(f"  Converged: {metrics['converged']}")
        
        return results
        
    except Exception as e:
        print(f"✗ Rank {rank} failed: {e}")
        return {
            'rank': rank,
            'success': False,
            'error': str(e)
        }


def calculate_decomposition_metrics(tensor, decomposition, model):
    """
    Calculate reconstruction quality and other metrics.
    """
    # Reconstruct tensor
    gene_factors = decomposition.factors[0]
    sample_factors = decomposition.factors[1]
    time_factors = decomposition.factors[2]
    weights = decomposition.weights
    
    rank = len(weights)
    reconstructed = np.zeros_like(tensor)
    
    for r in range(rank):
        outer_prod = np.outer(gene_factors[:, r], sample_factors[:, r])
        outer_prod = np.outer(outer_prod.flatten(), time_factors[:, r])
        outer_prod = outer_prod.reshape(tensor.shape)
        reconstructed += weights[r] * outer_prod
    
    # Calculate errors
    mse = np.mean((tensor - reconstructed) ** 2)
    rmse = np.sqrt(mse)
    
    # Relative error
    tensor_norm = np.linalg.norm(tensor)
    error_norm = np.linalg.norm(tensor - reconstructed)
    relative_error = error_norm / tensor_norm if tensor_norm > 0 else np.inf
    
    # Variance explained
    total_var = np.var(tensor)
    residual_var = np.var(tensor - reconstructed)
    variance_explained = 1 - (residual_var / total_var) if total_var > 0 else 0
    
    # Sparsity metrics
    gene_sparsity = np.mean(np.abs(gene_factors) < 0.1)
    sample_sparsity = np.mean(np.abs(sample_factors) < 0.1)
    time_sparsity = np.mean(np.abs(time_factors) < 0.1)
    
    # Component weight statistics
    weight_max = np.max(weights)
    weight_min = np.min(weights)
    weight_ratio = weight_max / weight_min if weight_min > 0 else np.inf
    weight_cv = np.std(weights) / np.mean(weights) if np.mean(weights) > 0 else np.inf
    
    # Convergence info
    # Note: Barnacle's SparseCP class does not have a converged_ attribute
    # We determine convergence by checking if iterations < max_iter
    n_iterations = len(model.loss_) if hasattr(model, 'loss_') and model.loss_ is not None else 0
    final_loss = model.loss_[-1] if hasattr(model, 'loss_') and model.loss_ is not None else np.nan
    # Converged if stopped before reaching max iterations (met tolerance criterion)
    converged = (n_iterations < model.n_iter_max) if n_iterations > 0 else False
    
    return {
        'mse': mse,
        'rmse': rmse,
        'relative_error': relative_error,
        'variance_explained': variance_explained,
        'gene_sparsity': gene_sparsity,
        'sample_sparsity': sample_sparsity,
        'time_sparsity': time_sparsity,
        'weight_max': weight_max,
        'weight_min': weight_min,
        'weight_ratio': weight_ratio,
        'weight_cv': weight_cv,
        'converged': converged,
        'n_iterations': n_iterations,
        'final_loss': final_loss
    }


def save_rank_results(results, output_dir, gene_ids, sample_labels, timepoint_labels, species_sample_map):
    """
    Save decomposition results for a specific rank.
    """
    rank = results['rank']
    rank_dir = os.path.join(output_dir, f'rank_{rank:02d}')
    os.makedirs(rank_dir, exist_ok=True)
    
    decomposition = results['decomposition']
    metrics = results['metrics']
    
    # Save factor matrices
    gene_factors = pd.DataFrame(
        decomposition.factors[0],
        index=gene_ids,
        columns=[f'Component_{i+1}' for i in range(rank)]
    )
    gene_factors.to_csv(os.path.join(rank_dir, 'gene_factors.csv'))
    
    sample_factors = pd.DataFrame(
        decomposition.factors[1],
        index=sample_labels,
        columns=[f'Component_{i+1}' for i in range(rank)]
    )
    sample_factors['Species'] = [species_sample_map[i]['species'] for i in range(len(sample_labels))]
    sample_factors['Sample_ID'] = [species_sample_map[i]['sample_id'] for i in range(len(sample_labels))]
    sample_factors.to_csv(os.path.join(rank_dir, 'sample_factors.csv'))
    
    time_factors = pd.DataFrame(
        decomposition.factors[2],
        index=timepoint_labels,
        columns=[f'Component_{i+1}' for i in range(rank)]
    )
    time_factors.to_csv(os.path.join(rank_dir, 'time_factors.csv'))
    
    # Save weights
    weights_df = pd.DataFrame({
        'Component': [f'Component_{i+1}' for i in range(rank)],
        'Weight': decomposition.weights
    })
    weights_df.to_csv(os.path.join(rank_dir, 'component_weights.csv'), index=False)
    
    # Save loss history
    if hasattr(results['model'], 'loss_') and results['model'].loss_ is not None:
        loss_df = pd.DataFrame({
            'Iteration': list(range(1, len(results['model'].loss_) + 1)),
            'Loss': results['model'].loss_
        })
        loss_df.to_csv(os.path.join(rank_dir, 'loss_history.csv'), index=False)
    
    # Save metadata
    metadata = {
        'rank': rank,
        'metrics': {k: float(v) if not isinstance(v, bool) else v for k, v in metrics.items()},
        'tensor_shape': list(decomposition.factors[0].shape[0:1]) + 
                       [decomposition.factors[1].shape[0]] + 
                       [decomposition.factors[2].shape[0]]
    }
    
    import json
    with open(os.path.join(rank_dir, 'metadata.json'), 'w') as f:
        json.dump(metadata, f, indent=2)
    
    print(f"  Saved results to: {rank_dir}")


print("Rank comparison helper functions loaded")
```

### Run barnacle

```{python barnacle-multiple-ranks, eval=TRUE}
# ===========================================
# CONFIGURE RANKS TO TEST
# ===========================================
# NOTE: This chunk creates a machine-readable log file
# (rank_comparison_log.json) that subsequent chunks will use.
# This allows you to skip re-running this chunk if you've
# already generated results.
# ===========================================
ranks_to_test = [5, 8, 10, 12, 15, 20, 25, 35, 45, 55, 65, 75]  # Modify this list as needed
# ===========================================

import os
import json
import datetime

print("="*60)
print("COMPARING MULTIPLE RANKS")
print("="*60)
print(f"Ranks to test: {ranks_to_test}")
print(f"This will run {len(ranks_to_test)} decompositions")
print(f"Estimated time: {len(ranks_to_test) * 2}-{len(ranks_to_test) * 5} minutes")
print("="*60)

# Run decompositions for each rank
all_rank_results = []

for test_rank in sorted(ranks_to_test):
    result = run_single_rank_decomposition(
        tensor=tensor_3d,
        rank=test_rank,
        lambdas=[0.1, 0.0, 0.1],
        n_iter_max=10000,
        random_state=41
    )
    all_rank_results.append(result)
    
    # Save individual results
    if result['success']:
        save_rank_results(
            result,
            output_dir,
            common_genes_list,
            sample_labels,
            [f'TP{tp}' for tp in common_timepoints],
            species_sample_map
        )

print(f"\n{'='*60}")
print(f"ALL RANKS COMPLETE")
print(f"{'='*60}")
print(f"Successful: {sum(1 for r in all_rank_results if r['success'])}/{len(all_rank_results)}")

# ===========================================
# SAVE COMPREHENSIVE LOG FILE
# ===========================================
import json
import datetime

log_file = os.path.join(output_dir, 'rank_comparison_log.json')

log_data = {
    'timestamp': datetime.datetime.now().isoformat(),
    'ranks_tested': ranks_to_test,
    'total_runs': len(all_rank_results),
    'successful_runs': sum(1 for r in all_rank_results if r['success']),
    'failed_runs': sum(1 for r in all_rank_results if not r['success']),
    'results': []
}

for result in all_rank_results:
    rank_info = {
        'rank': result['rank'],
        'success': result['success']
    }
    
    if result['success']:
        rank_info['metrics'] = {
            k: float(v) if isinstance(v, (int, float, np.number)) else bool(v) if isinstance(v, (bool, np.bool_)) else v
            for k, v in result['metrics'].items()
        }
        rank_info['output_dir'] = os.path.join(output_dir, f'rank_{result["rank"]:02d}')
        rank_info['has_loss_history'] = hasattr(result['model'], 'loss_') and result['model'].loss_ is not None
    else:
        rank_info['error'] = result.get('error', 'Unknown error')
    
    log_data['results'].append(rank_info)

with open(log_file, 'w') as f:
    json.dump(log_data, f, indent=2)

print(f"\nLog file saved to: {log_file}")
print("This log can be used by subsequent chunks to avoid re-running decompositions")
```

### Run barnacle with zero regularization (lambda = 0.0)

This chunk runs barnacle decomposition with all lambdas set to 0.0 (no regularization) to test how the model performs without sparsity constraints. Results are saved to a separate subdirectory for comparison.

```{python barnacle-zero-lambda, eval=TRUE}
# ===========================================
# CONFIGURE RANKS TO TEST - ZERO LAMBDA
# ===========================================
# This chunk tests decomposition with NO regularization
# All lambda values are set to 0.0
# ===========================================
ranks_to_test_zero_lambda = [10, 20, 35, 40]  # Specific ranks to test
# ===========================================

import os
import json
import datetime

# Create separate output directory for zero-lambda results
zero_lambda_output_dir = os.path.join(output_dir, 'zero_lambda_results')
os.makedirs(zero_lambda_output_dir, exist_ok=True)

print("="*60)
print("BARNACLE DECOMPOSITION - ZERO REGULARIZATION")
print("="*60)
print(f"Ranks to test: {ranks_to_test_zero_lambda}")
print(f"Lambda values: [0.0, 0.0, 0.0] (no regularization)")
print(f"Output directory: {zero_lambda_output_dir}")
print(f"This will run {len(ranks_to_test_zero_lambda)} decompositions")
print(f"Estimated time: {len(ranks_to_test_zero_lambda) * 2}-{len(ranks_to_test_zero_lambda) * 5} minutes")
print("="*60)

# Run decompositions for each rank with zero lambda
all_rank_results_zero_lambda = []

for test_rank in sorted(ranks_to_test_zero_lambda):
    result = run_single_rank_decomposition(
        tensor=tensor_3d,
        rank=test_rank,
        lambdas=[0.0, 0.0, 0.0],  # All lambdas set to 0.0
        n_iter_max=10000,
        random_state=41
    )
    all_rank_results_zero_lambda.append(result)
    
    # Save individual results
    if result['success']:
        save_rank_results(
            result,
            zero_lambda_output_dir,
            common_genes_list,
            sample_labels,
            [f'TP{tp}' for tp in common_timepoints],
            species_sample_map
        )

print(f"\n{'='*60}")
print(f"ZERO LAMBDA RANKS COMPLETE")
print(f"{'='*60}")
print(f"Successful: {sum(1 for r in all_rank_results_zero_lambda if r['success'])}/{len(all_rank_results_zero_lambda)}")

# ===========================================
# SAVE COMPREHENSIVE LOG FILE FOR ZERO LAMBDA
# ===========================================
import json
import datetime

log_file_zero_lambda = os.path.join(zero_lambda_output_dir, 'rank_comparison_log_zero_lambda.json')

log_data_zero_lambda = {
    'timestamp': datetime.datetime.now().isoformat(),
    'lambda_values': [0.0, 0.0, 0.0],
    'ranks_tested': ranks_to_test_zero_lambda,
    'total_runs': len(all_rank_results_zero_lambda),
    'successful_runs': sum(1 for r in all_rank_results_zero_lambda if r['success']),
    'failed_runs': sum(1 for r in all_rank_results_zero_lambda if not r['success']),
    'results': []
}

for result in all_rank_results_zero_lambda:
    rank_info = {
        'rank': result['rank'],
        'success': result['success']
    }
    
    if result['success']:
        rank_info['metrics'] = {
            k: float(v) if isinstance(v, (int, float, np.number)) else bool(v) if isinstance(v, (bool, np.bool_)) else v
            for k, v in result['metrics'].items()
        }
        rank_info['output_dir'] = os.path.join(zero_lambda_output_dir, f'rank_{result["rank"]:02d}')
        rank_info['has_loss_history'] = hasattr(result['model'], 'loss_') and result['model'].loss_ is not None
    else:
        rank_info['error'] = result.get('error', 'Unknown error')
    
    log_data_zero_lambda['results'].append(rank_info)

with open(log_file_zero_lambda, 'w') as f:
    json.dump(log_data_zero_lambda, f, indent=2)

print(f"\nLog file saved to: {log_file_zero_lambda}")
print(f"Results directory: {zero_lambda_output_dir}")
print("This log can be used to compare zero-lambda results with regularized decompositions")
```

### Plot Ranks vs Variance Explained (Zero Lambda Results)

This chunk creates a visualization showing the relationship between rank and variance explained for the zero-lambda decompositions, helping to identify the optimal rank based on the variance explained metric.

```{python plot-ranks-vs-variance-zero-lambda, eval=TRUE}
# ===========================================
# PLOT RANKS VS VARIANCE EXPLAINED - ZERO LAMBDA
# ===========================================
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

print(f"\n{'='*60}")
print("PLOTTING RANKS VS VARIANCE EXPLAINED (ZERO LAMBDA)")
print(f"{'='*60}")

# Extract data from zero-lambda results
ranks_zero_lambda = []
variance_explained_zero_lambda = []
successful_results = []

for result in all_rank_results_zero_lambda:
    if result['success'] and 'metrics' in result:
        ranks_zero_lambda.append(result['rank'])
        variance_explained_zero_lambda.append(result['metrics']['variance_explained'])
        successful_results.append(result)

if len(successful_results) == 0:
    print("WARNING: No successful decompositions found for zero-lambda results!")
    print("Cannot create variance plot.")
else:
    # Sort by rank for plotting
    sorted_data = sorted(zip(ranks_zero_lambda, variance_explained_zero_lambda))
    ranks_zero_lambda, variance_explained_zero_lambda = zip(*sorted_data)
    
    # Create the plot
    plt.figure(figsize=(10, 6))
    plt.plot(ranks_zero_lambda, [v*100 for v in variance_explained_zero_lambda], 
             'o-', linewidth=2.5, markersize=8, color='steelblue', markerfacecolor='lightblue', markeredgecolor='steelblue', markeredgewidth=2)
    
    # Customize the plot
    plt.xlabel('Rank', fontsize=12, fontweight='bold')
    plt.ylabel('Variance Explained (%)', fontsize=12, fontweight='bold')
    plt.title('Variance Explained vs Rank', fontsize=14, fontweight='bold')
    plt.grid(True, alpha=0.3)
    plt.xticks(ranks_zero_lambda)
    
    # Add value labels on points
    for rank, var_exp in zip(ranks_zero_lambda, variance_explained_zero_lambda):
        plt.annotate(f'{var_exp*100:.1f}%', 
                    (rank, var_exp*100), 
                    textcoords="offset points", 
                    xytext=(0,10), 
                    ha='center',
                    fontsize=10,
                    fontweight='bold')
    
    # Formatting
    plt.ylim(bottom=0)
    plt.tight_layout()
    
    # Save the plot
    plot_file = os.path.join(zero_lambda_output_dir, 'ranks_vs_variance_explained_zero_lambda.png')
    plt.savefig(plot_file, dpi=300, bbox_inches='tight')
    print(f"Plot saved to: {plot_file}")
    
    plt.show()
    
    # Print summary statistics
    print(f"\nSUMMARY - Zero Lambda Results:")
    print(f"{'='*50}")
    print(f"Number of successful decompositions: {len(successful_results)}")
    print(f"Ranks tested: {sorted(ranks_zero_lambda)}")
    print(f"Variance explained range: {min(variance_explained_zero_lambda)*100:.1f}% - {max(variance_explained_zero_lambda)*100:.1f}%")
    
    # Find best rank
    best_idx = np.argmax(variance_explained_zero_lambda)
    best_rank = ranks_zero_lambda[best_idx]
    best_variance = variance_explained_zero_lambda[best_idx]
    
    print(f"Best performing rank: {best_rank} (variance explained: {best_variance*100:.1f}%)")
    
    # Create a DataFrame for easy inspection
    results_df = pd.DataFrame({
        'Rank': ranks_zero_lambda,
        'Variance_Explained_Percent': [v*100 for v in variance_explained_zero_lambda]
    })
    
    print(f"\nDetailed Results:")
    print(results_df.to_string(index=False))

print(f"{'='*60}")
```

### Factor Match Score evaluation

This section evaluates each rank's decomposition quality using multiple Factor Match Score (FMS) metrics:

1. **TLViz FMS** (self-comparison): Sanity check only, not useful for rank selection
2. **Reconstruction FMS** (fast): Quick approximation based on tensor reconstruction correlation
3. **Synthetic FMS** (RECOMMENDED): Most rigorous - tests factor recovery from noisy data

The Synthetic FMS is the most reliable metric for determining optimal rank as it validates robustness.

```{python factor-match-score-evaluation, eval=TRUE}
# ===========================================
# FACTOR MATCH SCORE EVALUATION
# ===========================================
import numpy as np
import pandas as pd
import os

print(f"\n{'='*60}")
print("FACTOR MATCH SCORE (FMS) EVALUATION")
print(f"{'='*60}")
print("This analysis evaluates decomposition quality using multiple FMS metrics:")
print("  1. TLViz FMS: Self-comparison (sanity check only)")
print("  2. Reconstruction FMS: Fast approximation (correlation-based)")
print("  3. Synthetic FMS: RECOMMENDED (tests factor recovery robustness)")
print(f"{'='*60}")

def calculate_factor_match_score(tensor, decomposition, rank):
    """
    Calculate Factor Match Score using tlviz library.
    
    WARNING: This implementation compares the decomposition against itself,
    which will always yield FMS ≈ 1.0. This is only useful as a sanity check
    to verify the decomposition is internally consistent.
    
    For actual rank selection, use calculate_reconstruction_fms() instead,
    which compares the reconstruction against the original tensor.
    
    NOTE: May fail with scipy >= 1.14 due to _lazywhere removal.
    """
    try:
        # Import required libraries
        import tensorly as tl
        import tlviz.factor_tools
        
        # IMPORTANT: This creates a "reference" from the decomposition itself
        # This means we're comparing the decomposition to itself, which will
        # always give a perfect or near-perfect score (~1.0).
        # This is NOT a meaningful metric for rank selection!
        # It only validates that the decomposition object is internally consistent.
        reference_cp = tl.cp_tensor.CPTensor((decomposition.weights, decomposition.factors))
        
        # Calculate factor match score between decomposition and itself
        # Expected result: FMS ≈ 1.0 for all ranks (self-comparison)
        fms, perm = tlviz.factor_tools.factor_match_score(
            reference_cp,
            decomposition,
            return_permutation=True,
            allow_smaller_rank=True
        )
        
        return {
            'fms': float(fms),
            'permutation': perm,
            'success': True,
            'method': 'tlviz_self_consistency',
            'warning': 'Self-comparison: FMS ≈ 1.0 expected (not useful for rank selection)'
        }
        
    except ImportError as e:
        # Skip silently - this is expected with newer scipy versions
        error_msg = "SciPy/TLViz compatibility issue (expected with scipy >= 1.14)"
        if '_lazywhere' in str(e):
            error_msg = "SciPy version incompatibility (_lazywhere removed in scipy >= 1.14)"
        return {
            'fms': np.nan,
            'permutation': None,
            'success': False,
            'error': error_msg,
            'method': 'scipy_incompatibility'
        }
    except Exception as e:
        return {
            'fms': np.nan,
            'permutation': None,
            'success': False,
            'error': str(e),
            'method': 'failed'
        }

def calculate_proper_factor_match_score(tensor, decomposition, rank, noise_level=0.1, random_state=42):
    """
    RECOMMENDED METRIC: Synthetic validation approach for FMS.
    
    This creates a synthetic tensor from the decomposition, adds noise,
    re-decomposes it, then measures how well we recover the original factors.
    
    This is the most rigorous validation of decomposition quality because it tests
    the method's ability to recover known factors in the presence of noise.
    
    Higher scores indicate:
    - Better factor recovery from noisy data
    - More robust decomposition
    - Less sensitivity to noise
    - Better rank selection
    
    NOTE: Computationally expensive (requires additional decomposition per rank).
          
    """
    try:
        import tensorly as tl
        import tlviz.factor_tools
        from barnacle.decomposition import SparseCP
        
        # Step 1: Create synthetic "ground truth" tensor from current decomposition
        gene_factors = decomposition.factors[0]
        sample_factors = decomposition.factors[1]
        time_factors = decomposition.factors[2]
        weights = decomposition.weights
        
        # Reconstruct clean synthetic tensor
        synthetic_clean = np.zeros(tensor.shape)
        for r in range(len(weights)):
            outer_prod = np.outer(gene_factors[:, r], sample_factors[:, r])
            outer_prod = np.outer(outer_prod.flatten(), time_factors[:, r])
            outer_prod = outer_prod.reshape(tensor.shape)
            synthetic_clean += weights[r] * outer_prod
        
        # Step 2: Add realistic noise
        np.random.seed(random_state)
        noise_std = noise_level * np.std(synthetic_clean)
        noise = np.random.normal(0, noise_std, synthetic_clean.shape)
        synthetic_noisy = synthetic_clean + noise
        
        # Step 3: Decompose the noisy synthetic tensor
        # This tests: "Can we recover the factors from noisy data?"
        model_synthetic = SparseCP(
            rank=rank,
            lambdas=[0.1, 0.0, 0.1],
            nonneg_modes=[0],
            norm_constraint=True,
            init='random',
            tol=1e-5,
            n_iter_max=5000,  # Fewer iterations for speed
            random_state=random_state + 1,
            n_initializations=1
        )
        
        recovered_decomposition = model_synthetic.fit_transform(synthetic_noisy, verbose=0)
        
        # Step 4: Compare recovered factors to ground truth (original decomposition)
        # Create proper CPTensor objects
        ground_truth_cp = tl.cp_tensor.CPTensor((weights, [gene_factors, sample_factors, time_factors]))
        
        # MockDecomposition for recovered (same structure as before)
        class MockDecomposition:
            def __init__(self, weights, factors):
                self.weights = weights
                self.factors = factors
                self._tuple = (weights, factors)
            def __getitem__(self, index):
                return self._tuple[index]
            def __len__(self):
                return 2
        
        recovered_mock = MockDecomposition(recovered_decomposition.weights, recovered_decomposition.factors)
        
        # Calculate FMS between ground truth and recovered
        fms, perm = tlviz.factor_tools.factor_match_score(
            ground_truth_cp,
            recovered_mock,
            return_permutation=True,
            allow_smaller_rank=True
        )
        
        return {
            'fms': float(fms),
            'permutation': perm,
            'success': True,
            'method': 'synthetic_validation',
            'noise_level': noise_level,
            'recovery_quality': 'excellent' if fms > 0.9 else 'good' if fms > 0.7 else 'fair' if fms > 0.5 else 'poor'
        }
        
    except ImportError as e:
        error_msg = "TLViz not available (scipy incompatibility)"
        if '_lazywhere' in str(e):
            error_msg = "SciPy version incompatibility (_lazywhere removed in scipy >= 1.14)"
        return {
            'fms': np.nan,
            'success': False,
            'error': error_msg,
            'method': 'import_failed'
        }
    except Exception as e:
        return {
            'fms': np.nan,
            'success': False,
            'error': str(e),
            'method': 'failed'
        }

def calculate_reconstruction_fms(tensor, decomposition, rank):
    """
    Calculate FMS based on reconstruction quality.
    
    This compares the original tensor to its reconstruction from the decomposition,
    providing a measure of how well the decomposition captures the data.
    
    Higher scores indicate better reconstruction:
    - 1.0 = perfect reconstruction
    - 0.5 = no correlation
    - 0.0 = perfect negative correlation
    
    Note: This is a fast approximation. For more rigorous validation, use
    synthetic validation FMS (calculate_proper_factor_match_score).
    """
    try:
        # Create reconstructed tensor
        gene_factors = decomposition.factors[0]
        sample_factors = decomposition.factors[1] 
        time_factors = decomposition.factors[2]
        weights = decomposition.weights
        
        # Reconstruct tensor
        reconstructed = np.zeros_like(tensor)
        for r in range(len(weights)):
            outer_prod = np.outer(gene_factors[:, r], sample_factors[:, r])
            outer_prod = np.outer(outer_prod.flatten(), time_factors[:, r])
            outer_prod = outer_prod.reshape(tensor.shape)
            reconstructed += weights[r] * outer_prod
        
        # Calculate correlation-based FMS
        # Flatten tensors and calculate correlation
        orig_flat = tensor.flatten()
        recon_flat = reconstructed.flatten()
        
        # Remove NaN values for correlation calculation
        valid_mask = np.isfinite(orig_flat) & np.isfinite(recon_flat)
        orig_valid = orig_flat[valid_mask]
        recon_valid = recon_flat[valid_mask]
        
        if len(orig_valid) > 0:
            correlation = np.corrcoef(orig_valid, recon_valid)[0, 1]
            # Convert correlation to FMS-like score (0-1 scale)
            fms_score = (correlation + 1) / 2  # Map [-1,1] to [0,1]
        else:
            fms_score = 0.0
        
        return {
            'fms': float(fms_score),
            'correlation': float(correlation) if len(orig_valid) > 0 else np.nan,
            'success': True,
            'method': 'reconstruction_correlation',
            'valid_elements': len(orig_valid)
        }
        
    except Exception as e:
        print(f"  Error calculating reconstruction FMS for rank {rank}: {e}")
        return {
            'fms': np.nan,
            'success': False,
            'error': str(e),
            'method': 'failed'
        }

def evaluate_fms_for_rank(rank, output_dir, tensor, use_synthetic_validation=False):
    """
    Evaluate Factor Match Score for a specific rank.
    
    Parameters:
    -----------
    rank : int
        Rank to evaluate
    output_dir : str
        Directory containing rank results
    tensor : np.ndarray
        Original tensor
    use_synthetic_validation : bool, default=False
        If True, calculate synthetic validation FMS (RECOMMENDED but computationally expensive).
        If False, only calculate reconstruction FMS (fast approximation).
    """
    print(f"\n--- Evaluating FMS for Rank {rank} ---")
    
    rank_dir = os.path.join(output_dir, f'rank_{rank:02d}')
    
    # Check if rank directory exists
    if not os.path.exists(rank_dir):
        print(f"  Rank directory not found: {rank_dir}")
        return None
    
    # Load factor matrices
    try:
        gene_factors = pd.read_csv(os.path.join(rank_dir, 'gene_factors.csv'), index_col=0)
        sample_factors = pd.read_csv(os.path.join(rank_dir, 'sample_factors.csv'), index_col=0)
        time_factors = pd.read_csv(os.path.join(rank_dir, 'time_factors.csv'), index_col=0)
        weights_df = pd.read_csv(os.path.join(rank_dir, 'component_weights.csv'))
        
        # Remove metadata columns from sample_factors if present
        factor_cols = [col for col in sample_factors.columns if col.startswith('Component_')]
        sample_factor_values = sample_factors[factor_cols]
        
        print(f"  Loaded factor matrices - Genes: {gene_factors.shape}, Samples: {sample_factor_values.shape}, Time: {time_factors.shape}")
        
    except Exception as e:
        print(f"  Error loading factor matrices: {e}")
        return None
    
    # Create mock decomposition object that mimics tensorly CPTensor structure
    # TensorLy expects a tuple-like object: (weights, factors)
    # We'll create a simple class that supports both attribute and subscript access
    class MockDecomposition:
        def __init__(self, weights, factors):
            self.weights = weights
            self.factors = factors
            # Make it subscriptable for tensorly compatibility
            self._tuple = (weights, factors)
        
        def __getitem__(self, index):
            """Support subscript access like a tuple"""
            return self._tuple[index]
        
        def __len__(self):
            """Support len() for tuple-like behavior"""
            return 2
    
    weights = weights_df['Weight'].values
    factors = [gene_factors.values, sample_factor_values.values, time_factors.values]
    mock_decomposition = MockDecomposition(weights, factors)
    
    # Calculate FMS metrics
    # 1. TLViz FMS (self-comparison - sanity check only, not useful for rank selection)
    fms_result1 = calculate_factor_match_score(tensor, mock_decomposition, rank)
    
    # 2. Reconstruction FMS (fast approximation - useful for quick comparison)
    fms_result2 = calculate_reconstruction_fms(tensor, mock_decomposition, rank)
    
    # 3. Synthetic validation FMS (RECOMMENDED - most rigorous but computationally expensive)
    if use_synthetic_validation:
        print(f"  Running synthetic validation (this will take a long time)...")
        fms_result3 = calculate_proper_factor_match_score(tensor, mock_decomposition, rank)
    else:
        fms_result3 = {'success': False, 'method': 'skipped', 'fms': np.nan}
    
    # Combine results
    result = {
        'rank': rank,
        'fms_tlviz': fms_result1['fms'] if fms_result1['success'] else np.nan,
        'fms_reconstruction': fms_result2['fms'] if fms_result2['success'] else np.nan,
        'fms_synthetic': fms_result3['fms'] if fms_result3['success'] else np.nan,
        'reconstruction_correlation': fms_result2.get('correlation', np.nan),
        'tlviz_success': fms_result1['success'],
        'reconstruction_success': fms_result2['success'],
        'synthetic_success': fms_result3['success'],
        'tlviz_method': fms_result1.get('method', 'unknown'),
        'reconstruction_method': fms_result2.get('method', 'unknown'),
        'synthetic_method': fms_result3.get('method', 'unknown'),
        'valid_elements': fms_result2.get('valid_elements', 0)
    }
    
    # Print results
    # TLViz FMS note: This compares decomposition to itself (always ~1.0)
    if fms_result1['success']:
        warning_note = fms_result1.get('warning', '')
        if warning_note:
            print(f"  TLViz FMS: {fms_result1['fms']:.4f} ⚠️  (self-consistency check only)")
        else:
            print(f"  TLViz FMS: {fms_result1['fms']:.4f} ✅")
    elif 'scipy_incompatibility' not in fms_result1.get('method', ''):
        # Only print if it's a real error (not the expected scipy issue)
        print(f"  TLViz FMS: Failed - {fms_result1.get('error', 'unknown')} (method: {fms_result1.get('method', 'unknown')})")
    
    # Reconstruction FMS (fast approximation)
    if fms_result2['success']:
        print(f"  📊 RECONSTRUCTION FMS: {fms_result2['fms']:.4f} (correlation: {fms_result2.get('correlation', 'N/A'):.4f})")
        print(f"     ↳ Fast approximation for quick comparison")
        print(f"     ↳ Valid elements used: {fms_result2.get('valid_elements', 0):,}")
    else:
        print(f"  Reconstruction FMS: Failed ({fms_result2.get('error', 'unknown error')})")
    
    # Synthetic validation FMS (RECOMMENDED for rank selection)
    if use_synthetic_validation and fms_result3['success']:
        recovery = fms_result3.get('recovery_quality', 'unknown')
        noise = fms_result3.get('noise_level', 0.1)
        print(f"  ⭐ SYNTHETIC FMS: {fms_result3['fms']:.4f} ({recovery} recovery with {noise*100:.0f}% noise) ⭐")
        print(f"     ↳ RECOMMENDED metric for rank selection")
        print(f"     ↳ Tests factor recovery robustness with noise")
    elif use_synthetic_validation and not fms_result3['success']:
        print(f"  Synthetic FMS: Failed ({fms_result3.get('error', 'unknown error')})")
    
    # Save individual rank FMS results
    fms_file = os.path.join(rank_dir, 'factor_match_score.json')
    import json
    with open(fms_file, 'w') as f:
        json.dump({
            'rank': rank,
            'tlviz_fms': fms_result1,
            'reconstruction_fms': fms_result2,
            'synthetic_fms': fms_result3,
            'summary': result
        }, f, indent=2, default=str)
    
    print(f"  Saved FMS results to: factor_match_score.json")
    
    return result

# ===========================================
# EVALUATE FMS FOR ALL RANKS
# ===========================================

# Get output directory if not available
if 'output_dir' not in globals():
    output_dir = r.output_dir

# Get tensor if not available
if 'tensor_3d' not in globals():
    print("ERROR: tensor_3d not found in memory!")
    print("Please run the tensor creation chunks first.")
    tensor_3d = None

if tensor_3d is not None:
    # Load log file to get successful ranks
    log_file = os.path.join(output_dir, 'rank_comparison_log.json')
    
    if os.path.exists(log_file):
        import json
        with open(log_file, 'r') as f:
            log_data = json.load(f)
        
        successful_ranks = [r['rank'] for r in log_data['results'] if r['success']]
        print(f"Found {len(successful_ranks)} successful ranks to evaluate: {successful_ranks}")
        
    else:
        # Fallback: find rank directories
        import glob
        rank_dirs = sorted(glob.glob(os.path.join(output_dir, 'rank_*')))
        successful_ranks = []
        for rank_dir in rank_dirs:
            try:
                rank_num = int(os.path.basename(rank_dir).split('_')[1])
                successful_ranks.append(rank_num)
            except:
                pass
        print(f"Found {len(successful_ranks)} rank directories: {successful_ranks}")
    
    # ===========================================
    # CONFIGURATION: Synthetic Validation FMS
    # ===========================================
    # Set to True to enable the RECOMMENDED synthetic validation FMS metric.
    # This is the most rigorous method for rank selection as it tests factor
    # recovery from noisy data, providing the best indication of optimal rank.
    # 
    # COMPUTATION TIME: Adds ~1-2 minutes per rank
    # For 12 ranks, expect ~15-30 minutes additional runtime
    # 
    # If False, only reconstruction FMS will be calculated (faster but less rigorous)
    # ===========================================
    USE_SYNTHETIC_VALIDATION = True  # Change to False for faster (but less rigorous) analysis
    # ===========================================
    
    if USE_SYNTHETIC_VALIDATION:
        print("\n" + "="*60)
        print("⭐ SYNTHETIC VALIDATION ENABLED (RECOMMENDED)")
        print("="*60)
        print("This will test factor recovery robustness by:")
        print("  1. Creating synthetic tensor from each decomposition")
        print("  2. Adding noise to test robustness")
        print("  3. Re-decomposing to see if factors are recovered")
        print("  4. Comparing recovered factors to originals")
        print(f"\nExpected additional time: {len(successful_ranks) * 1.5:.0f}-{len(successful_ranks) * 2:.0f} minutes")
        print("This provides the most reliable metric for rank selection.")
        print("="*60 + "\n")
    
    # Evaluate FMS for each rank
    all_fms_results = []
    
    for rank in sorted(successful_ranks):
        fms_result = evaluate_fms_for_rank(rank, output_dir, tensor_3d, 
                                          use_synthetic_validation=USE_SYNTHETIC_VALIDATION)
        if fms_result is not None:
            all_fms_results.append(fms_result)
    
    # ===========================================
    # CREATE FMS COMPARISON TABLE AND PLOTS
    # ===========================================
    if len(all_fms_results) > 0:
        print(f"\n{'='*60}")
        print("FMS RESULTS SUMMARY")
        print(f"{'='*60}")
        
        # Create comparison dataframe
        fms_df = pd.DataFrame(all_fms_results)
        
        # Save comprehensive FMS results
        comparison_dir = os.path.join(output_dir, 'rank_comparison')
        os.makedirs(comparison_dir, exist_ok=True)
        
        fms_output_file = os.path.join(comparison_dir, 'factor_match_scores.csv')
        fms_df.to_csv(fms_output_file, index=False)
        print(f"Saved FMS comparison table to: {fms_output_file}")
        
        # Display results
        # Determine which columns to display based on what was calculated
        display_cols = ['rank', 'fms_tlviz', 'fms_reconstruction', 'reconstruction_correlation']
        if 'fms_synthetic' in fms_df.columns and fms_df['synthetic_success'].any():
            display_cols.append('fms_synthetic')
        
        print("\n" + "="*60)
        print("FACTOR MATCH SCORE RESULTS")
        print("="*60)
        print("⚠️  TLViz FMS: Self-comparison only (always ~1.0, not useful)")
        print("📊 Reconstruction FMS: Fast approximation (useful for quick comparison)")
        if 'fms_synthetic' in display_cols:
            print("⭐ SYNTHETIC FMS: RECOMMENDED for rank selection (tests robustness)")
        print("="*60)
        print(fms_df[display_cols].to_string(index=False, float_format='%.4f'))
        
        # Check if TLViz worked for any ranks
        tlviz_working = fms_df['tlviz_success'].any()
        if not tlviz_working:
            print("\n" + "="*60)
            print("NOTE: TLViz FMS unavailable (scipy version incompatibility)")
            print("="*60)
            print("This is expected with scipy >= 1.14.")
            print("Use Reconstruction FMS (fast) or Synthetic FMS (recommended).")
            print("="*60)
        elif tlviz_working and fms_df.loc[fms_df['tlviz_success'], 'fms_tlviz'].min() > 0.95:
            print("\n" + "="*60)
            print("⚠️  WARNING: TLViz FMS shows all values near 1.0")
            print("="*60)
            print("This indicates self-comparison (not useful for rank selection).")
            print("→ Use Synthetic FMS (recommended) or Reconstruction FMS (fast)")
            print("="*60)
        
        # Create visualization
        import matplotlib.pyplot as plt
        
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
        
        # Plot 1: Both FMS measures (but emphasize reconstruction if TLViz failed)
        valid_tlviz = fms_df['tlviz_success']
        valid_recon = fms_df['reconstruction_success']
        
        if valid_tlviz.any():
            ax1.plot(fms_df.loc[valid_tlviz, 'rank'], 
                    fms_df.loc[valid_tlviz, 'fms_tlviz'], 
                    'o-', linewidth=2, markersize=8, label='TLViz FMS', color='steelblue')
        
        if valid_recon.any():
            ax1.plot(fms_df.loc[valid_recon, 'rank'], 
                    fms_df.loc[valid_recon, 'fms_reconstruction'], 
                    's-', linewidth=2, markersize=8, label='Reconstruction FMS', color='coral')
        
        # Adjust title based on whether TLViz worked
        if not tlviz_working:
            title = 'Factor Match Score vs Rank\n(Reconstruction FMS: Fast Approximation)'
        else:
            title = 'Factor Match Score vs Rank\n(Higher = Better Factor Recovery)'
        
        ax1.set_xlabel('Rank (Number of Components)', fontsize=11)
        ax1.set_ylabel('Factor Match Score', fontsize=11)
        ax1.set_title(title, fontweight='bold')
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        ax1.set_ylim([0, 1.05])
        
        # Plot 2: Reconstruction correlation
        if valid_recon.any():
            ax2.plot(fms_df.loc[valid_recon, 'rank'], 
                    fms_df.loc[valid_recon, 'reconstruction_correlation'], 
                    'o-', linewidth=2, markersize=8, color='green')
            ax2.axhline(y=0.8, color='orange', linestyle='--', alpha=0.5, label='80% correlation')
            ax2.axhline(y=0.9, color='red', linestyle='--', alpha=0.5, label='90% correlation')
        
        ax2.set_xlabel('Rank (Number of Components)', fontsize=11)
        ax2.set_ylabel('Reconstruction Correlation', fontsize=11)
        ax2.set_title('Tensor Reconstruction Quality vs Rank\n(Correlation between original and reconstructed)', fontweight='bold')
        ax2.grid(True, alpha=0.3)
        ax2.legend()
        ax2.set_ylim([-1.05, 1.05])
        
        plt.tight_layout()
        fms_plot_file = os.path.join(comparison_dir, 'factor_match_scores_plot.png')
        plt.savefig(fms_plot_file, dpi=150, bbox_inches='tight')
        print(f"Saved FMS plot to: {fms_plot_file}")
        plt.close()
        
        # ===========================================
        # FMS-BASED RANK RECOMMENDATION
        # ===========================================
        print(f"\n{'='*60}")
        print("FMS-BASED RANK RECOMMENDATION")
        print(f"{'='*60}")
        
        # Find best ranks based on different criteria
        valid_df = fms_df[fms_df['reconstruction_success']]
        
        if len(valid_df) > 0:
            # Check for synthetic FMS results (RECOMMENDED metric)
            synthetic_valid = fms_df[fms_df['synthetic_success']]
            
            if len(synthetic_valid) > 0:
                # Synthetic FMS is available - use it as PRIMARY recommendation
                best_synthetic_idx = synthetic_valid['fms_synthetic'].idxmax()
                best_synthetic_rank = synthetic_valid.iloc[best_synthetic_idx]['rank']
                best_synthetic_score = synthetic_valid.iloc[best_synthetic_idx]['fms_synthetic']
                
                print(f"⭐ RECOMMENDED: Rank {int(best_synthetic_rank)} (Synthetic FMS: {best_synthetic_score:.4f})")
                print(f"   ↳ Based on factor recovery robustness with noise")
                
                # Show top 3 ranks by synthetic FMS
                top_synthetic = synthetic_valid.nlargest(3, 'fms_synthetic')[['rank', 'fms_synthetic']]
                print(f"\nTop 3 ranks by Synthetic FMS:")
                for idx, row in top_synthetic.iterrows():
                    print(f"  Rank {int(row['rank'])}: {row['fms_synthetic']:.4f}")
            else:
                print(f"⚠️  Synthetic FMS not calculated (set USE_SYNTHETIC_VALIDATION=True)")
                print(f"   ↳ Using Reconstruction FMS for recommendation (less rigorous)")
                
                # Fall back to reconstruction FMS
                best_recon_idx = valid_df['fms_reconstruction'].idxmax()
                best_recon_rank = valid_df.iloc[best_recon_idx]['rank']
                best_recon_score = valid_df.iloc[best_recon_idx]['fms_reconstruction']
                
                print(f"\n📊 Rank {int(best_recon_rank)} (Reconstruction FMS: {best_recon_score:.4f})")
                print(f"   ↳ Based on reconstruction correlation (fast approximation)")
            
            # Additional metrics for context
            print(f"\n--- Additional Metrics ---")
            
            # Best reconstruction correlation
            best_corr_idx = valid_df['reconstruction_correlation'].idxmax()
            best_corr_rank = valid_df.iloc[best_corr_idx]['rank']
            best_corr_score = valid_df.iloc[best_corr_idx]['reconstruction_correlation']
            print(f"Best Reconstruction Correlation: Rank {int(best_corr_rank)} ({best_corr_score:.4f})")
            
            # Check for TLViz results (informational only)
            tlviz_valid = fms_df[fms_df['tlviz_success']]
            if len(tlviz_valid) > 0:
                best_tlviz_idx = tlviz_valid['fms_tlviz'].idxmax()
                best_tlviz_rank = tlviz_valid.iloc[best_tlviz_idx]['rank']
                best_tlviz_score = tlviz_valid.iloc[best_tlviz_idx]['fms_tlviz']
                print(f"TLViz FMS (informational): Rank {int(best_tlviz_rank)} ({best_tlviz_score:.4f})")
            
            # Provide interpretation
            print(f"\n--- FMS Metrics Explained ---")
            if len(synthetic_valid) > 0:
                print(f"⭐ Synthetic FMS (RECOMMENDED):")
                print(f"   • Tests factor recovery from noisy data")
                print(f"   • Most rigorous validation of decomposition quality")
                print(f"   • Best metric for determining optimal rank")
                print(f"   • Values > 0.9 = excellent, > 0.7 = good, > 0.5 = fair")
            print(f"📊 Reconstruction FMS (fast approximation):")
            print(f"   • Based on correlation with original tensor")
            print(f"   • Useful for quick comparison")
            print(f"   • Less rigorous than synthetic validation")
            print(f"⚠️  TLViz FMS (self-comparison only):")
            print(f"   • Not useful for rank selection (always ~1.0)")
            print(f"   • Only validates internal consistency")
            
        else:
            print("No valid FMS results available for recommendation.")
    
    else:
        print("No FMS results to analyze.")

else:
    print("Skipping FMS evaluation - tensor not available.")

print(f"\n{'='*60}")
print("FACTOR MATCH SCORE EVALUATION COMPLETE")
print(f"{'='*60}")
if 'all_fms_results' in locals() and len(all_fms_results) > 0:
    working_recon = sum(1 for r in all_fms_results if r['reconstruction_success'])
    working_tlviz = sum(1 for r in all_fms_results if r['tlviz_success'])
    print(f"Summary: {working_recon}/{len(all_fms_results)} ranks have valid Reconstruction FMS")
    print(f"         {working_tlviz}/{len(all_fms_results)} ranks have valid TLViz FMS")
    if working_recon > 0:
        print(f"→ Use Reconstruction FMS and correlation for rank selection")
    else:
        print(f"→ No valid FMS results - check tensor and factor matrices")
```




# VISUALIZATION

## Plot rank comparisons

### Plot Synthetic FMS Elbow (Standalone)

This chunk can be run independently to generate the synthetic FMS elbow plot without re-running the expensive Factor Match Score evaluation. It loads the saved FMS results from the CSV file.

```{python plot-synthetic-fms-elbow-standalone, eval=TRUE}
# ===========================================
# STANDALONE SYNTHETIC FMS ELBOW PLOT
# ===========================================
# This chunk loads saved FMS results and creates the elbow plot
# Can be run independently without re-running FMS evaluation
# ===========================================

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Get output directory
output_dir = r.output_dir if 'r' in dir() else '../output/13.00-multiomics-barnacle'
comparison_dir = os.path.join(output_dir, 'rank_comparison')
fms_csv_file = os.path.join(comparison_dir, 'factor_match_scores.csv')

print(f"\n{'='*60}")
print("SYNTHETIC FMS ELBOW PLOT (STANDALONE)")
print(f"{'='*60}")
print(f"Loading FMS results from: {fms_csv_file}")

# Check if FMS results file exists
if not os.path.exists(fms_csv_file):
    print(f"\nERROR: FMS results file not found!")
    print(f"Expected location: {fms_csv_file}")
    print(f"\nPlease run the Factor Match Score evaluation chunk first.")
    print(f"It should create: {comparison_dir}/factor_match_scores.csv")
else:
    # Load FMS results
    fms_df = pd.read_csv(fms_csv_file)
    
    print(f"Loaded FMS data for {len(fms_df)} ranks")
    print(f"Columns available: {list(fms_df.columns)}")
    
    # Check if synthetic FMS data is available
    if 'fms_synthetic' not in fms_df.columns:
        print(f"\nERROR: No 'fms_synthetic' column found in FMS results!")
        print(f"Synthetic validation was not performed.")
        print(f"\nTo generate synthetic FMS data:")
        print(f"1. Set USE_SYNTHETIC_VALIDATION = True in the FMS evaluation chunk")
        print(f"2. Re-run the Factor Match Score evaluation chunk")
    else:
        # Check if any ranks have valid synthetic FMS
        valid_synthetic = fms_df['synthetic_success'].fillna(False)
        
        if not valid_synthetic.any():
            print(f"\nERROR: No valid synthetic FMS values found!")
            print(f"All synthetic validation attempts failed or were skipped.")
            print(f"\nTo generate synthetic FMS data:")
            print(f"1. Set USE_SYNTHETIC_VALIDATION = True in the FMS evaluation chunk")
            print(f"2. Re-run the Factor Match Score evaluation chunk")
        else:
            print(f"Found {valid_synthetic.sum()} ranks with valid synthetic FMS")
            
            # ===========================================
            # CREATE SYNTHETIC FMS ELBOW PLOT
            # ===========================================
            
            fig, ax = plt.subplots(1, 1, figsize=(10, 7))
            
            # Get data for plotting
            valid_tlviz = fms_df['tlviz_success'].fillna(False)
            valid_recon = fms_df['reconstruction_success'].fillna(False)
            
            # Synthetic FMS (RECOMMENDED - should show elbow)
            synthetic_data = fms_df.loc[valid_synthetic].sort_values('rank')
            ax.plot(synthetic_data['rank'], synthetic_data['fms_synthetic'], 
                    'o-', linewidth=3, markersize=10, label='⭐ Synthetic Validation FMS (RECOMMENDED)', 
                    color='#2ecc71', zorder=3)
            
            # Mark the peak (optimal rank)
            max_synthetic_idx = synthetic_data['fms_synthetic'].idxmax()
            optimal_rank = synthetic_data.loc[max_synthetic_idx, 'rank']
            optimal_score = synthetic_data.loc[max_synthetic_idx, 'fms_synthetic']
            ax.plot(optimal_rank, optimal_score, 'r*', markersize=20, 
                    label=f'★ Optimal: Rank {int(optimal_rank)} (FMS={optimal_score:.4f})',
                    zorder=4)
            
            # Reconstruction FMS (fast approximation - typically increases monotonically)
            if valid_recon.any():
                recon_data = fms_df.loc[valid_recon].sort_values('rank')
                ax.plot(recon_data['rank'], recon_data['fms_reconstruction'], 
                        's--', linewidth=2, markersize=8, label='📊 Reconstruction FMS (fast)', 
                        color='#3498db', alpha=0.7, zorder=2)
            
            # TLViz FMS (if available - self-comparison, typically near 1.0)
            if valid_tlviz.any():
                tlviz_data = fms_df.loc[valid_tlviz].sort_values('rank')
                ax.plot(tlviz_data['rank'], tlviz_data['fms_tlviz'], 
                        '^:', linewidth=1.5, markersize=6, label='⚠️  TLViz FMS (sanity check only)', 
                        color='#95a5a6', alpha=0.5, zorder=1)
            
            ax.set_xlabel('Rank (Number of Components)', fontsize=12, fontweight='bold')
            ax.set_ylabel('Factor Match Score', fontsize=12, fontweight='bold')
            ax.set_title('Factor Match Score Comparison: Finding Optimal Rank\n' +
                         'Synthetic FMS (RECOMMENDED): Peaks at optimal rank, then declines with overfitting',
                         fontsize=13, fontweight='bold')
            ax.grid(True, alpha=0.3, linestyle='--')
            ax.legend(loc='best', fontsize=10)
            ax.set_ylim([0, 1.05])
            
            # Add quality threshold lines
            ax.axhline(y=0.9, color='green', linestyle='--', alpha=0.3, linewidth=1)
            ax.text(ax.get_xlim()[1]*0.95, 0.91, 'Excellent (>0.9)', 
                    ha='right', va='bottom', fontsize=9, color='green', alpha=0.7)
            ax.axhline(y=0.7, color='orange', linestyle='--', alpha=0.3, linewidth=1)
            ax.text(ax.get_xlim()[1]*0.95, 0.71, 'Good (>0.7)', 
                    ha='right', va='bottom', fontsize=9, color='orange', alpha=0.7)
            ax.axhline(y=0.5, color='red', linestyle='--', alpha=0.3, linewidth=1)
            ax.text(ax.get_xlim()[1]*0.95, 0.51, 'Fair (>0.5)', 
                    ha='right', va='bottom', fontsize=9, color='red', alpha=0.7)
            
            plt.tight_layout()
            synthetic_plot_file = os.path.join(comparison_dir, 'synthetic_fms_elbow_plot.png')
            plt.savefig(synthetic_plot_file, dpi=150, bbox_inches='tight')
            print(f"\n✓ Saved Synthetic FMS elbow plot to: {synthetic_plot_file}")
            plt.close()
            
            # Print summary
            print(f"\n{'='*60}")
            print("SYNTHETIC FMS ELBOW ANALYSIS")
            print(f"{'='*60}")
            print(f"Optimal Rank: {int(optimal_rank)} (FMS = {optimal_score:.4f})")
            
            # Show all ranks with their scores
            print(f"\nComplete Synthetic FMS Results:")
            print(f"{'Rank':<8} {'Synthetic FMS':<15} {'Quality':<12} {'Status'}")
            print(f"{'-'*8} {'-'*15} {'-'*12} {'-'*20}")
            for idx, row in synthetic_data.iterrows():
                rank_val = int(row['rank'])
                fms_val = row['fms_synthetic']
                
                if fms_val > 0.9:
                    quality = "Excellent"
                elif fms_val > 0.7:
                    quality = "Good"
                elif fms_val > 0.5:
                    quality = "Fair"
                else:
                    quality = "Poor"
                
                status = "★ OPTIMAL" if rank_val == optimal_rank else ""
                print(f"{rank_val:<8} {fms_val:<15.4f} {quality:<12} {status}")
            
            print(f"\nInterpretation:")
            print(f"  • Synthetic FMS tests factor recovery from noisy data")
            print(f"  • Peak at rank {int(optimal_rank)} indicates best balance:")
            print(f"    - Sufficient complexity to capture patterns")
            print(f"    - Robust enough to avoid overfitting noise")
            print(f"  • Higher ranks show declining FMS (overfitting)")
            print(f"  • Lower ranks show lower FMS (underfitting)")
            
            # Additional statistics
            if len(synthetic_data) > 1:
                decline_after_optimal = after_optimal['fms_synthetic'].mean() if len(after_optimal) > 0 else np.nan
                improvement_to_optimal = optimal_score - before_optimal['fms_synthetic'].max() if len(before_optimal) > 0 else np.nan
                
                print(f"\nStatistics:")
                if not np.isnan(improvement_to_optimal):
                    print(f"  • Improvement to optimal: +{improvement_to_optimal:.4f}")
                if not np.isnan(decline_after_optimal):
                    print(f"  • Average FMS after optimal: {decline_after_optimal:.4f}")
                    print(f"  • Decline from peak: -{(optimal_score - decline_after_optimal):.4f}")
            
            print(f"\n{'='*60}")

print(f"\nSYNTHETIC FMS ELBOW PLOT GENERATION COMPLETE")
print(f"{'='*60}")
```

### 1. Comparison Plots

Located in `barnacle_factors/rank_comparison/`:

- **`elbow_plot.png`** - Shows variance explained vs rank (MOST IMPORTANT)
  - Look for the "elbow" where gains diminish
  - Marginal gains are annotated on the plot
  
- **`rank_comparison_summary.png`** - 6-panel overview
- **`loss_trajectories.png`** - Convergence speed comparison

### 2. Metrics Table

File: `barnacle_factors/rank_comparison/rank_comparison_metrics.csv`

Contains all quantitative metrics for comparison.

### 3. Individual Rank Results

Each rank gets its own directory (`rank_05/`, `rank_10/`, etc.) with full factor matrices.

```{python plot-rank-comparison, eval=TRUE}
# ===========================================
# This chunk reads from rank_comparison_log.json
# and does NOT require variables from previous chunks
# ===========================================
import matplotlib.pyplot as plt
import seaborn as sns
import json
import glob
import os
import pandas as pd
import numpy as np

print(f"\n{'='*60}")
print("Creating rank comparison visualizations")
print(f"{'='*60}")

# Check if output_dir exists, if not get it from R
if 'output_dir' not in globals():
    print("Getting output_dir from R environment...")
    output_dir = r.output_dir
    print(f"  output_dir = {output_dir}")

comparison_dir = os.path.join(output_dir, 'rank_comparison')
os.makedirs(comparison_dir, exist_ok=True)

# ===========================================
# LOAD RESULTS FROM LOG FILE
# ===========================================
log_file = os.path.join(output_dir, 'rank_comparison_log.json')

if os.path.exists(log_file):
    print(f"Loading results from log file: {log_file}")
    with open(log_file, 'r') as f:
        log_data = json.load(f)
    
    print(f"  Log timestamp: {log_data['timestamp']}")
    print(f"  Total runs: {log_data['total_runs']}")
    print(f"  Successful: {log_data['successful_runs']}")
    print(f"  Failed: {log_data['failed_runs']}")
    
    # Reconstruct all_rank_results from log
    all_rank_results = []
    for result_info in log_data['results']:
        if result_info['success']:
            all_rank_results.append({
                'rank': result_info['rank'],
                'metrics': result_info['metrics'],
                'success': True
            })
        else:
            all_rank_results.append({
                'rank': result_info['rank'],
                'success': False,
                'error': result_info.get('error', 'Unknown error')
            })
    
    print(f"  Loaded {len(all_rank_results)} results from log file")
    
elif 'all_rank_results' not in globals():
    # Fallback: Load from individual metadata files
    print("Log file not found. Loading from individual rank directories...")
    all_rank_results = []
    
    # Find all rank_* directories
    rank_dirs = sorted(glob.glob(os.path.join(output_dir, 'rank_*')))
    print(f"  Found {len(rank_dirs)} rank directories")
    
    for rank_dir in rank_dirs:
        metadata_file = os.path.join(rank_dir, 'metadata.json')
        if os.path.exists(metadata_file):
            with open(metadata_file, 'r') as f:
                metadata = json.load(f)
                all_rank_results.append({
                    'rank': metadata['rank'],
                    'metrics': metadata['metrics'],
                    'success': True
                })
    
    print(f"  Loaded {len(all_rank_results)} results from disk")
else:
    print("Using all_rank_results from memory")

# Extract successful results
successful = [r for r in all_rank_results if r['success']]
if len(successful) == 0:
    print("No successful decompositions to compare!")
else:
    ranks = [r['rank'] for r in successful]
    
    # Create comparison dataframe
    comparison_df = pd.DataFrame([
        {
            'Rank': r['rank'],
            **r['metrics']
        }
        for r in successful
    ])
    
    # Save comparison table
    comparison_df.to_csv(os.path.join(comparison_dir, 'rank_comparison_metrics.csv'), index=False)
    print(f"  Saved metrics table")
    
    # Create multi-panel comparison plot
    fig, axes = plt.subplots(2, 3, figsize=(18, 10))
    fig.suptitle('Rank Comparison: Key Metrics', fontsize=16, fontweight='bold')
    
    # 1. Variance Explained
    ax = axes[0, 0]
    ax.plot(ranks, comparison_df['variance_explained'], 'o-', linewidth=2, markersize=8, color='steelblue')
    ax.axhline(y=0.8, color='green', linestyle='--', alpha=0.5, label='80% threshold')
    ax.axhline(y=0.6, color='orange', linestyle='--', alpha=0.5, label='60% threshold')
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Variance Explained', fontsize=11)
    ax.set_title('Variance Explained vs Rank', fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.legend()
    ax.set_ylim([0, 1.05])
    
    # 2. Relative Error
    ax = axes[0, 1]
    ax.plot(ranks, comparison_df['relative_error'], 'o-', linewidth=2, markersize=8, color='coral')
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Relative Reconstruction Error', fontsize=11)
    ax.set_title('Reconstruction Error vs Rank', fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 3. Final Loss
    ax = axes[0, 2]
    ax.plot(ranks, comparison_df['final_loss'], 'o-', linewidth=2, markersize=8, color='purple')
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Final Loss', fontsize=11)
    ax.set_title('Final Loss vs Rank', fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 4. Gene Sparsity
    ax = axes[1, 0]
    ax.plot(ranks, comparison_df['gene_sparsity'], 'o-', linewidth=2, markersize=8, color='green')
    ax.axhline(y=0.9, color='red', linestyle='--', alpha=0.5, label='90% sparse')
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Gene Factor Sparsity', fontsize=11)
    ax.set_title('Gene Sparsity vs Rank', fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.legend()
    ax.set_ylim([0, 1.05])
    
    # 5. Weight Distribution (CV)
    ax = axes[1, 1]
    ax.plot(ranks, comparison_df['weight_cv'], 'o-', linewidth=2, markersize=8, color='darkorange')
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Weight Coefficient of Variation', fontsize=11)
    ax.set_title('Component Weight Variability vs Rank', fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 6. Iterations to Convergence
    ax = axes[1, 2]
    colors = ['green' if c else 'red' for c in comparison_df['converged']]
    ax.bar(ranks, comparison_df['n_iterations'], color=colors, alpha=0.7)
    ax.set_xlabel('Rank (Number of Components)', fontsize=11)
    ax.set_ylabel('Iterations to Convergence', fontsize=11)
    ax.set_title('Convergence Speed vs Rank\n(Green=Converged, Red=Max Iterations)', fontweight='bold')
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    output_path = os.path.join(comparison_dir, 'rank_comparison_summary.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    print(f"  Saved: rank_comparison_summary.png")
    plt.close()
    
    # Create elbow plot
    fig, ax = plt.subplots(figsize=(10, 6))
    
    var_exp = comparison_df['variance_explained'].values
    ax.plot(ranks, var_exp, 'o-', linewidth=3, markersize=10, color='steelblue', label='Variance Explained')
    
    # Calculate marginal gains
    if len(ranks) > 1:
        marginal_gains = np.diff(var_exp)
        for i in range(len(marginal_gains)):
            ax.annotate(f'+{marginal_gains[i]*100:.1f}%',
                       xy=(ranks[i+1], var_exp[i+1]),
                       xytext=(10, -10),
                       textcoords='offset points',
                       fontsize=9,
                       color='darkgreen',
                       bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.3))
    
    ax.axhline(y=0.8, color='green', linestyle='--', alpha=0.5, linewidth=2, label='80% threshold (excellent)')
    ax.axhline(y=0.6, color='orange', linestyle='--', alpha=0.5, linewidth=2, label='60% threshold (good)')
    ax.set_xlabel('Rank (Number of Components)', fontsize=13, fontweight='bold')
    ax.set_ylabel('Variance Explained', fontsize=13, fontweight='bold')
    ax.set_title('Elbow Plot: Finding Optimal Rank\n(Look for diminishing returns)', 
                fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.legend(fontsize=11)
    ax.set_ylim([0, 1.05])
    ax.set_xticks(ranks)
    
    plt.tight_layout()
    output_path = os.path.join(comparison_dir, 'elbow_plot.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    print(f"  Saved: elbow_plot.png")
    plt.close()
    
    # Loss trajectories
    fig, ax = plt.subplots(figsize=(10, 6))
    
    loss_plotted = False
    for r in successful:
        # Try to get loss from model object first
        if 'model' in r and hasattr(r['model'], 'loss_') and r['model'].loss_ is not None:
            loss = r['model'].loss_
            ax.plot(range(1, len(loss)+1), loss, linewidth=2, label=f'Rank {r["rank"]}', alpha=0.8)
            loss_plotted = True
        else:
            # Load from disk
            rank_dir = os.path.join(output_dir, f'rank_{r["rank"]:02d}')
            loss_file = os.path.join(rank_dir, 'loss_history.csv')
            if os.path.exists(loss_file):
                loss_df = pd.read_csv(loss_file)
                ax.plot(loss_df['Iteration'], loss_df['Loss'], linewidth=2, label=f'Rank {r["rank"]}', alpha=0.8)
                loss_plotted = True
    
    if loss_plotted:
        ax.set_xlabel('Iteration', fontsize=12)
        ax.set_ylabel('Loss', fontsize=12)
        ax.set_title('Loss Convergence Trajectories by Rank', fontsize=14, fontweight='bold')
        ax.grid(True, alpha=0.3)
        ax.legend(fontsize=10)
        ax.set_yscale('log')
        
        plt.tight_layout()
        output_path = os.path.join(comparison_dir, 'loss_trajectories.png')
        plt.savefig(output_path, dpi=150, bbox_inches='tight')
        print(f"  Saved: loss_trajectories.png")
    else:
        print(f"  Skipped loss_trajectories.png (no loss history available)")
    
    plt.close()
    
    print(f"\nAll comparison plots saved to: {comparison_dir}")
```



# INDIVIDUAL RANK VISUALIZATIONS

## Detailed plots for Rank 10

This section creates comprehensive visualizations for rank 10 (the optimal rank determined from previous analysis).

```{python individual-rank-visualizations, eval=TRUE}
# ===========================================
# CREATE DETAILED VISUALIZATIONS FOR RANK 10
# ===========================================
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

print(f"\n{'='*60}")
print("CREATING RANK 10 VISUALIZATIONS")
print(f"{'='*60}")

# Get output directory
if 'output_dir' not in globals():
    output_dir = r.output_dir

# Set rank to analyze
SELECTED_RANK = 10
TOP_GENES_HEATMAP = 50  # Number of top genes to show in heatmap

# Process only rank 10
rank = SELECTED_RANK
print(f"\n--- Rank {rank} ---")
rank_dir = os.path.join(output_dir, f'rank_{rank:02d}')

# Load factor matrices
gene_factors_file = os.path.join(rank_dir, 'gene_factors.csv')
time_factors_file = os.path.join(rank_dir, 'time_factors.csv')
weights_file = os.path.join(rank_dir, 'component_weights.csv')

if not all([os.path.exists(f) for f in [gene_factors_file, time_factors_file, weights_file]]):
    print(f"ERROR: Missing factor files for rank {rank}")
    print(f"Please ensure barnacle decomposition has been run for rank {rank}")
else:
    gene_factors = pd.read_csv(gene_factors_file, index_col=0)
    time_factors = pd.read_csv(time_factors_file, index_col=0)
    weights = pd.read_csv(weights_file)
    
    # ===========================================
    # PLOT 1: GENE HEATMAPS (Multiple selection criteria)
    # ===========================================
    print(f"  Creating gene-component heatmaps with different selection criteria...")
    
    # ===========================================
    # 1A: Sum of absolute loadings (broadly active genes)
    # ===========================================
    print(f"    1A: Sum of absolute loadings (broadly active genes)...")
    gene_importance_sum = gene_factors.abs().sum(axis=1)
    top_genes_sum = gene_importance_sum.nlargest(TOP_GENES_HEATMAP).index
    top_gene_factors_sum = gene_factors.loc[top_genes_sum]
    
    # Save data to CSV
    csv_path = os.path.join(rank_dir, f'gene_heatmap_1A_sum_top{TOP_GENES_HEATMAP}_data.csv')
    top_gene_factors_sum.to_csv(csv_path)
    print(f"      Saved data: gene_heatmap_1A_sum_top{TOP_GENES_HEATMAP}_data.csv")
    
    fig, ax = plt.subplots(figsize=(max(12, rank*0.5), 10))
    sns.heatmap(
        top_gene_factors_sum,
        cmap='RdBu_r',
        center=0,
        cbar_kws={'label': 'Component Loading'},
        xticklabels=True,
        yticklabels=True,
        ax=ax
    )
    ax.set_xlabel('Components', fontsize=12, fontweight='bold')
    ax.set_ylabel('Gene ID', fontsize=12, fontweight='bold')
    ax.set_title(f'Rank {rank}: Top {TOP_GENES_HEATMAP} Genes by Sum of Absolute Loadings\n' + 
                f'(Genes that participate across multiple components)',
                fontsize=13, fontweight='bold', pad=15)
    plt.tight_layout()
    output_path = os.path.join(rank_dir, f'gene_heatmap_1A_sum_top{TOP_GENES_HEATMAP}.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"      Saved plot: gene_heatmap_1A_sum_top{TOP_GENES_HEATMAP}.png")
    
    # ===========================================
    # 1B: Maximum absolute loading (component-specific genes)
    # ===========================================
    print(f"    1B: Maximum absolute loading (component-specific genes)...")
    gene_importance_max = gene_factors.abs().max(axis=1)
    top_genes_max = gene_importance_max.nlargest(TOP_GENES_HEATMAP).index
    top_gene_factors_max = gene_factors.loc[top_genes_max]
    
    # Save data to CSV
    csv_path = os.path.join(rank_dir, f'gene_heatmap_1B_max_top{TOP_GENES_HEATMAP}_data.csv')
    top_gene_factors_max.to_csv(csv_path)
    print(f"      Saved data: gene_heatmap_1B_max_top{TOP_GENES_HEATMAP}_data.csv")
    
    fig, ax = plt.subplots(figsize=(max(12, rank*0.5), 10))
    sns.heatmap(
        top_gene_factors_max,
        cmap='RdBu_r',
        center=0,
        cbar_kws={'label': 'Component Loading'},
        xticklabels=True,
        yticklabels=True,
        ax=ax
    )
    ax.set_xlabel('Components', fontsize=12, fontweight='bold')
    ax.set_ylabel('Gene ID', fontsize=12, fontweight='bold')
    ax.set_title(f'Rank {rank}: Top {TOP_GENES_HEATMAP} Genes by Maximum Loading\n' + 
                f'(Genes highly specific to individual components)',
                fontsize=13, fontweight='bold', pad=15)
    plt.tight_layout()
    output_path = os.path.join(rank_dir, f'gene_heatmap_1B_max_top{TOP_GENES_HEATMAP}.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"      Saved plot: gene_heatmap_1B_max_top{TOP_GENES_HEATMAP}.png")
    
    # ===========================================
    # 1C: Mean absolute loading (consistently active genes)
    # ===========================================
    print(f"    1C: Mean absolute loading (consistently active genes)...")
    gene_importance_mean = gene_factors.abs().mean(axis=1)
    top_genes_mean = gene_importance_mean.nlargest(TOP_GENES_HEATMAP).index
    top_gene_factors_mean = gene_factors.loc[top_genes_mean]
    
    # Save data to CSV
    csv_path = os.path.join(rank_dir, f'gene_heatmap_1C_mean_top{TOP_GENES_HEATMAP}_data.csv')
    top_gene_factors_mean.to_csv(csv_path)
    print(f"      Saved data: gene_heatmap_1C_mean_top{TOP_GENES_HEATMAP}_data.csv")
    
    fig, ax = plt.subplots(figsize=(max(12, rank*0.5), 10))
    sns.heatmap(
        top_gene_factors_mean,
        cmap='RdBu_r',
        center=0,
        cbar_kws={'label': 'Component Loading'},
        xticklabels=True,
        yticklabels=True,
        ax=ax
    )
    ax.set_xlabel('Components', fontsize=12, fontweight='bold')
    ax.set_ylabel('Gene ID', fontsize=12, fontweight='bold')
    ax.set_title(f'Rank {rank}: Top {TOP_GENES_HEATMAP} Genes by Mean Absolute Loading\n' + 
                f'(Genes with consistent moderate activity across components)',
                fontsize=13, fontweight='bold', pad=15)
    plt.tight_layout()
    output_path = os.path.join(rank_dir, f'gene_heatmap_1C_mean_top{TOP_GENES_HEATMAP}.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"      Saved plot: gene_heatmap_1C_mean_top{TOP_GENES_HEATMAP}.png")
    
    # ===========================================
    # 1D: Variance of loadings (genes with high specificity)
    # ===========================================
    print(f"    1D: Variance of loadings (component-specific markers)...")
    gene_variance = gene_factors.var(axis=1)
    top_genes_var = gene_variance.nlargest(TOP_GENES_HEATMAP).index
    top_gene_factors_var = gene_factors.loc[top_genes_var]
    
    # Save data to CSV
    csv_path = os.path.join(rank_dir, f'gene_heatmap_1D_variance_top{TOP_GENES_HEATMAP}_data.csv')
    top_gene_factors_var.to_csv(csv_path)
    print(f"      Saved data: gene_heatmap_1D_variance_top{TOP_GENES_HEATMAP}_data.csv")
    
    fig, ax = plt.subplots(figsize=(max(12, rank*0.5), 10))
    sns.heatmap(
        top_gene_factors_var,
        cmap='RdBu_r',
        center=0,
        cbar_kws={'label': 'Component Loading'},
        xticklabels=True,
        yticklabels=True,
        ax=ax
    )
    ax.set_xlabel('Components', fontsize=12, fontweight='bold')
    ax.set_ylabel('Gene ID', fontsize=12, fontweight='bold')
    ax.set_title(f'Rank {rank}: Top {TOP_GENES_HEATMAP} Genes by Loading Variance\n' + 
                f'(Genes with high specificity - strong in some components, weak in others)',
                fontsize=13, fontweight='bold', pad=15)
    plt.tight_layout()
    output_path = os.path.join(rank_dir, f'gene_heatmap_1D_variance_top{TOP_GENES_HEATMAP}.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"      Saved plot: gene_heatmap_1D_variance_top{TOP_GENES_HEATMAP}.png")
    
    # ===========================================
    # Summary comparison
    # ===========================================
    print(f"\n    Gene selection summary:")
    print(f"      Sum method:      Favors genes active across many components")
    print(f"      Max method:      Favors genes strongly specific to one component")
    print(f"      Mean method:     Favors genes with consistent moderate activity")
    print(f"      Variance method: Favors genes with high component specificity")
    
    # Check overlap between methods
    overlap_sum_max = len(set(top_genes_sum) & set(top_genes_max))
    overlap_sum_mean = len(set(top_genes_sum) & set(top_genes_mean))
    overlap_max_var = len(set(top_genes_max) & set(top_genes_var))
    
    print(f"\n    Gene overlap between methods:")
    print(f"      Sum & Max:      {overlap_sum_max}/{TOP_GENES_HEATMAP} genes in common")
    print(f"      Sum & Mean:     {overlap_sum_mean}/{TOP_GENES_HEATMAP} genes in common")
    print(f"      Max & Variance: {overlap_max_var}/{TOP_GENES_HEATMAP} genes in common")
    
    # ===========================================
    # EXPORT TOP 500 GENES PER COMPONENT
    # ===========================================
    print(f"\n    Exporting top 500 genes per component...")
    TOP_GENES_PER_COMPONENT = 500
    
    for component in gene_factors.columns:
        # Get top genes by absolute loading for this component
        component_loadings = gene_factors[component].abs().sort_values(ascending=False)
        top_genes_for_component = component_loadings.head(TOP_GENES_PER_COMPONENT)
        
        # Create dataframe with gene ID and loading value (both absolute and signed)
        top_genes_df = pd.DataFrame({
            'gene_id': top_genes_for_component.index,
            'absolute_loading': top_genes_for_component.values,
            'signed_loading': gene_factors.loc[top_genes_for_component.index, component].values
        })
        
        # Save to CSV
        csv_filename = f'top_{TOP_GENES_PER_COMPONENT}_genes_{component}.csv'
        csv_path = os.path.join(rank_dir, csv_filename)
        top_genes_df.to_csv(csv_path, index=False)
        
    print(f"      Saved top {TOP_GENES_PER_COMPONENT} genes for each of {len(gene_factors.columns)} components")
    print(f"      Files: top_{TOP_GENES_PER_COMPONENT}_genes_Component_X.csv")
    
    # ===========================================
    # PLOT 2: TIME COMPONENT LINE PLOT
    # ===========================================
    print(f"  Creating time-component line plot...")
    
    fig, ax = plt.subplots(figsize=(12, 8))
    
    # Define different marker shapes for each component
    markers = ['o', 's', '^', 'D', 'v', '<', '>', 'p', '*', 'h', '+', 'x', '8', 'H', 'd', '|', '_']
    
    # Plot each component
    for i, col in enumerate(time_factors.columns):
        marker_style = markers[i % len(markers)]  # Cycle through markers if more components than markers
        ax.plot(time_factors.index, time_factors[col], 
               marker=marker_style, linewidth=2, markersize=8,
               label=f'{col}', alpha=0.8)
    
    ax.set_xlabel('Timepoint', fontsize=12, fontweight='bold')
    ax.set_ylabel('Component Loading', fontsize=12, fontweight='bold')
    ax.set_title(f'Rank {rank}: Component Loadings Across Timepoints',
                fontsize=13, fontweight='bold', pad=15)
    ax.grid(True, alpha=0.3)
    ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=9)
    ax.axhline(y=0, color='black', linestyle='-', linewidth=0.5, alpha=0.3)
    
    plt.tight_layout()
    output_path = os.path.join(rank_dir, 'time_component_lineplot.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"    Saved: time_component_lineplot.png")
    
    # ===========================================
    # PLOT 3: SAMPLE HEATMAP (Samples x Components)
    # ===========================================
    print(f"  Creating sample-component heatmap...")
    
    # Load sample factors with metadata
    sample_factors_file = os.path.join(rank_dir, 'sample_factors.csv')
    sample_factors = pd.read_csv(sample_factors_file, index_col=0)
    
    # Extract only the component columns
    component_cols = [col for col in sample_factors.columns if col.startswith('Component_')]
    sample_component_data = sample_factors[component_cols]
    
    # Create labels that include species and sample info
    sample_labels = [f"{row['Species']}_{row['Sample_ID']}" for idx, row in sample_factors.iterrows()]
    
    # Add sample labels as index for clustering
    sample_component_data.index = sample_labels
    
    # Create clustered heatmap using seaborn's clustermap
    # This performs hierarchical clustering on samples (rows) based on their component loadings
    print(f"    Performing hierarchical clustering on samples...")
    
    g = sns.clustermap(
        sample_component_data,
        cmap='RdBu_r',
        center=0,
        figsize=(max(12, rank*0.5), max(10, len(sample_labels)*0.3)),
        cbar_kws={'label': 'Component Loading'},
        row_cluster=True,   # Cluster samples (rows)
        col_cluster=False,  # Don't cluster components (columns) - keep original order
        xticklabels=True,
        yticklabels=True,
        dendrogram_ratio=0.15,
        cbar_pos=(0.02, 0.8, 0.03, 0.15),
        linewidths=0,
        method='average',  # Linkage method for hierarchical clustering
        metric='euclidean'  # Distance metric
    )
    
    # Adjust labels and title
    g.ax_heatmap.set_xlabel('Components', fontsize=12, fontweight='bold')
    g.ax_heatmap.set_ylabel('Sample (Species_SampleID)', fontsize=12, fontweight='bold')
    g.fig.suptitle(f'Rank {rank}: Sample Loadings Across Components\n' + 
                   f'(Hierarchically clustered by component loading similarity)',
                   fontsize=13, fontweight='bold', y=0.98)
    
    plt.tight_layout()
    output_path = os.path.join(rank_dir, 'sample_component_heatmap.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"    Saved: sample_component_heatmap.png")
    
    # ===========================================
    # PLOT 4: COMPONENT CORRELATION HEATMAP
    # ===========================================
    print(f"  Creating component correlation heatmap...")
    
    # Calculate correlation between components based on gene loadings
    component_corr = gene_factors.corr()
    
    fig, ax = plt.subplots(figsize=(max(10, rank*0.4), max(8, rank*0.35)))
    
    sns.heatmap(
        component_corr,
        annot=True if rank <= 20 else False,
        fmt='.2f' if rank <= 20 else '.1f',
        cmap='coolwarm',
        center=0,
        square=True,
        cbar_kws={'label': 'Correlation'},
        ax=ax
    )
    
    ax.set_title(f'Rank {rank}: Component Correlation Matrix\n' +
                '(Based on gene loadings - low correlation = good separation)',
                fontsize=13, fontweight='bold', pad=15)
    
    plt.tight_layout()
    output_path = os.path.join(rank_dir, 'component_correlation_heatmap.png')
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"    Saved: component_correlation_heatmap.png")
    
    print(f"\n{'='*60}")
    print("RANK 10 VISUALIZATIONS COMPLETE")
    print(f"{'='*60}")
    print(f"\nThe following plots were created in rank_{rank:02d}/ directory:")
    print(f"  1A. gene_heatmap_1A_sum_top{TOP_GENES_HEATMAP}.png - Broadly active genes (sum)")
    print(f"  1B. gene_heatmap_1B_max_top{TOP_GENES_HEATMAP}.png - Component-specific genes (max)")
    print(f"  1C. gene_heatmap_1C_mean_top{TOP_GENES_HEATMAP}.png - Consistently active genes (mean)")
    print(f"  1D. gene_heatmap_1D_variance_top{TOP_GENES_HEATMAP}.png - High specificity genes (variance)")
    print(f"  2. time_component_lineplot.png - Component patterns over time")
    print(f"  3. sample_component_heatmap.png - Individual samples across components")
    print(f"  4. component_correlation_heatmap.png - Component independence check")
```





## Close log file

```{python close-log, eval=TRUE}
from datetime import datetime

print()
print("="*60)
print(f"ANALYSIS COMPLETE")
print(f"Finished: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
print(f"All output saved to: {output_dir}")
print(f"Log file: {log_file}")
print("="*60)

# Restore stdout and close log file
if hasattr(sys.stdout, 'log'):
    sys.stdout.log.close()
    sys.stdout = sys.stdout.terminal
```

# GOSLIM ANNOTATION ANALYSIS

This section annotates the top 500 genes for components 5, 6, 7, and 9 with GOslim Biological Process terms from the rank_10 directory.

**Note:** Only Biological Process (BP) GOslim terms are included in this analysis. Cellular Component and Molecular Function GOslim terms are excluded.

## Setup and load annotations

```{r goslim-setup, eval=TRUE}
library(tidyverse)
library(ggplot2)

# Define paths
rank_10_dir <- file.path(output_dir, "rank_10")
goslim_output_dir <- file.path(output_dir, "rank_10", "goslim_analysis")
dir.create(goslim_output_dir, showWarnings = FALSE, recursive = TRUE)

# Path to goslim_generic.obo file
goslim_obo_file <- "../output/12-ortho-annot/run_20250831_172744/goslim_generic.obo"

# Load the ortholog groups annotation file
ortho_annot <- read.csv(ortholog_groups_file, stringsAsFactors = FALSE)

cat(strrep("=", 60), "\n")
cat("GOSLIM ANNOTATION SETUP\n")
cat(strrep("=", 60), "\n\n")
cat("Rank 10 directory:", rank_10_dir, "\n")
cat("GOslim output directory:", goslim_output_dir, "\n")
cat("Loaded ortholog annotations:", nrow(ortho_annot), "ortholog groups\n\n")
```

## Extract Biological Process GOslim IDs

```{r extract-bp-goslims, eval=TRUE}
# Function to parse OBO file and extract Biological Process GOslim IDs
extract_bp_goslim_ids <- function(obo_file) {
  if (!file.exists(obo_file)) {
    cat("WARNING: GOslim OBO file not found:", obo_file, "\n")
    cat("Proceeding without BP filtering (all GOslims will be included)\n")
    return(NULL)
  }
  
  # Read OBO file
  obo_lines <- readLines(obo_file)
  
  # Parse OBO file
  bp_goslim_ids <- c()
  current_term <- list(id = NULL, namespace = NULL, is_goslim = FALSE)
  
  for (line in obo_lines) {
    line <- trimws(line)
    
    if (line == "[Term]") {
      # Save previous term if it's a BP GOslim
      if (!is.null(current_term$id) && 
          current_term$namespace == "biological_process" && 
          current_term$is_goslim) {
        bp_goslim_ids <- c(bp_goslim_ids, current_term$id)
      }
      # Reset for new term
      current_term <- list(id = NULL, namespace = NULL, is_goslim = FALSE)
      
    } else if (grepl("^id:", line)) {
      current_term$id <- sub("^id:\\s*", "", line)
      
    } else if (grepl("^namespace:", line)) {
      current_term$namespace <- sub("^namespace:\\s*", "", line)
      
    } else if (grepl("^subset:\\s*goslim_generic", line)) {
      current_term$is_goslim <- TRUE
    }
  }
  
  # Don't forget the last term
  if (!is.null(current_term$id) && 
      current_term$namespace == "biological_process" && 
      current_term$is_goslim) {
    bp_goslim_ids <- c(bp_goslim_ids, current_term$id)
  }
  
  return(unique(bp_goslim_ids))
}

# Extract BP GOslim IDs
bp_goslim_ids <- extract_bp_goslim_ids(goslim_obo_file)

if (!is.null(bp_goslim_ids)) {
  cat("Extracted Biological Process GOslim terms:", length(bp_goslim_ids), "\n")
  cat("First 10 BP GOslim IDs:", paste(head(bp_goslim_ids, 10), collapse = ", "), "\n\n")
} else {
  cat("NOTE: BP filtering disabled - all GOslim terms will be used\n\n")
}
```

## Load top 500 genes for components 5, 6, 7, and 9

```{r load-top-genes, eval=TRUE}
# Components to analyze
components_to_analyze <- c(5, 6, 7, 9)

# Load top 500 genes for each component
top_genes_list <- list()

cat(strrep("=", 60), "\n")
cat("LOADING TOP 500 GENES\n")
cat(strrep("=", 60), "\n\n")

for (comp_num in components_to_analyze) {
  file_path <- file.path(rank_10_dir, 
                         paste0("top_500_genes_Component_", comp_num, ".csv"))
  
  if (!file.exists(file_path)) {
    cat("WARNING: File not found:", file_path, "\n")
    next
  }
  
  # Load the top genes file
  top_genes_df <- read.csv(file_path, stringsAsFactors = FALSE)
  
  # Rename gene_id column to group_id for consistency with annotation file
  if ("gene_id" %in% colnames(top_genes_df)) {
    top_genes_df <- top_genes_df %>%
      rename(group_id = gene_id)
  }
  
  top_genes_list[[paste0("Component_", comp_num)]] <- top_genes_df
  
  cat("Component", comp_num, ":\n")
  cat("  Loaded", nrow(top_genes_df), "genes\n")
  cat("  Loading range:", 
      round(min(top_genes_df$absolute_loading), 4), "to", 
      round(max(top_genes_df$absolute_loading), 4), "\n")
  cat("  Top 3 genes:", paste(head(top_genes_df$group_id, 3), collapse = ", "), "\n\n")
}

cat("Total components loaded:", length(top_genes_list), "\n\n")
```

## Annotate with GOslim terms

```{r annotate-goslim, eval=TRUE}
# Function to parse GOslim IDs and names, filtering for BP terms only
parse_goslim <- function(goslim_ids_str, goslim_names_str, bp_filter_ids = NULL) {
  if (is.na(goslim_ids_str) || goslim_ids_str == "" || 
      is.na(goslim_names_str) || goslim_names_str == "") {
    return(NULL)
  }
  
  # Split GOslim IDs and names
  ids <- strsplit(goslim_ids_str, ";\\s*")[[1]]
  names <- strsplit(goslim_names_str, ";\\s*")[[1]]
  
  # Clean up whitespace
  ids <- trimws(ids)
  names <- trimws(names)
  
  # Return data frame if lengths match
  if (length(ids) == length(names) && length(ids) > 0) {
    result_df <- data.frame(
      goslim_id = ids,
      goslim_name = names,
      stringsAsFactors = FALSE
    )
    
    # Filter for Biological Process GOslims if bp_filter_ids is provided
    if (!is.null(bp_filter_ids)) {
      result_df <- result_df %>%
        filter(goslim_id %in% bp_filter_ids)
      
      # Return NULL if no BP terms remain after filtering
      if (nrow(result_df) == 0) {
        return(NULL)
      }
    }
    
    return(result_df)
  }
  
  return(NULL)
}

# Annotate each component's top genes
annotated_components <- list()

cat(strrep("=", 60), "\n")
cat("ANNOTATING WITH GOSLIM TERMS (Biological Process only)\n")
cat(strrep("=", 60), "\n\n")

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(top_genes_list)) {
    cat("Component", comp_num, ": Skipping (not loaded)\n\n")
    next
  }
  
  top_genes_df <- top_genes_list[[comp_col]]
  
  cat("Component", comp_num, ":\n")
  
  # Merge with annotation data
  annotated_df <- top_genes_df %>%
    left_join(ortho_annot %>% 
                select(group_id, goslim_ids, goslim_names),
              by = "group_id")
  
  # Expand GOslim terms (one row per GOslim term)
  # Filter for Biological Process terms only
  goslim_expanded <- list()
  
  for (i in seq_len(nrow(annotated_df))) {
    row <- annotated_df[i, ]
    # Pass bp_goslim_ids to filter for BP terms only
    goslim_data <- parse_goslim(row$goslim_ids, row$goslim_names, bp_goslim_ids)
    
    if (!is.null(goslim_data)) {
      goslim_data$group_id <- row$group_id
      goslim_data$loading <- row$absolute_loading
      goslim_expanded[[i]] <- goslim_data
    }
  }
  
  # Combine all GOslim annotations
  if (length(goslim_expanded) > 0) {
    goslim_all <- bind_rows(goslim_expanded)
    
    # Count occurrences of each GOslim term
    goslim_counts <- goslim_all %>%
      group_by(goslim_id, goslim_name) %>%
      summarize(
        count = n(),
        mean_loading = mean(loading),
        .groups = "drop"
      ) %>%
      arrange(desc(count))
    
    annotated_components[[comp_col]] <- list(
      top_genes = annotated_df,
      goslim_counts = goslim_counts
    )
    
    cat("  Genes with BP GOslim:", 
        length(unique(goslim_all$group_id)), "/", nrow(top_genes_df), "\n")
    cat("  Unique BP GOslim terms:", nrow(goslim_counts), "\n")
    cat("  Top 5 BP GOslim terms:\n")
    top5 <- head(goslim_counts, 5)
    for (j in seq_len(nrow(top5))) {
      cat("    ", j, ". ", top5$goslim_name[j], " (", top5$goslim_id[j], ", n=", top5$count[j], ")\n", sep = "")
    }
  } else {
    cat("  WARNING: No BP GOslim annotations found!\n")
    annotated_components[[comp_col]] <- list(
      top_genes = annotated_df,
      goslim_counts = data.frame()
    )
  }
  cat("\n")
}
```

## Save GOslim annotation tables

```{r save-goslim-tables, eval=TRUE}
cat(strrep("=", 60), "\n")
cat("SAVING GOSLIM TABLES (Biological Process only)\n")
cat(strrep("=", 60), "\n\n")

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(annotated_components)) {
    next
  }
  
  # Save top genes with annotations
  top_genes_file <- file.path(goslim_output_dir, 
                               paste0("component_", comp_num, "_top500_genes_goslim_BP_annotated.csv"))
  write.csv(annotated_components[[comp_col]]$top_genes, 
            top_genes_file, 
            row.names = FALSE)
  cat("Saved:", basename(top_genes_file), "\n")
  
  # Save GOslim counts
  if (nrow(annotated_components[[comp_col]]$goslim_counts) > 0) {
    goslim_counts_file <- file.path(goslim_output_dir, 
                                    paste0("component_", comp_num, "_goslim_BP_counts.csv"))
    write.csv(annotated_components[[comp_col]]$goslim_counts, 
              goslim_counts_file, 
              row.names = FALSE)
    cat("Saved:", basename(goslim_counts_file), "\n")
  }
}

cat("\nAll tables saved to:", goslim_output_dir, "\n\n")
```

## Create GOslim bar charts

```{r plot-goslim-barcharts, eval=TRUE, fig.width=10, fig.height=8}
cat(strrep("=", 60), "\n")
cat("CREATING GOSLIM BAR CHARTS (Biological Process only)\n")
cat(strrep("=", 60), "\n\n")

# Color palette for components
component_colors <- c(
  "5" = "#E64B35",
  "6" = "#4DBBD5", 
  "7" = "#00A087",
  "9" = "#F39B7F"
)

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(annotated_components)) {
    next
  }
  
  goslim_counts <- annotated_components[[comp_col]]$goslim_counts
  
  if (nrow(goslim_counts) == 0) {
    cat("Component", comp_num, ": No BP GOslim data to plot\n")
    next
  }
  
  # Take top 20 GOslim terms for visualization
  top_goslim <- goslim_counts %>%
    head(20) %>%
    mutate(goslim_name = reorder(goslim_name, count))
  
  # Get total number of genes for subtitle
  top_genes <- annotated_components[[comp_col]]$top_genes
  total_genes <- nrow(top_genes)
  
  # Create bar chart
  p <- ggplot(top_goslim, aes(x = goslim_name, y = count)) +
    geom_bar(stat = "identity", fill = component_colors[as.character(comp_num)], 
             alpha = 0.8) +
    geom_text(aes(label = count), hjust = -0.2, size = 3.5) +
    coord_flip() +
    labs(
      title = paste0("Component ", comp_num, ": Top 20 GOslim Biological Process Terms"),
      subtitle = paste0("Top 500 genes (genes can have multiple BP GOslim annotations)"),
      x = "GOslim Biological Process Term",
      y = "Number of Ortholog Groups"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      plot.subtitle = element_text(size = 11, color = "gray30"),
      axis.text.y = element_text(size = 10),
      axis.text.x = element_text(size = 10),
      axis.title = element_text(size = 11, face = "bold"),
      panel.grid.major.y = element_blank(),
      panel.grid.minor = element_blank()
    ) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.15)))
  
  # Save plot
  plot_file <- file.path(goslim_output_dir, 
                         paste0("component_", comp_num, "_goslim_BP_barchart.png"))
  ggsave(plot_file, p, width = 10, height = 8, dpi = 300)
  cat("Saved:", basename(plot_file), "\n")
}

cat("\nAll bar charts saved to:", goslim_output_dir, "\n\n")
```

## Create GOslim percentage bar charts

```{r plot-goslim-percentage-barcharts, eval=TRUE, fig.width=10, fig.height=8}
cat(strrep("=", 60), "\n")
cat("CREATING GOSLIM PERCENTAGE BAR CHARTS (Biological Process only)\n")
cat(strrep("=", 60), "\n\n")

# Color palette for components
component_colors <- c(
  "5" = "#E64B35",
  "6" = "#4DBBD5", 
  "7" = "#00A087",
  "9" = "#F39B7F"
)

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(annotated_components)) {
    next
  }
  
  goslim_counts <- annotated_components[[comp_col]]$goslim_counts
  
  if (nrow(goslim_counts) == 0) {
    cat("Component", comp_num, ": No BP GOslim data to plot\n")
    next
  }
  
  # Calculate total genes with BP GOslim annotations
  # Note: We're using the full top 500 genes as the denominator
  top_genes <- annotated_components[[comp_col]]$top_genes
  total_genes <- nrow(top_genes)
  
  # Calculate percentage for each GOslim term
  # (genes can have multiple GOslim terms, so percentages may sum to >100%)
  goslim_percentages <- goslim_counts %>%
    mutate(percentage = (count / total_genes) * 100) %>%
    head(20) %>%
    mutate(goslim_name = reorder(goslim_name, percentage))
  
  # Create percentage bar chart
  p <- ggplot(goslim_percentages, aes(x = goslim_name, y = percentage)) +
    geom_bar(stat = "identity", fill = component_colors[as.character(comp_num)], 
             alpha = 0.8) +
    geom_text(aes(label = sprintf("%.1f%%", percentage)), hjust = -0.2, size = 3.5) +
    coord_flip() +
    labs(
      title = paste0("Component ", comp_num, ": Top 20 GOslim Biological Process Terms (%)"),
      subtitle = paste0("Percentage of 500 genes with each BP GOslim (genes can have multiple annotations)"),
      x = "GOslim Biological Process Term",
      y = "Percentage of Ortholog Groups (%)"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      plot.subtitle = element_text(size = 11, color = "gray30"),
      axis.text.y = element_text(size = 10),
      axis.text.x = element_text(size = 10),
      axis.title = element_text(size = 11, face = "bold"),
      panel.grid.major.y = element_blank(),
      panel.grid.minor = element_blank()
    ) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.15)))
  
  # Save plot
  plot_file <- file.path(goslim_output_dir, 
                         paste0("component_", comp_num, "_goslim_BP_barchart_percentage.png"))
  ggsave(plot_file, p, width = 10, height = 8, dpi = 300)
  cat("Saved:", basename(plot_file), "\n")
  
  # Also save the percentage data
  percentage_file <- file.path(goslim_output_dir, 
                               paste0("component_", comp_num, "_goslim_BP_percentages.csv"))
  goslim_percentages_full <- goslim_counts %>%
    mutate(percentage = (count / total_genes) * 100) %>%
    arrange(desc(percentage))
  write.csv(goslim_percentages_full, percentage_file, row.names = FALSE)
  cat("Saved:", basename(percentage_file), "\n\n")
}

cat("All percentage bar charts saved to:", goslim_output_dir, "\n\n")
```

## Create comprehensive GOslim percentage bar charts (all 72 BP terms)

```{r plot-goslim-all-percentage-barcharts, eval=TRUE, fig.width=12, fig.height=22}
cat(strrep("=", 60), "\n")
cat("CREATING COMPREHENSIVE GOSLIM PERCENTAGE BAR CHARTS (All 72 BP GOslim terms)\n")
cat(strrep("=", 60), "\n\n")

# Color palette for components
component_colors <- c(
  "5" = "#E64B35",
  "6" = "#4DBBD5", 
  "7" = "#00A087",
  "9" = "#F39B7F"
)

# Create a reference dataframe with ALL 72 BP GOslim terms
# Use the bp_goslim_ids that were already extracted, and get names from OBO file
extract_bp_goslim_names <- function(obo_file, bp_ids) {
  if (!file.exists(obo_file)) {
    return(data.frame(goslim_id = bp_ids, goslim_name = bp_ids, stringsAsFactors = FALSE))
  }
  
  obo_lines <- readLines(obo_file)
  id_to_name <- list()
  current_id <- NULL
  current_name <- NULL
  
  for (line in obo_lines) {
    line <- trimws(line)
    
    if (line == "[Term]") {
      # Save previous term
      if (!is.null(current_id) && !is.null(current_name)) {
        id_to_name[[current_id]] <- current_name
      }
      current_id <- NULL
      current_name <- NULL
      
    } else if (grepl("^id:", line)) {
      current_id <- sub("^id:\\s*", "", line)
      
    } else if (grepl("^name:", line)) {
      current_name <- sub("^name:\\s*", "", line)
    }
  }
  
  # Don't forget the last term
  if (!is.null(current_id) && !is.null(current_name)) {
    id_to_name[[current_id]] <- current_name
  }
  
  # Create dataframe with all BP GOslim IDs and their names
  data.frame(
    goslim_id = bp_ids,
    goslim_name = sapply(bp_ids, function(id) {
      if (id %in% names(id_to_name)) {
        id_to_name[[id]]
      } else {
        id  # Use ID as name if name not found
      }
    }, USE.NAMES = FALSE),
    stringsAsFactors = FALSE
  )
}

# Get names for all 72 BP GOslim IDs
all_bp_goslim_ref <- extract_bp_goslim_names(goslim_obo_file, bp_goslim_ids)

cat("BP GOslim reference table created with", nrow(all_bp_goslim_ref), "terms\n\n")

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(annotated_components)) {
    next
  }
  
  goslim_counts <- annotated_components[[comp_col]]$goslim_counts
  
  # Get total number of genes
  top_genes <- annotated_components[[comp_col]]$top_genes
  total_genes <- nrow(top_genes)
  
  # Create complete dataframe with ALL 72 BP GOslim terms
  # Join with actual counts, filling missing terms with 0
  goslim_all_percentages <- all_bp_goslim_ref %>%
    left_join(goslim_counts, by = c("goslim_id", "goslim_name")) %>%
    mutate(
      count = replace_na(count, 0),
      percentage = (count / total_genes) * 100
    ) %>%
    arrange(desc(percentage)) %>%
    mutate(goslim_name = reorder(goslim_name, percentage))
  
  # Count how many terms have representation
  terms_with_data <- sum(goslim_all_percentages$count > 0)
  
  # Create comprehensive percentage bar chart with all 72 terms
  p <- ggplot(goslim_all_percentages, aes(x = goslim_name, y = percentage)) +
    geom_bar(stat = "identity", fill = component_colors[as.character(comp_num)], 
             alpha = 0.8) +
    geom_text(data = subset(goslim_all_percentages, percentage > 0),
              aes(label = sprintf("%.1f%%", percentage)), hjust = -0.2, size = 2.5) +
    coord_flip() +
    labs(
      title = paste0("Component ", comp_num, ": All 72 GOslim Biological Process Terms (%)"),
      subtitle = paste0(terms_with_data, " of 72 BP GOslim terms found in top 500 genes (genes can have multiple annotations)"),
      x = "GOslim Biological Process Term",
      y = "Percentage of Ortholog Groups (%)"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      plot.subtitle = element_text(size = 11, color = "gray30"),
      axis.text.y = element_text(size = 8),
      axis.text.x = element_text(size = 10),
      axis.title = element_text(size = 11, face = "bold"),
      panel.grid.major.y = element_blank(),
      panel.grid.minor = element_blank()
    ) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.15)))
  
  # Save plot with fixed height for 72 terms
  plot_file <- file.path(goslim_output_dir, 
                         paste0("component_", comp_num, "_goslim_BP_barchart_all_percentage.png"))
  ggsave(plot_file, p, width = 12, height = 22, dpi = 300)
  cat("Saved:", basename(plot_file), "(", terms_with_data, "of 72 BP GOslim terms represented)\n")
}

cat("\nAll comprehensive percentage bar charts saved to:", goslim_output_dir, "\n\n")
```

## Summary statistics

```{r goslim-summary, eval=TRUE}
cat(strrep("=", 60), "\n")
cat("GOSLIM ANNOTATION SUMMARY (Biological Process only)\n")
cat(strrep("=", 60), "\n\n")

summary_data <- list()

for (comp_num in components_to_analyze) {
  comp_col <- paste0("Component_", comp_num)
  
  if (!comp_col %in% names(annotated_components)) {
    next
  }
  
  top_genes <- annotated_components[[comp_col]]$top_genes
  goslim_counts <- annotated_components[[comp_col]]$goslim_counts
  
  genes_with_goslim <- sum(!is.na(top_genes$goslim_ids) & top_genes$goslim_ids != "")
  
  if (nrow(goslim_counts) > 0) {
    top_term <- goslim_counts$goslim_name[1]
    top_term_id <- goslim_counts$goslim_id[1]
    top_count <- goslim_counts$count[1]
  } else {
    top_term <- NA
    top_term_id <- NA
    top_count <- NA
  }
  
  summary_data[[comp_col]] <- data.frame(
    Component = comp_num,
    Total_Genes = nrow(top_genes),
    Genes_With_BP_GOslim = length(unique(top_genes$group_id[!is.na(top_genes$goslim_ids) & top_genes$goslim_ids != ""])),
    Percent_Annotated = round(genes_with_goslim / nrow(top_genes) * 100, 1),
    Unique_BP_GOslim_Terms = nrow(goslim_counts),
    Top_BP_GOslim_Term = top_term,
    Top_BP_GOslim_ID = top_term_id,
    Top_BP_GOslim_Count = top_count,
    stringsAsFactors = FALSE
  )
}

summary_df <- bind_rows(summary_data)

print(summary_df)

# Save summary table
summary_file <- file.path(goslim_output_dir, "goslim_BP_annotation_summary.csv")
write.csv(summary_df, summary_file, row.names = FALSE)
cat("\nSaved summary table:", basename(summary_file), "\n")

cat("\n", strrep("=", 60), "\n")
cat("GOSLIM ANNOTATION ANALYSIS COMPLETE\n")
cat("Output directory:", goslim_output_dir, "\n")
cat(strrep("=", 60), "\n\n")

cat("Output files created:\n")
cat("  - component_[5,6,7,9]_top500_genes_goslim_BP_annotated.csv\n")
cat("  - component_[5,6,7,9]_goslim_BP_counts.csv\n")
cat("  - component_[5,6,7,9]_goslim_BP_percentages.csv\n")
cat("  - component_[5,6,7,9]_goslim_BP_barchart.png (counts)\n")
cat("  - component_[5,6,7,9]_goslim_BP_barchart_percentage.png (percentages)\n")
cat("  - goslim_BP_annotation_summary.csv\n")
cat("\nNote: Only Biological Process GOslim terms are included in this analysis.\n\n")
```

# SYSTEM INFO

```{r system-info, eval=TRUE}
sessionInfo()
```

# REFERENCES
