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
@@ -22,8 +22,8 @@
 from micall.core.project_config import ProjectConfig
 from micall.utils.alignment_wrapper import align_nucs
 from micall.utils.contig_stitcher_contigs import Contig, GenotypedContig, AlignedContig
-from micall.utils.contig_stitcher_context import StitcherContext
-import micall.utils.contig_stitcher_events as events
+from micall.utils.contig_stitcher_context import ReferencefullStitcherContext
+import micall.utils.referencefull_contig_stitcher_events as events
 from micall.data.landmark_reader import LandmarkReader
 
 
@@ -401,7 +401,7 @@ def build_coverage_figure(genome_coverage_csv, blast_csv=None, use_concordance=F
 
 
 def plot_stitcher_coverage(logs: Iterable[events.EventType], genome_coverage_svg_path: str):
-    with StitcherContext.stage():
+    with ReferencefullStitcherContext.stage():
         f = build_stitcher_figure(logs)
         f.show(w=970).save_svg(genome_coverage_svg_path, context=draw.Context(invert_y=True))
         return f
@@ -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,9 +88,14 @@
     "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/utils/fasta_to_fastq.py",
     "micall/utils/append_primers.py",
     "micall/utils/randomize_fastq.py",
+    "micall/utils/calculate_kmer_frequencies.py",
+    "micall/utils/fastq_to_fasta.py",
     "micall/utils/analyze_kive_batches/analyze_kive_batches.py",
 ]
 
@@ -120,10 +125,10 @@ def execute_module_as_main(module_name: str, arguments: Sequence[str]) -> int:
 
 
 def get_version() -> str:
-    if __package__ is None:
-        return "development"
-    else:
+    try:
         return str(version(__package__))
+    except BaseException:
+        return "development"
 
 
 def get_parser() -> argparse.ArgumentParser: