Reference-less contig stitching

micall/core/plot_contigs.py

@@ -1,5 +1,5 @@
 import typing
-from typing import Dict, Tuple, List, Set, Iterable, NoReturn
+from typing import Dict, Tuple, List, Set, Iterable, NoReturn, Sequence
 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter, FileType
 from collections import Counter, defaultdict
 from csv import DictReader
@@ -555,11 +555,11 @@ def hit_to_insertions(contig: GenotypedContig, hit: CigarHit):
         yield CigarHit.from_default_alignment(q_st=hit.q_ei + 1, q_ei=len(contig.seq) - 1,
                                               r_st=hit.r_ei + 1, r_ei=hit.r_ei)
 
-    def hits_to_insertions(contig: GenotypedContig, hits: List[CigarHit]):
+    def hits_to_insertions(contig: GenotypedContig, hits: Iterable[CigarHit]):
         for hit in hits:
             yield from hit_to_insertions(contig, hit)
 
-    def record_initial_hit(contig: GenotypedContig, hits: List[CigarHit]):
+    def record_initial_hit(contig: GenotypedContig, hits: Sequence[CigarHit]):
         insertions = [gap for gap in hits_to_insertions(contig, hits)]
         unaligned_map[contig.id] = insertions
 

micall/drivers/sample.py

@@ -9,7 +9,6 @@
 from micall.core.aln2counts import aln2counts
 from micall.core.amplicon_finder import write_merge_lengths_plot, merge_for_entropy
 from micall.core.cascade_report import CascadeReport
-from micall.core.contig_stitcher import contig_stitcher
 from micall.core.coverage_plots import coverage_plot, concordance_plot
 from micall.core.plot_contigs import plot_genome_coverage
 from micall.core.prelim_map import prelim_map
@@ -21,6 +20,10 @@
 from micall.g2p.fastq_g2p import fastq_g2p, DEFAULT_MIN_COUNT, MIN_VALID, MIN_VALID_PERCENT
 from micall.utils.driver_utils import makedirs
 from micall.utils.fasta_to_csv import fasta_to_csv
+from micall.utils.csv_to_fasta import csv_to_fasta, NoContigsInCSV
+from micall.utils.referencefull_contig_stitcher import referencefull_contig_stitcher
+from micall.utils.referenceless_contig_stitcher import referenceless_contig_stitcher
+from micall.utils.cat import cat as concatenate_files
 from contextlib import contextmanager
 
 logger = logging.getLogger(__name__)
@@ -150,7 +153,7 @@ def get_default_path(self, output_name):
         if self.scratch_path is None:
             raise AttributeError(
                 'Unknown output {} and no scratch path.'.format(output_name))
-        for extension in ('csv', 'fastq', 'pdf', 'svg', 'png'):
+        for extension in ('csv', 'fastq', 'pdf', 'svg', 'png', 'fasta'):
             if output_name.endswith('_'+extension):
                 file_name = output_name[:-(len(extension)+1)] + '.' + extension
                 break
@@ -426,18 +429,30 @@ def run_denovo(self, excluded_seeds):
                    merged_contigs_csv,
                    )
 
+        with open(self.merged_contigs_csv, 'r'):
+            try:
+                csv_to_fasta(self.merged_contigs_csv, Path(self.merged_contigs_fasta))
+            except NoContigsInCSV:
+                Path(self.merged_contigs_fasta).touch()
+
+        concatenate_files(inputs=[self.unstitched_contigs_fasta,
+                                  self.merged_contigs_fasta],
+                          output=self.combined_contigs_fasta)
+
+        with open(self.combined_contigs_fasta, 'r') as combined_contigs_fasta, \
+             open(self.stitched_contigs_fasta, 'w') as stitched_contigs_fasta:
+            referenceless_contig_stitcher(combined_contigs_fasta, stitched_contigs_fasta)
+
         with open(self.unstitched_contigs_csv, 'w') as unstitched_contigs_csv, \
-             open(self.merged_contigs_csv, 'r') as merged_contigs_csv, \
              open(self.blast_csv, 'w') as blast_csv:
