---
title: "22.20-Apul-multiomics-MOFA2"
author: "Sam White"
date: "2025-09-15"
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 script performs multi-omics data integration and analysis for *A. pulchra* using the MOFA2 framework ([MOFA2 documentation](https://biofam.github.io/MOFA2/)) .

**Software:**
- MOFA2: Multi-Omics Factor Analysis v2, a statistical framework for comprehensive integration of multi-modal data.

**Input files:**

Data files downloaded from [*A.pulchra* Expression Count GitHub wiki](https://github.com/urol-e5/timeseries_molecular/wiki/02%E2%80%90Expression-Count-Matrices#a-pulchra).

-   Genes: [apul-gene_count_matrix.csv](https://gannet.fish.washington.edu/gitrepos/urol-e5/timeseries_molecular/D-Apul/output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv)
-   Transcripts: [apul-transcript_count_matrix.csv](https://gannet.fish.washington.edu/gitrepos/urol-e5/timeseries_molecular/D-Apul/output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-transcript_count_matrix.csv)
-   lncRNA: [lncRNA_counts.clean.filtered.txt](https://gannet.fish.washington.edu/v1_web/owlshell/bu-github/timeseries_molecular/D-Apul/output/31.5-Apul-lncRNA-discovery/lncRNA_counts.clean.filtered.txt)
-   miRNA: [Apul_miRNA_counts_formatted_updated.txt](../data/Apul_miRNA_counts_formatted_updated.txt)
-   Methylated CpGs: [merged-WGBS-CpG-counts_filtered.csv](https://gannet.fish.washington.edu/metacarcinus/E5/20250903_meth_Apul/merged-WGBS-CpG-counts_filtered.csv)

**Output files:**

NOTE: Due to the large file sizes of the output files, they are not uploaded to GitHub.

They are available here: https://gannet.fish.washington.edu/gitrepos/urol-e5/timeseries_molecular/D-Apul/output/22.20-Apul-multiomics-MOFA2/

Creating these files (MOFA2 model HDF5 and RDS files) requires running the MOFA2 training steps in this script, which can take several hours depending on the number of factors and computational resources. It will be faster to download these output files from the above link and then run this script.

- Output files are written to: `../output/22.20-Apul-multiomics-MOFA2/`
    - MOFA2 model HDF5 files: `MOFA2_model_factors-<N>.hdf5` (where `<N>` is the number of factors, e.g. 10-15)
    - MOFA2 model RDS files: `MOFA2_model_factors-<N>.rds` (where `<N>` is the number of factors, e.g. 10-15)

**Script workflow:**
- Loads and harmonizes multi-omics input data (gene, transcript, lncRNA, miRNA, methylation counts)
- Preprocesses data (removes zero-variance features, harmonizes sample columns)
- Trains MOFA2 models for a range of factor numbers (10-15) or loads pre-trained models from disk
- Performs downstream analysis and visualization of MOFA2 results (variance explained, factor characterization, feature weights)

# Setup

## Libraries

```{r setup, include=FALSE}
library(knitr)
library(data.table)
library(MOFA2)
library(ggplot2)
library(tidyverse)
library(ggpubr)
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = FALSE,        # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  comment = ""         # Prevents appending '##' to beginning of lines in code output
)
```

## Variables

```{r variables, eval=TRUE}
# Output directory
output_dir <- "../output/22.20-Apul-multiomics-MOFA2"

# Model seed
model_seed <- 42

# INPUT FILES
data_files <- list(
  gene_counts = "../data/apul-gene_count_matrix.csv",
  transcript_counts = "../data/apul-transcript_count_matrix.csv",
  lncRNA_counts = "../data/lncRNA_counts.clean.filtered.txt",
  miRNA_counts = "../data/Apul_miRNA_counts_formatted_updated.txt",
  methylation_counts = "../data/merged-WGBS-CpG-counts_filtered_n20.csv"
)

# Covariate: Time point
# Will be extracted from sample names after harmonization
```

# Data Structure Overview

```{r data-structure, eval=TRUE}
# Show structure and preview of each data file in the data folder using tidyverse/dplyr

for (name in names(data_files)) {
  cat("\n\n===== ", name, " =====\n")
  file <- data_files[[name]]
  dat <- fread(file)
  print(glimpse(dat))
  print(dat %>% slice_head(n = 3))
}
```

# DATA PREP FOR MOFA2

## View data structures

Primarily converts data to tibbles for easier viewing/compatability with tidyverse.

```{r view-data-structures, eval=TRUE}
# Prepare data for MOFA2: ensure features are rows, samples are columns using dplyr/tidyverse
lncRNA_raw <- fread(data_files$lncRNA_counts) %>% as_tibble()
lncRNA_counts <- lncRNA_raw %>% select(-(1:6))
lncRNA_counts <- lncRNA_counts %>% mutate(Geneid = lncRNA_raw$Geneid) %>% relocate(Geneid)

gene_counts <- fread(data_files$gene_counts) %>% as_tibble()

transcript_counts <- fread(data_files$transcript_counts) %>% as_tibble()

miRNA_counts <- fread(data_files$miRNA_counts) %>% as_tibble()

methylation_counts <- fread(data_files$methylation_counts) %>% as_tibble()

# Show structure of processed data
cat("\n\n===== lncRNA_counts (processed) =====\n")
print(glimpse(lncRNA_counts))
print(lncRNA_counts %>% slice_head(n = 3))

cat("\n\n===== gene_counts (original) =====\n")
print(glimpse(gene_counts))
print(gene_counts %>% slice_head(n = 3))

cat("\n\n===== transcript_counts (original) =====\n")
print(glimpse(transcript_counts))
print(transcript_counts %>% slice_head(n = 3))

cat("\n\n===== miRNA_counts (original) =====\n")
print(glimpse(miRNA_counts))
print(miRNA_counts %>% slice_head(n = 3))

cat("\n\n===== methylation_counts (original) =====\n")
print(glimpse(methylation_counts))
print(methylation_counts %>% slice_head(n = 3))
```

## Compare Column Names Across Input Data

```{r compare-column-names, eval=TRUE}
# List of all data frames
input_dfs <- list(
  gene_counts = gene_counts,
  transcript_counts = transcript_counts,
  lncRNA_counts = lncRNA_counts,
  miRNA_counts = miRNA_counts,
  methylation_counts = methylation_counts
)

# Get column names for each (excluding feature/ID column)
colnames_list <- lapply(input_dfs, function(df) colnames(df)[-1])

# Use the first as reference
ref_names <- colnames_list[[1]]
ref_file <- names(colnames_list)[1]

for (i in seq_along(colnames_list)) {
  this_names <- colnames_list[[i]]
  this_file <- names(colnames_list)[i]
  if (!identical(ref_names, this_names)) {
    cat("\nFile with different column names:", this_file, "\n")
    cat("Columns in", this_file, "not in", ref_file, ":\n")
    print(setdiff(this_names, ref_names))
    cat("Columns in", ref_file, "not in", this_file, ":\n")
    print(setdiff(ref_names, this_names))
  }
}
cat("\nColumn name comparison complete.\n")
```

## Harmonize Columns Across Data Frames

```{r harmonize-columns, eval=TRUE}
# Use the intersection of all sample columns (excluding feature/ID column)
colnames_list <- lapply(list(
  gene_counts = gene_counts,
  transcript_counts = transcript_counts,
  lncRNA_counts = lncRNA_counts,
  miRNA_counts = miRNA_counts,
  methylation_counts = methylation_counts
), function(df) colnames(df)[-1])

common_cols <- Reduce(intersect, colnames_list)

# Function to remove extra columns and print what is removed
remove_extra_cols <- function(df, df_name) {
  orig_cols <- colnames(df)[-1]
  extra <- setdiff(orig_cols, common_cols)
  if (length(extra) > 0) {
    cat("\nRemoving columns from", df_name, ":\n")
    print(extra)
  }
  # Keep only the ID column and common columns
  df <- df[, c(colnames(df)[1], common_cols)]
  return(df)
}

gene_counts <- remove_extra_cols(gene_counts, "gene_counts")
transcript_counts <- remove_extra_cols(transcript_counts, "transcript_counts")
lncRNA_counts <- remove_extra_cols(lncRNA_counts, "lncRNA_counts")
miRNA_counts <- remove_extra_cols(miRNA_counts, "miRNA_counts")
methylation_counts <- remove_extra_cols(methylation_counts, "methylation_counts")


# Print final structure
cat("\nFinal structure after harmonizing columns:\n")
cat("\ngene_counts:\n")
print(str(gene_counts))
cat("\ntranscript_counts:\n")
print(str(transcript_counts))
cat("\nlncRNA_counts:\n")
print(str(lncRNA_counts))
cat("\nmiRNA_counts:\n")
print(str(miRNA_counts))
cat("\nmethylation_counts:\n")
print(str(methylation_counts))

# Create sample metadata with time point covariate
sample_names <- colnames(gene_counts)[-1]
time_point <- stringr::str_extract(sample_names, "TP[1-4]")
sample_metadata <- data.frame(
  sample = sample_names,
  time_point = time_point,
  stringsAsFactors = FALSE
)
cat("\nSample metadata with time point covariate:\n")
print(head(sample_metadata))
```

## Remove Zero-Variance Features

MOFA2 model training recommends removing features with no variance.

```{r remove-zero-variance-features, eval=TRUE}
remove_zero_var <- function(df) {
  feature_data <- df[ , -1, drop=FALSE]
  # Compute variance for each row (feature)
  vars <- apply(feature_data, 1, function(x) var(as.numeric(x), na.rm=TRUE))
  zero_var_idx <- which(vars == 0 | is.na(vars))
  keep_idx <- which(vars > 0 & !is.na(vars))
  removed <- df[zero_var_idx, , drop=FALSE]
  filtered <- df[keep_idx, , drop=FALSE]
  list(filtered=filtered, removed=removed)
}

# Apply to each data frame
gene_counts_zero <- remove_zero_var(gene_counts)
transcript_counts_zero <- remove_zero_var(transcript_counts)
lncRNA_counts_zero <- remove_zero_var(lncRNA_counts)
miRNA_counts_zero <- remove_zero_var(miRNA_counts)
methylation_counts_zero <- remove_zero_var(methylation_counts)

# Extract filtered and removed data frames
gene_counts_nzv <- gene_counts_zero$filtered
gene_counts_zero_removed <- gene_counts_zero$removed
transcript_counts_nzv <- transcript_counts_zero$filtered
transcript_counts_zero_removed <- transcript_counts_zero$removed
lncRNA_counts_nzv <- lncRNA_counts_zero$filtered
lncRNA_counts_zero_removed <- lncRNA_counts_zero$removed
miRNA_counts_nzv <- miRNA_counts_zero$filtered
miRNA_counts_zero_removed <- miRNA_counts_zero$removed
methylation_counts_nzv <- methylation_counts_zero$filtered
methylation_counts_zero_removed <- methylation_counts_zero$removed

# Print structure of all created data frames
cat("\nFiltered data frames (zero-variance features removed):\n")
cat("\ngene_counts_nzv:\n")
print(str(gene_counts_nzv))
cat("\ntranscript_counts_nzv:\n")
print(str(transcript_counts_nzv))
cat("\nlncRNA_counts_nzv:\n")
print(str(lncRNA_counts_nzv))
cat("\nmiRNA_counts_nzv:\n")
print(str(miRNA_counts_nzv))
cat("\nmethylation_counts_nzv:\n")
print(str(methylation_counts_nzv))

cat("\nData frames of removed zero-variance features:\n")
cat("\ngene_counts_zero_removed:\n")
print(str(gene_counts_zero_removed))
cat("\ntranscript_counts_zero_removed:\n")
print(str(transcript_counts_zero_removed))
cat("\nlncRNA_counts_zero_removed:\n")
print(str(lncRNA_counts_zero_removed))
cat("\nmiRNA_counts_zero_removed:\n")
print(str(miRNA_counts_zero_removed))
cat("\nmethylation_counts_zero_removed:\n")
print(str(methylation_counts_zero_removed))
```

## Create Matrices After Zero-Variance Removal

```{r create-matrices, eval=TRUE}

# Create matrices from filtered data frames
gene_mat_nzv <- gene_counts_nzv %>% 
  column_to_rownames(var = colnames(gene_counts_nzv)[1]) %>% 
  as.matrix()

transcript_mat_nzv <- transcript_counts_nzv %>% 
  column_to_rownames(var = colnames(transcript_counts_nzv)[1]) %>% 
  as.matrix()

lncRNA_mat_nzv <- lncRNA_counts_nzv %>% 
  column_to_rownames(var = "Geneid") %>% 
  as.matrix()

miRNA_mat_nzv <- miRNA_counts_nzv %>% 
  column_to_rownames(var = "Name") %>% 
  as.matrix()

methylation_mat_nzv <- methylation_counts_nzv %>% 
  column_to_rownames(var = colnames(methylation_counts_nzv)[1]) %>% 
  as.matrix()

# Combine into a named list for MOFA2
data_list_nzv <- list(
  gene = gene_mat_nzv,
  transcript = transcript_mat_nzv,
  lncRNA = lncRNA_mat_nzv,
  miRNA = miRNA_mat_nzv,
  methylation = methylation_mat_nzv
)
# Show structure of the new list
str(data_list_nzv, max.level = 1)
cat("\nList of matrices after zero-variance feature removal ready for MOFA2 object creation.\n")
```

## Extract sample names
```{r get-sample-names, eval=TRUE}
# Get sample names (columns) from the harmonized data
sample_names <- colnames(gene_counts_nzv)[-1]

# Extract time point from sample names using regex
time_point <- stringr::str_extract(sample_names, "TP[1-4]")

# Create metadata data frame
sample_metadata <- data.frame(
  sample = sample_names,
  time_point = time_point,
  stringsAsFactors = FALSE
)

# Preview
head(sample_metadata)
```

# CREATE MOFA2 OBJECT


```{r mofa2-create-object, eval=TRUE}
# Create a MOFA object from the list of matrices
MOFAobject <- create_mofa(data_list_nzv)

# Add sample metadata (with time point) to MOFA object
samples_metadata(MOFAobject) <- sample_metadata

# Show summary of the MOFA object
print(MOFAobject)
cat("\nMOFA object created.\n")
```

## Plot data overview

```{r plot-data-overview, eval=TRUE}
plot_data_overview(MOFAobject)
```

# DEFINE MOFA2 OPTIONS

## Data options

```{r, data-options, eval=TRUE}
data_opts <- get_default_data_options(MOFAobject)
data_opts
```

## Model options

```{r, model-options, eval=TRUE}
model_opts <- get_default_model_options(MOFAobject)
model_opts$num_factors <- 10

model_opts
```

## Training options

```{r training-options, eval=TRUE}
train_opts <- get_default_training_options(MOFAobject)
train_opts$convergence_mode <- "fast"
train_opts$seed <- model_seed

train_opts
```

# TRAIN MOFA2 OBJECT

A range from 10 - 15 factors is recommended by the developers.

## Load saved MOFA objects from disk

If this has been previously run, this will load existing, trained models.

```{r load-mofa-objects, eval=TRUE}
# Create a new list to store loaded MOFA objects
mofa_objects <- list()
for (nf in 10:15) {
  rds_file <- file.path(output_dir, paste0("MOFA2_model_factors-", nf, ".rds"))
  obj_name <- paste0("MOFAobject_factors", nf)
  if (file.exists(rds_file)) {
    mofa_objects[[obj_name]] <- readRDS(rds_file)
  } else {
    warning(paste("File not found:", rds_file))
  }
}
```

## Train MOFA objects for different numbers of factors

If this is the first time running this, this will train the models.

```{r train-mofa-multi-factors, eval=FALSE}
# List to store MOFA objects
mofa_objects <- list()

data_opts <- get_default_data_options(MOFAobject)

for (nf in 10:15) {
  # Set model options for this number of factors
  model_opts_i <- get_default_model_options(MOFAobject)
  model_opts_i$num_factors <- nf
  train_opts_i <- get_default_training_options(MOFAobject)
  train_opts_i$convergence_mode <- "fast"
  train_opts_i$seed <- model_seed

  # Prepare MOFA object (no covariate)
  MOFAobject_i <- prepare_mofa(MOFAobject,
    data_options = data_opts,
    model_options = model_opts_i,
    training_options = train_opts_i
  )

  # Output filenames and object name
  model_hdf5_file_i <- file.path(output_dir, paste0("MOFA2_model_factors-", nf, ".hdf5"))
  model_rds_file_i <- file.path(output_dir, paste0("MOFA2_model_factors-", nf, ".rds"))
  mofa_obj_name <- paste0("MOFAobject_factors", nf)

  # Train and save
  MOFAobject_i <- run_mofa(MOFAobject_i, outfile = model_hdf5_file_i, use_basilisk = TRUE)
  saveRDS(MOFAobject_i, model_rds_file_i)

  # Store in list with descriptive name
  mofa_objects[[mofa_obj_name]] <- MOFAobject_i
  cat("\nTrained and saved MOFA object with", nf, "factors.\n")
}
```

# PLOTS

## Factor correlation

```{r plot-factor-correlation, eval=TRUE}
for (obj_name in names(mofa_objects)) {
  cat("\nFactor correlation for", obj_name, ":\n")
  print(plot_factor_cor(mofa_objects[[obj_name]]))
}
```

## Variance decomposition

```{r plot-variance-decomposition-all, eval=TRUE}
for (obj_name in names(mofa_objects)) {
  cat("\nVariance explained for", obj_name, ":\n")
  print(plot_variance_explained(mofa_objects[[obj_name]], max_r2=0.1))
}
```

Looks like the MOFA object with Factor count of `15` is the most informative, or at least provides some of the best resulotion. It shows a variance range in methylation (i.e. variation percentage isn't pegged at the max/min colors across all factors), as well as some possible effects in:

- Factor 1 showing variation across all views, so possibly important.

- Factor 8 showing co-variation in only methylation and miRNA, suggesting some level of importance.

- Factor 9 showing co-variation in only methylation lncRNA, suggesting some level of importance.

# FACTOR 1 CHARACTERIZATION


## Plot variance explained per view

```{r plot-variance-explained-view, eval=TRUE}
cat("\nVariance explained per view for MOFAobject_factors15:\n")
print(plot_variance_explained(mofa_objects[["MOFAobject_factors15"]], plot_total = T)[[2]])
```

## Print percent variance explained per view as a table

```{r print-variance-explained-table, eval=TRUE}
var_list <- plot_variance_explained(mofa_objects[["MOFAobject_factors15"]], plot_total = TRUE)
# Extract the ggplot object
var_gg <- var_list[[2]]
# Get the data frame from the ggplot object
var_df <- var_gg$data
print(var_df)
```

## Plot Factor 1 values

```{r plot-factor1-factor-values, eval=TRUE}
plot_factor(mofa_objects[["MOFAobject_factors15"]], 
  factors = 1, 
  color_by = "Factor1"
)
```

## Plot Factor 1 weights

### Genes

```{r plot-factor1-gene-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor1-top-gene-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-gene-weights, eval=TRUE}
# Print top 10 gene features (by absolute weight) and their actual weights for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 1)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```


#### Plot heatmap

```{r plot-factor1-heatmap-gene, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 1,
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 genes

```{r plot-factor1-top-gene-factor1, eval=TRUE}
# Programmatically get top gene by absolute weight for factor 1

# Get top 10 genes by absolute weight for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 1)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_gene <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_genes <- names(top_features)
cat("Top 10 genes by absolute weight for factor 1:\n")
print(top_genes)

# Create a separate plot for each top gene
for (gene in top_genes) {
  cat("\nPlotting factor 1 colored by time point:", gene, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 1,
      color_by = gene,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot gene by transcript scatter

```{r plot-factor1-top-gene-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_gene
) + labs(y="Transcript expression")
```

#### Plot gene by lncRNA scatter

```{r plot-factor1-top-gene-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_gene
) + labs(y="lncRNA expression")
```

#### Plot gene by miRNA scatter

```{r plot-factor1-top-gene-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_gene
) + labs(y="miRNA expression")
```

#### Plot gene by methylation scatter

```{r plot-factor1-top-gene-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_gene
) + labs(y="Methylated CpG counts")
```

### Transcripts

```{r plot-factor1-transcripts-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor1-top-transcripts-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-transcripts-weights, eval=TRUE}
# Print top 10 transcripts features (by absolute weight) and their actual weights for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 1)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```


#### Plot heatmap

```{r plot-factor1-heatmap-transcripts, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 1,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 transcripts

```{r plot-factor1-top-transcripts, eval=TRUE}
# Programmatically get top transcripts by absolute weight for factor 1

# Get top 10 transcripts by absolute weight for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 1)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_transcript <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_transcripts <- names(top_features)
cat("Top 10 transcripts by absolute weight for factor 1:\n")
print(top_transcripts)

# Create a separate plot for each top transcript
for (transcript in top_transcripts) {
  cat("\nPlotting factor 1 colored by transcript:", transcript, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 1, 
      color_by = transcript,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot transcript by gene scatter

```{r plot-factor1-top-transcript-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_transcript
) + labs(y="Transcript expression")
```

#### Plot transcript by lncRNA scatter

```{r plot-factor1-top-transcript-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_transcript
) + labs(y="lncRNA expression")
```

#### Plot transcript by miRNA scatter

```{r plot-factor1-top-transcript-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_transcript
) + labs(y="miRNA expression")
```

#### Plot transcript by methylation scatter

```{r plot-factor1-top-transcript-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_transcript
) + labs(y="Methylated CpG counts")
```

### lncRNA

#### Plot weights

```{r plot-factor1-lncRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA", 
  factor = 1, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor1-top-lncRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "lncRNA",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-lncRNA-weights, eval=TRUE}
# Print top 10 lncRNA features (by absolute weight) and their actual weights for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 1)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor1-heatmap-lncRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 1,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top lncRNA

```{r plot-factor1-top-lncRNA, eval=TRUE}
# Programmatically get top lncRNA by absolute weight for factor 1

# Get top 10 lncRNAs by absolute weight for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 1)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_lncRNA <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_lncRNAs <- names(top_features)
cat("Top 10 lncRNAs by absolute weight for factor 1:\n")
print(top_lncRNAs)

# Create a separate plot for each top lncRNA
for (lncRNA in top_lncRNAs) {
  cat("\nPlotting factor 1 colored by lncRNA:", lncRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 1,
      color_by = lncRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot lncRNA by gene scatter

```{r plot-factor1-top-lncRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 1,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot lncRNA by transcript scatter

```{r plot-factor1-top-lncRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 1,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot lncRNA by miRNA scatter

```{r plot-factor1-top-lncRNA-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 1,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```

#### Plot lncRNA by methylation scatter

```{r plot-factor1-top-lncRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 1,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### miRNA

#### Plot weights

```{r plot-factor1-miRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA", 
  factor = 1, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor1-top-miRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "miRNA",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-miRNA-weights, eval=TRUE}
# Print top 10 miRNA features (by absolute weight) and their actual weights for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 1)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor1-heatmap-miRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 1,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top miRNA

```{r plot-factor1-top-miRNA, eval=TRUE}
# Programmatically get top miRNA by absolute weight for factor 1

# Get top 10 miRNAs by absolute weight for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 1)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_miRNA <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_miRNAs <- names(top_features)
cat("Top 10 miRNAs by absolute weight for factor 1:\n")
print(top_miRNAs)

# Create a separate plot for each top miRNA
for (miRNA in top_miRNAs) {
  cat("\nPlotting factor 1 colored by miRNA:", miRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 1,
      color_by = miRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot miRNA by gene scatter

```{r plot-factor1-top-miRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot miRNA by transcript scatter

```{r plot-factor1-top-miRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot miRNA by lncRNA scatter

```{r plot-factor1-top-miRNA-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot miRNA by methylation scatter

```{r plot-factor1-top-miRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### Methylation

#### Plot weights

```{r plot-factor1-methylation-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation", 
  factor = 1, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor1-top-methylation-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "methylation",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-methylation-weights, eval=TRUE}
# Print top 10 methylation features (by absolute weight) and their actual weights for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 1)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor1-heatmap-methylation, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 1,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top methylation

```{r plot-factor1-top-methylation, eval=TRUE}
# Programmatically get top methylation by absolute weight for factor 1

# Get top 10 methylation features by absolute weight for factor 1
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 1)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_methylations <- names(top_features)

top_methylation <- names(sort(abs(weights_vec), decreasing = TRUE))[1]

cat("Top 10 methylation features by absolute weight for factor 1:\n")
print(top_methylations)

# Create a separate plot for each top methylation feature
for (meth in top_methylations) {
  cat("\nPlotting factor 1 colored by methylation feature:", meth, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 1, 
      color_by = meth,
      group_by = "time_point",
      dodge = TRUE
    )
  )
}
```

#### Plot methylation by gene scatter

```{r plot-factor1-top-methylation-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_methylation,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot methylation by transcript scatter

```{r plot-factor1-top-methylation-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_methylation,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot methylation by lncRNA scatter

```{r plot-factor1-top-methylation-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_methylation,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot methylation by miRNA scatter

```{r plot-factor1-top-methylation-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 1,  
  features = 12,
  sign = "positive",
  color_by = top_methylation,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```



# FACTOR 8 CHARACTERIZATION

## Plot Factor 8 values

```{r plot-factor8-factor-values, eval=TRUE}
plot_factor(mofa_objects[["MOFAobject_factors15"]], 
  factors = 8, 
  color_by = "Factor8"
)
```

## Plot Factor 8 weights

### Genes

```{r plot-factor8-gene-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor8-top-gene-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-gene-weights-factor8, eval=TRUE}
# Print top 10 gene features (by absolute weight) and their actual weights for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 8)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor8-heatmap-gene, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 8,
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 genes

```{r plot-factor8-top-gene-factor8, eval=TRUE}
# Programmatically get top gene by absolute weight for factor 8

# Get top 10 genes by absolute weight for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 8)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_gene_f8 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_genes_f8 <- names(top_features)
cat("Top 10 genes by absolute weight for factor 8:\n")
print(top_genes_f8)

# Create a separate plot for each top gene
for (gene in top_genes_f8) {
  cat("\nPlotting factor 8 colored by gene:", gene, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 8,
      color_by = gene,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot gene by transcript scatter

```{r plot-factor8-top-gene-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f8
) + labs(y="Transcript expression")
```

#### Plot gene by lncRNA scatter

```{r plot-factor8-top-gene-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f8
) + labs(y="lncRNA expression")
```

#### Plot gene by miRNA scatter

```{r plot-factor8-top-gene-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f8
) + labs(y="miRNA expression")
```

#### Plot gene by methylation scatter

```{r plot-factor8-top-gene-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f8
) + labs(y="Methylated CpG counts")
```

### Transcripts

```{r plot-factor8-transcripts-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor8-top-transcripts-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-transcripts-weights-factor8, eval=TRUE}
# Print top 10 transcripts features (by absolute weight) and their actual weights for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 8)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor8-heatmap-transcripts, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 8,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 transcripts

```{r plot-factor8-top-transcripts, eval=TRUE}
# Programmatically get top transcripts by absolute weight for factor 8

# Get top 10 transcripts by absolute weight for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 8)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_transcript_f8 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_transcripts_f8 <- names(top_features)
cat("Top 10 transcripts by absolute weight for factor 8:\n")
print(top_transcripts_f8)

# Create a separate plot for each top transcript
for (transcript in top_transcripts_f8) {
  cat("\nPlotting factor 8 colored by transcript:", transcript, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 8, 
      color_by = transcript,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot transcript by gene scatter

```{r plot-factor8-top-transcript-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f8
) + labs(y="Gene expression")
```

#### Plot transcript by lncRNA scatter

```{r plot-factor8-top-transcript-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f8
) + labs(y="lncRNA expression")
```

#### Plot transcript by miRNA scatter

```{r plot-factor8-top-transcript-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f8
) + labs(y="miRNA expression")
```

#### Plot transcript by methylation scatter

```{r plot-factor8-top-transcript-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f8
) + labs(y="Methylated CpG counts")
```

### lncRNA

#### Plot weights

```{r plot-factor8-lncRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA", 
  factor = 8, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor8-top-lncRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "lncRNA",
 factor = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-lncRNA-weights-factor8, eval=TRUE}
# Print top 10 lncRNA features (by absolute weight) and their actual weights for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 8)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor8-heatmap-lncRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 8,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top lncRNA

```{r plot-factor8-top-lncRNA, eval=TRUE}
# Programmatically get top lncRNA by absolute weight for factor 8

# Get top 10 lncRNAs by absolute weight for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 8)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_lncRNA_f8 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_lncRNAs_f8 <- names(top_features)
cat("Top 10 lncRNAs by absolute weight for factor 8:\n")
print(top_lncRNAs_f8)

# Create a separate plot for each top lncRNA
for (lncRNA in top_lncRNAs_f8) {
  cat("\nPlotting factor 8 colored by lncRNA:", lncRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 8,
      color_by = lncRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot lncRNA by gene scatter

```{r plot-factor8-top-lncRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 8,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA_f8,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot lncRNA by transcript scatter

```{r plot-factor8-top-lncRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 8,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA_f8,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot lncRNA by miRNA scatter

```{r plot-factor8-top-lncRNA-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 8,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA_f8,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```

#### Plot lncRNA by methylation scatter

```{r plot-factor8-top-lncRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 8,  
  features = 12,
  sign = "negative",
  color_by = top_lncRNA_f8,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### miRNA

#### Plot weights

```{r plot-factor8-miRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA", 
  factor = 8, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor8-top-miRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "miRNA",
 factor = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-miRNA-weights-factor8, eval=TRUE}
# Print top 10 miRNA features (by absolute weight) and their actual weights for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 8)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor8-heatmap-miRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 8,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top miRNA

```{r plot-factor8-top-miRNA, eval=TRUE}
# Programmatically get top miRNA by absolute weight for factor 8

# Get top 10 miRNAs by absolute weight for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 8)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_miRNA_f8 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_miRNAs_f8 <- names(top_features)
cat("Top 10 miRNAs by absolute weight for factor 8:\n")
print(top_miRNAs_f8)

# Create a separate plot for each top miRNA
for (miRNA in top_miRNAs_f8) {
  cat("\nPlotting factor 8 colored by miRNA:", miRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 8,
      color_by = miRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot miRNA by gene scatter

```{r plot-factor8-top-miRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f8,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot miRNA by transcript scatter

```{r plot-factor8-top-miRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f8,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot miRNA by lncRNA scatter

```{r plot-factor8-top-miRNA-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f8,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot miRNA by methylation scatter

```{r plot-factor8-top-miRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f8,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### Methylation

#### Plot weights

```{r plot-factor8-methylation-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation", 
  factor = 8, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor8-top-methylation-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "methylation",
 factor = 8,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-methylation-weights-factor8, eval=TRUE}
# Print top 10 methylation features (by absolute weight) and their actual weights for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 8)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor8-heatmap-methylation, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 8,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top methylation

```{r plot-factor8-top-methylation, eval=TRUE}
# Programmatically get top methylation by absolute weight for factor 8

# Get top 10 methylation features by absolute weight for factor 8
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 8)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_methylations_f8 <- names(top_features)

top_methylation_f8 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]

cat("Top 10 methylation features by absolute weight for factor 8:\n")
print(top_methylations_f8)

# Create a separate plot for each top methylation feature
for (meth in top_methylations_f8) {
  cat("\nPlotting factor 8 colored by methylation feature:", meth, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 8, 
      color_by = meth,
      group_by = "time_point",
      dodge = TRUE
    )
  )
}
```

#### Plot methylation by gene scatter

```{r plot-factor8-top-methylation-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f8,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot methylation by transcript scatter

```{r plot-factor8-top-methylation-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f8,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot methylation by lncRNA scatter

```{r plot-factor8-top-methylation-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f8,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot methylation by miRNA scatter

```{r plot-factor8-top-methylation-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 8,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f8,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```

# FACTOR 9 CHARACTERIZATION


## Plot Factor 9 values

```{r plot-factor9-factor-values, eval=TRUE}
plot_factor(mofa_objects[["MOFAobject_factors15"]], 
  factors = 9, 
  color_by = "Factor9"
)
```

## Plot Factor 9 weights

### Genes

```{r plot-factor9-gene-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor9-top-gene-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "gene",
 factors = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-gene-weights-factor9, eval=TRUE}
# Print top 10 gene features (by absolute weight) and their actual weights for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 9)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor9-heatmap-gene, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 9,
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 genes

```{r plot-factor9-top-gene-factor9, eval=TRUE}
# Programmatically get top gene by absolute weight for factor 9

# Get top 10 genes by absolute weight for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "gene", factor = 9)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_gene_f9 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_genes_f9 <- names(top_features)
cat("Top 10 genes by absolute weight for factor 9:\n")
print(top_genes_f9)

# Create a separate plot for each top gene
for (gene in top_genes_f9) {
  cat("\nPlotting factor 9 colored by gene:", gene, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 9,
      color_by = gene,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot gene by transcript scatter

```{r plot-factor9-top-gene-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f9
) + labs(y="Transcript expression")
```

#### Plot gene by lncRNA scatter

```{r plot-factor9-top-gene-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f9
) + labs(y="lncRNA expression")
```

#### Plot gene by miRNA scatter

```{r plot-factor9-top-gene-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f9
) + labs(y="miRNA expression")
```

#### Plot gene by methylation scatter

```{r plot-factor9-top-gene-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_gene_f9
) + labs(y="Methylated CpG counts")
```

### Transcripts

```{r plot-factor9-transcripts-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Plot top weights

```{r plot-factor9-top-transcripts-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "transcript",
 factor = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-transcripts-weights-factor9, eval=TRUE}
# Print top 10 transcripts features (by absolute weight) and their actual weights for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 9)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor9-heatmap-transcripts, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 9,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top 10 transcripts

```{r plot-factor9-top-transcripts, eval=TRUE}
# Programmatically get top transcripts by absolute weight for factor 9

# Get top 10 transcripts by absolute weight for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "transcript", factor = 9)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_transcript_f9 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_transcripts_f9 <- names(top_features)
cat("Top 10 transcripts by absolute weight for factor 9:\n")
print(top_transcripts_f9)

# Create a separate plot for each top transcript
for (transcript in top_transcripts_f9) {
  cat("\nPlotting factor 9 colored by transcript:", transcript, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 9, 
      color_by = transcript,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot transcript by gene scatter

```{r plot-factor9-top-transcript-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f9
) + labs(y="Gene expression")
```

#### Plot transcript by lncRNA scatter

```{r plot-factor9-top-transcript-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f9
) + labs(y="lncRNA expression")
```

#### Plot transcript by miRNA scatter

```{r plot-factor9-top-transcript-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f9
) + labs(y="miRNA expression")
```

#### Plot transcript by methylation scatter

```{r plot-factor9-top-transcript-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_transcript_f9
) + labs(y="Methylated CpG counts")
```

### lncRNA

#### Plot weights

```{r plot-factor9-lncRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA", 
  factor = 9, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor9-top-lncRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "lncRNA",
 factor = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-lncRNA-weights-factor9, eval=TRUE}
# Print top 10 lncRNA features (by absolute weight) and their actual weights for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 9)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor9-heatmap-lncRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 9,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top lncRNA

```{r plot-factor9-top-lncRNA, eval=TRUE}
# Programmatically get top lncRNA by absolute weight for factor 9

# Get top 10 lncRNAs by absolute weight for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "lncRNA", factor = 9)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_lncRNA_f9 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_lncRNAs_f9 <- names(top_features)
cat("Top 10 lncRNAs by absolute weight for factor 9:\n")
print(top_lncRNAs_f9)

# Create a separate plot for each top lncRNA
for (lncRNA in top_lncRNAs_f9) {
  cat("\nPlotting factor 9 colored by lncRNA:", lncRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 9,
      color_by = lncRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot lncRNA by gene scatter

```{r plot-factor9-top-lncRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_lncRNA_f9,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot lncRNA by transcript scatter

```{r plot-factor9-top-lncRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_lncRNA_f9,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot lncRNA by miRNA scatter

```{r plot-factor9-top-lncRNA-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_lncRNA_f9,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```

#### Plot lncRNA by methylation scatter

```{r plot-factor9-top-lncRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_lncRNA_f9,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### miRNA

#### Plot weights

```{r plot-factor9-miRNA-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA", 
  factor = 9, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor9-top-miRNA-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "miRNA",
 factor = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-miRNA-weights-factor9, eval=TRUE}
# Print top 10 miRNA features (by absolute weight) and their actual weights for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 9)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor9-heatmap-miRNA, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 9,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top miRNA

```{r plot-factor9-top-miRNA, eval=TRUE}
# Programmatically get top miRNA by absolute weight for factor 9

# Get top 10 miRNAs by absolute weight for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "miRNA", factor = 9)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_miRNA_f9 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]
top_miRNAs_f9 <- names(top_features)
cat("Top 10 miRNAs by absolute weight for factor 9:\n")
print(top_miRNAs_f9)

# Create a separate plot for each top miRNA
for (miRNA in top_miRNAs_f9) {
  cat("\nPlotting factor 9 colored by miRNA:", miRNA, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 9,
      color_by = miRNA,
      group_by = "time_point",
      add_violin = TRUE,
      dodge = TRUE
    )
  )
}
```

#### Plot miRNA by gene scatter

```{r plot-factor9-top-miRNA-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f9,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot miRNA by transcript scatter

```{r plot-factor9-top-miRNA-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f9,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot miRNA by lncRNA scatter

```{r plot-factor9-top-miRNA-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f9,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot miRNA by methylation scatter

```{r plot-factor9-top-miRNA-by-methylation, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_miRNA_f9,
  shape_by = "time_point"
) + labs(y="Methylated CpG counts")
```

### Methylation

#### Plot weights

```{r plot-factor9-methylation-weights, eval=TRUE}
plot_weights(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation", 
  factor = 9, 
  nfeatures = 10
)
```

#### Plot top weights

```{r plot-factor9-top-methylation-weights, eval=TRUE}
plot_top_weights(mofa_objects[["MOFAobject_factors15"]],
 view = "methylation",
 factor = 9,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

#### Print top 10 weights

```{r print-top-methylation-weights-factor9, eval=TRUE}
# Print top 10 methylation features (by absolute weight) and their actual weights for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 9)

# Extract the matrix from the list
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_feature_names <- names(top_features)
top_feature_weights <- weights_vec[top_feature_names]
top_df <- data.frame(Feature = top_feature_names, Weight = top_feature_weights)
print(top_df)
```

#### Plot heatmap

```{r plot-factor9-heatmap-methylation, eval=TRUE}
plot_data_heatmap(mofa_objects[["MOFAobject_factors15"]], 
  view = "methylation",
  factor = 9,  
  features = 25,
  denoise = TRUE,
  cluster_rows = FALSE, cluster_cols = TRUE,
  show_rownames = TRUE, show_colnames = TRUE,
  scale = "row"
)
```

#### Plot top methylation

```{r plot-factor9-top-methylation, eval=TRUE}
# Programmatically get top methylation by absolute weight for factor 9

# Get top 10 methylation features by absolute weight for factor 9
weights <- get_weights(mofa_objects[["MOFAobject_factors15"]], view = "methylation", factor = 9)
weights_vec <- weights[[1]][,1]
names(weights_vec) <- rownames(weights[[1]])
top_n <- 10
top_features <- sort(abs(weights_vec), decreasing = TRUE)[1:top_n]
top_methylations_f9 <- names(top_features)

top_methylation_f9 <- names(sort(abs(weights_vec), decreasing = TRUE))[1]

cat("Top 10 methylation features by absolute weight for factor 9:\n")
print(top_methylations_f9)

# Create a separate plot for each top methylation feature
for (meth in top_methylations_f9) {
  cat("\nPlotting factor 9 colored by methylation feature:", meth, "\n")
  print(
    plot_factor(mofa_objects[["MOFAobject_factors15"]], 
      factors = 9, 
      color_by = meth,
      group_by = "time_point",
      dodge = TRUE
    )
  )
}
```

#### Plot methylation by gene scatter

```{r plot-factor9-top-methylation-by-gene, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "gene",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f9,
  shape_by = "time_point"
) + labs(y="Gene expression")
```

#### Plot methylation by transcript scatter

```{r plot-factor9-top-methylation-by-transcript, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "transcript",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f9,
  shape_by = "time_point"
) + labs(y="Transcript expression")
```

#### Plot methylation by lncRNA scatter

```{r plot-factor9-top-methylation-by-lncRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "lncRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f9,
  shape_by = "time_point"
) + labs(y="lncRNA expression")
```

#### Plot methylation by miRNA scatter

```{r plot-factor9-top-methylation-by-miRNA, eval=TRUE}
plot_data_scatter(mofa_objects[["MOFAobject_factors15"]], 
  view = "miRNA",
  factor = 9,  
  features = 12,
  sign = "positive",
  color_by = top_methylation_f9,
  shape_by = "time_point"
) + labs(y="miRNA expression")
```

# SESSION INFO
```{r session-info, eval=TRUE}
sessionInfo()
```

# REFERENCES