## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
    echo = TRUE, warning = FALSE, message = FALSE,
    fig.width = 7, fig.height = 5,
    fig.crop = FALSE   ## avoids magick dependency during R CMD check
)

## ----eval=FALSE---------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("CGRphylo2")

## ----eval=FALSE---------------------------------------------------------------
# 
# if (!require("devtools", quietly = TRUE)) {
#     install.packages("devtools")
# }
# devtools::install_github("amarinderthind/CGRphylo2")

## ----load_packages------------------------------------------------------------
library(CGRphylo2)
library(ape)

## ----load_data----------------------------------------------------------------
library(seqinr)
fasta_path <- system.file("extdata", "example_sequences.fasta",
                           package = "CGRphylo2")
fastafile  <- seqinr::read.fasta(fasta_path, seqtype = "DNA",
                                  as.string = TRUE, set.attributes = FALSE)
N_filter       <- 50
fasta_filtered <- filter_N(fastafile, N_filter)

cat("Sequences retained:", length(fasta_filtered), "\n")
for (nm in names(fasta_filtered)) {
    cat(nm, ": length =", nchar(fasta_filtered[[nm]]), "bp\n")
}

## ----metadata, fig.cap="Sequence lengths (left) and GC content (right) for the bundled SARS-CoV-2 sequences."----
meta <- create_meta(fasta_filtered, N_filter = 50)
print(meta)

par(mfrow = c(1, 2), mar = c(4, 6, 3, 1))

dotchart(meta$length,
    labels = meta$name, xlab = "Sequence length (bp)",
    pch = 21, bg = "green", pt.cex = 1.4, cex = 0.9,
    main = "Sequence Lengths"
)
dotchart(meta$GC_content,
    labels = meta$name, xlab = "GC content (%)",
    pch = 21, bg = "steelblue", pt.cex = 1.4, cex = 0.9,
    main = "GC Content"
)
par(mfrow = c(1, 1))

## ----calculate_frequencies----------------------------------------------------
library(BiocParallel)

k_mer    <- 4   ## use 6 for real genome-scale data (virus) ~30k bp
len_trim <- min(nchar(unlist(fasta_filtered)))

## SerialParam() required in vignettes — R CMD check forbids >2 parallel processes
freq_matrices <- parallelCGR(fasta_filtered,
    k_mer    = k_mer,
    len_trim = len_trim,
    BPPARAM  = SerialParam()
)

cat("Sequences processed:", length(freq_matrices), "\n")
cat("K-mers per matrix  :", nrow(freq_matrices[[1]]),
    paste0("(4^", k_mer, ")"), "\n")

## ----cgr_plots, fig.cap="CGR plots for the first two sequences, showing the k-mer density patterns characteristic of each genome."----
par(mfrow = c(1, 2), mar = c(2, 2, 3, 1))

cgr1 <- cgrplot(1)
plot(cgr1[, 1], cgr1[, 2],
    main = paste("CGR:", names(fasta_filtered)[1]),
    xlab = "", ylab = "", cex = 0.4, pch = 4, col = "steelblue"
)

cgr2 <- cgrplot(2)
plot(cgr2[, 1], cgr2[, 2],
    main = paste("CGR:", names(fasta_filtered)[2]),
    xlab = "", ylab = "", cex = 0.4, pch = 4, col = "tomato"
)
par(mfrow = c(1, 1))

## ----calculate_distances------------------------------------------------------
 
distance_matrix <- calculateDistanceMatrix(freq_matrices,
    distance_type = "Euclidean"
)
print(round(distance_matrix, 4))

## ----heatmap, fig.cap="Pairwise Euclidean distance matrix. Lighter colours indicate smaller distances (more similar sequences)."----
col_pal <- colorRampPalette(c("white", "steelblue", "darkblue"))(50)
n       <- nrow(distance_matrix)

image(seq_len(n), seq_len(n), distance_matrix,
    col  = col_pal, xaxt = "n", yaxt = "n",
    xlab = "", ylab = "",
    main = "Pairwise Distance Heatmap"
)
axis(1, at = seq_len(n), labels = rownames(distance_matrix),
    las = 2, cex.axis = 0.9)
axis(2, at = seq_len(n), labels = colnames(distance_matrix),
    las = 1, cex.axis = 0.9)