-            fasta_to_csv(Path(self.unstitched_contigs_fasta),
+            fasta_to_csv(Path(self.stitched_contigs_fasta),
                          unstitched_contigs_csv,
-                         merged_contigs_csv,
                          blast_csv=blast_csv,
                          )
 
         with open(self.unstitched_contigs_csv, 'r') as unstitched_contigs_csv, \
              open(self.contigs_csv, 'w') as contigs_csv:
-            contig_stitcher(unstitched_contigs_csv, contigs_csv, self.stitcher_plot_svg)
+            referencefull_contig_stitcher(unstitched_contigs_csv, contigs_csv, self.stitcher_plot_svg)
 
         logger.info('Running remap on %s.', self)
         if self.debug_remap:

micall/main.py

@@ -88,6 +88,9 @@
     "micall/monitor/update_qai.py",
     "micall/monitor/micall_watcher.py",
     "micall/tcr/igblast.py",
+    "micall/utils/find_maximum_overlap.py",
+    "micall/utils/csv_to_fasta.py",
+    "micall/utils/cat.py",
 ]
 
 

micall/tests/test_contig_stitcher.py

@@ -7,27 +7,27 @@
 
 from aligntools import CigarActions, CigarHit, Cigar
 
-import micall.core.contig_stitcher as stitcher
-from micall.core.contig_stitcher import (
+import micall.utils.referencefull_contig_stitcher as stitcher
+from micall.utils.referencefull_contig_stitcher import (
     split_contigs_with_gaps,
     stitch_contigs,
     GenotypedContig,
     merge_intervals,
     find_covered_contig,
     stitch_consensus,
-    calculate_concordance,
+    calculate_concordance_norm,
     align_all_to_reference,
-    disambiguate_concordance,
     lstrip,
     rstrip,
 )
+from micall.utils.overlap_stitcher import disambiguate_concordance
 from micall.core.plot_contigs import plot_stitcher_coverage
 from micall.tests.utils import mock_align_consensus, MockAlignment, fixed_random_seed
 from micall.tests.test_fasta_to_csv import check_hcv_db, DEFAULT_DATABASE  # activates the fixture
 from micall.tests.test_remap import load_projects  # activates the "projects" fixture
 
 
-logging.getLogger("micall.core.contig_stitcher").setLevel(logging.DEBUG)
+logging.getLogger("micall.utils.referencefull_contig_stitcher").setLevel(logging.DEBUG)
 logging.getLogger("micall.core.plot_contigs").setLevel(logging.DEBUG)
 
 
@@ -39,7 +39,7 @@
 
 @pytest.fixture()
 def exact_aligner(monkeypatch):
-    monkeypatch.setattr("micall.core.contig_stitcher.align_consensus", mock_align_consensus)
+    monkeypatch.setattr("micall.utils.referencefull_contig_stitcher.align_consensus", mock_align_consensus)
 
 
 @pytest.fixture
@@ -875,7 +875,7 @@ def test_stitching_contig_with_small_covered_gap(exact_aligner, visualizer):
     assert len(visualizer().elements) > len(contigs)
 
     assert all(x.seq == lstrip(rstrip(x)).seq for x in results)
-    assert {contig.seq for contig in contigs} == {contig.seq for contig in results}
+    assert {contig.seq for contig in contigs} != {contig.seq for contig in results}
 
 
 def test_stitching_partial_align(exact_aligner, visualizer):
@@ -1395,7 +1395,7 @@ def mock_align(reference_seq: str, consensus: str) -> Tuple[List[MockAlignment],
         algorithm = 'mock'
         return (alignments, algorithm)
 
-    monkeypatch.setattr("micall.core.contig_stitcher.align_consensus", mock_align)
+    monkeypatch.setattr("micall.utils.referencefull_contig_stitcher.align_consensus", mock_align)
 
     ref = 'A' * 700
     seq = 'C' * 600
@@ -1457,10 +1457,11 @@ def test_correct_stitching_of_one_normal_and_one_unknown(exact_aligner, visualiz
 
 
 def test_main_invocation(exact_aligner, tmp_path, hcv_db):
+    from micall.core.contig_stitcher import main
     pwd = os.path.dirname(__file__)
     contigs = os.path.join(pwd, "data", "exact_parts_contigs.csv")
     stitched_contigs = os.path.join(tmp_path, "stitched.csv")
-    stitcher.main([contigs, stitched_contigs])
+    main(['with-references', contigs, stitched_contigs])
 
     assert os.path.exists(contigs)
     assert os.path.exists(stitched_contigs)
@@ -1479,11 +1480,12 @@ def test_main_invocation(exact_aligner, tmp_path, hcv_db):
 
 
 def test_visualizer_simple(exact_aligner, tmp_path, hcv_db):
+    from micall.core.contig_stitcher import main
     pwd = os.path.dirname(__file__)
     contigs = os.path.join(pwd, "data", "exact_parts_contigs.csv")
     stitched_contigs = os.path.join(tmp_path, "stitched.csv")
     plot = os.path.join(tmp_path, "exact_parts_contigs.plot.svg")
-    stitcher.main([contigs, stitched_contigs, "--debug", "--plot", plot])
+    main(['with-references', contigs, stitched_contigs, "--debug", "--plot", plot])
 
     assert os.path.exists(contigs)
     assert os.path.exists(stitched_contigs)
@@ -1626,17 +1628,16 @@ def test_merge_intervals(intervals, expected):
     assert merge_intervals(intervals) == expected
 
 
-@dataclass
-class TestMockAlignment:
-    r_st: int
-    r_ei: int
-
-
 class MockAlignedContig:
+    @dataclass
+    class TestMockAlignment:
+        r_st: int
+        r_ei: int
+
     def __init__(self, ref_name, group_ref, r_st, r_ei, name="contig"):
         self.ref_name = ref_name
         self.group_ref = group_ref
-        self.alignment = TestMockAlignment(r_st, r_ei)
+        self.alignment = MockAlignedContig.TestMockAlignment(r_st, r_ei)
         self.name = name
         self.id = id(self)
 
@@ -1744,12 +1745,12 @@ def test_find_covered(contigs, expected_covered_name):
 
 def test_concordance_same_length_inputs():
     with pytest.raises(ValueError):
-        calculate_concordance("abc", "ab")
+        calculate_concordance_norm("abc", "ab")
 
 
 def test_concordance_completely_different_strings():
-    result = calculate_concordance("a" * 30, "b" * 30)
-    assert all(n == 0 for n in result)
+    result = calculate_concordance_norm("a" * 30, "b" * 30)
+    assert any(n > 0 for n in result)
 
 
 def generate_random_string_pair(length):
@@ -1761,19 +1762,21 @@ def generate_random_string_pair(length):
 @pytest.mark.parametrize(
     "left, right, expected",
     [
-        ("aaaaa", "aaaaa", [0.6, 0.68, 0.7, 0.68, 0.6]),
-        ("abcdd", "abcdd", [0.6, 0.68, 0.7, 0.68, 0.6]),
-        ("aaaaaaaa", "baaaaaab", [0.3, 0.62, 0.71, 0.75, 0.75, 0.71, 0.62, 0.3]),
-        ("aaaaaaaa", "aaaaaaab", [0.64, 0.73, 0.79, 0.8, 0.79, 0.73, 0.64, 0.31]),
-        ("aaaaaaaa", "aaaaaaab", [0.64, 0.73, 0.79, 0.8, 0.79, 0.73, 0.64, 0.31]),
-        ("aaaaaaaa", "aaaaabbb", [0.6, 0.68, 0.7, 0.68, 0.6, 0.29, 0.19, 0.13]),
-        ("aaaaaaaa", "aaabbaaa", [0.56, 0.63, 0.62, 0.39, 0.39, 0.62, 0.63, 0.56]),
-        ("aaaaa", "bbbbb", [0] * 5),
+        # ("aaaaa", "aaaaa", [0.0, 0.78, 1.0, 0.78, 0.0]),
+        # ("abcdd", "abcdd", [0.0, 0.78, 1.0, 0.78, 0.0]),
+        ("aaaaaaaa", "baaaaaab", [0.0, 0.95, 0.99, 1.0, 1.0, 0.99, 0.95, 0.0]),
+        ("aaaaaaaa", "aaaaaaab", [0.94, 0.98, 0.99, 1.0, 0.99, 0.98, 0.94, 0.0]),
+        ("aaaaaaaa", "aaaaabbb", [0.96, 0.99, 1.0, 0.99, 0.96, 0.02, 0.0, 0.02]),
+        ("aaaaaaaaaaa", "aaaaabaaaaa", [0.96, 0.99, 1.0, 0.99, 0.96, 0.0, 0.96, 0.99, 1.0, 0.99, 0.96]),
+        ("aaaaaaaaaaa", "aaaabbbaaaa", [0.98, 1.0, 1.0, 0.98, 0.02, 0.0, 0.02, 0.98, 1.0, 1.0, 0.98]),
+        ("aaaaaaaaaaa", "aaabbbbbaaa", [0.98, 1.0, 0.98, 0.04, 0.01, 0.0, 0.01, 0.04, 0.98, 1.0, 0.98]),
+        ("aaaaaaaa", "aaabbaaa", [0.98, 1.0, 0.98, 0.0, 0.0, 0.98, 1.0, 0.98]),
+        ("aaaaa", "bbbbb", [1.0, 0.22, 0.0, 0.22, 1.0]),
         ("", "", []),
     ],
 )
 def test_concordance_simple(left, right, expected):
-    result = [round(float(x), 2) for x in calculate_concordance(left, right)]
+    result = [round(float(x), 2) for x in calculate_concordance_norm(left, right)]
     assert result == expected
 
 
@@ -1790,8 +1793,8 @@ def test_concordance_simple(left, right, expected):
         (
             "a" * 128,
             "b" * 48 + "a" * 15 + "ab" + "a" * 15 + "b" * 48,
-            48 + 16 // 2,
-        ),  # two peaks - choosing 1nd
+            48 + 15 // 2,
+        ),  # two peaks - choosing 1st
         (
             "a" * 128,
             "b" * 48 + "a" * 15 + "ba" + "a" * 15 + "b" * 48,
@@ -1805,7 +1808,7 @@ def test_concordance_simple(left, right, expected):
     ],
 )
 def test_concordance_simple_index(left, right, expected):
-    concordance = calculate_concordance(left, right)
+    concordance = calculate_concordance_norm(left, right)
     concordance_d = list(disambiguate_concordance(concordance))
     index = max(range(len(concordance)), key=lambda i: concordance_d[i])
     if abs(index - expected) > 1:
@@ -1823,7 +1826,7 @@ def generate_test_cases(num_cases):
 
 @pytest.mark.parametrize("left, right", concordance_cases)
 def test_concordance_output_range(left, right):
-    result = calculate_concordance(left, right)
+    result = calculate_concordance_norm(left, right)
     assert all(
         0 <= n <= 1 for n in result
     ), "All values in result should be between 0 and 1"
@@ -1845,8 +1848,8 @@ def test_concordance_higher_if_more_matches_added(left, right):
         + right[insert_position + len(matching_sequence):]
     )
 
