diff --git a/docs/design/referencefull_stitcher.md b/docs/design/referencefull_stitcher.md new file mode 100644 index 000000000..0aa631e8d --- /dev/null +++ b/docs/design/referencefull_stitcher.md @@ -0,0 +1,762 @@ +--- +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 +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 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 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 referencefull stitcher: + +```sh +micall contig_stitcher with-references --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 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 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, +such as the inferred reference genome's name. + + + +# Operational procedure + +To clarify operations of the referencefull 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 + referencefull stitcher. +- The **final output** refers to the contents of the only output CSV + file produced by the referencefull stitcher. + +## Principles + +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: + +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 referencefull 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 referencefull 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 referencefull 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 2. + +#### 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 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. + +## 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 with-references "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 referencefull stitcher's operation. Understanding them +allows users to better interpret the results and apply post-processing +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 +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 referencefull 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. diff --git a/docs/design/referenceless_stitcher.md b/docs/design/referenceless_stitcher.md new file mode 100644 index 000000000..d5ba0f92e --- /dev/null +++ b/docs/design/referenceless_stitcher.md @@ -0,0 +1,500 @@ +--- +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 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 +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 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 + 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. 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. + +### 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 — 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 +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. diff --git a/docs/design/stitcher.md b/docs/design/stitcher.md index 52e14e930..757de3a4b 100644 --- a/docs/design/stitcher.md +++ b/docs/design/stitcher.md @@ -2,743 +2,121 @@ 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`. +MiCall has two different ways to perform that refinement. -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` +## Referencefull stitcher -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]`. +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. -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. +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. Its additional inputs are +the reference sequence and, optionally, remap read counts — not the +raw-read junction validation used on the referenceless path. - ``` - Aligned sequences: +This is useful and intentional. - A: --ATAC - B: AAATAC +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. - Concordance: - 0.1 0.2 0.3 0.8 0.8 0.3 +## Referenceless stitcher - Based on the concordance, cut between the positions: - A: --AT|AC - B: AAAT|AC - ``` +The referenceless stitcher +(`micall/utils/referenceless_contig_stitcher.py`, +`micall contig_stitcher without-references`) +exists for a different purpose. -6. **Segment and Combine**: - - Cut the sequences at the determined cut points. - - Combine sequence parts: `G[GGCC][--AT][AC][C]GGG`. +Its goal is: -**Result:** +> Improve the initial de novo assembly without using a reference +> sequence or reference-derived structural assumptions. -![overlaping example result illustration](stitcher_rule_2_result.svg) +It relies only on evidence intrinsic to the sample, such as: -- The new sequence `G[GGC--ATACC]GGG` spans positions 10 to 20 on Reference X, - representing the most accurate combined sequence. +* the contig sequences themselves; +* sequence overlap between contigs; +* short-read evidence; +* other sample-derived linkage evidence actually available to the + implementation. -#### Rationale +Reference independence is a deliberate property of this output, not +an implementation deficiency. -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: +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. -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**. +## Why both exist -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: +Neither stitcher is simply "better" than the other. They answer +different questions: -**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: +```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. +The referencefull result can be more complete. The referenceless +result has stronger reference-independence semantics. These are +complementary algorithms — different assembly interpretations of +the same contigs. -- **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: +Conceptually — not the literal pipeline dataflow — MiCall has two +refinement semantics: ```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 reference-guided assembly + | [+ optional remap counts]) + | + +-- 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. +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. -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 +interpretation: -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)