## ----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