From cca58677d3f566b564e515af4f880c74c0e02712 Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:02:00 -0700 Subject: [PATCH 1/9] Add referencefull_stitcher.md as copy of stitcher.md --- docs/design/referencefull_stitcher.md | 744 ++++++++++++++++++++++++++ 1 file changed, 744 insertions(+) create mode 100644 docs/design/referencefull_stitcher.md diff --git a/docs/design/referencefull_stitcher.md b/docs/design/referencefull_stitcher.md new file mode 100644 index 000000000..52e14e930 --- /dev/null +++ b/docs/design/referencefull_stitcher.md @@ -0,0 +1,744 @@ +--- +title: Contig Stitching in MiCall +--- + +DeNovo assembly does not invariably translate input reads into a +single contiguous sequence akin to a genomic consensus. Typically, +errors in input data lead to fragmented sequences — referred to as +contigs — which furthermore may overlap, thus encoding the same region +of a genome more than once. Assembling a unified consensus sequence +necessitates the systematic arrangement of these contigs while +addressing discrepancies within overlapping regions. That is the +Stitcher's function. + +# Structure + +The Stitcher is a specialized component within the MiCall system. It +is designed to operate as an independent module which processes the +assembled contigs, generally derived from DeNovo assembler outputs, +and produce a singular, coherent sequence. + +## Modular Aspect + +The Stitcher maintains a distinct and autonomous role within +MiCall. Its implementation is fully isolated to the +`contig_stitcher*.py` files within the MiCall's source code. The +stitcher module can be run as a CLI script, separately from the rest +of the pipeline. The following command runs the Stitcher: + +```sh +micall contig_stitcher --help +``` + + + + + + + + + + +## Interaction + +Stitching is initiated either as a pipeline step in MiCall, or as a +command line call given above. In each case: + +**Input:** The Stitcher receives a single input file in CSV +format. This file contains 1 or more contigs that are the outcomes of +the previous assembly step, together with associated reference genome +information. These contigs are essentially segments of DNA +sequences. They can vary significantly in length. + +**Output:** The sole output from the Stitcher is a CSV +file. This file holds the stitched sequences -- longer or fully +continuous sequences that represent the genomic consensus formed by +merging the initial fragmented contigs, and additional metadata, +such as the inferred reference genome's name. + + + +# Operational procedure + +To clarify operations of the Stitcher, the subsequent section +introduces a vocabulary that is necessary for a precise description. + +## Definitions + +- An **input nucleotide** refers to a nucleotide of an initial + assembly contig sequence. +- A **reference nucleotide** refers to a nucleotide of a reference + genome sequence. +- A **non-conflicting nucleotide** is a **reference nucleotide** that + has at most one candidate **input nucleotide**. +- A **non-ambiguous nucleotide** is an **input nucleotide**, which has + a clear positioning with respect to all **input nucleotides** of all + other contigs associated with the same reference genome. In + particular, all **conflicting nucleotides** are **ambiguous + nucleotides** because they do not have a clear positioning with + respect to their competing **conflicting nucleotide**. +- An **overlap** is a continuos segement of **conflicting + nucleotides**. +- **Multidirectional alignment** is a property of a contig such that: + 1. the contig has aligned in multiple parts. + 2. some parts have been aligned to the forward strand, and some to + the reverse strand of the reference genome. +- **Cross-alignment** is a property of a contig such that: + 1. the contig has aligned in multiple parts. + 2. the contig-order of the aligned parts does not agree with the + reference-order of the aligned parts. +- A **non-aligned contig** is a contig that has been assinged a + reference sequence, but did not align to it. +- An **invalid contig** is a contig with **multidirectional + alignment**. +- A **stitched consensus** is a **valid contig** in the output of the + Stitcher. +- The **final output** refers to the contents of the only output CSV + file produced by the Stitcher. + +## Principles + +The reason the Stitcher operates effectively is due to its utilization +of reference genomes as additional source of truth. More precisely, +the Stitcher integrates two sets of data: + +1. Sequences generated by the initial assembly. +2. Sequences of reference genomes to which assembled contigs get aligned. + +We will say that 1. is the assembler's data, and 2. is aligner's. + +The core belief is that a reference genome can be used to enhance the +quality of and resolve conflicts within initial assembly contigs. + +In applying this approach, the Stitcher is guided by the following principles: + +### Principle of Scale-Dependent Credibility + +The reliability of sequence alignments increases as the length of the +aligned segment increases. +Therefore: + +- **Micro Scale**: For shorter segments, assembler's findings are more + reliable, because of expected abundance of small, local mutations + not present in the reference genome. + +- **Macro Scale**: For longer segments, the aligner's interpretations + are prioritized. The exponential decrease in alignment errors with + increased sequence length makes long alignments particularly + trustworthy. + +### Principle of Length Prioritization + +A longer contig typically arises from a greater number of reads +spanning a larger genomic region. While this does not imply more reads +per individual position, it suggests that the initial set of reads has +successfully assembled over a more extensive sequence, reflecting a +broader and more robust dataset. Moreover, aligning a longer sequence +to the reference genome is statistically less probable, compared to a +shorter sequence. This means that a successful alignment of a longer +contig to the reference genome provides further confidence in its +accuracy. + +Therefore in scenarios where multiple contigs cover the same region of +the reference genome, longer contigs are prioritized over shorter +ones. + +### Ambiguity Omission Principle + +To mitigate the potential propagation of uncertainties, any data that +lacks a definitive, unambiguous position within the reference genome +should be entirely excluded. This approach acknowledges that absolute +certainty in complex genomic datasets is often unattainable, and tries +to establish a reasonable default. + +## Regulations + +Guided by the previously outlined principles, +several precise regulations governing the Stitcher can be extracted: + +1. For every reference genome, at most one **stitched consensus** + must result. +2. No **ambiguous, non-conflicting nucleotide** + shall be included into the **final output**. +3. Every **non-conflicting-** and **non-ambiguous-** nucleotide + pertaining to a **valid contig** is required to be included in the + **stitched consensus** for the associated reference genome. +4. The relative positions of **non-conflicting-** and + **non-ambiguous-** nucleotides must be preserved in the **final output**. +5. All nucleotides present in the **final output** must exclusively + originate from the initial assembly data. + +## Setup + +The setup process for the Stitcher ensures that each contig is +properly aligned and prepared for the stitching process. The steps are +as follows: + +1. **Align Contigs**: Align each contig to its corresponding reference + genome to approximate their positions within a global reference + framework, allowing for spatial comparison between different contigs. + +2. **Split Multi-Alignment Contigs**: Split contigs that align to + multiple distinct parts of the reference genome into separate + segments. + +3. **Handle Reverse Complement**: Reverse complement contigs that + align to the reverse strand of the reference genome to ensure all + sequences are oriented in the same direction. + +4. **Sort Contigs**: Arrange the contigs based on their starting + positions along the reference genome. + +5. **Group by Reference**: Group contigs such that all contigs + associated with the same reference genome are processed together. + +These setup steps perform minimal alteration to the original contigs +and are primarily guided by straightforward, logical +considerations. Therefore, they do not require extensive +rationalization compared to the subsequent rules. + +## Rules of operation + +Stitching is an iterative process, governed by the following rules: + +### Rule 1: Merge Non-Overlapping Contigs + +1. **Verify Non-Overlap**: Ensure that the end of the first contig is + less or equal to the start of the second contig according to their + positions on the reference genome. + +2. **Delete adjacent non-aligned parts**: Filter out any non-aligned + nucleotides positioned after the first contig's aligned part and + before the second contig's aligned part. + +3. **Concatenate Sequences**: Directly join the end of the first + contig to the start of the second contig. + +#### Example: + +**Input:** + +![non overlaping example input illustration](stitcher_rule_1_input.svg) + +- Contig 1: Sequence = `GG[ATGCCC]AA`, aligned to Referece X at + position 10, with first two and last two nucleotides not aligned. +- Contig 2: Sequence = `AC[TTAG]TA`, aligned to Referece X at position + 30, with first two and last two nucleotides not aligned. + +**Procedure:** +- Verify that Contig 1 ends before Contig 2 begins. +- Delete non-aligned nucleotides resulting in Contig 1 = `GG[ATGCCC]` and Contig 2 = `[TTAG]TA`. +- Concatenate Contig 1 and Contig 2 to form `GG[ATGCCC][TTAG]TA`. + +**Result:** + +![non overlaping example result illustration](stitcher_rule_1_result.svg) + +- The new sequence, `GG[ATGCCCTTAG]TA`, spans positions 10 to 34 on the reference genome. + +#### Rationale + +There isn't many alternative actions available to us in these circumstances. +This enables us to consider all of them: + +1. **Leaving contigs as separate**: + + Separate contigs would result in multiple consensus outputs for one genome. + Thus it fails to comply with **regulation 1**. + +2. **Omitting the strip step**: + + Note that the adjacent non-aligned nucleotides of the two sequences + are **ambiguous, non-conflicting nucleotides**. Therefore, leaving + them in place violates **regulation 2**. + +3. **Introducing additional modifications**: + + Since given contigs do not overlap, every nucleotide in them is **non-conflicting**. + Additionally, we have stripped all the **ambiguous nucleotides**. + Therefore, all modifications that can be introduced + would either violate **regulation 3**, **regulation 4** or **regulation 5**. + +### Rule 2: Merge Overlapping Contigs + +1. **Verify Overlap**: Check if the ending position of the first + contig is greater than the starting position of the second contig. + +2. **Delete adjacent non-aligned parts**: Filter out any non-aligned + nucleotides positioned after the first contig's aligned part and + before the second contig's aligned part. + +3. **Align Overlapping Regions**: + - Extract the sequences from the overlapping region in both + contigs. + - Use a global alignment method to align these overlapping + sub-sequences. + +4. **Calculate Concordance Scores**: + - Compute concordance scores for each position within the + overlapping region. Importantly, the concordance calculation is + done purely between the aligned overlapping subsequences of the + contigs, with no regard to the reference genome sequence. The + concordance score represents how well the nucleotides from the + two contigs match at each position. + - The score is calculated using a sliding average approach, + emphasizing regions with high sequence agreement. + +5. **Determine Optimal Cut Point**: + - Identify the cut point based on the concordance scores such that + the it lies in the middle of regions with the highest + concordance. + - This means making cuts as far away from disagreeing nucleotides + as possible. + +6. **Segment and Combine**: + - Segment the overlapping sequences at the determined cut point. + - Concatenate the non-overlapping parts of the contigs with the + segmented parts from the overlapping region. + +#### Example + +**Input:** + +![overlaping example input illustration](stitcher_rule_2_input.svg) + +- Contig 1: Sequence = `G[GGCC A--TAC]T T`, aligned to Reference X from positions 10 to 19. +- Contig 2: Sequence = `--CCAC[AAATAC C]GGG`, aligned to Reference X from positions 14 to 20. + +**Procedure:** + +1. **Verify Overlap**: + - Contig 1 ends at position 19, and Contig 2 starts at position 14 + (both on Reference X), resulting in an overlap from positions 14 + to 19. + +2. **Delete adjacent non-aligned parts**: Contig 1 is right-stripped + to become `G[GGCCA--TAC]`, contig B is left-stripped to become + `[AAATACC]GGG`. + +3. **Align Overlapping Regions**: + - The overlaping sequence is `A--TAC` from contig A, and `AAATAC` + from contig B. + - Align them globally to produce the following alignments: `--ATAC` + and `AAATAC` + +4. **Calculate Concordance**: + - Calculate concordance scores for positions 15 to 20, considering + only the overlap between the two aligned sequences. + - Approximate concordance: `[0.1, 0.2, 0.3, 0.8, 0.8, 0.3]`. + +5. **Determine Cut Point**: + - Use the computed concordance scores to identify the cut point. + - In this example, the highest concordance scores are around + positions with the score 0.9, so choose it as the cut point. + + ``` + Aligned sequences: + + A: --ATAC + B: AAATAC + + Concordance: + 0.1 0.2 0.3 0.8 0.8 0.3 + + Based on the concordance, cut between the positions: + A: --AT|AC + B: AAAT|AC + ``` + +6. **Segment and Combine**: + - Cut the sequences at the determined cut points. + - Combine sequence parts: `G[GGCC][--AT][AC][C]GGG`. + +**Result:** + +![overlaping example result illustration](stitcher_rule_2_result.svg) + +- The new sequence `G[GGC--ATACC]GGG` spans positions 10 to 20 on Reference X, + representing the most accurate combined sequence. + +#### Rationale + +This rule is similar to Rule 1, but deals with overlapping +regions. When contigs overlap, there is a need to choose a cut point +due to: + +1. **Aligner Constraints**: The aligner constrains the size of the + overlapping sequence (by the **Principle of Scale-Dependent + Credibility**), making it impossible to keep both versions of the + overlapping region simultaneously. +2. **Small scale adjustments**: Overlaps are usually small enough that + assembler data is the highest quality data we have for the + nucleotide positions within it. Thus interleaving segments from + both contigs would again violate the **Principle of Scale Dependent + Credibility**. + +We base the choice on concordance +scores, which measure the degree of agreement between the overlapping +sequences of the two contigs. We look for the highest concordance +because: + +**Choice of Cut Point**: +- If a cut point is chosen where concordance is lower than the + maximum, it implies that in the neighboring region around the cut + point, either to the left or right, there will almost certainly be + some incorrect nucleotides due to disagreement between the contigs. +- Conversely, if the concordance is high at the chosen cut point, the + neighboring region is similar between the two contigs. The selected + extensions (left of the cut point from the left contig and right of + the cut point from the right contig) are longer than the alternative + from the conflicting contig, ensuring greater trust in these regions + based on their length (by the **Principle of Length Prioritization**). + +While this method of choosing the cut point based on concordance +scores aligns with the Principles, we acknowledge that there might be +other ways to determine the optimal cut point. However, given the +complexity of overlapping regions and the necessity to preserve +relative ordering, this concordance-based approach is the best we have +identified so far. + +### Rule 3: Split Contigs with Gaps Which Are Covered by Other Contigs + +1. **Identify Large Gaps**: + - For each contig, identify regions within its alignment to the + reference genome that lack coverage, i.e., gaps. Both small gaps + resulting from sequencing errors and large gaps are recognized. + - Significant gaps are determined based on a pre-defined + threshold. In the context of HIV genome analysis, a gap size of + greater than 21 nucleotides is considered significant due to + common RNA secondary structure phenomena. + +2. **Verify Coverage by Other Contigs**: + - For each identified significant gap, check if other contigs span + or cover this gap. Specifically, check if other contigs have + aligned reference coordinates that overlap with the coordinates + of the gap. + +3. **Split Contig at Gap Midpoint**: + - If a significant gap is covered by another contig, split the + contig containing the gap into two separate contigs at the + midpoint of the gap. + - Left-trim the new right contig segment and right-trim the new + left contig segment to remove ambiguity from their ends. + +4. **Update Contig List**: + - Replace the original contig with its two new segments in the list + of contigs. + +#### Example + +**Input:** + +![gap example input illustration](stitcher_rule_3_input.svg) + +- Contig 1: Sequence = `AGC[TTAC---------------------GGCACATATCATA]CTA`, + aligned to Reference X from positions 10 to 48. +- Contig 2: Sequence = `G[TGAC-----GGACG-TCGTCG--TACGATCAG]G`, + aligned to Reference X from positions 8 to 40. + +**Procedure:** + +1. **Identify Large Gaps**: + - Contig 1 has a significant gap between positions 14 and 35. + +2. **Verify Coverage by Other Contigs**: + - Contig 2 covers the gap region from positions 8 to 40. + +3. **Split Contig at Gap Midpoint**: + - Split Contig 1 into two parts at the midpoint of the gap (i.e., position 24). + This creates two new contigs: + - Contig 1a: Sequence = `AGC[TTAC----------]`, + aligned to Reference X from positions 10 to 24. + - Contig 1b: Sequence = `[-----------GGCACATATCATA]CTA`, + aligned to Reference X from positions 25 to 48. + - Trim the new segments: + - Contig 1a becomes `AGC[TTAC]`. + - Contig 1b becomes `[GGCACATATCATA]CTA`. + +4. **Update Contig List**: + - Discard the original Contig 1 and add Contig 1a and Contig 1b to + the list of contigs. + +**Result:** + +![gap example result illustration](stitcher_rule_3_result.svg) + +- Modified list of contigs now includes Contig 2, Contig 11, and Contig 12. + +#### Rationale + +The decision to split contigs at large gaps covered by other contigs +is grounded in the **Principle of Scale-Dependent +Credibility**. Assemblers can occasionally join sequence fragments +incorrectly if the end of one segment appears similar to the start of +another. Relying on the aligner's macro-scale credibility helps +identify these erroneous joins. Large gaps within a contig are +suspicious and suggest potential assembler errors, whereas small gaps +are generally due to sequencing errors or micro-scale mutations and do +not warrant splitting. By leveraging the aligner's high reliability on +a macro scale, we can effectively pinpoint these errors. If other +contigs cover large gaps, it confirms the aligner's indication that +the assembly might have joined unrelated segments. Splitting contigs +at the midpoint of significant gaps ensures that only those segments +supported by both the assembler's micro-scale data and the aligner's +macro-scale alignment are included in the final stitched consensus. + +The threshold for considering a gap significant is set at 21 +nucleotides. This value was chosen because it correlates with the +average pitch of the RNA helix, which reflects how reverse +transcription periodic deletions are structured around 21 nucleotides +in HIV sequences. Choosing this cutoff recognizes that deletions of +approximately this length are a common feature due to RNA secondary +structures and should not automatically warrant a split. This way, we +avoid splitting on every small gap, which is expected given the nature +of micro-scale mutations, but effectively identify and act on larger, +suspect gaps indicative of potential assembler errors. + +### Rule 4: Discard Contigs That Are Fully Covered By Other Contigs + +1. **Identify Covered Contigs**: + - For each contig in the input set, calculate its aligned interval on the reference genome. + - Identify intervals (regions) that are completely covered by input contigs. + +2. **Compare Intervals**: + - Assess the intervals of each contig to find any contig that falls entirely within the span of other contig intervals. + These are the contigs that are fully covered by others. + +3. **Discard Fully Covered Contigs**: + - Once identified, remove the covered contigs. + +#### Example + +**Input:** + +![covered example input illustration](stitcher_rule_4_input.svg) + +- Contig 1: Sequence = `A[ATCGA]GCT`, aligned to Reference X from positions 10 to 15. +- Contig 2: Sequence = `C[TAGTTG]A`, aligned to Reference X from positions 14 to 19. +- Contig 3: Sequence = `G[CGTACC]G`, aligned to Reference X from positions 12 to 17. + +**Procedure:** + +1. **Identify Covered Contigs**: + - Calculate the intervals: + - Contig 1: `[10-15]` + - Contig 2: `[14-19]` + - Contig 3: `[12-17]` + +2. **Compare Intervals**: + - Assess intervals and find Contig 3: `[12-17]` is completely within the intervals `[10-15]` of Contig 1 and `[14-19]` of Contig 2. + +3. **Discard Fully Covered Contigs**: + - Remove Contig 3 from the analysis. + +**Result:** + +![covered example result illustration](stitcher_rule_4_result.svg) + +- Unchanged remaining contigs Contig 1 and Contig 3. + +#### Rationale + +The underlying idea for this rule is founded on the two following principles: + +1. **Principle of Length Prioritization**: longer contigs are + inherently more reliable. + +2. **Principle of Scale-Dependent Credibility**: Fully covered contigs + might introduce small-scale inconsistencies that the longer + contig can resolve more credibly, given the enhanced reliability + associated with its length and alignment. + +Moreover, keeping all contigs would violate **Regulation 1**. + +--- + +**Note**: rules apply to contigs that are in the same group. + +# Diagnostics + +The Stitcher includes diagnostic tools to ensure transparency and +correctness throughout the stitching process. Two primary methods are +used for diagnostics: visualizer plots and traditional log +files. These tools help users understand and verify the decisions made +by the Stitcher during the stitching process. + +## The Optional Visualizer Tool + +The visualizer can be enabled through the `--plot` flag when running +the Stitcher executable. Running the Stitcher with this flag will +produce an SVG file that visualizes the stitching process, helping to +confirm and debug the Stitcher's operations. + +To use the visualizer, run the Stitcher with an additional argument +specifying the path to the output plot file. Here's an example of how +to stitch contigs and retrieve a visualizer plot: + +```sh +PYTHONPATH="/path/to/micall/repository" python3 -m micall.core.contig_stitcher "contigs.csv" "stitched_contigs.csv" --plot "visualized.svg" +``` + +**Command Line Arguments:** + +- `contigs.csv`: Input file in CSV format containing assembled + contigs and related information. +- `stitched_contigs.csv`: Output CSV file that will contain the + stitched contigs. +- `--plot visualized.svg`: The optional argument to generate a visual + representation of the stitching process, saved as `visualized.svg`. + +### Understanding the Output + +In practice, a visualizer plot might look something like this: + +![practical visualizer plot](stitcher_practical_plot.svg) + +From such a diagram, you can gain insights into the following aspects +of the stitching process: + +- **Reference genome**: The best matching reference genome for this + group of contigs was determined to be `HIV1-A1-RW-KF716472`. + +- **Dropped Contigs**: Contigs that were dropped due to being fully + covered by other contigs, as per Rule 4. In the example plot: + - Contigs 2, 4, 7, 8, and 6 were dropped. + +- **Split Contigs**: Contigs split at large gaps covered by other + contigs, according to Rule 3. The resulting parts are shown as + individual segments. + - Contig 1 was split around Contig 3, producing segments labeled as + 1.1 and 1.3. + +- **Joined Contigs**: Contigs that were merged due to overlap: + - Contigs 1 and 3, which were joined as per Rule 2, with + **ambiguous, non-conflicting** nucleotides discarded, shown as + segments labeled 1.2 and 3.1. + +- **Unaligned Contigs**: Contigs that failed to align to the reference + genome during the alignment step of the setup. + - Contig 5 failed to align. + +- **Contigs without a Reference**: Contigs for which a reference + genome could not be determined during the reference detection step + of the setup. + - Contigs 9 and 10 failed to determine a reference genome. + +Understanding these basics will help to interpret other scenarios +displayed by the visualizer plot. + +## Traditional Logs + +In addition to visual tools, the Stitcher produces traditional log +files that provide textual details of the stitching process. These +logs are crucial for debugging and understanding the sequence of +operations performed by the Stitcher. The verbosity of logs can be +adjusted using command-line options (`--verbose`, `--debug`, `--quiet`). + +Here is an example of typical log entries: + +```text +DEBUG:micall.core.contig_stitcher:Introduced contig 'contig.00001' (seq = TA...CA) of ref 'HIV1-C-BR-JX140663-seed', group_ref HIV1-A1-RW-KF716472-seed (seq = GA...AC), and length 7719. +DEBUG:micall.core.contig_stitcher:Introduced contig 'contig.00002' (seq = CG...AG) of ref 'HIV1-A1-RW-KF716472-seed', group_ref HIV1-A1-RW-KF716472-seed (seq = GA...AC), and length 1634. +... +DEBUG:micall.core.contig_stitcher:Contig 'contig.00006' produced 1 aligner hits. After connecting them, the number became 1. +DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00006' re-aligned as (5) at 7M...3D@[8,1433]->[7461,8946]. +DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00007' aligned at 76M...3D@[0,732]->[5536,6277]. +DEBUG:micall.core.contig_stitcher:Contig 'contig.00007' produced 1 aligner hits. After connecting them, the number became 1. +DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00007' re-aligned as (6) at 76M...3D@[0,732]->[5536,6277]. +... +DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 3D@[790,789]->[8280,8282]. +DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 19D@[1324,1323]->[8817,8835]. +DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 2D@[1354,1353]->[8866,8867]. +... +DEBUG:micall.core.contig_stitcher:Created contigs (8) at 24M...1I@[14,3864]->[0,4558] and (9) at 708D...92I@[3865,7691]->[4559,9032] by cutting (1) at 24M...1I@[14,7691]->[0,9032] at cut point = 4558.5. +DEBUG:micall.core.contig_stitcher:Doing rstrip of (8) at 24M...1I@[14,3864]->[0,4558] (len 7719) resulted in (10) at 24M...1I@[14,3864]->[0,3850] (len 3865). +DEBUG:micall.core.contig_stitcher:Doing lstrip of (9) at 708D...92I@[3865,7691]->[4559,9032] (len 7719) resulted in (11) at 14M...1I@[0,3734]->[5267,9032] (len 3762). +DEBUG:micall.core.contig_stitcher:Split contig (1) at 24M...1I@[14,7691]->[0,9032] around its gap at [3864, 3863]->[3851, 5266]. Left part: (10) at 24M...1I@[14,3864]->[0,3850], right part: (11) at 14M...1I@[0,3734]->[5267,9032]. +... +DEBUG:micall.core.contig_stitcher:Created a frankenstein (34) at 24M...1I@[14,4185]->[0,4171] (len 4186) from [(26) at 24M...1I@[14,3041]->[0,3027] (len 3042), (28) at 271M2D3M2I395M@[0,670]->[3028,3698] (len 671), (30) at 152M@[0,151]->[3699,3850] (len 152), (31) at 321M@[0,320]->[3851,4171] (len 321)]. +DEBUG:micall.core.plot_contigs:Contig name (26) is displayed as '1.1'. +DEBUG:micall.core.plot_contigs:Contig name (36) is displayed as '1.3'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00002' is displayed as '2'. +DEBUG:micall.core.plot_contigs:Contig name (2) is displayed as '2'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00003' is displayed as '3'. +DEBUG:micall.core.plot_contigs:Contig name (31) is displayed as '3.2'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00004' is displayed as '4'. +DEBUG:micall.core.plot_contigs:Contig name (4) is displayed as '4'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00005' is displayed as '5'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00006' is displayed as '6'. +DEBUG:micall.core.plot_contigs:Contig name (5) is displayed as '6'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00007' is displayed as '7'. +DEBUG:micall.core.plot_contigs:Contig name (6) is displayed as '7'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00008' is displayed as '8'. +DEBUG:micall.core.plot_contigs:Contig name (7) is displayed as '8'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00009' is displayed as '9'. +DEBUG:micall.core.plot_contigs:Contig name 'contig.00010' is displayed as '10'. +``` + +The following points illustrate how these logs can facilitate +understanding the stitching process: + +- **Contig Introduction**: Provides details about the contigs + introduced for stitching. + - `Introduced contig 'contig.00001'...` + +- **Alignment Details**: Shows the alignment results for each contig. + - `Part 0 of contig 'contig.00006' re-aligned as (5) at + 7M...3D@[8,1433]->[7461,8946].` + +- **Gap Handling**: Indicates which gaps were ignored as + insignificant. + - `Ignored insignificant gap of (5), 3D@[790,789]->[8280,8282].` + +- **Splitting and Merging Contigs**: Documents the splitting of + contigs at identified gaps and merging of overlapping segments. + - `Split contig (1) at 24M...1I@[14,7691]->[0,9032]...` + - `Created a frankenstein (34) at 24M...1I@[14,4185]->[0,4171]...` + +- **Visualizer Compatibility**: The visualizer diagrams are produced + exclusively from these logs, ensuring compatibility and consistency + between the logs and visual output. + +# Limitations + +Following limitations stem from the choice of principles and various +assumptions that guide the Stitcher's operation. Understanding them +allows users to better interpret the results and apply post-processing +steps to mitigate potential issues. + +One of the critical challenges is the handling of ambiguous +nucleotides. The Stitcher's **Ambiguity Omission Principle**, which +aims to avoid propagating uncertainties, might lead to the exclusion +of significant sequence data, resulting in the loss of potentially +valuable variations or mutations. + +Moreover, the calculation of concordance in overlapping regions +assumes that local concordance is the best indicator of the correct +sequence. This approach may not fully account for complex genomic +rearrangements or context outside the overlap, potentially +compromising the accuracy of the stitched sequence. + +The predefined threshold for significant gaps, based on specific +assumptions about RNA secondary structures of organisms like HIV, +might not generalize well to other organisms or genomic regions. This +can lead to over-splitting or under-splitting contigs, further +fragmenting the consensus sequence. + +Additionally, The Stitcher’s principle of scale-dependent credibility +might overlook important small-scale variations, such as single +nucleotide polymorphisms (SNPs) or small indels, especially if they +are lost in longer contigs deemed more reliable. + +Another critical limitation arises in the context of pipelines dealing +with proviral sequences. The Stitcher might attempt to "fix" sequences +that are inherently "broken", such as those that are scrambled, +contain long deletions, or exhibit hypermutation. In such cases, the +tool's corrective measures may not be desirable, as they risk +introducing inaccuracies. This limitation makes the Stitcher +unsuitable for certain pipelines where the integrity of such broken +sequences should be preserved without alteration. + +Finally, the handling of multidirectional and cross-alignments may +fall short when addressing complex genomic rearrangements, such as +translocations or inversions, potentially resulting in misalignments +and stitching errors in the consensus sequence. From 6cebb034ffbb34e62595dcc817f1cd12333db92c Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:03:14 -0700 Subject: [PATCH 2/9] Scope referencefull doc and fix stale CLI and typo --- docs/design/referencefull_stitcher.md | 78 ++++++++++++++++----------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/docs/design/referencefull_stitcher.md b/docs/design/referencefull_stitcher.md index 52e14e930..0aa631e8d 100644 --- a/docs/design/referencefull_stitcher.md +++ b/docs/design/referencefull_stitcher.md @@ -1,7 +1,16 @@ --- -title: Contig Stitching in MiCall +title: Referencefull Contig Stitching in MiCall --- +This document describes the **referencefull** contig stitcher +(`micall/utils/referencefull_contig_stitcher.py`, +invoked as `micall contig_stitcher with-references`). +It is one of two stitching algorithms in MiCall; see +[Contig Stitching in MiCall](stitcher.md) for the overview and +[Referenceless stitcher](referenceless_stitcher.md) for the other +algorithm. Unless noted otherwise, "the stitcher" below means the +referencefull stitcher. + DeNovo assembly does not invariably translate input reads into a single contiguous sequence akin to a genomic consensus. Typically, errors in input data lead to fragmented sequences — referred to as @@ -13,21 +22,26 @@ Stitcher's function. # Structure -The Stitcher is a specialized component within the MiCall system. It +The referencefull stitcher is a specialized component within the +MiCall system. It is designed to operate as an independent module which processes the assembled contigs, generally derived from DeNovo assembler outputs, and produce a singular, coherent sequence. ## Modular Aspect -The Stitcher maintains a distinct and autonomous role within -MiCall. Its implementation is fully isolated to the -`contig_stitcher*.py` files within the MiCall's source code. The +The referencefull stitcher maintains a distinct and autonomous role within +MiCall. Its implementation lives primarily in +`micall/utils/referencefull_contig_stitcher.py` (with supporting +types in `micall/utils/contig_stitcher_contigs.py`, events in +`micall/utils/referencefull_contig_stitcher_events.py`, context in +`micall/utils/contig_stitcher_context.py`, and entry point in +`micall/core/contig_stitcher.py`). The stitcher module can be run as a CLI script, separately from the rest -of the pipeline. The following command runs the Stitcher: +of the pipeline. The following command runs the referencefull stitcher: ```sh -micall contig_stitcher --help +micall contig_stitcher with-references --help ``` @@ -44,13 +58,13 @@ micall contig_stitcher --help Stitching is initiated either as a pipeline step in MiCall, or as a command line call given above. In each case: -**Input:** The Stitcher receives a single input file in CSV +**Input:** The referencefull stitcher receives a single input file in CSV format. This file contains 1 or more contigs that are the outcomes of the previous assembly step, together with associated reference genome information. These contigs are essentially segments of DNA sequences. They can vary significantly in length. -**Output:** The sole output from the Stitcher is a CSV +**Output:** The sole output from the referencefull stitcher is a CSV file. This file holds the stitched sequences -- longer or fully continuous sequences that represent the genomic consensus formed by merging the initial fragmented contigs, and additional metadata, @@ -60,7 +74,7 @@ such as the inferred reference genome's name. # Operational procedure -To clarify operations of the Stitcher, the subsequent section +To clarify operations of the referencefull stitcher, the subsequent section introduces a vocabulary that is necessary for a precise description. ## Definitions @@ -92,15 +106,15 @@ introduces a vocabulary that is necessary for a precise description. - An **invalid contig** is a contig with **multidirectional alignment**. - A **stitched consensus** is a **valid contig** in the output of the - Stitcher. + referencefull stitcher. - The **final output** refers to the contents of the only output CSV - file produced by the Stitcher. + file produced by the referencefull stitcher. ## Principles -The reason the Stitcher operates effectively is due to its utilization +The reason the referencefull stitcher operates effectively is due to its utilization of reference genomes as additional source of truth. More precisely, -the Stitcher integrates two sets of data: +the stitcher integrates two sets of data: 1. Sequences generated by the initial assembly. 2. Sequences of reference genomes to which assembled contigs get aligned. @@ -110,7 +124,7 @@ We will say that 1. is the assembler's data, and 2. is aligner's. The core belief is that a reference genome can be used to enhance the quality of and resolve conflicts within initial assembly contigs. -In applying this approach, the Stitcher is guided by the following principles: +In applying this approach, the referencefull stitcher is guided by the following principles: ### Principle of Scale-Dependent Credibility @@ -154,7 +168,7 @@ to establish a reasonable default. ## Regulations Guided by the previously outlined principles, -several precise regulations governing the Stitcher can be extracted: +several precise regulations governing the referencefull stitcher can be extracted: 1. For every reference genome, at most one **stitched consensus** must result. @@ -170,7 +184,7 @@ several precise regulations governing the Stitcher can be extracted: ## Setup -The setup process for the Stitcher ensures that each contig is +The setup process for the referencefull stitcher ensures that each contig is properly aligned and prepared for the stitching process. The steps are as follows: @@ -535,7 +549,7 @@ suspect gaps indicative of potential assembler errors. ![covered example result illustration](stitcher_rule_4_result.svg) -- Unchanged remaining contigs Contig 1 and Contig 3. +- Unchanged remaining contigs Contig 1 and Contig 2. #### Rationale @@ -557,25 +571,25 @@ Moreover, keeping all contigs would violate **Regulation 1**. # Diagnostics -The Stitcher includes diagnostic tools to ensure transparency and +The referencefull stitcher includes diagnostic tools to ensure transparency and correctness throughout the stitching process. Two primary methods are used for diagnostics: visualizer plots and traditional log files. These tools help users understand and verify the decisions made -by the Stitcher during the stitching process. +by the stitcher during the stitching process. ## The Optional Visualizer Tool The visualizer can be enabled through the `--plot` flag when running -the Stitcher executable. Running the Stitcher with this flag will +the stitcher executable. Running the stitcher with this flag will produce an SVG file that visualizes the stitching process, helping to -confirm and debug the Stitcher's operations. +confirm and debug the stitcher's operations. -To use the visualizer, run the Stitcher with an additional argument +To use the visualizer, run the stitcher with an additional argument specifying the path to the output plot file. Here's an example of how to stitch contigs and retrieve a visualizer plot: ```sh -PYTHONPATH="/path/to/micall/repository" python3 -m micall.core.contig_stitcher "contigs.csv" "stitched_contigs.csv" --plot "visualized.svg" +PYTHONPATH="/path/to/micall/repository" python3 -m micall.core.contig_stitcher with-references "contigs.csv" "stitched_contigs.csv" --plot "visualized.svg" ``` **Command Line Arguments:** @@ -702,12 +716,16 @@ understanding the stitching process: # Limitations Following limitations stem from the choice of principles and various -assumptions that guide the Stitcher's operation. Understanding them +assumptions that guide the referencefull stitcher's operation. Understanding them allows users to better interpret the results and apply post-processing -steps to mitigate potential issues. +steps to mitigate potential issues. Where structural fidelity matters +more than contiguity, consider the +[referenceless stitcher](referenceless_stitcher.md), which deliberately +avoids reference-derived ordering; see [Contig Stitching in +MiCall](stitcher.md) for the comparison. One of the critical challenges is the handling of ambiguous -nucleotides. The Stitcher's **Ambiguity Omission Principle**, which +nucleotides. The stitcher's **Ambiguity Omission Principle**, which aims to avoid propagating uncertainties, might lead to the exclusion of significant sequence data, resulting in the loss of potentially valuable variations or mutations. @@ -724,17 +742,17 @@ might not generalize well to other organisms or genomic regions. This can lead to over-splitting or under-splitting contigs, further fragmenting the consensus sequence. -Additionally, The Stitcher’s principle of scale-dependent credibility +Additionally, The stitcher’s principle of scale-dependent credibility might overlook important small-scale variations, such as single nucleotide polymorphisms (SNPs) or small indels, especially if they are lost in longer contigs deemed more reliable. Another critical limitation arises in the context of pipelines dealing -with proviral sequences. The Stitcher might attempt to "fix" sequences +with proviral sequences. The stitcher might attempt to "fix" sequences that are inherently "broken", such as those that are scrambled, contain long deletions, or exhibit hypermutation. In such cases, the tool's corrective measures may not be desirable, as they risk -introducing inaccuracies. This limitation makes the Stitcher +introducing inaccuracies. This limitation makes the referencefull stitcher unsuitable for certain pipelines where the integrity of such broken sequences should be preserved without alteration. From 97718e39526e706fbbf3102178df15558541198e Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:03:38 -0700 Subject: [PATCH 3/9] Rewrite stitcher.md as two-stitcher overview wrapper --- docs/design/stitcher.md | 787 ++++------------------------------------ 1 file changed, 77 insertions(+), 710 deletions(-) diff --git a/docs/design/stitcher.md b/docs/design/stitcher.md index 52e14e930..8422efae0 100644 --- a/docs/design/stitcher.md +++ b/docs/design/stitcher.md @@ -2,743 +2,110 @@ title: Contig Stitching in MiCall --- -DeNovo assembly does not invariably translate input reads into a -single contiguous sequence akin to a genomic consensus. Typically, -errors in input data lead to fragmented sequences — referred to as -contigs — which furthermore may overlap, thus encoding the same region -of a genome more than once. Assembling a unified consensus sequence -necessitates the systematic arrangement of these contigs while -addressing discrepancies within overlapping regions. That is the -Stitcher's function. +MiCall performs post-assembly stitching after de novo assembly. +This document explains why there are two stitchers and where each +detailed design lives: -# Structure +* [Referencefull stitcher](referencefull_stitcher.md) +* [Referenceless stitcher](referenceless_stitcher.md) -The Stitcher is a specialized component within the MiCall system. It -is designed to operate as an independent module which processes the -assembled contigs, generally derived from DeNovo assembler outputs, -and produce a singular, coherent sequence. +## The shared problem -## Modular Aspect +De novo assemblers such as IVA or Haploflow do not always turn input +reads into a single contiguous sequence. They may return multiple +contigs: fragmented sequences that can represent adjacent or +overlapping portions of the same biological sequence, sometimes +encoding the same region more than once. -The Stitcher maintains a distinct and autonomous role within -MiCall. Its implementation is fully isolated to the -`contig_stitcher*.py` files within the MiCall's source code. The -stitcher module can be run as a CLI script, separately from the rest -of the pipeline. The following command runs the Stitcher: +A post-assembly stitching/refinement step can therefore improve the +assembly by systematically arranging those contigs and resolving +discrepancies in overlapping regions. -```sh -micall contig_stitcher --help -``` - - - - - - - - - - -## Interaction - -Stitching is initiated either as a pipeline step in MiCall, or as a -command line call given above. In each case: - -**Input:** The Stitcher receives a single input file in CSV -format. This file contains 1 or more contigs that are the outcomes of -the previous assembly step, together with associated reference genome -information. These contigs are essentially segments of DNA -sequences. They can vary significantly in length. - -**Output:** The sole output from the Stitcher is a CSV -file. This file holds the stitched sequences -- longer or fully -continuous sequences that represent the genomic consensus formed by -merging the initial fragmented contigs, and additional metadata, -such as the inferred reference genome's name. - - - -# Operational procedure - -To clarify operations of the Stitcher, the subsequent section -introduces a vocabulary that is necessary for a precise description. - -## Definitions - -- An **input nucleotide** refers to a nucleotide of an initial - assembly contig sequence. -- A **reference nucleotide** refers to a nucleotide of a reference - genome sequence. -- A **non-conflicting nucleotide** is a **reference nucleotide** that - has at most one candidate **input nucleotide**. -- A **non-ambiguous nucleotide** is an **input nucleotide**, which has - a clear positioning with respect to all **input nucleotides** of all - other contigs associated with the same reference genome. In - particular, all **conflicting nucleotides** are **ambiguous - nucleotides** because they do not have a clear positioning with - respect to their competing **conflicting nucleotide**. -- An **overlap** is a continuos segement of **conflicting - nucleotides**. -- **Multidirectional alignment** is a property of a contig such that: - 1. the contig has aligned in multiple parts. - 2. some parts have been aligned to the forward strand, and some to - the reverse strand of the reference genome. -- **Cross-alignment** is a property of a contig such that: - 1. the contig has aligned in multiple parts. - 2. the contig-order of the aligned parts does not agree with the - reference-order of the aligned parts. -- A **non-aligned contig** is a contig that has been assinged a - reference sequence, but did not align to it. -- An **invalid contig** is a contig with **multidirectional - alignment**. -- A **stitched consensus** is a **valid contig** in the output of the - Stitcher. -- The **final output** refers to the contents of the only output CSV - file produced by the Stitcher. - -## Principles - -The reason the Stitcher operates effectively is due to its utilization -of reference genomes as additional source of truth. More precisely, -the Stitcher integrates two sets of data: - -1. Sequences generated by the initial assembly. -2. Sequences of reference genomes to which assembled contigs get aligned. - -We will say that 1. is the assembler's data, and 2. is aligner's. - -The core belief is that a reference genome can be used to enhance the -quality of and resolve conflicts within initial assembly contigs. - -In applying this approach, the Stitcher is guided by the following principles: - -### Principle of Scale-Dependent Credibility - -The reliability of sequence alignments increases as the length of the -aligned segment increases. -Therefore: - -- **Micro Scale**: For shorter segments, assembler's findings are more - reliable, because of expected abundance of small, local mutations - not present in the reference genome. - -- **Macro Scale**: For longer segments, the aligner's interpretations - are prioritized. The exponential decrease in alignment errors with - increased sequence length makes long alignments particularly - trustworthy. - -### Principle of Length Prioritization - -A longer contig typically arises from a greater number of reads -spanning a larger genomic region. While this does not imply more reads -per individual position, it suggests that the initial set of reads has -successfully assembled over a more extensive sequence, reflecting a -broader and more robust dataset. Moreover, aligning a longer sequence -to the reference genome is statistically less probable, compared to a -shorter sequence. This means that a successful alignment of a longer -contig to the reference genome provides further confidence in its -accuracy. - -Therefore in scenarios where multiple contigs cover the same region of -the reference genome, longer contigs are prioritized over shorter -ones. - -### Ambiguity Omission Principle - -To mitigate the potential propagation of uncertainties, any data that -lacks a definitive, unambiguous position within the reference genome -should be entirely excluded. This approach acknowledges that absolute -certainty in complex genomic datasets is often unattainable, and tries -to establish a reasonable default. - -## Regulations - -Guided by the previously outlined principles, -several precise regulations governing the Stitcher can be extracted: - -1. For every reference genome, at most one **stitched consensus** - must result. -2. No **ambiguous, non-conflicting nucleotide** - shall be included into the **final output**. -3. Every **non-conflicting-** and **non-ambiguous-** nucleotide - pertaining to a **valid contig** is required to be included in the - **stitched consensus** for the associated reference genome. -4. The relative positions of **non-conflicting-** and - **non-ambiguous-** nucleotides must be preserved in the **final output**. -5. All nucleotides present in the **final output** must exclusively - originate from the initial assembly data. - -## Setup - -The setup process for the Stitcher ensures that each contig is -properly aligned and prepared for the stitching process. The steps are -as follows: - -1. **Align Contigs**: Align each contig to its corresponding reference - genome to approximate their positions within a global reference - framework, allowing for spatial comparison between different contigs. - -2. **Split Multi-Alignment Contigs**: Split contigs that align to - multiple distinct parts of the reference genome into separate - segments. - -3. **Handle Reverse Complement**: Reverse complement contigs that - align to the reverse strand of the reference genome to ensure all - sequences are oriented in the same direction. - -4. **Sort Contigs**: Arrange the contigs based on their starting - positions along the reference genome. - -5. **Group by Reference**: Group contigs such that all contigs - associated with the same reference genome are processed together. - -These setup steps perform minimal alteration to the original contigs -and are primarily guided by straightforward, logical -considerations. Therefore, they do not require extensive -rationalization compared to the subsequent rules. - -## Rules of operation - -Stitching is an iterative process, governed by the following rules: - -### Rule 1: Merge Non-Overlapping Contigs - -1. **Verify Non-Overlap**: Ensure that the end of the first contig is - less or equal to the start of the second contig according to their - positions on the reference genome. - -2. **Delete adjacent non-aligned parts**: Filter out any non-aligned - nucleotides positioned after the first contig's aligned part and - before the second contig's aligned part. - -3. **Concatenate Sequences**: Directly join the end of the first - contig to the start of the second contig. - -#### Example: - -**Input:** - -![non overlaping example input illustration](stitcher_rule_1_input.svg) - -- Contig 1: Sequence = `GG[ATGCCC]AA`, aligned to Referece X at - position 10, with first two and last two nucleotides not aligned. -- Contig 2: Sequence = `AC[TTAG]TA`, aligned to Referece X at position - 30, with first two and last two nucleotides not aligned. - -**Procedure:** -- Verify that Contig 1 ends before Contig 2 begins. -- Delete non-aligned nucleotides resulting in Contig 1 = `GG[ATGCCC]` and Contig 2 = `[TTAG]TA`. -- Concatenate Contig 1 and Contig 2 to form `GG[ATGCCC][TTAG]TA`. - -**Result:** - -![non overlaping example result illustration](stitcher_rule_1_result.svg) - -- The new sequence, `GG[ATGCCCTTAG]TA`, spans positions 10 to 34 on the reference genome. - -#### Rationale - -There isn't many alternative actions available to us in these circumstances. -This enables us to consider all of them: - -1. **Leaving contigs as separate**: - - Separate contigs would result in multiple consensus outputs for one genome. - Thus it fails to comply with **regulation 1**. - -2. **Omitting the strip step**: - - Note that the adjacent non-aligned nucleotides of the two sequences - are **ambiguous, non-conflicting nucleotides**. Therefore, leaving - them in place violates **regulation 2**. - -3. **Introducing additional modifications**: - - Since given contigs do not overlap, every nucleotide in them is **non-conflicting**. - Additionally, we have stripped all the **ambiguous nucleotides**. - Therefore, all modifications that can be introduced - would either violate **regulation 3**, **regulation 4** or **regulation 5**. - -### Rule 2: Merge Overlapping Contigs - -1. **Verify Overlap**: Check if the ending position of the first - contig is greater than the starting position of the second contig. - -2. **Delete adjacent non-aligned parts**: Filter out any non-aligned - nucleotides positioned after the first contig's aligned part and - before the second contig's aligned part. - -3. **Align Overlapping Regions**: - - Extract the sequences from the overlapping region in both - contigs. - - Use a global alignment method to align these overlapping - sub-sequences. - -4. **Calculate Concordance Scores**: - - Compute concordance scores for each position within the - overlapping region. Importantly, the concordance calculation is - done purely between the aligned overlapping subsequences of the - contigs, with no regard to the reference genome sequence. The - concordance score represents how well the nucleotides from the - two contigs match at each position. - - The score is calculated using a sliding average approach, - emphasizing regions with high sequence agreement. - -5. **Determine Optimal Cut Point**: - - Identify the cut point based on the concordance scores such that - the it lies in the middle of regions with the highest - concordance. - - This means making cuts as far away from disagreeing nucleotides - as possible. - -6. **Segment and Combine**: - - Segment the overlapping sequences at the determined cut point. - - Concatenate the non-overlapping parts of the contigs with the - segmented parts from the overlapping region. - -#### Example - -**Input:** - -![overlaping example input illustration](stitcher_rule_2_input.svg) - -- Contig 1: Sequence = `G[GGCC A--TAC]T T`, aligned to Reference X from positions 10 to 19. -- Contig 2: Sequence = `--CCAC[AAATAC C]GGG`, aligned to Reference X from positions 14 to 20. - -**Procedure:** - -1. **Verify Overlap**: - - Contig 1 ends at position 19, and Contig 2 starts at position 14 - (both on Reference X), resulting in an overlap from positions 14 - to 19. - -2. **Delete adjacent non-aligned parts**: Contig 1 is right-stripped - to become `G[GGCCA--TAC]`, contig B is left-stripped to become - `[AAATACC]GGG`. - -3. **Align Overlapping Regions**: - - The overlaping sequence is `A--TAC` from contig A, and `AAATAC` - from contig B. - - Align them globally to produce the following alignments: `--ATAC` - and `AAATAC` +MiCall has two different ways to perform that refinement. -4. **Calculate Concordance**: - - Calculate concordance scores for positions 15 to 20, considering - only the overlap between the two aligned sequences. - - Approximate concordance: `[0.1, 0.2, 0.3, 0.8, 0.8, 0.3]`. +## Referencefull stitcher -5. **Determine Cut Point**: - - Use the computed concordance scores to identify the cut point. - - In this example, the highest concordance scores are around - positions with the score 0.9, so choose it as the cut point. +The referencefull stitcher +(`micall/utils/referencefull_contig_stitcher.py`, +`micall contig_stitcher with-references`) +is allowed to use a reference sequence and positions relative to that +reference. - ``` - Aligned sequences: +That gives it valuable information unavailable from the contigs +alone. In particular, it can reason about ordering and adjacency +using the reference coordinate system and can often produce +substantially more complete assemblies. - A: --ATAC - B: AAATAC +This is useful and intentional. - Concordance: - 0.1 0.2 0.3 0.8 0.8 0.3 +However, its output is therefore reference-guided. If the biological +sequence genuinely differs structurally from the reference — for +example through a large deletion, inversion, rearrangement, +duplication, or other scrambled structure — reference-derived +ordering can potentially obscure or normalize that structure. - Based on the concordance, cut between the positions: - A: --AT|AC - B: AAAT|AC - ``` +## Referenceless stitcher -6. **Segment and Combine**: - - Cut the sequences at the determined cut points. - - Combine sequence parts: `G[GGCC][--AT][AC][C]GGG`. +The referenceless stitcher +(`micall/utils/referenceless_contig_stitcher.py`, +`micall contig_stitcher without-references`) +exists for a different purpose. -**Result:** +Its goal is: -![overlaping example result illustration](stitcher_rule_2_result.svg) +> Improve the initial de novo assembly without using a reference +> sequence or reference-derived structural assumptions. -- The new sequence `G[GGC--ATACC]GGG` spans positions 10 to 20 on Reference X, - representing the most accurate combined sequence. +It relies only on evidence intrinsic to the sample, such as: -#### Rationale +* the contig sequences themselves; +* sequence overlap between contigs; +* short-read evidence; +* other sample-derived linkage evidence actually available to the + implementation. -This rule is similar to Rule 1, but deals with overlapping -regions. When contigs overlap, there is a need to choose a cut point -due to: +Reference independence is a deliberate property of this output, not +an implementation deficiency. -1. **Aligner Constraints**: The aligner constrains the size of the - overlapping sequence (by the **Principle of Scale-Dependent - Credibility**), making it impossible to keep both versions of the - overlapping region simultaneously. -2. **Small scale adjustments**: Overlaps are usually small enough that - assembler data is the highest quality data we have for the - nucleotide positions within it. Thus interleaving segments from - both contigs would again violate the **Principle of Scale Dependent - Credibility**. +The important use case is that a user may want an improved de novo +assembly while still preserving unusual biological structure exactly +as supported by the sample. Note that referenceless refinement is a +step applied after de novo assembly; "unstitched" output is not +automatically reference-free in the same sense. -We base the choice on concordance -scores, which measure the degree of agreement between the overlapping -sequences of the two contigs. We look for the highest concordance -because: +## Why both exist -**Choice of Cut Point**: -- If a cut point is chosen where concordance is lower than the - maximum, it implies that in the neighboring region around the cut - point, either to the left or right, there will almost certainly be - some incorrect nucleotides due to disagreement between the contigs. -- Conversely, if the concordance is high at the chosen cut point, the - neighboring region is similar between the two contigs. The selected - extensions (left of the cut point from the left contig and right of - the cut point from the right contig) are longer than the alternative - from the conflicting contig, ensuring greater trust in these regions - based on their length (by the **Principle of Length Prioritization**). +Neither stitcher is simply "better" than the other. They answer +different questions: -While this method of choosing the cut point based on concordance -scores aligns with the Principles, we acknowledge that there might be -other ways to determine the optimal cut point. However, given the -complexity of overlapping regions and the necessity to preserve -relative ordering, this concordance-based approach is the best we have -identified so far. - -### Rule 3: Split Contigs with Gaps Which Are Covered by Other Contigs - -1. **Identify Large Gaps**: - - For each contig, identify regions within its alignment to the - reference genome that lack coverage, i.e., gaps. Both small gaps - resulting from sequencing errors and large gaps are recognized. - - Significant gaps are determined based on a pre-defined - threshold. In the context of HIV genome analysis, a gap size of - greater than 21 nucleotides is considered significant due to - common RNA secondary structure phenomena. - -2. **Verify Coverage by Other Contigs**: - - For each identified significant gap, check if other contigs span - or cover this gap. Specifically, check if other contigs have - aligned reference coordinates that overlap with the coordinates - of the gap. - -3. **Split Contig at Gap Midpoint**: - - If a significant gap is covered by another contig, split the - contig containing the gap into two separate contigs at the - midpoint of the gap. - - Left-trim the new right contig segment and right-trim the new - left contig segment to remove ambiguity from their ends. - -4. **Update Contig List**: - - Replace the original contig with its two new segments in the list - of contigs. - -#### Example - -**Input:** - -![gap example input illustration](stitcher_rule_3_input.svg) - -- Contig 1: Sequence = `AGC[TTAC---------------------GGCACATATCATA]CTA`, - aligned to Reference X from positions 10 to 48. -- Contig 2: Sequence = `G[TGAC-----GGACG-TCGTCG--TACGATCAG]G`, - aligned to Reference X from positions 8 to 40. - -**Procedure:** - -1. **Identify Large Gaps**: - - Contig 1 has a significant gap between positions 14 and 35. - -2. **Verify Coverage by Other Contigs**: - - Contig 2 covers the gap region from positions 8 to 40. - -3. **Split Contig at Gap Midpoint**: - - Split Contig 1 into two parts at the midpoint of the gap (i.e., position 24). - This creates two new contigs: - - Contig 1a: Sequence = `AGC[TTAC----------]`, - aligned to Reference X from positions 10 to 24. - - Contig 1b: Sequence = `[-----------GGCACATATCATA]CTA`, - aligned to Reference X from positions 25 to 48. - - Trim the new segments: - - Contig 1a becomes `AGC[TTAC]`. - - Contig 1b becomes `[GGCACATATCATA]CTA`. - -4. **Update Contig List**: - - Discard the original Contig 1 and add Contig 1a and Contig 1b to - the list of contigs. - -**Result:** - -![gap example result illustration](stitcher_rule_3_result.svg) - -- Modified list of contigs now includes Contig 2, Contig 11, and Contig 12. - -#### Rationale - -The decision to split contigs at large gaps covered by other contigs -is grounded in the **Principle of Scale-Dependent -Credibility**. Assemblers can occasionally join sequence fragments -incorrectly if the end of one segment appears similar to the start of -another. Relying on the aligner's macro-scale credibility helps -identify these erroneous joins. Large gaps within a contig are -suspicious and suggest potential assembler errors, whereas small gaps -are generally due to sequencing errors or micro-scale mutations and do -not warrant splitting. By leveraging the aligner's high reliability on -a macro scale, we can effectively pinpoint these errors. If other -contigs cover large gaps, it confirms the aligner's indication that -the assembly might have joined unrelated segments. Splitting contigs -at the midpoint of significant gaps ensures that only those segments -supported by both the assembler's micro-scale data and the aligner's -macro-scale alignment are included in the final stitched consensus. - -The threshold for considering a gap significant is set at 21 -nucleotides. This value was chosen because it correlates with the -average pitch of the RNA helix, which reflects how reverse -transcription periodic deletions are structured around 21 nucleotides -in HIV sequences. Choosing this cutoff recognizes that deletions of -approximately this length are a common feature due to RNA secondary -structures and should not automatically warrant a split. This way, we -avoid splitting on every small gap, which is expected given the nature -of micro-scale mutations, but effectively identify and act on larger, -suspect gaps indicative of potential assembler errors. - -### Rule 4: Discard Contigs That Are Fully Covered By Other Contigs - -1. **Identify Covered Contigs**: - - For each contig in the input set, calculate its aligned interval on the reference genome. - - Identify intervals (regions) that are completely covered by input contigs. - -2. **Compare Intervals**: - - Assess the intervals of each contig to find any contig that falls entirely within the span of other contig intervals. - These are the contigs that are fully covered by others. - -3. **Discard Fully Covered Contigs**: - - Once identified, remove the covered contigs. - -#### Example - -**Input:** - -![covered example input illustration](stitcher_rule_4_input.svg) - -- Contig 1: Sequence = `A[ATCGA]GCT`, aligned to Reference X from positions 10 to 15. -- Contig 2: Sequence = `C[TAGTTG]A`, aligned to Reference X from positions 14 to 19. -- Contig 3: Sequence = `G[CGTACC]G`, aligned to Reference X from positions 12 to 17. - -**Procedure:** - -1. **Identify Covered Contigs**: - - Calculate the intervals: - - Contig 1: `[10-15]` - - Contig 2: `[14-19]` - - Contig 3: `[12-17]` - -2. **Compare Intervals**: - - Assess intervals and find Contig 3: `[12-17]` is completely within the intervals `[10-15]` of Contig 1 and `[14-19]` of Contig 2. - -3. **Discard Fully Covered Contigs**: - - Remove Contig 3 from the analysis. - -**Result:** - -![covered example result illustration](stitcher_rule_4_result.svg) - -- Unchanged remaining contigs Contig 1 and Contig 3. - -#### Rationale - -The underlying idea for this rule is founded on the two following principles: - -1. **Principle of Length Prioritization**: longer contigs are - inherently more reliable. - -2. **Principle of Scale-Dependent Credibility**: Fully covered contigs - might introduce small-scale inconsistencies that the longer - contig can resolve more credibly, given the enhanced reliability - associated with its length and alignment. - -Moreover, keeping all contigs would violate **Regulation 1**. - ---- - -**Note**: rules apply to contigs that are in the same group. - -# Diagnostics - -The Stitcher includes diagnostic tools to ensure transparency and -correctness throughout the stitching process. Two primary methods are -used for diagnostics: visualizer plots and traditional log -files. These tools help users understand and verify the decisions made -by the Stitcher during the stitching process. - -## The Optional Visualizer Tool - -The visualizer can be enabled through the `--plot` flag when running -the Stitcher executable. Running the Stitcher with this flag will -produce an SVG file that visualizes the stitching process, helping to -confirm and debug the Stitcher's operations. - -To use the visualizer, run the Stitcher with an additional argument -specifying the path to the output plot file. Here's an example of how -to stitch contigs and retrieve a visualizer plot: +```text +referencefull: + What assembly can we obtain when reference-derived structure + is allowed as evidence? -```sh -PYTHONPATH="/path/to/micall/repository" python3 -m micall.core.contig_stitcher "contigs.csv" "stitched_contigs.csv" --plot "visualized.svg" +referenceless: + What improvement over the de novo assembly can be justified + from sample-intrinsic evidence alone? ``` -**Command Line Arguments:** - -- `contigs.csv`: Input file in CSV format containing assembled - contigs and related information. -- `stitched_contigs.csv`: Output CSV file that will contain the - stitched contigs. -- `--plot visualized.svg`: The optional argument to generate a visual - representation of the stitching process, saved as `visualized.svg`. - -### Understanding the Output - -In practice, a visualizer plot might look something like this: - -![practical visualizer plot](stitcher_practical_plot.svg) - -From such a diagram, you can gain insights into the following aspects -of the stitching process: - -- **Reference genome**: The best matching reference genome for this - group of contigs was determined to be `HIV1-A1-RW-KF716472`. - -- **Dropped Contigs**: Contigs that were dropped due to being fully - covered by other contigs, as per Rule 4. In the example plot: - - Contigs 2, 4, 7, 8, and 6 were dropped. - -- **Split Contigs**: Contigs split at large gaps covered by other - contigs, according to Rule 3. The resulting parts are shown as - individual segments. - - Contig 1 was split around Contig 3, producing segments labeled as - 1.1 and 1.3. +The referencefull result can be more complete. The referenceless +result has stronger reference-independence semantics. These are +complementary products. -- **Joined Contigs**: Contigs that were merged due to overlap: - - Contigs 1 and 3, which were joined as per Rule 2, with - **ambiguous, non-conflicting** nucleotides discarded, shown as - segments labeled 1.2 and 3.1. - -- **Unaligned Contigs**: Contigs that failed to align to the reference - genome during the alignment step of the setup. - - Contig 5 failed to align. - -- **Contigs without a Reference**: Contigs for which a reference - genome could not be determined during the reference detection step - of the setup. - - Contigs 9 and 10 failed to determine a reference genome. - -Understanding these basics will help to interpret other scenarios -displayed by the visualizer plot. - -## Traditional Logs - -In addition to visual tools, the Stitcher produces traditional log -files that provide textual details of the stitching process. These -logs are crucial for debugging and understanding the sequence of -operations performed by the Stitcher. The verbosity of logs can be -adjusted using command-line options (`--verbose`, `--debug`, `--quiet`). - -Here is an example of typical log entries: +Conceptually: ```text -DEBUG:micall.core.contig_stitcher:Introduced contig 'contig.00001' (seq = TA...CA) of ref 'HIV1-C-BR-JX140663-seed', group_ref HIV1-A1-RW-KF716472-seed (seq = GA...AC), and length 7719. -DEBUG:micall.core.contig_stitcher:Introduced contig 'contig.00002' (seq = CG...AG) of ref 'HIV1-A1-RW-KF716472-seed', group_ref HIV1-A1-RW-KF716472-seed (seq = GA...AC), and length 1634. -... -DEBUG:micall.core.contig_stitcher:Contig 'contig.00006' produced 1 aligner hits. After connecting them, the number became 1. -DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00006' re-aligned as (5) at 7M...3D@[8,1433]->[7461,8946]. -DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00007' aligned at 76M...3D@[0,732]->[5536,6277]. -DEBUG:micall.core.contig_stitcher:Contig 'contig.00007' produced 1 aligner hits. After connecting them, the number became 1. -DEBUG:micall.core.contig_stitcher:Part 0 of contig 'contig.00007' re-aligned as (6) at 76M...3D@[0,732]->[5536,6277]. -... -DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 3D@[790,789]->[8280,8282]. -DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 19D@[1324,1323]->[8817,8835]. -DEBUG:micall.core.contig_stitcher:Ignored insignificant gap of (5), 2D@[1354,1353]->[8866,8867]. -... -DEBUG:micall.core.contig_stitcher:Created contigs (8) at 24M...1I@[14,3864]->[0,4558] and (9) at 708D...92I@[3865,7691]->[4559,9032] by cutting (1) at 24M...1I@[14,7691]->[0,9032] at cut point = 4558.5. -DEBUG:micall.core.contig_stitcher:Doing rstrip of (8) at 24M...1I@[14,3864]->[0,4558] (len 7719) resulted in (10) at 24M...1I@[14,3864]->[0,3850] (len 3865). -DEBUG:micall.core.contig_stitcher:Doing lstrip of (9) at 708D...92I@[3865,7691]->[4559,9032] (len 7719) resulted in (11) at 14M...1I@[0,3734]->[5267,9032] (len 3762). -DEBUG:micall.core.contig_stitcher:Split contig (1) at 24M...1I@[14,7691]->[0,9032] around its gap at [3864, 3863]->[3851, 5266]. Left part: (10) at 24M...1I@[14,3864]->[0,3850], right part: (11) at 14M...1I@[0,3734]->[5267,9032]. -... -DEBUG:micall.core.contig_stitcher:Created a frankenstein (34) at 24M...1I@[14,4185]->[0,4171] (len 4186) from [(26) at 24M...1I@[14,3041]->[0,3027] (len 3042), (28) at 271M2D3M2I395M@[0,670]->[3028,3698] (len 671), (30) at 152M@[0,151]->[3699,3850] (len 152), (31) at 321M@[0,320]->[3851,4171] (len 321)]. -DEBUG:micall.core.plot_contigs:Contig name (26) is displayed as '1.1'. -DEBUG:micall.core.plot_contigs:Contig name (36) is displayed as '1.3'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00002' is displayed as '2'. -DEBUG:micall.core.plot_contigs:Contig name (2) is displayed as '2'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00003' is displayed as '3'. -DEBUG:micall.core.plot_contigs:Contig name (31) is displayed as '3.2'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00004' is displayed as '4'. -DEBUG:micall.core.plot_contigs:Contig name (4) is displayed as '4'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00005' is displayed as '5'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00006' is displayed as '6'. -DEBUG:micall.core.plot_contigs:Contig name (5) is displayed as '6'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00007' is displayed as '7'. -DEBUG:micall.core.plot_contigs:Contig name (6) is displayed as '7'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00008' is displayed as '8'. -DEBUG:micall.core.plot_contigs:Contig name (7) is displayed as '8'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00009' is displayed as '9'. -DEBUG:micall.core.plot_contigs:Contig name 'contig.00010' is displayed as '10'. +de novo contigs (e.g. IVA, Haploflow) + | + +-- referencefull refinement -------> more complete, + | (reference + contigs + reads) reference-guided assembly + | + +-- referenceless refinement --------> improved assembly with + (contigs + reads, reference-independent + no reference structure) structure preserved ``` -The following points illustrate how these logs can facilitate -understanding the stitching process: - -- **Contig Introduction**: Provides details about the contigs - introduced for stitching. - - `Introduced contig 'contig.00001'...` - -- **Alignment Details**: Shows the alignment results for each contig. - - `Part 0 of contig 'contig.00006' re-aligned as (5) at - 7M...3D@[8,1433]->[7461,8946].` - -- **Gap Handling**: Indicates which gaps were ignored as - insignificant. - - `Ignored insignificant gap of (5), 3D@[790,789]->[8280,8282].` - -- **Splitting and Merging Contigs**: Documents the splitting of - contigs at identified gaps and merging of overlapping segments. - - `Split contig (1) at 24M...1I@[14,7691]->[0,9032]...` - - `Created a frankenstein (34) at 24M...1I@[14,4185]->[0,4171]...` - -- **Visualizer Compatibility**: The visualizer diagrams are produced - exclusively from these logs, ensuring compatibility and consistency - between the logs and visual output. - -# Limitations - -Following limitations stem from the choice of principles and various -assumptions that guide the Stitcher's operation. Understanding them -allows users to better interpret the results and apply post-processing -steps to mitigate potential issues. - -One of the critical challenges is the handling of ambiguous -nucleotides. The Stitcher's **Ambiguity Omission Principle**, which -aims to avoid propagating uncertainties, might lead to the exclusion -of significant sequence data, resulting in the loss of potentially -valuable variations or mutations. - -Moreover, the calculation of concordance in overlapping regions -assumes that local concordance is the best indicator of the correct -sequence. This approach may not fully account for complex genomic -rearrangements or context outside the overlap, potentially -compromising the accuracy of the stitched sequence. - -The predefined threshold for significant gaps, based on specific -assumptions about RNA secondary structures of organisms like HIV, -might not generalize well to other organisms or genomic regions. This -can lead to over-splitting or under-splitting contigs, further -fragmenting the consensus sequence. - -Additionally, The Stitcher’s principle of scale-dependent credibility -might overlook important small-scale variations, such as single -nucleotide polymorphisms (SNPs) or small indels, especially if they -are lost in longer contigs deemed more reliable. - -Another critical limitation arises in the context of pipelines dealing -with proviral sequences. The Stitcher might attempt to "fix" sequences -that are inherently "broken", such as those that are scrambled, -contain long deletions, or exhibit hypermutation. In such cases, the -tool's corrective measures may not be desirable, as they risk -introducing inaccuracies. This limitation makes the Stitcher -unsuitable for certain pipelines where the integrity of such broken -sequences should be preserved without alteration. +See the detailed designs for the algorithm behind each product: -Finally, the handling of multidirectional and cross-alignments may -fall short when addressing complex genomic rearrangements, such as -translocations or inversions, potentially resulting in misalignments -and stitching errors in the consensus sequence. +* [Referencefull stitcher](referencefull_stitcher.md) +* [Referenceless stitcher](referenceless_stitcher.md) From ea28fb59dd736c2ff299f15a0f202b77f2e4fd62 Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:04:42 -0700 Subject: [PATCH 4/9] Add referenceless stitcher design document --- docs/design/referenceless_stitcher.md | 477 ++++++++++++++++++++++++++ 1 file changed, 477 insertions(+) create mode 100644 docs/design/referenceless_stitcher.md diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md new file mode 100644 index 000000000..2db871744 --- /dev/null +++ b/docs/design/referenceless_stitcher.md @@ -0,0 +1,477 @@ +--- +title: Referenceless Contig Stitching in MiCall +--- + +This document describes the **referenceless** contig stitcher +(`micall/utils/referenceless_contig_stitcher.py`, +invoked as `micall contig_stitcher without-references`). +It is one of two stitching algorithms in MiCall; see +[Contig Stitching in MiCall](stitcher.md) for the overview and +[Referencefull stitcher](referencefull_stitcher.md) for the other +algorithm. + +Unless noted otherwise, "the stitcher" below means the referenceless +stitcher. + +## 1. Design objective + +The referenceless stitcher is a post-de-novo-assembly refinement step +that attempts to combine contigs only when the relationship can be +supported without reference-derived structural information. + +Reference independence is intentional. + +The stitcher must not decide that contigs belong together merely +because: + +* they align near each other on a reference; +* they have the expected reference ordering; +* they have the expected reference orientation; +* joining them would make the result look more like a canonical + genome; +* a subtype or reference label suggests that they should form one + genome. + +The implementation may carry metadata around (contig names, read +counts where available, per-run caches), but reference-derived +biological expectations must not be the evidence used to establish a +join. + +In the current code this restriction is structural: the +referenceless path takes FASTA contigs +(`Contig` / `ContigWithAligner` in +`micall/utils/contig_stitcher_contigs.py` and +`micall/utils/referenceless_contig_with_aligner.py`), never a +reference sequence or reference coordinates. Ordering, overlap +windows, alignments, and scores are all computed from the contig +sequences and from short-read evidence. The standard denovo pipeline +(`micall/drivers/sample.py`) feeds this stitcher the combined +assembler FASTA and writes a stitched FASTA; the referencefull CSV +fields (`ref`, `group_ref`, `match`) do not exist on this path. + +## 2. Why this constraint exists + +A sample can genuinely contain structure that differs from the +canonical reference, including things such as: + +* large deletions; +* inversions; +* rearrangements; +* duplications; +* recombinant or otherwise noncanonical structure. + +Such structure may itself be biologically important. + +A reference-guided assembly can sometimes improve contiguity by +imposing reference-derived order, but for analyses where structural +fidelity matters this can also be undesirable: the output may look +complete while silently normalizing away the unusual structure. + +The referenceless stitcher therefore deliberately asks a narrower +question: + +> What additional assembly structure is supported by the sample +> itself? + +When evidence is insufficient, leaving contigs separate can be +preferable to inventing a relationship. + +## 3. Error asymmetry + +A false negative generally means: + +```text +two truly related contigs remain separate +``` + +The original sequence evidence remains visible. A downstream user or +tool can still see both pieces. + +A false positive means: + +```text +two contigs that should remain separate are fused +``` + +That can destroy or obscure biological structure. The fused sequence +asserts a junction the sample did not support. + +Therefore the algorithm is intentionally conservative. Several +safeguards below — the minimum-agreement score, the independent +shared-k-mer check, the perfect-match rule for contained contigs, +and read validation around proposed joins — all raise the bar for +accepting a merge rather than lowering it. + +## 4. Inputs, outputs, and overall flow + +**Input:** a FASTA file of de novo contigs. Each record becomes a +`ContigWithAligner` (sequence plus cached aligner views; `reads_count` +is currently `None` on this file path). + +**Output:** a FASTA file of refined contigs. Some outputs combine +several input contigs; others pass through unchanged when no +supported join was found. + +The top-level flow (`stitch_consensus` in +`micall/utils/referenceless_contig_stitcher.py`) has two phases: + +1. **Overlap-path stitching** (`stitch_consensus_overlaps`): + iteratively select the most probable compatible path through the + remaining contigs, emit its merged sequence, remove its members + from consideration, and repeat. +2. **Greedy pairwise cleanup** (`o2_loop` / `try_combine_1`): try + every unordered pair once per round and merge the first acceptable + pair found, repeating until no acceptable pair remains. + +All per-run state lives in `ReferencelessStitcherContext` +(`micall/utils/contig_stitcher_context.py`): overlap, k-mer, +alignment, cutoff, and read-evidence caches plus read-validation +parameters. This keeps repeated pairwise checks cheap without +changing the acceptance rules. + +## 5. Evidence the stitcher uses + +A candidate join must survive every applicable check below. The +checks are conjunctive safeguards, not alternative theories of +relatedness: + +* **Terminal overlap placement** — a coarse convolution estimate of + where two contigs would sit relative to each other, reduced to a + terminal overlap window (section 6). +* **Overlap alignment** — a global pairwise alignment of the two + overlap windows (`align_queries` in + `micall/utils/overlap_stitcher.py`). +* **Overlap scoring** — a statistical score of the alignment that + must clear a minimum-agreement threshold derived from + `MIN_MATCHES = 40` (section 6). +* **Shared k-mers** — an independent exact-match requirement + (`KMER_SIZE = 30`) that rejects statistically plausible overlaps + with no shared exact sequence (section 7). +* **Containment handling** — a separate perfect-match rule when one + contig is fully covered by another (section 8). +* **Raw-read support** — an independent check that the proposed + junction is crossed by sample reads (section 9 and + [Read-Supported Join Validation](../specs/referenceless-stitcher-read-information-handling.md)). +* **Path competition** — individually plausible edges compete for + membership in a bounded set of candidate paths; only winners + survive (section 10). +* **Concordance-based construction** — the merged sequence itself is + cut where local agreement is strongest (section 11). + +No step consults a reference genome, reference coordinates, or +expected gene order. + +## 6. Overlap discovery and scoring + +### 6.1 Coarse placement + +The stitcher first needs a hypothesis for *where* two contigs +overlap. `find_maximum_overlap` (in +`micall/utils/referenceless_contig_with_aligner.py`, built on +`micall/utils/find_maximum_overlap.py` and +`micall/utils/overlap_stitcher.py`) answers this with a fast +convolution: + +* each contig is expanded into per-symbol indicator vectors; +* each vector is smoothed with an exponential drop-off + (`exp_dropoff_array`, factor 8), so near-misses still contribute + weakly and the estimate tolerates small local disagreements; +* cross-correlating the softened vectors across all shifts yields an + expected-match profile; +* each shift is scored with the same statistical overlap model used + later (`calculate_overlap_score`), and the best shift becomes the + candidate placement. + +A non-positive best value means "no convincing overlap": the pair is +abandoned (`shift == 0` in `get_overlap`). Otherwise the shift is +converted to a terminal overlap window (`Overlap(shift, size)` in +`micall/utils/referenceless_contig_stitcher_overlap.py`; +`compute_overlap_size`, `normalize_orientation`, +`initial_overlap_windows`). + +The model here is deliberately coarse. It proposes a window worth +aligning; it does not itself accept a merge. + +### 6.2 End-aware anchoring and cutoffs + +Before aligning, the stitcher trims the problem to the part of each +contig it is willing to trust. `map_overlap` queries lightweight +`mappy`-backed views of a contig under a stitching relation: + +* `"left"` — anchor at the left end (keep the earliest start); +* `"right"` — anchor at the right end (keep the latest end); +* `"cover"` — unconstrained mapping (used when one contig may be + fully covered). + +End anchoring is implemented with synthetic homogeneous padding +(`ForwardAligner` / `ReversedAligner`) so the underlying mapper +respects the chosen edge rather than sliding to an interior repeat. + +The returned anchors become cutoffs +(`compute_overlap_cutoffs` / `find_overlap_cutoffs`, with the +`cutoffs_left_*` / `cutoffs_right_*` helpers) delimiting the overlap +region to align and score. A theoretical upper bound +(`find_max_overlap_length`) can additionally narrow the contig view +presented to the aligner when the required score cannot use the full +length. Cutoffs are cached per contig pair; because they are +monotonic in the acceptance threshold, a cutoff computed for a lower +threshold remains valid for a higher one. + +### 6.3 Alignment and concordance + +The trimmed overlap windows are globally aligned with Biopython's +`PairwiseAligner` in global mode with penalized end gaps +(`align_queries`). Matches, mismatches, and indels in that alignment +are the evidence for or against the join. + +From the alignment the stitcher derives two things: + +* a **concordance** profile (`calculate_concordance`): a sliding + average of per-position agreement, accumulated forward and + backward with a square-root weighting so sustained runs of matches + score higher than isolated matches. The eventual merge point is + chosen where concordance is strongest + (`sort_concordance_indexes`), with ties broken toward the middle + of the overlap so cuts stay far from disagreements + (`merge_by_concordance`; see section 11). +* an **overlap score** (`calculate_overlap_score` with + `score_alignment`): a z-like rarity score over a four-letter + alphabet, generalized for correlated genomic sequence with an + exponent (`alpha = -0.60`, i.e. standard deviation growing as + `L^0.8` rather than `sqrt(L)`). Higher means more unexpected under + the null model and therefore stronger evidence. The scored length + includes a small bonus (`+1` for ordinary overlaps, `+2` for + covering overlaps) expressing that the overlap is flanked by + non-matching context. + +### 6.4 Minimum agreement and fast rejection + +A minimum amount of sequence agreement is required. The raw +threshold is set by `MIN_MATCHES = 40`: +`ACCEPTABLE_STITCHING_SCORE` is the transformed score of an overlap +just above that size, and every candidate edge must ultimately reach +at least the pool's minimum acceptable score (section 10). + +To avoid wasted alignments, `try_combine_contigs` / +`precheck_and_prepare_overlap` applies optimistic upper bounds first: +if even a perfect overlap of the available lengths +(`max_possible_overlap_score`) or of the discovered window +(`optimistic_overlap_score`) cannot reach the needed score, the pair +is rejected before alignment. The transformed score +(`calculate_referenceless_overlap_score`) monotonically amplifies the +raw score (a `999 + (999 * base)^2` shaping) and keeps genuine +scores far from the `SCORE_EPSILON = 1` sentinel used for +covered-contig bookkeeping, so scoring and containment signalling can +never be confused. + +## 7. Shared-k-mer requirement + +Statistical similarity alone is not always sufficient evidence of a +meaningful overlap: repeats, low-complexity sequence, and smoothed +convolution estimates can all produce plausible-looking scores for +unrelated contigs. + +The stitcher therefore applies an independent shared-k-mer check +(`get_kmers`, `does_share_kmers`, `get_overlap`): + +* every contig yields the set of its exact k-mers with + `KMER_SIZE = 30`; +* if both contigs are at least 30 bases long and their k-mer sets + are disjoint, the pair is rejected before any alignment, however + good its statistical score would have been; +* k-mer sets are cached per sequence in the stitching context. + +Requiring shared exact sequence provides additional specificity: a +genuine terminal overlap of sufficient length should normally share +at least one 30-mer, while coincidental similarity often shares +none. + +Its scope is deliberately narrow: + +* contigs shorter than 30 bases are exempt (they cannot contain a + full k-mer to share); +* sharing a k-mer does **not** by itself prove an overlap — it only + permits the statistical and read checks to proceed; +* failing to share a k-mer rejects the candidate merge but does not + prove the contigs are biologically unrelated. + +## 8. Covered-contig handling + +When one contig is fully covered by another +(`calculate_covered`: one sequence length is at most the overlap +size), the stitcher does not perform a normal concordance merge. +Instead it applies a strict perfect-match rule in +`try_combine_contigs`: + +* the covered sequence and the corresponding window of the larger + contig are aligned; +* if every base of the overlap matches (`number_of_matches == + overlap.size`), the larger contig is kept and the smaller one is + recorded as contained (returned with `SCORE_EPSILON` and + `covered_input` marking which side was absorbed); +* any mismatch means no merge at all: the pair is rejected. + +The conservative rationale is: + +> An imperfect contained contig may represent error or redundancy, +> but it may also encode real variation. Without sufficient +> evidence, silently absorbing it would destroy that uncertainty. + +This does not claim that all imperfect contained contigs are +biologically important. The point is that the algorithm intentionally +refuses to assume that they are not. Exact duplicates collapse +safely; near-duplicates are left alone for downstream analysis +rather than fused on statistical grounds. + +Containment is tracked separately from path membership +(`ContigsPath.contigs_ids` versus `contains_contigs_ids`), so a +perfectly covered contig is remembered as explained without +contributing a second copy of its sequence to the merged result. + +## 9. Read support + +Even a join that passes overlap, k-mer, and containment checks still +proposes a new junction — a sequence that neither input contig +contained on its own. The raw sample reads provide independent +evidence about whether that junction is supported. + +At a high level (`check_merged_sequence_support` and its cached +caller in `try_combine_contigs`): + +* the candidate merged contig and its join boundary (`join_boundary` + from `merge_by_concordance`) define a cut position; +* the stitcher requires exact placements of sample reads that + strictly cross the cut, plus exact coverage of every base in a + read-length-sized window centred on the cut; +* placements are canonicalized (`min(seq, reverse_complement(seq))`) + so either strand counts, weighted by FASTQ multiplicity, with each + valid placement contributing; +* if support is below `minimum_read_depth`, the merge is rejected + (emitting `ReadSupportRejected` in debug2). + +Full contracts — cut-spanning definitions, window geometry, +counting model including the accepted placement-times-multiplicity +overcounting tradeoff, disabled states (`read_index is None` or +`minimum_read_depth == 0` accepts; enabled-but-empty `{}` rejects), +CLI flags (`--fastq1` / `--fastq2`, `--minimum-read-depth`, +`--read-length`), and pipeline defaults (enabled with trimmed FASTQs +at depth 1 in `micall/drivers/sample.py`) — belong to the +implementation spec and are not repeated here. See: + +* [Read-Supported Join Validation](../specs/referenceless-stitcher-read-information-handling.md) + +## 10. Path and candidate competition + +The stitcher does not accept every individually plausible overlap +independently. Candidate relationships compete, and compatible joins +form paths. + +The mechanism (`ContigsPath` in +`micall/utils/referenceless_contig_path.py`, `Pool` in +`micall/utils/referenceless_contig_stitcher_pool.py`, +`calculate_all_paths` / `extend_by_1` / `calc_extension`): + +* every remaining contig starts as a singleton seed path with + `SCORE_NOTHING = 0`, seeds sorted longest-first; +* each cycle tries to extend every retained path with every + remaining contig via `try_combine_contigs`, scoring extensions by + summing edge scores (`combine_scores`); +* a bounded `Pool` (a `SortedRing` plus sequence deduplication) + keeps only the best paths: same merged sequence keeps only its + highest score, and the pool's minimum acceptable score only rises, + pruning progressively weaker extensions; +* capacity per cycle is set by + `intrapolate_number_of_alternatives` (`999 / max(1, n - 2)`, + clamped to `[1, 999]`), bounding total work while still exploring + alternatives when few contigs remain; +* the best surviving path (`find_most_probable_path`) is emitted, + its members (including contained ones) are removed from the + remaining set, and the loop repeats; +* if the best path is a singleton, the stitcher gives up on further + path extension (`GiveUp`) and emits the rest unchanged; +* the later `o2_loop` performs a final greedy pairwise pass for + leftovers. + +The important consequence is: + +> Evidence is evaluated locally at candidate edges, while a final +> multi-contig result can be produced through a chain of supported +> relationships. + +A final component therefore asserts a chain of pairwise-supported +joins, not an all-pairs guarantee about every member. Two contigs at +opposite ends of an emitted component were never directly compared; +they are joined because each link in the chain cleared the +thresholds. + +## 11. How merged sequence is constructed + +When a non-covering pair is accepted, +`merge_by_concordance` builds the output from the global alignment +of the two overlap windows: + +* the alignment's concordance profile selects the best split index; +* the left part of the left alignment and the right part of the + right alignment (dashes removed) become the overlap contribution; +* outer remainders (`left_remainder`, `right_remainder`) are + prepended and appended unchanged; +* the boundary between the left-derived and right-derived overlap + chunks is recorded as `join_boundary` for read validation. + +The merged contig sums input read counts only when both are +available; on the file path both are currently `None`, so the result +carries `None`. The merged sequence then participates in further +extension rounds as an ordinary contig. + +## 12. Limitations + +Reference independence does not mean that the algorithm can always +recover biological truth. + +If two distinct molecules contain a long, highly similar region and +the sample-intrinsic evidence available to the algorithm does not +phase that region to distinguishing sequence — no conclusive overlap +placement, no shared k-mer that anchors the true junction, no +cut-spanning reads — the relationship may be fundamentally ambiguous +to the stitcher. It will leave the contigs separate, even if a +reference would have suggested an order. + +That restriction is about the evidence the implementation is allowed +to use, not an absolute claim about all reference-independent or all +short-read methods. Additional sample-intrinsic evidence such as +longer reads or stronger linkage could, in principle, resolve such +ambiguities without using a reference. The defining restriction is: + +> Do not use external reference-derived structural assumptions to +> resolve the ambiguity. + +Other limits follow from the conservative design: + +* exact read matching undercounts true support when reads carry + errors or variation relative to the contigs; see the spec for the + accepted tradeoffs; +* in repetitive sequence, one read may contribute at multiple + placements, inflating support counts without creating support out + of nothing; +* the minimum-agreement threshold (`MIN_MATCHES = 40`) and k-mer + size (30) will miss true short overlaps — a deliberate price for + specificity; +* greedy path selection and the bounded pool can in principle prefer + a locally strong chain over a globally better one; capacity tuning + bounds the search rather than guaranteeing optimality. + +## 13. Relationship to referencefull + +The [referencefull stitcher](referencefull_stitcher.md) is allowed +to use information the referenceless stitcher deliberately excludes: +a reference sequence, reference coordinates, and reference-derived +ordering and adjacency. + +Referencefull may therefore intentionally resolve cases that +referenceless leaves unresolved — for example, placing two +non-overlapping contigs in reference order, or bridging a gap with +no sample-supported overlap. That is not automatically a failure of +either algorithm. One trades structural caution for completeness; +the other trades completeness for reference independence. See +[Contig Stitching in MiCall](stitcher.md) for guidance on which +question each output answers. From a63aa09ad06340ca7eed3f40bbc8ae4c86cba367 Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:33:33 -0700 Subject: [PATCH 5/9] Fix overlap-score exponent: code divides by L^0.60 --- docs/design/referenceless_stitcher.md | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md index 2db871744..681c6bf3c 100644 --- a/docs/design/referenceless_stitcher.md +++ b/docs/design/referenceless_stitcher.md @@ -236,10 +236,15 @@ From the alignment the stitcher derives two things: (`merge_by_concordance`; see section 11). * an **overlap score** (`calculate_overlap_score` with `score_alignment`): a z-like rarity score over a four-letter - alphabet, generalized for correlated genomic sequence with an - exponent (`alpha = -0.60`, i.e. standard deviation growing as - `L^0.8` rather than `sqrt(L)`). Higher means more unexpected under - the null model and therefore stronger evidence. The scored length + alphabet. The centered match excess `(4*M - L)` is scaled by + `L^-0.60` — that is, the implementation divides by `L^0.60` + rather than by the independent-match `sqrt(L)`, penalizing long + overlaps more strongly. (The function's own comments motivate + this with a correlated-match model whose standard deviation grows + as `L^a` with an empirical `a ≈ 0.8`, but the constant the code + actually uses is `0.60`; this document describes the computation, + not the comment.) Higher means more unexpected under the null + model and therefore stronger evidence. The scored length includes a small bonus (`+1` for ordinary overlaps, `+2` for covering overlaps) expressing that the overlap is flanked by non-matching context. From 6084dd2e673fa46faf486be6ea28d0b5ac768ee6 Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:33:41 -0700 Subject: [PATCH 6/9] Tighten k-mer gate wording to shared exact 30-mer --- docs/design/referenceless_stitcher.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md index 681c6bf3c..6727369ae 100644 --- a/docs/design/referenceless_stitcher.md +++ b/docs/design/referenceless_stitcher.md @@ -146,7 +146,7 @@ relatedness: `MIN_MATCHES = 40` (section 6). * **Shared k-mers** — an independent exact-match requirement (`KMER_SIZE = 30`) that rejects statistically plausible overlaps - with no shared exact sequence (section 7). + with no shared exact 30-mer (section 7). * **Containment handling** — a separate perfect-match rule when one contig is fully covered by another (section 8). * **Raw-read support** — an independent check that the proposed From de962d287915c651d785fc602631fb89de3454aa Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:33:48 -0700 Subject: [PATCH 7/9] Soften false-positive claim to local evidence vs phasing --- docs/design/referenceless_stitcher.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md index 6727369ae..b45ff2e8b 100644 --- a/docs/design/referenceless_stitcher.md +++ b/docs/design/referenceless_stitcher.md @@ -94,7 +94,12 @@ two contigs that should remain separate are fused ``` That can destroy or obscure biological structure. The fused sequence -asserts a junction the sample did not support. +asserts that both sides came from the same biological molecule — a +claim that local evidence alone cannot always establish. Two +distinct molecules can share a long, highly similar region, so a +proposed junction may be locally supported by sample reads yet still +join sequences that were never adjacent in any single molecule. +Local evidence is necessary but not sufficient for phasing. Therefore the algorithm is intentionally conservative. Several safeguards below — the minimum-agreement score, the independent From 16d6adaf134a3efb5cf7876b60ee838bf9f1fc7d Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:33:54 -0700 Subject: [PATCH 8/9] Rewrite limitations around local evidence vs phasing --- docs/design/referenceless_stitcher.md | 31 +++++++++++++++++++-------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md index b45ff2e8b..d5ba0f92e 100644 --- a/docs/design/referenceless_stitcher.md +++ b/docs/design/referenceless_stitcher.md @@ -436,15 +436,28 @@ extension rounds as an ordinary contig. ## 12. Limitations Reference independence does not mean that the algorithm can always -recover biological truth. - -If two distinct molecules contain a long, highly similar region and -the sample-intrinsic evidence available to the algorithm does not -phase that region to distinguishing sequence — no conclusive overlap -placement, no shared k-mer that anchors the true junction, no -cut-spanning reads — the relationship may be fundamentally ambiguous -to the stitcher. It will leave the contigs separate, even if a -reference would have suggested an order. +recover biological truth — in either direction. + +**False separation.** If two truly related contigs leave too little +distinguishing evidence — no conclusive overlap placement, no shared +k-mer anchoring the true junction, no cut-spanning reads — the +relationship is ambiguous to the stitcher and it leaves the contigs +separate, even if a reference would have suggested an order. + +**False joining is the more important ceiling.** Two biologically +distinct molecules (different haplotypes, repeat copies, or +recombinant forms) can share a long, highly similar region. In that +case every *local* check the stitcher performs — overlap alignment +score, a shared 30-mer, exact read placements crossing the chosen +cut plus local window coverage — can look convincing while still +failing to establish that the two sides came from the same molecule. +The read check validates a local junction hypothesis: it asks +whether sample reads exactly match the sequence around the proposed +cut. It does not phase the flanking sequence to a single haplotype. +Short reads falling entirely inside the shared region are consistent +with either origin; only linkage reaching into distinguishing +sequence (or longer reads) could resolve which molecule each side +belongs to — without any reference. That restriction is about the evidence the implementation is allowed to use, not an absolute claim about all reference-independent or all From 62a4dc3f79c3056b1f04f3ee96a0e06cc6b4111d Mon Sep 17 00:00:00 2001 From: Vitaliy Mysak Date: Wed, 16 Sep 2026 16:34:07 -0700 Subject: [PATCH 9/9] Make wrapper diagram conceptual and fix evidence labels --- docs/design/stitcher.md | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/docs/design/stitcher.md b/docs/design/stitcher.md index 8422efae0..757de3a4b 100644 --- a/docs/design/stitcher.md +++ b/docs/design/stitcher.md @@ -34,7 +34,9 @@ reference. That gives it valuable information unavailable from the contigs alone. In particular, it can reason about ordering and adjacency using the reference coordinate system and can often produce -substantially more complete assemblies. +substantially more complete assemblies. Its additional inputs are +the reference sequence and, optionally, remap read counts — not the +raw-read junction validation used on the referenceless path. This is useful and intentional. @@ -90,22 +92,31 @@ referenceless: The referencefull result can be more complete. The referenceless result has stronger reference-independence semantics. These are -complementary products. +complementary algorithms — different assembly interpretations of +the same contigs. -Conceptually: +Conceptually — not the literal pipeline dataflow — MiCall has two +refinement semantics: ```text de novo contigs (e.g. IVA, Haploflow) | +-- referencefull refinement -------> more complete, - | (reference + contigs + reads) reference-guided assembly + | (reference + contigs reference-guided assembly + | [+ optional remap counts]) | +-- referenceless refinement --------> improved assembly with (contigs + reads, reference-independent no reference structure) structure preserved ``` -See the detailed designs for the algorithm behind each product: +In the current denovo pipeline (`micall/drivers/sample.py`), the +referencefull stitcher in fact runs on the referenceless stitcher's +output, so the two branches above describe the refinement semantics +rather than two independent pipeline paths. + +See the detailed designs for the algorithm behind each +interpretation: * [Referencefull stitcher](referencefull_stitcher.md) * [Referenceless stitcher](referenceless_stitcher.md)