## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)

## -----------------------------------------------------------------------------
# library(amplican)
# 
# config <- "config.csv"
# fastq_folder <- "./fastq"
# results_folder <- "./results"

## -----------------------------------------------------------------------------
# amplicanPipeline(
#   config,
#   fastq_folder,
#   results_folder = "./results_quicktest",
#   sample = 10000,       # only read 10000 reads per FASTQ file
#   seed = 42,            # reproducible sampling
#   knit_reports = TRUE   # knit reports so you can inspect them
# )

## -----------------------------------------------------------------------------
# amplicanPipeline(
#   config,
#   fastq_folder,
#   results_folder,
#   knit_reports = TRUE
# )

## -----------------------------------------------------------------------------
# amplicanPipeline(
#   config,
#   fastq_folder,
#   results_folder = "./results_custom",
#   primer_mismatch = 0,         # strict primer matching
#   fastqfiles = 0.5,            # one primer sufficient (default)
#   cut_buffer = 0,              # no buffer (whole amplicon uppercased)
#   normalize = c("ID"),         # disable cross-sample normalization
#   PRIMER_DIMER = 20,           # stricter primer dimer detection
#   donor_strict = TRUE,         # strict HDR detection
#   donor_mismatch = 0,          # zero tolerance for donor noise
#   use_parallel = TRUE,         # enable multicore
#   knit_reports = TRUE
# )

## -----------------------------------------------------------------------------
# library(BiocParallel)
# register(MulticoreParam(workers = 8), default = TRUE)
# 
# amplicanPipeline(
#   config,
#   fastq_folder,
#   results_folder,
#   use_parallel = TRUE
# )

## -----------------------------------------------------------------------------
# library(data.table)
# cfg <- fread("results/config_summary.csv")
# 
# # Quick overview: editing rates per experiment
# cfg[, .(ID, Reads_Filtered, Reads_Edited,
#         percent_edited = round(Reads_Edited / Reads_Filtered * 100, 1))]

## -----------------------------------------------------------------------------
# bdf <- fread("results/barcode_reads_filters.csv")
# bdf[, .(Barcode, read_count, bad_base_quality, bad_average_quality,
#         bad_alphabet, filtered_read_count, assigned_reads, unassigned_reads)]

## -----------------------------------------------------------------------------
# # Run with primer_mismatch = 0
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_pm0",
#                  primer_mismatch = 0, sample = 10000, seed = 42,
#                  knit_reports = FALSE)
# 
# # Run with primer_mismatch = 2
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_pm2",
#                  primer_mismatch = 2, sample = 10000, seed = 42,
#                  knit_reports = FALSE)
# 
# # Compare
# pm0 <- fread("./results_pm0/barcode_reads_filters.csv")
# pm2 <- fread("./results_pm2/barcode_reads_filters.csv")
# cat("pm=0 assigned:", sum(pm0$assigned_reads), "\n")
# cat("pm=2 assigned:", sum(pm2$assigned_reads), "\n")

## -----------------------------------------------------------------------------
# # NULL skips normalization completely
# amplicanPipeline(config, fastq_folder, results_folder,
#                  normalize = NULL)

## -----------------------------------------------------------------------------
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_with_norm",
#                  normalize = c("guideRNA", "Group"),
#                  sample = 10000, seed = 42, knit_reports = TRUE)
# 
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_no_norm",
#                  normalize = NULL,
#                  sample = 10000, seed = 42, knit_reports = TRUE)

## -----------------------------------------------------------------------------
# with_norm <- fread("./results_with_norm/config_summary.csv")
# no_norm   <- fread("./results_no_norm/config_summary.csv")
# 
# comparison <- merge(
#   with_norm[, .(ID, Edited_norm = Reads_Edited)],
#   no_norm[, .(ID, Edited_no_norm = Reads_Edited)],
#   by = "ID")
# comparison[, diff_pct := round((Edited_norm - Edited_no_norm) /
#                                  Edited_no_norm * 100, 1)]
# comparison

## -----------------------------------------------------------------------------
# aln <- readRDS("results/alignments/AlignmentsExperimentSet.rds")
# lookupAlignment(aln, ID = "ID_1")  # view the most frequent alignment

## -----------------------------------------------------------------------------
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_hdr_permissive",
#                  donor_strict = FALSE,
#                  donor_mismatch = 3,
#                  sample = 10000, seed = 42, knit_reports = TRUE)

## -----------------------------------------------------------------------------
# amplicanPipeline(config, fastq_folder,
#                  results_folder = "./results_hdr_strict",
#                  donor_strict = TRUE,
#                  donor_mismatch = 0,
#                  sample = 10000, seed = 42, knit_reports = TRUE)

## -----------------------------------------------------------------------------
# permissive <- fread("./results_hdr_permissive/config_summary.csv")
# strict     <- fread("./results_hdr_strict/config_summary.csv")
# 
# hdr_comparison <- merge(
#   permissive[, .(ID, HDR_perm = HDR, Reads_Filtered)],
#   strict[, .(ID, HDR_strict = HDR)],
#   by = "ID")
# 
# hdr_comparison[, HDR_perm_pct := round(HDR_perm / Reads_Filtered * 100, 2)]
# hdr_comparison[, HDR_strict_pct := round(HDR_strict / Reads_Filtered * 100, 2)]
# hdr_comparison[, .(ID, HDR_perm_pct, HDR_strict_pct)]

## -----------------------------------------------------------------------------
# cfg <- fread("results/config_summary.csv")
# cfg[, HDR_per := round(HDR * 100 / Reads_Filtered, 2)]
# fwrite(cfg, "results/config_summary.csv")

## -----------------------------------------------------------------------------
# # Examples from a real project:
# # results_pm0_dsT_dm0_cb0/         - pm=0, donor_strict=T, dm=0, cb=0
# # results_pm0_dsT_dm0_cb0_noNorm/  - same, normalization disabled
# # results_pm1_dsT_dm3_cb0/         - pm=1, donor_strict=T, dm=3, cb=0

## -----------------------------------------------------------------------------
# # Load all config_summary.csv files from different runs
# runs <- list(
#   pm0 = fread("./results_pm0/config_summary.csv"),
#   pm2 = fread("./results_pm2/config_summary.csv")
# )
# 
# # Compare editing rates
# comparison <- data.table(
#   ID = runs$pm0$ID,
#   Edited_pm0 = runs$pm0$Reads_Edited,
#   Edited_pm2 = runs$pm2$Reads_Edited,
#   Filtered_pm0 = runs$pm0$Reads_Filtered,
#   Filtered_pm2 = runs$pm2$Reads_Filtered
# )
# comparison