-    old_conc = calculate_concordance(left, right)
-    new_conc = calculate_concordance(new_left, new_right)
+    old_conc = calculate_concordance_norm(left, right)
+    new_conc = calculate_concordance_norm(new_left, new_right)
     old_average = sum(old_conc) / len(old_conc)
     new_average = sum(new_conc) / len(new_conc)
     assert old_average <= new_average
@@ -1868,7 +1871,7 @@ def test_concordance_higher_in_matching_areas(left, right):
         + right[insert_position + len(matching_sequence):]
     )
 
-    concordance_scores = calculate_concordance(new_left, new_right)
+    concordance_scores = calculate_concordance_norm(new_left, new_right)
 
     # Check concordance in the matching area
     matching_area_concordance = concordance_scores[

micall/tests/test_contig_stitcher_fuzz.py

@@ -1,15 +1,15 @@
 import pytest
 import json
 import os
-from micall.core.contig_stitcher import (
+from micall.utils.referencefull_contig_stitcher import (
     GenotypedContig,
     AlignedContig,
     stitch_consensus,
     stitch_contigs,
     drop_completely_covered,
     StitcherContext,
 )
-import micall.core.contig_stitcher as stitcher
+import micall.utils.referencefull_contig_stitcher as stitcher
 from micall.core.plot_contigs import build_stitcher_figure
 from aligntools import CigarHit, Cigar, CigarActions
 from typing import Dict, List
@@ -18,7 +18,7 @@
 
 @pytest.fixture
 def no_aligner(monkeypatch):
-    monkeypatch.setattr("micall.core.contig_stitcher.align_to_reference", lambda x: [x])
+    monkeypatch.setattr("micall.utils.referencefull_contig_stitcher.align_to_reference", lambda x: [x])
 
 
 @pytest.fixture(autouse=True)

micall/tests/test_csv_to_fasta.py

@@ -0,0 +1,64 @@
+import sys
+import pytest
+from pathlib import Path
+from micall.utils.csv_to_fasta import cli, main
+
+
+def test_normal_input(tmpdir):
+    contigs_fasta = Path(tmpdir) / "contigs.fasta"
+    contigs_csv = Path(tmpdir) / "contigs.csv"
+
+    contigs_csv.write_text("""\
+ref,match,group_ref,contig
+HCV-1a,1.0,HCV-1a,TCACCAGGACAGCGGGTTGAATTCCTCGTGCAAGCGTGGAA
+HCV-1a,0.75,HCV-1a,CATCACATAGGAGACAGGGCTCCAGGACTGCACCATGCTCGTGTGTGGCGACGAC
+""")
+
+    expected_contigs_fasta = """\
+>1
+TCACCAGGACAGCGGGTTGAATTCCTCGTGCAAGCGTGGAA
+>2
+CATCACATAGGAGACAGGGCTCCAGGACTGCACCATGCTCGTGTGTGGCGACGAC
+"""
+
+    assert main([str(contigs_csv), str(contigs_fasta)]) == 0
+    assert expected_contigs_fasta == contigs_fasta.read_text()
+
+
+def test_empty(tmpdir):
+    contigs_fasta = Path(tmpdir) / "contigs.fasta"
+    contigs_csv = Path(tmpdir) / "contigs.csv"
+
+    contigs_csv.write_text("""\
+ref,match,group_ref,contig
+""")
+
+    expected_contigs_fasta = ""
+    assert main([str(contigs_csv), str(contigs_fasta)]) == 0
+    assert expected_contigs_fasta == contigs_fasta.read_text()
+
+
+def test_main_invocation(tmpdir, monkeypatch):
+    contigs_fasta = Path(tmpdir) / "contigs.fasta"
+    contigs_csv = Path(tmpdir) / "contigs.csv"
+
+    contigs_csv.write_text("""\
+ref,match,group_ref,contig
+HCV-1a,1.0,HCV-1a,TCACCAGGACAGCGGGTTGAATTCCTCGTGCAAGCGTGGAA
+HCV-1a,0.75,HCV-1a,CATCACATAGGAGACAGGGCTCCAGGACTGCACCATGCTCGTGTGTGGCGACGAC
+""")
+
+    expected_contigs_fasta = """\
+>1
+TCACCAGGACAGCGGGTTGAATTCCTCGTGCAAGCGTGGAA
+>2
+CATCACATAGGAGACAGGGCTCCAGGACTGCACCATGCTCGTGTGTGGCGACGAC
+"""
+
+    argv = ['program', str(contigs_csv), str(contigs_fasta)]
+    monkeypatch.setattr(sys, 'argv', argv)
+
+    with pytest.raises(SystemExit):
+        cli()
+
+    assert expected_contigs_fasta == contigs_fasta.read_text()