for (i in seq_len(n)) {
    for (j in seq_len(n)) {
        text(i, j, round(distance_matrix[i, j], 3), cex = 0.75)
    }
}

## ----build_tree, fig.cap="Neighbour-joining tree built from CGR Euclidean distances for the bundled SARS-CoV-2 sequences."----
nj_tree <- nj(as.dist(distance_matrix))

plot(nj_tree,
    main = "Neighbour-Joining Tree (CGR distances)",
    cex  = 0.9
)

## ----export_results-----------------------------------------------------------
 
mega_file   <- tempfile(fileext = ".meg")
phylip_file <- tempfile(fileext = ".phy")

saveMegaDistance(mega_file, distance_matrix)
savePhylipDistance(phylip_file, distance_matrix, mode = "relaxed")

cat("--- MEGA file (first 6 lines) ---\n")
writeLines(readLines(mega_file, n = 6))

cat("\n--- PHYLIP file (first 6 lines) ---\n")
writeLines(readLines(phylip_file, n = 6))

## ----biostrings_integration---------------------------------------------------
library(Biostrings)
dna <- DNAStringSet(c(
    seq1 = "ATCGATCGATCGATCG",
    seq2 = "GCTAGCTAGCTAGCTA"
))
seqs <- from_DNAStringSet(dna)
cat("Sequences loaded:", length(seqs), "\n")

## ----detailed_load_data-------------------------------------------------------
library(seqinr)
fasta_path     <- system.file("extdata", "example_sequences.fasta",
                               package = "CGRphylo2")
fastafile      <- seqinr::read.fasta(fasta_path, seqtype = "DNA",
                                      as.string = TRUE, set.attributes = FALSE)
N_filter       <- 50
fasta_filtered <- filter_N(fastafile, N_filter)
meta           <- create_meta(fasta_filtered, N_filter)
len_trim       <- min(meta$length)
k_mer          <- 6

## ----distance_methods, fig.cap="NJ trees using three distance metrics."-------
## Recompute frequency matrices for this section
freq_matrices <- parallelCGR(fasta_filtered,
    k_mer    = k_mer,
    len_trim = len_trim,
    BPPARAM  = SerialParam()
)

par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
for (dtype in c("Euclidean", "Manhattan", "S_Euclidean")) {
    dist_d <- calculateDistanceMatrix(freq_matrices, distance_type = dtype)
    tree_d <- nj(as.dist(dist_d))
    plot(tree_d, main = dtype, cex = 0.8)
}
par(mfrow = c(1, 1))

## ----kmer_comparison, fig.cap="NJ trees at k = 2, 3, and 4."------------------
par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
for (k in c(2, 3, 4)) {
    freq_k <- parallelCGR(fasta_filtered, k_mer = k, len_trim = len_trim,
                          BPPARAM = SerialParam())
    dist_k <- calculateDistanceMatrix(freq_k, distance_type = "Euclidean")
    tree_k <- nj(as.dist(dist_k))
    plot(tree_k, main = paste0("k = ", k), cex = 0.8)
}
par(mfrow = c(1, 1))

## ----parallel_processing, eval=TRUE-------------------------------------------
freq_matrices <- parallelCGR(fasta_filtered,
    k_mer    = 6,
    len_trim = len_trim,
    BPPARAM  = SerialParam()
)
cat("Matrices computed:", length(freq_matrices), "\n")

## ----memory_optimization, eval=FALSE------------------------------------------
# batch_size        <- 100
# n_sequences       <- length(fasta_filtered)
# n_batches         <- ceiling(n_sequences / batch_size)
# all_freq_matrices <- list()
# 
# for (i in seq_len(n_batches)) {
#     start_idx         <- (i - 1) * batch_size + 1
#     end_idx           <- min(i * batch_size, n_sequences)
#     batch_seqs        <- fasta_filtered[start_idx:end_idx]
#     batch_freq        <- parallelCGR(batch_seqs, k_mer = 6, len_trim = len_trim,
#                                      BPPARAM = BiocParallel::bpparam())
#     all_freq_matrices <- c(all_freq_matrices, batch_freq)
# }

## ----session_info-------------------------------------------------------------
sessionInfo()