--- title: "ampliCan Code Explained" author: "Kornel Labun" date: "`r BiocStyle::doc_date()`" package: "`r pkg_ver('amplican')`" output: BiocStyle::html_document vignette: > %\VignetteIndexEntry{code explained} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- This vignette is a plain-English walkthrough of what the amplican code does, step by step, for readers who do not read code. Each section describes what happens inside the pipeline, which settings matter most, and which behaviours may be surprising. File and function names are mentioned so you can locate the relevant code, but no code is shown. # The big picture `amplican` takes raw sequencing reads (FASTQ files) and a configuration file describing your experiments, and produces a set of HTML reports with plots and tables summarising CRISPR editing outcomes. The entry point is a single function called `amplicanPipeline`, defined in `amplican.R`. Internally it calls a sequence of specialised functions, each living in its own file: 1. **Read and validate** the configuration file (`helpers_warnings.R`). 2. **Read and quality-filter** the FASTQ reads (`helpers_alignment.R`, `helpers_filters.R`). 3. **Assign** each read to an experiment by finding its primers (`helpers_alignment.R`). 4. **Align** each read to its amplicon using a global alignment algorithm (`helpers_alignment.R`). 5. **Detect HDR** events if a donor template is provided (`helpers_alignment.R`). 6. **Extract** individual mutations (deletions, insertions, mismatches) from the alignments (`AlignmentsExperimentSet-class.R`, `helpers_general.R`). 7. **Filter** out artifacts: primer-overlap deletions, primer dimers, and off-target reads (`amplicanFilter.R`, `helpers_filters.R`). 8. **Build consensus** between forward and reverse reads (`amplicanSummarize.R`). 9. **Normalize** against control samples to remove background noise (`amplicanNormalize.R`). 10. **Summarize** editing rates, frameshifts, and cut rates (`amplicanSummarize.R`). 11. **Generate reports** as knitted HTML pages (`amplicanReport.R`). A conceptual map of this flow is shown in the overview vignette. Throughout the pipeline, intermediate results are saved to disk using a temporary-file-and-rename pattern: each file is first written with a `.temp` extension and then atomically renamed to its final name. This means that if the process crashes midway, you will not see half-written files -- only `.temp` files that are cleaned up on the next run. # Step 1: Setup and validation **File:** `amplican.R` (`amplicanPipe` function) **Helper files:** `helpers_directory.R`, `helpers_warnings.R` ## What happens When you call `amplicanPipeline`, the first thing it does is resolve all file paths to absolute paths and check that the output folder is writable (via `checkFileWriteAccess` in `helpers_directory.R`). If it is not, the pipeline stops immediately with an error. Next, it handles the `continue` setting. When `continue` is `TRUE` (the default), the pipeline looks for any leftover `.temp` files from a previous crashed run and deletes them, but otherwise tries to reuse existing results to avoid redundant computation. When `continue` is `FALSE`, the entire results folder is deleted and recreated from scratch. The configuration file is then read and validated by three functions in `helpers_warnings.R`: - **`checkConfigFile`** verifies that there are no duplicate experiment IDs, that every combination of barcode, forward primer, and reverse primer is unique, that no primer contains numeric characters (a common copy-paste error), and that all referenced FASTQ files actually exist and are readable. - **`checkPrimers`** searches for the forward primer and the reverse complement of the reverse primer inside each amplicon sequence, and records their positions. If a primer cannot be found in its amplicon, the pipeline stops. - **`checkTarget`** checks whether each guideRNA sequence can be found inside its amplicon. If it cannot, a warning message is printed, but the pipeline does **not** stop -- it simply marks that guide as not found and continues. This is intentional, because you might be working with guides that have mismatches or unusual architectures. ## Key options and their impact - **`continue`** (default: TRUE): Set to FALSE when you want a completely fresh run and are sure you do not need any previous intermediate files. - **`fastqfiles`** (default: 0.5): Controls how forward and reverse reads are used. A value of 0 means both forward and reverse primers must be found in each read pair. A value of 0.5 (the default) means only one of the two primers needs to match -- this is useful when, for example, reverse reads have been primer-trimmed by the sequencing facility. A value of 1 uses only forward reads; 2 uses only reverse reads. This setting flows through the entire pipeline and affects primer matching, consensus, and plotting. - **`sample`** and **`seed`** (defaults: 0): When `sample` is set above zero, only that many reads are randomly sampled from each FASTQ file, rather than reading the full file. This is useful for quick testing or parameter tuning. The `seed` makes the sampling reproducible. ## Quirks - The guideRNA not being found in the amplicon is a warning, not an error. The analysis proceeds, but downstream plots and cut-site detection may not work correctly for that experiment. - If `continue` is TRUE and a previous run completed alignment but not reporting, the pipeline will skip alignment entirely and jump ahead, which can be confusing if you changed alignment parameters expecting them to take effect. # Step 2: Reading and quality filtering **Files:** `helpers_alignment.R` (`makeAlignment` function), `helpers_filters.R` (`goodBaseQuality`, `goodAvgQuality`, `alphabetQuality`) ## What happens For each barcode (a group of experiments sharing the same sequencing pool), the `makeAlignment` function opens the corresponding FASTQ file(s) using a streaming reader. This means reads are loaded in manageable chunks rather than all at once, which keeps memory usage low even for very large files. Each chunk of reads passes through three quality filters: 1. **`goodBaseQuality`** checks that every single nucleotide in a read has a quality score at or above the `min_quality` threshold. If even one base is below the threshold, the entire read is rejected. 2. **`goodAvgQuality`** checks that the average quality score across all nucleotides in a read meets the `average_quality` threshold. A read might pass the per-base filter but still have an overall low average. 3. **`alphabetQuality`** (only active when `filter_n` is TRUE) rejects any read containing the letter "N", which represents a nucleotide that the sequencer could not confidently identify. The function keeps a running tally of how many reads were rejected for each reason, and these counts end up in the `barcode_reads_filters.csv` output file. ## Key options and their impact - **`average_quality`** (default: 30 in the pipeline): Lowering this value keeps more reads but may include low-quality data that produces spurious mismatches. Raising it makes filtering stricter. The scale depends on your sequencing platform (Phred+33 for Illumina). - **`min_quality`** (default: **0** in `amplicanPipeline`, but **20** in the standalone `amplicanAlign` function): This is the single most important quality parameter. The pipeline default of 0 is very permissive -- it effectively disables per-base quality filtering. If you are seeing a lot of background mismatches in your reports, raising this to 15-20 can make a large difference by removing reads with individual low-quality bases. See the quirks section below for why this default differs. - **`filter_n`** (default: FALSE): Turn this on if your sequencing tends to produce many N calls. Reads with N bases cannot be meaningfully aligned. - **`batch_size`** (default: 10 million in the pipeline): This controls how many reads are loaded into memory at once. It has no effect on results -- it is purely a memory management parameter. Lower it if you are running on a machine with limited RAM. ## Quirks - The default `min_quality` differs between `amplicanPipeline` (0) and the standalone `amplicanAlign` (20). This can be confusing. The pipeline defaults were chosen to be permissive so that fewer reads are lost, but for many datasets, a value of 15-20 for `min_quality` gives cleaner results. - All quality filtering happens **before** alignment. Reads that fail quality checks are never aligned and never appear in any downstream analysis. # Step 3: Collapsing identical reads **File:** `helpers_alignment.R` (`makeAlignment` function) ## What happens After quality filtering, each chunk of surviving reads is examined for duplicates. The function creates a table where each row is a unique combination of forward and reverse read sequences, along with a count of how many times that exact pair appeared. When all chunks have been processed, the per-chunk tables are merged and duplicate pairs across chunks are summed. This step is critical for performance: if the same read sequence appears ten thousand times, it is aligned only once, and the count is carried forward. All downstream statistics (editing rates, frameshift counts, etc.) weight each unique read by this count. ## Key options This step has no adjustable parameters. It is always performed. ## Quirks - "Unique reads" (the number after collapsing) is a useful heterogeneity metric. A very low number of unique reads relative to total reads suggests either no editing (most reads are wild-type) or very uniform editing. A very high number of unique reads can indicate sequencing errors or mosaic editing. - The collapsing treats forward and reverse reads as a pair. Two reads with the same forward sequence but different reverse sequences count as two separate unique reads. # Step 4: Assigning reads to experiments **Files:** `helpers_alignment.R` (`locate_pr_start` and `makeAlignment`), `helpers_general.R` (`comb_along`) ## What happens Each unique read must be assigned to one of the experiments defined in the configuration file. Assignment is done by searching for the experiment's primers within the read. The `locate_pr_start` function searches for each primer inside each read using an overlap alignment (a variant of sequence alignment that allows the primer to be partially present at the beginning of the read). This handles the common situation where reads are slightly truncated at the 5' end and do not start exactly at the primer. The function accepts a certain number of mismatches (controlled by `primer_mismatch`) and requires a minimum overlap of approximately half the primer length. For each experiment (identified by its ID), the function checks whether the forward primer is found in the forward read and/or the reverse primer is found in the reverse read. Which primers must match depends on the `fastqfiles` setting: - `fastqfiles = 0`: both forward and reverse primers must be found. - `fastqfiles = 0.5` (default): either the forward or the reverse primer is sufficient. - `fastqfiles = 1`: only the forward primer is checked (reverse reads not used). - `fastqfiles = 2`: only the reverse primer is checked (forward reads not used). Reads whose primers match a given experiment are assigned to it. Reads that do not match any experiment's primers are labelled "unassigned" and saved to `unassigned_reads.csv`. Their forward and reverse sequences can later be aligned against each other in the barcode report to diagnose what went wrong. ## Key options and their impact - **`primer_mismatch`** (default: **2** in the pipeline, **0** in standalone `amplicanAlign`): This controls how many mismatches are tolerated when matching primers to reads. With the default of 2, the pipeline allows up to 2 positions in the primer to differ from the read, which helps recover reads with sequencing errors in the primer region. Setting this to 0 is stricter and will increase the number of unassigned reads, but eliminates any risk of mis-assigning reads to the wrong experiment due to similar primers. If you see a high unassigned rate, try increasing this value; if you suspect cross-contamination between experiments with similar primers, decrease it. - **`fastqfiles`**: See Step 1. The most common reason to change this is when your reverse reads have been trimmed of primers by the sequencing facility (use 0.5 or 1), or when you only have forward reads (use 1, and leave Reverse_Primer and Reverse_Reads_File empty in the config). ## Quirks - The primer search allows partial primer matches at the start of the read (5' truncation). This means a read that begins partway through the primer can still be assigned. This is handled by the `min_overlap` parameter inside `locate_pr_start`, which defaults to half the primer length. - Primer matching uses the cached result for identical primer sequences across experiments sharing the same barcode, so it is efficient. - The default values for `primer_mismatch` differ between `amplicanPipeline` (2) and `amplicanAlign` (0), for the same reason as `min_quality` -- the pipeline is designed to be permissive by default. # Step 5: Aligning reads to amplicons **File:** `helpers_alignment.R` (`makeAlignment` function) ## What happens Once reads are assigned to an experiment, each read is aligned against the amplicon sequence for that experiment. Before alignment, the read is trimmed so that it starts at the position where the forward primer was found (or, for reverse reads, the reverse primer position). The amplicon is likewise trimmed to the region between the two primers. This ensures that both the read and the amplicon start and end at comparable positions, which improves alignment quality. The alignment itself uses the Needleman-Wunsch global alignment algorithm, as implemented in the `pwalign` (formerly `Biostrings`) package. This algorithm finds the optimal way to line up two sequences by inserting gaps where needed. It uses a scoring matrix that rewards matching bases and penalises mismatches, plus penalties for opening and extending gaps. The alignment results -- the aligned read, the aligned amplicon, and a score representing alignment quality -- are stored in the `AlignmentsExperimentSet` object (described in Step 8). ## Key options and their impact - **`scoring_matrix`** (default: match = +5, mismatch = -4): This matrix determines how the alignment algorithm scores each position. Higher match rewards and more negative mismatch penalties make the alignment stricter about calling mismatches. The default works well for most CRISPR amplicon data. You would rarely change this unless you have a specific reason, such as aligning reads with unusual error profiles. - **`gap_opening`** (default: 25): This is the penalty for starting a new gap in the alignment. A high value (like 25) means the algorithm will strongly prefer to call mismatches rather than open a gap. **This is the single most impactful alignment parameter**: lowering it will cause the aligner to call more deletions and insertions (since gaps become cheaper), while raising it will suppress indels in favour of mismatches. For CRISPR analysis where indels are the primary signal, the default of 25 is a good balance. - **`gap_extension`** (default: 0): This is the additional penalty for each position added to an existing gap. A value of 0 means that a long gap costs the same as a short gap (only the opening penalty matters). This is appropriate for CRISPR data where deletions can be large but should not be penalised more for being longer. ## Quirks - The alignment is performed against the amplicon sub-region between the primers, not the full amplicon. This means that any sequence outside the primer region in the read is ignored. If your amplicon is very long and reads do not span the full length, the alignment covers only the overlapping portion. - Alignment scores are stored with each event and are later used by the consensus algorithm (Step 10) to decide between conflicting forward and reverse reads. # Step 6: HDR detection **File:** `helpers_alignment.R` (`is_hdr` and `is_hdr_strict` functions) ## What happens If you have filled in the **Donor** column in your configuration file (indicating a homology-directed repair template), amplican attempts to identify which reads represent successful HDR events. Two different algorithms are available: **Permissive mode** (`is_hdr`, used when `donor_strict = FALSE`): The donor template is first aligned against the amplicon to identify the positions where they differ. Each read is then aligned against the donor. A read is classified as HDR if two conditions are met: (1) it aligns to the donor at least as well as it aligns to the amplicon, and (2) the number of small events (single-base mismatches or indels) overlapping the donor-specific positions does not exceed the `donor_mismatch` threshold. This mode is fast and tolerant of sequencing noise. **Strict mode** (`is_hdr_strict`, used when `donor_strict = TRUE`): Instead of comparing alignment scores, the strict algorithm checks whether every single event that distinguishes the donor from the amplicon is present in the read, matched exactly by position, type, and sequence. A read must contain all donor-specific events to be counted as HDR. Additional events elsewhere in the read (noise) do not disqualify it, unless there are more than `donor_mismatch` extra events within the cut-site window. ## Key options and their impact - **`donor_strict`** (default: FALSE): Switch to TRUE if you need very precise HDR calling. The strict algorithm is more conservative -- it will not count a read as HDR unless it carries every expected donor change. This reduces false positives but may miss reads with partial HDR or sequencing errors in the donor region. The strict algorithm is also significantly slower because it requires event-by-event matching. - **`donor_mismatch`** (default: 3): In permissive mode, this is the maximum number of single-base discrepancies allowed within the donor region. Lower values are stricter; setting to 0 requires a perfect match to the donor in that region (not recommended due to typical sequencing error rates of ~0.1-1%). In strict mode with a finite value, this controls how many extra consensus events (beyond the required donor events) are tolerated within the cut-site window before a read is rejected. - **`cut_buffer`** (default: 5): In strict mode, this defines the window around the cut site within which extra events are counted toward the `donor_mismatch` threshold. ## Quirks - In permissive mode, `donor_mismatch` only counts events that fall within the donor-vs-amplicon event coordinate window. Mismatches or indels elsewhere in the read are invisible to this threshold and never disqualify a read. This is by design: the goal is to assess whether the donor-specific region looks like HDR, not whether the entire read is perfect. - A read that aligns equally well to both the donor and the amplicon (a perfect tie in scores) is treated as a candidate in permissive mode (the comparison uses "greater than or equal"). - If the donor and amplicon are identical (no differing positions), no HDR calling is possible and the function returns without changes. # Step 7: The AlignmentsExperimentSet container **File:** `AlignmentsExperimentSet-class.R` ## What happens All alignment results are stored in a structured data container called `AlignmentsExperimentSet`. This is an S4 class (a formal R object type with strict validation). It holds: - **Forward read alignments**: for each experiment ID, the aligned forward reads and their alignment scores. - **Reverse read alignments**: the same for reverse reads. - **Read types**: for each read, whether it was classified as HDR (TRUE) or not (FALSE). - **Read counts**: how many times each unique read sequence was observed (remember the collapsing step). - **Unassigned reads**: reads that could not be matched to any experiment, saved for diagnostic purposes. - **Experiment data**: the configuration file, extended with statistics that accumulate during the pipeline. - **Barcode data**: per-barcode statistics such as total reads, filtered reads, and unique reads. The class includes a validation method that enforces consistency: every experiment ID must be unique, the number of forward alignments must match the number of reverse alignments, read counts must not contain zeros, and so on. If any of these conditions are violated, an error is reported. The object supports subsetting (selecting individual experiments by index or name), combining multiple objects together, and a `lookupAlignment` method that prints a single alignment in a human-readable, BLAST-like format. The `writeAlignments` method can output all alignments to a file in FASTA or plain-text format. ## Key options - **`write_alignments_format`** (default: "None"): Set to "fasta", "txt", or both (by passing a vector like `c("fasta", "txt")`) to save human-readable alignments to disk. This is useful for debugging or for sharing alignments with collaborators. The default "None" skips this output to save disk space and time. ## Quirks - When you subset an `AlignmentsExperimentSet` (for example, by selecting one experiment), the barcode-level and unassigned-read data are dropped. This is because there is no reliable way to know which barcode-level statistics belong to which experiment after subsetting. - In chunked processing mode (used by the pipeline for large datasets), results are saved per-barcode as temporary RDS files inside the `temp` subfolder. These are later merged. If you delete these files, the pipeline will need to redo the alignments for those barcodes. - The `extractEvents` method (described next) is noted in the code comments as being slow "despite vectorization." If you have a multi-core machine, pass `use_parallel = TRUE` to speed it up. # Step 8: Extracting events **Files:** `AlignmentsExperimentSet-class.R` (`extractEvents` method), `helpers_general.R` (`getEventInfo`, `getEvents`, `defGR`) ## What happens After alignment, each read is represented as a pair of aligned strings (the read and the amplicon, with gap characters inserted). The `extractEvents` method converts these aligned strings into a table of individual events. For each alignment, the `getEvents` function scans the aligned strings character by character and identifies three types of events: - **Deletions**: bases present in the amplicon but missing from the read (gaps in the read). - **Insertions**: bases present in the read but not in the amplicon (gaps in the amplicon). - **Mismatches**: positions where both the read and amplicon have a base, but they differ. Each event is recorded with its start position, end position, width, the original amplicon base(s), the replacement read base(s), the type (deletion, insertion, or mismatch), a read identifier, the alignment score, and the strand (plus for forward reads, minus for reverse reads). The function also adds "artificial" deletion events at the ends of reads that do not span the full amplicon -- these represent the fact that the read ended early, and they are later used by the filtering step. The `getEventInfo` function is a wrapper that handles the coordinate shifts needed because the alignment may not start at position 1 of the amplicon (the amplicon was trimmed to the primer region). It also handles the asymmetry between forward and reverse strands: for forward reads, it checks whether the read extends to the end of the amplicon; for reverse reads, it checks whether the read starts at the beginning. ## Key options - **`use_parallel`** (default: FALSE): Set to TRUE if you have registered a multi-core back-end (via `BiocParallel::register`). This can significantly speed up event extraction for large datasets, as each experiment is processed independently. ## Quirks - Insertion coordinates are adjusted by a cumulative sum of preceding insertion widths. This is necessary because insertions are not present in the amplicon coordinate system, so their positions must be shifted to avoid overlapping with each other. - Multi-base mismatches (where several consecutive bases differ) are split into individual single-base events. This ensures that each mismatch event has a width of 1, which simplifies downstream processing and plotting. - The strand convention is important: "+" always means the event came from a forward read, and "-" means it came from a reverse read. This is true even for amplicons with `Direction = 1` (where the guide is reverse-complemented), because `flipRanges` (called later) handles the coordinate reversal. # Step 9: Filtering artifacts **Files:** `amplicanFilter.R`, `helpers_filters.R` (`findEOP`, `findPD`, `findLQR`) ## What happens Not all events extracted from alignments represent real biological mutations. Three types of artifacts are filtered out: **Events Overlapping Primers (EOP)** -- detected by `findEOP`: When a read is shorter than the amplicon and ends prematurely, the alignment may show what looks like a large deletion covering the end of the amplicon. But this is not a real deletion -- it is simply that the read ran out of sequence. The `findEOP` function identifies any deletion that touches or extends past the primer regions and marks it for removal. This also catches deletions that start before the forward primer ends. This filter operates in both absolute coordinates (before shifting) and relative coordinates (after shifting), and automatically adjusts based on whether the events have been shifted yet. **PRIMER DIMERS** -- detected by `findPD`: Primer dimers are PCR artifacts where the primers themselves amplify each other without any amplicon template, producing very short products. In sequencing, these appear as reads where the alignment shows a gap covering most of the amplicon. The `findPD` function calculates a cutoff: amplicon length minus the combined primer lengths minus a buffer. Any deletion wider than this cutoff is flagged as a primer dimer. All events from reads containing such a deletion are then removed. **Low-quality / off-target reads** -- detected by `findLQR`: Some reads may align to the amplicon but carry an unusually high number of events and a low alignment score, suggesting they originated from a different genomic locus (off-target) or are otherwise poor quality. The `findLQR` function uses a clustering approach (CLARA, a fast variant of k-medoids) on two dimensions: the normalized number of events per read and the normalized alignment score. It tries dividing reads into 2 or 3 clusters and uses the silhouette criterion to decide which is better. If 3 clusters fit better, the cluster with high event counts and low scores (top-left in a score-vs-events plot) is marked for removal. This filter only activates when there are at least 1000 events for the experiment -- below that threshold, it does nothing. After filtering, the `Reads_Filtered` column in the config summary is calculated as total reads minus PRIMER_DIMER reads minus low-score reads. ## Key options and their impact - **`PRIMER_DIMER`** (default: 30): This is the buffer size added to the primer lengths when computing the primer dimer cutoff. A larger value makes the filter more lenient (fewer reads classified as primer dimers); a smaller value makes it more aggressive. The default of 30 means that a read must have a deletion covering more than `amplicon_length - primer_lengths - 30` bases to be called a primer dimer. If you are seeing too many reads classified as primer dimers, increase this value. - **`event_filter`** (default: TRUE): Set to FALSE to disable the off-target detection (`findLQR`). This might be useful if you have a very unusual editing pattern that is being incorrectly flagged, or if you want to see all reads regardless of quality. In most cases, leave it on. ## Quirks - The `findLQR` off-target filter requires at least 1000 events per experiment to activate. For small datasets or pilot experiments, it silently does nothing. This means that the same dataset analyzed with different amounts of data could produce different filtering behavior at the boundary. - The clustering in `findLQR` uses CLARA with a sample size capped at 1000. This means it is an approximation -- results may vary slightly between runs if the random sampling differs. However, the silhouette-based decision between 2 and 3 clusters tends to be stable. - Events overlapping primers are filtered even if they are insertions or mismatches near the primer boundary, not just deletions. The exact rule differs by event type: for insertions, the start position is checked against the reverse primer position; for other events, the end position is checked. - The EOP filter correctly handles both absolute coordinates (events before shifting to relative positions) and relative coordinates (events after shifting). It detects which coordinate system is in use by checking for negative values. # Step 10: Consensus from both strands **File:** `amplicanSummarize.R` (`amplicanConsensus` function) ## What happens When you use both forward and reverse reads (i.e., `fastqfiles` is 0 or 0.5), each read pair provides two independent observations of the same DNA molecule. Ideally, both reads should agree on what mutations are present. The `amplicanConsensus` function compares events from the forward and reverse reads and decides which events to trust. The consensus logic follows these rules: 1. **Full agreement**: If a forward event and a reverse event have the same position, length, type, and replacement sequence, the event is confirmed. The forward read's version is kept as the representative. 2. **Partial disagreement**: If both reads have events in the same region but they differ in details (e.g., same start but different end), the event from whichever read has the higher alignment score is kept, **but only if both reads overlap the expected cut site**. If neither overlaps the cut site, both events are rejected. 3. **One strand broken**: If one read has an event overlapping a primer (indicating it is truncated or unreliable in that region) but the other read covers that region cleanly, the clean read's event is trusted. 4. **Single-strand-only events** (controlled by `promiscuous_consensus`): When an event appears on only one strand and there is no conflicting information from the other strand, the `promiscuous_consensus` setting determines what happens. When TRUE (the default), these single-strand events are accepted. When FALSE, they are rejected unless the other strand is broken in that region (in which case the intact strand is trusted). The output is a boolean vector (TRUE/FALSE for each event) indicating which events represent the consensus. ## Key options and their impact - **`promiscuous_consensus`** (default: TRUE): This is the key parameter for consensus behavior. When TRUE, indels detected on only one strand are kept. This is more sensitive but may include some artifacts. When FALSE (strict mode), single-strand indels are rejected unless the other strand is broken. Use FALSE when you want high-confidence events only, and TRUE when you do not want to miss any real editing events. The trade-off is sensitivity versus precision. - This step only applies when `fastqfiles <= 0.5` (i.e., both forward and reverse reads are used). When using only forward or only reverse reads, all events are automatically marked as consensus (TRUE), because there is no second strand to compare against. ## Quirks - The consensus rules are also documented with a visual diagram in the overview vignette, which may be easier to follow than the text description. - The `amplicanConsensus` function requires that `amplicanOverlap` has already been run (it needs to know which events overlap the cut site to apply rule 2). In the pipeline, overlap is always computed before consensus. - When both strands are "broken" (both have events overlapping primers in the cut-site region), the events from both strands are removed. This is a conservative choice: if neither strand can be trusted in that region, no event is reported. # Step 11: Cut-site overlap detection **File:** `amplicanSummarize.R` (`amplicanOverlap` function) ## What happens The `amplicanOverlap` function determines which events are likely caused by CRISPR editing, as opposed to being random sequencing errors or natural variants. It does this by checking whether each event overlaps the expected cut site. The expected cut site is specified by you in the configuration file: any UPPER CASE letters in the amplicon sequence are treated as the target region (typically the PAM sequence or the guide RNA binding site). The function detects these UPPER CASE groups using the `upperGroups` helper, which finds runs of capital letters in the amplicon string. Each detected cut-site region is then expanded by `cut_buffer` bases in both directions to create a window. Any event whose position overlaps this window is marked as overlapping the cut site (TRUE); all others are marked as not overlapping (FALSE). If an amplicon has no UPPER CASE letters at all, the function issues a warning and uses the entire amplicon as the cut-site region. This means all events will be marked as overlapping, which effectively disables the spatial filtering. ## Key options and their impact - **`cut_buffer`** (default: 5): This controls how forgiving the spatial filter is. With a buffer of 5, events up to 5 bases outside the uppercase region are still counted as overlapping. Increasing the buffer captures events from slightly imprecise cutting; decreasing it makes the filter stricter. This parameter also flows through to reporting, where it determines the highlighted region in plots. Note that CRISPR cut sites are typically a few bases upstream of the PAM, so the buffer helps account for this biological variability. ## Quirks - If you forget to specify UPPER CASE letters in your amplicon, every event will be treated as overlapping the cut site, and the spatial filtering that separates real edits from background noise will not work. This is one of the most common mistakes. - The cut site can consist of multiple disjoint UPPER CASE regions (for example, if you have multiple guide RNAs targeting the same amplicon). All regions are expanded by the buffer and checked independently. - The overlap check uses "any" overlap: even a single-base overlap between the event and the expanded cut-site window counts. # Step 12: Shifting to relative coordinates **File:** `helpers_general.R` (`amplicanMap`, `flipRanges`, `upperGroups`) ## What happens After filtering, events are re-encoded so that their coordinates are relative to the expected cut site rather than being absolute positions within the amplicon. This is done by the `amplicanMap` function. The "zero point" is defined as the first UPPER CASE letter in the amplicon. All event positions are shifted so that this zero point becomes position 0. Events that were upstream of the cut site get negative coordinates; events downstream get positive coordinates. This relative coordinate system is what makes it possible to create metaplots that overlay data from different amplicons of different lengths. For amplicons with `Direction = 1` (meaning the guide RNA needed to be reverse-complemented to match the amplicon), an additional transformation is applied by `flipRanges`: the strand labels are swapped (+ becomes -, and vice versa), coordinates are mirrored (reversed end-to-end relative to the amplicon length), and the original and replacement base sequences are reverse-complemented. This ensures that all events are reported from the perspective of the forward strand, regardless of the original guide direction. ## Key options This step has no adjustable parameters. ## Quirks - Events from amplicons that have **no UPPER CASE letters** are **filtered out** entirely during this step, with a warning. This means that if you forget to mark the cut site in your amplicon, you will lose all events. Do not use `amplicanMap` if you have no expected cut sites specified, and do not use the metaplot functions (which rely on relative coordinates) either. - The `flipRanges` function modifies the data in place for efficiency. This means that after this step, the strand column no longer reflects which read (forward or reverse) produced the event for Direction = 1 amplicons -- it has been flipped to maintain the forward-strand perspective. # Step 13: Normalization against controls **File:** `amplicanNormalize.R` ## What happens If your configuration file has any rows marked as controls (`Control = TRUE` or `Control = 1`), the `amplicanNormalize` function removes background events from the treatment (non-control) samples. The process works as follows: 1. Events from control samples are grouped by the columns specified in the `normalize` parameter (by default, guideRNA and Group). This means that a control is only used to normalize treatments that share the same guideRNA and Group. 2. For each group, the frequency of each unique event (defined by its position, type, width, original base, and replacement base) is calculated relative to the total number of filtered reads in that control group. 3. Events in the control group whose frequency falls below `min_freq` are treated as sequencing errors and discarded -- they are too rare to be considered real background. 4. For each treatment event, the function checks whether an exact match exists among the surviving control events (same position, type, width, original base, replacement base, guideRNA, and Group). If it does, the treatment event is removed, because it is indistinguishable from background. 5. All control events are kept (they are not removed). This step handles the common situation where the sequenced strain has natural sequence differences from the reference amplicon (SNPs), or where there is a low level of index hopping (cross-contamination between sequencing pools). Events containing the replacement base "N" are always removed before normalization, as they represent noise from poor sequencing. ## Key options and their impact - **`min_freq`** (default: 0.01, i.e., 1%): This is the single most important normalization parameter. It sets the frequency threshold below which control events are considered noise and not used for filtering. Lowering it (e.g., to 0.001) makes normalization more aggressive -- more events will be removed from treatments. Raising it (e.g., to 0.1) makes normalization more conservative. The default of 1% works well for most experiments with good sequencing depth. If you suspect index hopping (where reads from one sample contaminate another), raising this to 0.03-0.15 can help, because index hopping typically introduces low-frequency spurious events. You can inspect the mismatch plots in the control reports to see what the background noise level looks like and choose an appropriate threshold. - **`normalize`** (default: c("guideRNA", "Group")): This specifies which columns to use for matching controls to treatments. The default means that only events from the same guideRNA in the same Group are normalized against each other. You can change this to use fewer columns (e.g., `c("guideRNA")` to normalize by guide only, ignoring Group) or more columns. Pass `NULL` to skip normalization entirely even when controls are present. Pass an empty vector `c()` to normalize globally with no stratification (all controls averaged together for all treatments). - **`cut_buffer`** affects only the strict HDR detection path, not normalization directly. ## Quirks - Normalization only activates when the Control column has at least one TRUE/1 value. If all values are FALSE/0, the function prints a warning and returns the events unchanged. - The frequency calculation in controls uses `Reads_Filtered` (total reads after filtering), not `Reads` (total reads before filtering). This means that primer dimers and low-quality reads are not counted in the denominator. - Events in the treatment that match control events are removed **entirely** -- not just their frequency reduced. This is an all-or-nothing filter: if an event looks like background in the control, it is removed from the treatment, regardless of how frequent it is in the treatment. # Step 14: Summarization **File:** `amplicanSummarize.R` (`amplicanSummarize` function) ## What happens The `amplicanSummarize` function computes per-experiment statistics from the filtered, normalized, and consensus-confirmed events. It adds five new columns to the configuration table: - **HDR**: the number of reads classified as homology-directed repair events (only if a donor was provided). - **Reads_Del**: the number of reads containing at least one deletion that overlaps the expected cut site. - **Reads_In**: the number of reads containing at least one insertion that overlaps the expected cut site. - **Reads_Edited**: the number of reads with any edit (deletion, insertion, or HDR) overlapping the cut site. This is not necessarily the sum of Reads_Del and Reads_In, because a single read can have both a deletion and an insertion -- it is counted once in each category but only once in Reads_Edited. - **Reads_Frameshifted**: the number of reads where the net indel length (total insertions minus total deletions) is not a multiple of 3. A net indel length that is not divisible by 3 shifts the reading frame, which typically disrupts protein function. The calculation proceeds in two passes. First, for each individual read, the function determines whether it has any deletions, insertions, HDR status, or frameshift. Then, these per-read results are aggregated to the experiment level by summing the read counts (remembering that each unique read carries a weight equal to how many times it was observed). Only events that are both consensus (confirmed by forward-reverse agreement, or accepted via promiscuous mode) and overlapping the cut site are included in these counts. This ensures that the statistics reflect high-confidence CRISPR-induced editing. ## Key options This step has no direct parameters, but its results depend on all upstream choices: the quality filters (which reads survive), the alignment parameters (which events are called), the consensus mode (which events are confirmed), the cut-site buffer (which events overlap), and the normalization (which events are removed as background). ## Quirks - `Reads_Edited` can be less than `Reads_Del + Reads_In` because a read with both a deletion and an insertion is counted in both individual categories but only once in the combined category. This is explained in the FAQ vignette. - The frameshift calculation uses the **net** indel length across all events in a read. A read with a 4-base deletion and a 1-base insertion has a net indel of -3, which is divisible by 3 and therefore does **not** count as a frameshift. This is biologically correct: the net effect on the reading frame is what matters. - Experiments with no events at all get zeros for all five columns. # Step 15: Reporting **Files:** `amplicanReport.R`, `helpers_rmd.R` ## What happens The final step generates a set of reports in R Markdown format. R Markdown files are templates that combine text, R code, and plotting instructions. When "knitted" (rendered), they produce polished HTML pages with tables and figures. The `amplicanReport` function creates six types of report: 1. **ID report**: One section per experiment ID. Shows cut rates, deletions, insertions, mismatches, and variants for each individual experiment. This is the most detailed level. 2. **Amplicon report**: Aggregates data by unique amplicon sequence. Useful when multiple IDs share the same amplicon and you want to see the combined editing picture. 3. **Barcode report**: Grouped by barcode (sequencing pool). Shows read-level statistics (total reads, filtered reads, assigned/unassigned) and the top unassigned reads aligned against each other for diagnostics. 4. **Group report**: Grouped by the user-defined Group column. Allows comparison across experimental conditions, technicians, or any other grouping. 5. **Guide report**: Aggregates by unique guideRNA. Shows boxplots of cut rates across all experiments using the same guide, helping assess guide consistency. 6. **Index report**: A summary page with links to all other reports and overall statistics. Each report is generated from a template stored inside the package (in the `rmarkdown/templates` directory). The function injects the paths to the result files and relevant parameters into each template, then optionally knits them to HTML. ## Key options and their impact - **`knit_reports`** (default: TRUE): When TRUE, all reports are automatically rendered to HTML. This is the most time-consuming step because the reports contain high-quality vector figures. When FALSE, the Rmd template files are created but not rendered -- you can then open them in RStudio, customise the text, colours, or layout, and knit them manually. - **`cut_buffer`** (default: 5): Passed through to the report templates and used to highlight the expected cut-site region in plots. - **`top`** (default: 5): Controls how many of the most frequent unassigned reads are shown in the barcode report. Each unassigned read is aligned forward-against-reverse to help diagnose why it was not assigned (e.g., wrong primers, primer dimers, contamination). ## Quirks - The reports are always regenerated from scratch -- even if you only changed a plotting colour, the entire reports folder is deleted and recreated. This is why you might want to set `knit_reports = FALSE`, customise the Rmd files, and knit them yourself. - The Rmd templates are fully editable R code. They demonstrate how to load the result CSV files, filter events, and create each type of plot. You can use them as starting points for your own custom analysis. - Report rendering can take significant time (minutes to tens of minutes) depending on dataset size and the number of experiments. This is due to the high-resolution figures (especially the arch plots for deletions and cuts). # Importing external alignments **File:** `helpers_general.R` (`cigarsToEvents`, `pairToEvents`) ## What happens If you prefer to use an external alignment tool instead of amplican's built-in aligner, two functions allow you to import external alignment results and convert them into the same event table format that the pipeline uses internally. Once imported, the events can be filtered, normalised, summarised, and plotted using all the same downstream functions. **`cigarsToEvents`** accepts alignment results in CIGAR string format -- the standard way that alignment tools like BWA, minimap2, and Bowtie2 report alignments in SAM/BAM files. You provide the CIGAR strings, the alignment start positions, the query (read) sequences, the reference (amplicon) sequences, read IDs, mapping quality scores, sequence names, strands, and read counts. The function parses the CIGAR operations (M, I, D, X, N) and converts them into insertion, deletion, and mismatch events in the GRanges format. This function requires the `GenomicAlignments` package to be installed. **`pairToEvents`** accepts alignment results in the "pair" output format produced by EMBOSS needle (specifically `needleall`). This is a human-readable format that shows aligned sequences with match/mismatch/gap indicators. You provide the path to a pair-format file, and the function parses it into the same event table. It reads the alignment parameters (gap open, gap extend) and scores from the file header. Only the `needleall` "pair" format is supported -- other EMBOSS output formats will produce an error. ## Key options These functions have no tunable parameters beyond the input data. They are designed to faithfully convert external alignment representations into amplican's internal format. ## Quirks - `cigarsToEvents` handles extended CIGAR operations including "X" (mismatch) and "N" (skip/ref_skip). If your aligner does not produce "X" operations (some older aligners use "M" for both matches and mismatches), you will not get mismatch events -- only insertions and deletions. You would need to derive mismatches separately (e.g., from MD tags in the BAM file). - `pairToEvents` assumes that reads are not collapsed (each read has a count of 1). If you have pre-collapsed your reads, you would need to adjust the counts manually after import. - `pairToEvents` explicitly checks that the input file was produced by `needleall` with the `-aformat3 pair` option. Files from `needle` or other programs, or with different output formats, will be rejected with an error. - After importing with either function, you will need to run the downstream steps (filtering, consensus, normalization, summarization) manually, since these import functions bypass the `amplicanPipeline` wrapper. # Plotting and visualization **Files:** `helpers_plots.R`, `ggforce_bezier.R`, `src/bezier.cpp` ## What happens amplican provides a rich set of plotting functions, all built on top of the `ggplot2` graphics library. These functions read the event tables produced by the pipeline and create publication-quality figures. The main plot types are: - **`plot_deletions`**: Draws each deletion as a smooth arch (bezier curve) over the amplicon. The height of the arch represents the deletion's frequency. Arches are coloured by whether the deletion overlaps the expected cut site (pink for overlapping, blue for non-overlapping). This is the most distinctive amplican plot type. - **`plot_cuts`**: Overlays deletion arches from multiple experiments on top of each other, showing the distribution of cut sites across the amplicon. Useful for comparing editing outcomes. - **`plot_insertions`**: Represents each insertion as a triangle over the amplicon, with height proportional to frequency. - **`plot_mismatches`**: Shows mismatches as stacked bars at each position along the amplicon, coloured by the replacement nucleotide (red for A, blue for C, yellow for G, green for T). The top panel shows forward reads; the bottom panel shows reverse reads. - **`plot_variants`**: Shows the most frequent mutated sequences (variants) with their frequency, frameshift status, and codon-level information. Each variant is annotated with the amino acid changes it would cause. - **`plot_heterogeneity`**: Shows the cumulative share of reads accounted for by each unique sequence, sorted from most to least frequent. A steep curve indicates low heterogeneity (a few dominant sequences); a flat curve indicates high heterogeneity (many different sequences at similar frequencies). - **`waffle`**: Creates a grid-of-squares chart (like a pie chart but with discrete squares), used in the index report to show the proportion of good vs. filtered reads. There are also "meta" versions of some plots (`metaplot_deletions`, `metaplot_insertions`, `metaplot_mismatches`) that aggregate data across multiple amplicons by using relative coordinates. These are used in the group, guide, and amplicon reports. The smooth arches in deletion and cut plots are rendered using custom bezier curve mathematics implemented in C++ (in `src/bezier.cpp`) and exposed to R through Rcpp (in `ggforce_bezier.R`). The C++ code computes quadratic and cubic bezier paths -- the same mathematical curves used in vector graphics software. The arch height is set to approximately 4/3 of the deletion frequency, following a geometric formula that makes the arch shape visually proportional. ## Key options - Most plotting functions accept a `cut_buffer` parameter that determines which region is highlighted as the expected cut site. - Plotting functions are designed to work with the event tables as produced by the pipeline. If you modify the tables, ensure that the expected columns (seqnames, start, end, width, type, counts, strand, overlaps, consensus, replacement, originally) are still present. ## Quirks - The plots use a consistent colour scheme: A = red, C = blue, G = yellow, T = green for nucleotides; pink for cut-site-overlapping events, blue for non-overlapping. Ambiguous IUPAC codes (M, R, W, S, Y, K, etc.) and N are shown in grey. - The arch plots can become visually crowded when there are many different deletions. In the reports, the figure height is automatically calculated based on the number of elements to display (via the `plot_height` helper), so the HTML pages may be very tall. - The plotting functions are exported and documented, so you can call them directly on your own data subsets to create custom figures outside the report framework. # Quick reference: parameter summary | Parameter | Default | Where | Effect when increased | Effect when decreased | |---|---|---|---|---| | `average_quality` | 30 | Quality filter | Fewer reads pass (stricter) | More reads pass (looser) | | `min_quality` | 0 | Quality filter | Fewer reads pass (stricter) | More reads pass (looser) | | `filter_n` | FALSE | Quality filter | Reads with N bases removed | No N filtering | | `primer_mismatch` | 2 | Read assignment | More reads assigned (risk misassignment) | Fewer reads assigned (more unassigned) | | `gap_opening` | 25 | Alignment | Fewer gaps, more mismatches | More gaps, fewer mismatches | | `PRIMER_DIMER` | 30 | Artifact filter | Fewer reads called primer dimer | More reads called primer dimer | | `cut_buffer` | 5 | Overlap detection | More events count as cut-site | Fewer events count as cut-site | | `promiscuous_consensus` | TRUE | Consensus | Single-strand indels kept | Single-strand indels rejected | | `min_freq` | 0.01 | Normalization | More control events used (more aggressive removal) | Fewer control events used (gentler removal) | | `donor_mismatch` | 3 | HDR detection | More reads called HDR (looser) | Fewer reads called HDR (stricter) | | `event_filter` | TRUE | Artifact filter | Off-target detection active | Off-target detection disabled | # Quick reference: pipeline default differences The standalone `amplicanAlign` function has different defaults from `amplicanPipeline` for several parameters. This is by design: the pipeline is tuned for convenience and permissiveness, while the standalone function gives you more control. | Parameter | `amplicanPipeline` | `amplicanAlign` | |---|---|---| | `min_quality` | 0 | 20 | | `primer_mismatch` | 2 | 0 | | `batch_size` | 10,000,000 | 1,000,000 | If you call `amplicanAlign` directly and then feed its results into the rest of the pipeline manually, be aware that these stricter defaults may produce different filtering behavior than the all-in-one `amplicanPipeline`.