class AlignmentScoringScheme(object):
    """
    This class holds an alignment scoring scheme.
    """

    def __init__(self, scheme_string):
        scheme_parts = scheme_string.split(',')

        # Default scoring scheme
        self.match = 3
        self.mismatch = -6
        self.gap_open = -5
        self.gap_extend = -2

        if len(scheme_parts) == 4:
            self.match = int(scheme_parts[0])
            self.mismatch = int(scheme_parts[1])
            self.gap_open = int(scheme_parts[2])
            self.gap_extend = int(scheme_parts[3])

    def __repr__(self):
        return str(self.match) + ',' + str(self.mismatch) + ',' + str(self.gap_open) + ',' + \
               str(self.gap_extend) class AssemblyGraph(object):
    """
    This class holds an assembly graph in GFA with segments and links (and contig paths too if the
    graph was made by SPAdes).
    """

    def __init__(self, filename, overlap, insert_size_mean=250, insert_size_deviation=50):
        self.segments = {}  # Dict of unsigned segment number -> segment
        self.forward_links = {}  # Dict of signed segment number -> list of signed segment numbers
        self.reverse_links = {}  # Dict of signed segment number <- list of signed segment numbers
        self.copy_depths = {}  # Dict of unsigned segment number -> list of copy depths
        self.manual_multiplicity = {}  # Dict of unsigned segment number -> multiplicity
        self.paths = {}  # Dict of path name -> list of signed segment numbers
        self.overlap = overlap
        self.insert_size_mean = insert_size_mean
        self.insert_size_deviation = insert_size_deviation

        self.load_from_gfa(filename)
        if not overlap:
            self.overlap = get_overlap_from_gfa_link(filename)

    def load_from_gfa(self, filename):
        """
        Loads a Graph from a GFA file. It does not load any GFA file, but makes some restrictions:
        1) The segment names must be integers.
        2) The depths should be stored in a dp tag.
        3) All link overlaps are the same (equal to the graph overlap value).
        """
        # Load in the segments.
        with open(filename, 'rt') as gfa_file:
            for line in gfa_file:
                if line.startswith('S'):
                    line_parts = line.strip().split('\t')
                    num = int(line_parts[1])
                    depth = 1.0
                    for part in line_parts:
                        if part.lower().startswith('dp:'):
                            depth = float(part[5:])
                        if part.lower().startswith('ml:'):
                            self.manual_multiplicity[num] = int(part[5:])
                    sequence = line_parts[2]
                    self.segments[num] = Segment(num, depth, sequence, True)
                    self.segments[num].build_other_sequence_if_necessary()
                if line.startswith('i'):
                    line_parts = line.strip().split('\t')
                    try:
                        self.insert_size_mean = float(line_parts[1])
                        self.insert_size_deviation = float(line_parts[2])
                    except ValueError:
                        pass

        # Load in the links.
        with open(filename, 'rt') as gfa_file:
            for line in gfa_file:
                if line.startswith('L'):
                    line_parts = line.strip().split('\t')
                    start = signed_string_to_int(line_parts[1] + line_parts[2])
                    end = signed_string_to_int(line_parts[3] + line_parts[4])
                    if start not in self.forward_links:
                        self.forward_links[start] = [end]
                    else:
                        self.forward_links[start].append(end)
            self.forward_links = build_rc_links_if_necessary(self.forward_links)
            self.reverse_links = build_reverse_links(self.forward_links)
        self.sort_link_order()

        # Load in the paths.
        with open(filename, 'rt') as gfa_file:
            for line in gfa_file:
                if line.startswith('P'):
                    line_parts = line.strip().split('\t')
                    path_name = line_parts[1]
                    segments = [signed_string_to_int(x) for x in line_parts[2].split(',')]
                    if len(segments) > 1:
                        self.paths[path_name] = segments

    def get_median_read_depth(self, segment_list=None):
        """
        Returns the assembly graph's median read depth (by base).  Optionally, a list of segments
        can be given, in which case only those segments are used for the calculation.
        """
        if not segment_list:
            segment_list = self.segments.values()
        sorted_segments = sorted(segment_list, key=lambda x: x.depth)
        total_length = 0
        for segment in sorted_segments:
            total_length += segment.get_length_no_overlap(self.overlap)
        halfway_length = total_length // 2
        length_so_far = 0
        for segment in sorted_segments:
            length_so_far += segment.get_length_no_overlap(self.overlap)
            if length_so_far >= halfway_length:
                return segment.depth
        return 0.0

    def get_single_copy_depth(self):
        """
        Determines the single copy read depth for the graph. It uses the median depth (by base)
        using the 10 longest segments in the graph.
        """
        ten_longest_contigs = sorted(self.segments.values(), reverse=True,
                                     key=lambda x: x.get_length())[:10]
        median_depth = self.get_median_read_depth(ten_longest_contigs)
        log.log('Median depth of 10 longest contigs: ' + float_to_str(median_depth, 2), 2)
        log.log('', 2)
        return median_depth

    def reassign_read_depths(self):
        """
        This function looks for segments which have an unoriginal read depth. If they are connected
        to segments with original read depths, these neighbours will be used to reassign depth.
        """
        while True:
            for seg_num, segment in self.segments.items():
                if not segment.original_depth:

                    new_depth_downstream = None
                    downstream_seg_nums = self.get_exclusive_outputs(seg_num)
                    if downstream_seg_nums:
                        downstream_segs = [self.segments[abs(x)] for x in downstream_seg_nums]
                        if all(x.original_depth for x in downstream_segs):
                            new_depth_downstream = sum(x.depth for x in downstream_segs)

                    new_depth_upstream = None
                    upstream_seg_nums = self.get_exclusive_inputs(seg_num)
                    if upstream_seg_nums:
                        upstream_segs = [self.segments[abs(x)] for x in upstream_seg_nums]
                        if all(x.original_depth for x in upstream_segs):
                            new_depth_upstream = sum(x.depth for x in upstream_segs)

                    # If both an upstream and downstream depth is available, use the mean of the
                    # two. If only one is available, just use that.
                    if new_depth_downstream and new_depth_upstream:
                        new_depth = (new_depth_downstream + new_depth_upstream) / 2.0
                    elif new_depth_downstream:
                        new_depth = new_depth_downstream
                    elif new_depth_upstream:
                        new_depth = new_depth_upstream
                    else:
                        new_depth = None
                    if new_depth:
                        segment.depth = new_depth
                        segment.original_depth = True
                        break
            else:
                break

    def normalise_read_depths(self):
        """
        For every segment in the graph, divide its depth by the median depth of the 10 longest
        contigs. This makes segments with the median depth have a depth of 1, segments with more
        than the median a depth of greater than 1 and segments with less than the median a depth of
        less than 1.
        """
        ten_longest_contigs = sorted(self.segments.values(), reverse=True,
                                     key=lambda x: x.get_length())[:10]
        median_depth = self.get_median_read_depth(ten_longest_contigs)
        if median_depth == 0.0:
            return
        for segment in self.segments.values():
            segment.depth /= median_depth

    def get_total_length(self):
        """
        Returns the sum of all segment sequence lengths.
        """
        return sum([x.get_length() for x in self.segments.values()])

    def get_total_length_no_overlaps(self):
        """
        Returns the sum of all segment sequence lengths, subtracting the overlap size from each
        segment.
        """
        return sum([x.get_length_no_overlap(self.overlap) for x in self.segments.values()])

    def total_dead_end_count(self):
        """
        Returns the total number of dead ends in the assembly graph.
        """
        dead_ends = 0
        for seg_num in self.segments:
            dead_ends += self.dead_end_count(seg_num)
        return dead_ends

    def dead_end_count(self, seg_num):
        """
        Returns the number of dead ends for one segment: 0, 1 or 2.
        """
        dead_ends = 0
        if seg_num not in self.forward_links or not self.forward_links[seg_num]:
            dead_ends += 1
        if seg_num not in self.reverse_links or not self.reverse_links[seg_num]:
            dead_ends += 1
        return dead_ends

    def save_to_fasta(self, filename, newline=False, min_length=1, verbosity=1, silent=False):
        """
        Saves whole graph (only forward sequences) to a FASTA file.
        """
        if not silent:
            log.log(('\n' if newline else '') + 'Saving ' + filename, verbosity)
        circular_seg_nums = self.completed_circular_replicons()
        with open(filename, 'w') as fasta:
            sorted_segments = sorted(self.segments.values(), key=lambda x: x.number)
            for segment in sorted_segments:
                if segment.get_length() >= min_length:
                    fasta.write(segment.get_fasta_name_and_description_line(circular_seg_nums))
                    fasta.write(add_line_breaks_to_sequence(segment.forward_sequence))

    @staticmethod
    def save_specific_segments_to_fasta(filename, segments, silent=False):
        """
        Saves single copy segments (only forward sequences) to a FASTA file.
        """
        if not silent:
            log.log('Saving ' + filename)
        with open(filename, 'w') as fasta:
            sorted_segments = sorted(segments, key=lambda x: x.number)
            for segment in sorted_segments:
                fasta.write('>' + str(segment.number) + '\n')
                fasta.write(add_line_breaks_to_sequence(segment.forward_sequence))

    def save_to_gfa(self, filename, verbosity=1, save_copy_depth_info=False,
                    save_seg_type_info=False, newline=False, include_insert_size=False):
        """
        Saves whole graph to a GFA file.
        """
        log.log(('\n' if newline else '') + 'Saving ' + filename, verbosity)
        with open(filename, 'w') as gfa:
            sorted_segments = sorted(self.segments.values(), key=lambda x: x.number)
            for segment in sorted_segments:
                segment_line = segment.gfa_segment_line()
                segment_colour, label = '', ''
                if save_copy_depth_info and segment.number in self.copy_depths:
                    segment_colour = self.get_copy_number_colour(segment)
                    label = self.get_depth_string(segment)
                if save_seg_type_info and segment.bridge is not None:
                    segment_colour = 'pink'
                    label = segment.get_seg_type_label()
                if segment_colour or label:
                    segment_line = segment_line[:-1]  # Remove newline
                    segment_line += '\tLB:z:' + label.replace('\n', '\\n')
                    segment_line += '\tCL:z:' + segment_colour
                    segment_line += '\n'
                gfa.write(segment_line)
            gfa.write(self.get_all_gfa_link_lines())
            paths = sorted(self.paths.items())
            overlap_cigar = str(self.overlap) + 'M'
            for path_name, segment_list in paths:
                gfa.write('P\t' + path_name + '\t')
                gfa.write(','.join([int_to_signed_string(x) for x in segment_list]))
                gfa.write('\t')
                gfa.write(','.join([overlap_cigar] * (len(segment_list) - 1)))
                gfa.write('\n')
            if include_insert_size and self.insert_size_mean is not None and \
                    self.insert_size_deviation is not None:
                gfa.write('i\t')
                gfa.write(str(self.insert_size_mean))
                gfa.write('\t')
                gfa.write(str(self.insert_size_deviation))
                gfa.write('\n')

    def get_all_gfa_link_lines(self):
        """
        Returns a string of the link component of the GFA file for this graph.
        """
        gfa_link_lines = []
        for start, ends in self.forward_links.items():
            for end in ends:
                if is_link_positive(start, end):
                    gfa_link_lines.append(self.gfa_link_line(start, end))
        return ''.join(gfa_link_lines)

    def filter_by_read_depth(self, relative_depth_cutoff):
        """
        This function removes segments from the graph based on a relative depth cutoff. Segments
        are considered below the cutoff if they are less than the cutoff for the entire graph or
        less than the cutoff for their connected component.
        To be removed, one of the following must also be true:
          1) the segment has at least one dead end
          2) the segment is part of a connected component where all of the segments are below the
             whole graph cutoff
          3) deleting the segment would not create any dead ends
        """
        segment_nums_to_remove = []
        total_length_removed = 0
        ten_longest_contigs = sorted(self.segments.values(), reverse=True,
                                     key=lambda x: x.get_length())[:10]
        whole_graph_cutoff = self.get_median_read_depth(ten_longest_contigs) * relative_depth_cutoff
        connected_components = self.get_connected_components()
        for component in connected_components:
            component_segs = [self.segments[x] for x in component]
            component_cutoff = self.get_median_read_depth(component_segs) * relative_depth_cutoff
            for seg_num in component:
                segment = self.segments[seg_num]
                if segment.depth < whole_graph_cutoff or segment.depth < component_cutoff:
                    if self.dead_end_count(seg_num) > 0 or \
                            self.all_segments_below_depth(component, whole_graph_cutoff) or \
                            self.dead_end_change_if_deleted(seg_num) <= 0:
                        segment_nums_to_remove.append(seg_num)
                        total_length_removed += segment.get_length()
        self.remove_segments(segment_nums_to_remove)
        return len(segment_nums_to_remove), total_length_removed

    def filter_homopolymer_loops(self):
        """
        A common feature in SPAdes graphs is a small piece of the graph (often just one segment)
        which has nothing but one base.  Filter these out.
        """
        segment_nums_to_remove = []
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            component_segments = [self.segments[x] for x in component_nums]
            if all_segments_are_one_base(component_segments):
                segment_nums_to_remove += component_nums
        self.remove_segments(segment_nums_to_remove)
        if segment_nums_to_remove:
            log.log('Removed homopolymer loops:', 3)
            log.log_number_list(segment_nums_to_remove, 3)

    def choose_largest_component(self):
        """
        Special logic: throw out all of the graph's connected components except for the largest one.
        """
        largest_component_length = None
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            component_segments = [self.segments[x] for x in component_nums]
            component_length = sum(x.get_length() for x in component_segments)
            if largest_component_length is None or component_length > largest_component_length:
                largest_component_length = component_length
        segment_nums_to_remove = []
        for component_nums in connected_components:
            component_segments = [self.segments[x] for x in component_nums]
            component_length = sum(x.get_length() for x in component_segments)
            if component_length < largest_component_length:
                segment_nums_to_remove += component_nums
        self.remove_segments(segment_nums_to_remove)
        if segment_nums_to_remove:
            log.log('\nRemoved not-largest components:', 3)
            log.log_number_list(segment_nums_to_remove, 3)

    def remove_segments(self, nums_to_remove):
        """
        This function deletes all segments in the nums_to_remove list, along with their links. It
        also deletes any paths which contain those segments.
        """
        for num_to_remove in nums_to_remove:
            assert num_to_remove >= 0  # this function takes positive segment numbers only
            if num_to_remove in self.segments:
                seg_to_remove = self.segments[num_to_remove]

                # If the segment being removed is a bridge, and if that bridge's application took
                # depth away from other segments, then we now need to give that depth back.
                if seg_to_remove.bridge and seg_to_remove.bridge.segments_reduced_depth:
                    for num, depth, copy_depth in seg_to_remove.bridge.segments_reduced_depth:
                        if num in self.segments:
                            restore_depth_seg = self.segments[num]
                            restore_depth_seg.depth += depth
                            if copy_depth and num in self.copy_depths:
                                self.copy_depths[num].append(copy_depth)
                # Now actually delete the segment.
                del self.segments[num_to_remove]

        # Delete the copy depths for deleted segments.
        for num in nums_to_remove:
            if num in self.copy_depths:
                del self.copy_depths[num]

        # Remove links for the deleted segments.
        links_to_remove = set()
        for num_to_remove in nums_to_remove:
            for down_seg in self.get_downstream_seg_nums(num_to_remove):
                links_to_remove.add((num_to_remove, down_seg))
            for up_seg in self.get_upstream_seg_nums(num_to_remove):
                links_to_remove.add((up_seg, num_to_remove))
        for link in links_to_remove:
            self.remove_link(link[0], link[1])

        self.remove_segments_from_paths(nums_to_remove)

    def remove_segments_from_paths(self, seg_nums):
        """
        Deletes the given segment numbers (regardless of sign) from the paths. If this results in
        an invalid path, then the whole path is deleted.
        """
        fixed_paths = {}
        for path_name, path in self.paths.items():
            fixed_path = [x for x in path if x not in seg_nums and -x not in seg_nums]
            if len(fixed_path) > 1 and self.is_path_valid(fixed_path):
                fixed_paths[path_name] = fixed_path
        self.paths = fixed_paths

    def is_path_valid(self, path):
        """
        Returns True/False based on whether or not the given path is in the graph.
        """
        for i, seg_2 in enumerate(path):
            if i > 0:
                seg_1 = path[i - 1]
                if seg_1 not in self.forward_links or seg_2 not in self.forward_links[seg_1]:
                    return False
        return True

    def remove_small_components(self, min_component_size):
        """
        Remove small graph components, but only if they do not contain any bridges. The idea is
        to clean up parts of the graph that were orphaned by the bridging process. But if they
        contain a bridge, then they are more likely to be genuine and we keep them.
        """
        segment_nums_to_remove = []
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            component_segments = [self.segments[x] for x in component_nums]
            component_length = sum(x.get_length() for x in component_segments)
            if component_length >= min_component_size:
                continue
            if any(x.bridge is not None for x in component_segments):
                continue
            segment_nums_to_remove += component_nums
        self.remove_segments(segment_nums_to_remove)
        if segment_nums_to_remove:
            log.log('\nRemoved small components:', 2)
            log.log_number_list(segment_nums_to_remove, 2)

    def remove_small_dead_ends(self, min_dead_end_size):
        """
        Remove small segments which are graph dead-ends. This is just to tidy things up a bit
        before the final merge.
        """
        removed_segments = []
        while True:
            for seg_num, segment in self.segments.items():
                if segment.get_length() >= min_dead_end_size:
                    continue
                if self.dead_end_change_if_deleted(seg_num) < 0:
                    self.remove_segments([seg_num])
                    removed_segments.append(seg_num)
                    break
            else:
                break
        if removed_segments:
            log.log('\nRemoved small dead ends:', 2)
            log.log_number_list(removed_segments, 2)

    def merge_all_possible(self, anchor_segments, bridging_mode):
        """
        This function merges segments which are in a simple, unbranching path.
        """
        if anchor_segments is not None:
            anchor_seg_nums = set(x.number for x in anchor_segments)
        else:
            anchor_seg_nums = None
        while True:
            # Sort the segment numbers first so we apply the merging in a consistent order.
            seg_nums = sorted(list(self.segments.keys()))
            for num in seg_nums:
                path = self.get_simple_path(num, anchor_seg_nums, bridging_mode)
                assert len(path) > 0
                if len(path) > 1:
                    self.merge_simple_path(path)
                    break
            else:
                break
        self.renumber_segments()

    def merge_simple_path(self, merge_path):
        """
        Merges the path into a single segment and adjusts any graph paths as necessary.
        """
        start = merge_path[0]
        end = merge_path[-1]
        mean_depth, original_depth = self.get_mean_path_depth(merge_path)

        # Make sure this is indeed a simple path.
        for i in range(len(merge_path) - 1):
            s_1 = merge_path[i]
            s_2 = merge_path[i+1]
            if [s_2] != self.forward_links[s_1]:
                raise BadPath(str(merge_path) + ' is not a simple path')

        new_seg_num = self.get_next_available_seg_number()
        merged_forward_seq = self.get_path_sequence(merge_path)
        new_seg = Segment(new_seg_num, mean_depth, merged_forward_seq, True,
                          original_depth=original_depth)
        new_seg.build_other_sequence_if_necessary()

        # Save some info that we'll need, and then delete the old segments.
        paths_copy = self.paths.copy()
        outgoing_links = []
        if end in self.forward_links:
            outgoing_links = list(self.forward_links[end])
        incoming_links = []
        if start in self.reverse_links:
            incoming_links = list(self.reverse_links[start])
        outgoing_links = find_replace_one_val_in_list(outgoing_links, start, new_seg_num)
        outgoing_links = find_replace_one_val_in_list(outgoing_links, -end, -new_seg_num)
        incoming_links = find_replace_one_val_in_list(incoming_links, end, new_seg_num)
        incoming_links = find_replace_one_val_in_list(incoming_links, -start, -new_seg_num)

        self.remove_segments([abs(x) for x in merge_path])

        # Add the new segment to the graph and give it the links from its source segments.
        self.segments[new_seg_num] = new_seg
        for link in outgoing_links:
            self.add_link(new_seg_num, link)
        for link in incoming_links:
            self.add_link(link, new_seg_num)

        # Merge the segments in any paths.
        flipped_merge_path = [-x for x in reversed(merge_path)]
        for path_name in paths_copy:
            paths_copy[path_name] = find_replace_in_list(paths_copy[path_name], merge_path,
                                                         [new_seg_num])
            paths_copy[path_name] = find_replace_in_list(paths_copy[path_name], flipped_merge_path,
                                                         [-new_seg_num])

        # If any paths still contain the original segments, then split those paths into pieces,
        # removing the original segments.
        new_paths = {}
        for path_name, path_segments in paths_copy.items():
            split_paths = split_path_multiple(path_segments, merge_path + flipped_merge_path)
            if len(split_paths) == 1:
                new_paths[path_name] = split_paths[0]
            elif len(split_paths) > 1:
                for i, path in enumerate(split_paths):
                    new_paths[path_name + '_' + str(i + 1)] = path
        self.paths = new_paths

        return new_seg_num

    def get_mean_path_depth(self, path):
        """
        Returns the mean depth for the path. If any segments in the path are bridges, their depth
        isn't counted because bridges got their depth from the segments they are bridging, so to
        count them would be to count that depth twice.
        """
        non_bridge_seg_nums = [abs(x) for x in path if self.segments[abs(x)].bridge is None]

        # If possible, we'd like to only use the depth from segments which haven't had their depth
        # altered by being used in bridges. But if none are available (i.e. all segments have been
        # used in bridges), then we go ahead and use them anyway.
        original_depth_seg_nums = [x for x in non_bridge_seg_nums
                                   if self.segments[x].original_depth]
        if original_depth_seg_nums:
            segs_nums_for_depth = original_depth_seg_nums
            original_depth = True
        else:
            segs_nums_for_depth = non_bridge_seg_nums
            original_depth = False

        depths = [self.segments[x].depth for x in segs_nums_for_depth]
        lengths = [self.segments[x].get_length() - self.overlap for x in segs_nums_for_depth]
        if sum(lengths) > 0.0:
            new_depth = weighted_average_list(depths, lengths)
        else:
            new_depth = 1.0
        return new_depth, original_depth

    def add_link(self, start, end):
        """
        Adds a link to the graph in all necessary ways: forward and reverse, and for reverse
        complements too.
        """
        if start not in self.forward_links:
            self.forward_links[start] = []
        if end not in self.forward_links[start]:
            self.forward_links[start].append(end)

        if end not in self.reverse_links:
            self.reverse_links[end] = []
        if start not in self.reverse_links[end]:
            self.reverse_links[end].append(start)

        if -start not in self.reverse_links:
            self.reverse_links[-start] = []
        if -end not in self.reverse_links[-start]:
            self.reverse_links[-start].append(-end)

        if -end not in self.forward_links:
            self.forward_links[-end] = []
        if -start not in self.forward_links[-end]:
            self.forward_links[-end].append(-start)

    def remove_link(self, start, end):
        """
        Removes a link from the graph in all necessary ways: forward and reverse, and for reverse
        complements too.
        """
        if start in self.forward_links:
            try:
                self.forward_links[start].remove(end)
            except ValueError:
                pass
            if len(self.forward_links[start]) == 0:
                del self.forward_links[start]
        if -end in self.forward_links:
            try:
                self.forward_links[-end].remove(-start)
            except ValueError:
                pass
            if len(self.forward_links[-end]) == 0:
                del self.forward_links[-end]
        if end in self.reverse_links:
            try:
                self.reverse_links[end].remove(start)
            except ValueError:
                pass
            if len(self.reverse_links[end]) == 0:
                del self.reverse_links[end]
        if -start in self.reverse_links:
            try:
                self.reverse_links[-start].remove(-end)
            except ValueError:
                pass
            if len(self.reverse_links[-start]) == 0:
                del self.reverse_links[-start]

    def seq_from_signed_seg_num(self, signed_num):
        """
        Returns the forwards or reverse sequence of a segment, if the number is next_positive or
        negative, respectively. Assumes the segment number is in the graph.
        """
        if signed_num > 0:
            return self.segments[signed_num].forward_sequence
        else:
            return self.segments[-signed_num].reverse_sequence

    def get_connected_components(self):
        """
        Returns a list of lists, where each inner list is the segment numbers of one connected
        component of the graph.
        E.g. [[1, 2], [3, 4, 5]] would mean that segments 1 and 2 are in a connected component
        and segments 3, 4 and 5 are in another connected component.
        """
        visited = set()
        components = []
        for v in self.segments:
            if v not in visited:
                component = []
                q = deque()
                q.append(v)
                visited.add(v)
                while q:
                    w = q.popleft()
                    component.append(w)
                    connected_segments = self.get_connected_segments(w)
                    for k in connected_segments:
                        if k not in visited:
                            visited.add(k)
                            q.append(k)
                components.append(sorted(component))

        # Sort (just for consistency from one run to the next)
        return sorted(components)

    def get_connected_segments(self, segment_num):
        """
        Given a segment number, this function returns a list of all other segment numbers for
        segments that are directly connected.
        It only returns positive numbers (i.e. is not strand-specific).
        """
        connected_segments = set()
        if segment_num in self.forward_links:
            downstream_segments = self.forward_links[segment_num]
            connected_segments.update([abs(x) for x in downstream_segments])
        if segment_num in self.reverse_links:
            upstream_segments = self.reverse_links[segment_num]
            connected_segments.update([abs(x) for x in upstream_segments])
        return list(connected_segments)

    def all_segments_below_depth(self, segment_nums, cutoff):
        """
        Returns true if all segments in the list are below the depth cutoff.
        """
        for num in segment_nums:
            if self.segments[num].depth >= cutoff:
                return False
        return True

    def get_n_segment_length(self, n_percent):
        """
        Returns the length for which segments that length and longer make up >= n% of the total
        bases.  E.g. if n = 50, this function returns the N50.  n must be from 0 to 100.
        """
        total_length = self.get_total_length()
        target_length = total_length * (n_percent / 100.0)
        sorted_segments = sorted(self.segments.values(), key=lambda x: x.get_length(), reverse=True)
        length_so_far = 0
        for segment in sorted_segments:
            seg_length = segment.get_length()
            length_so_far += seg_length
            if length_so_far >= target_length:
                return seg_length
        return 0

    def gfa_link_line(self, start, end):
        """
        Returns an entire L line for GFA output, including the newline.
        """
        l_line = 'L\t'
        l_line += str(abs(start)) + '\t'
        l_line += get_sign_string(start) + '\t'
        l_line += str(abs(end)) + '\t'
        l_line += get_sign_string(end) + '\t'
        l_line += str(self.overlap) + 'M\n'
        return l_line

    def get_exclusive_inputs(self, segment_number):
        """
        This function finds all segments which lead into the given segment.  If those segments
        do not lead into any other segments, then this function returns them in a list.
        Specifically, this function returns a list of unsigned numbers.
        """
        if segment_number not in self.reverse_links:
            return []
        return [abs(x) for x in self.reverse_links[segment_number] if
                self.lead_exclusively_to(x, segment_number)]

    def get_exclusive_outputs(self, segment_number):
        """
        Does the same thing as get_exclusive_inputs, but in the other direction.
        """
        if segment_number not in self.forward_links:
            return []
        return [abs(x) for x in self.forward_links[segment_number] if
                self.lead_exclusively_from(x, segment_number)]

    def get_exclusive_inputs_signed(self, segment_number):
        """
        This function finds all segments which lead into the given segment.  If those segments
        do not lead into any other segments, then this function returns them in a list.
        Specifically, this function returns a list of signed numbers.
        """
        if segment_number not in self.reverse_links:
            return []
        return [x for x in self.reverse_links[segment_number] if
                self.lead_exclusively_to(x, segment_number)]

    def get_exclusive_outputs_signed(self, segment_number):
        """
        Does the same thing as get_exclusive_inputs_signed, but in the other direction.
        """
        if segment_number not in self.forward_links:
            return []
        return [x for x in self.forward_links[segment_number] if
                self.lead_exclusively_from(x, segment_number)]

    def lead_exclusively_to(self, segment_num_1, segment_num_2):
        """
        Returns whether or not the first segment leads to and only to the second segment.
        """
        if segment_num_1 not in self.forward_links:
            return False
        return self.forward_links[segment_num_1] == [segment_num_2]

    def lead_exclusively_from(self, segment_num_1, segment_num_2):
        """
        Returns whether the first segment has inputs from and only from the second segment.
        """
        if segment_num_1 not in self.reverse_links:
            return False
        return self.reverse_links[segment_num_1] == [segment_num_2]

    def dead_end_change_if_deleted(self, seg_num):
        """
        Returns the change in graph dead end count if this segment was deleted. 0 means no change,
        positive values mean an increase in dead ends, negative values mean a decrease.
        """
        potential_dead_ends = 0
        for downstream_segment in self.get_downstream_seg_nums(seg_num):
            if len(self.reverse_links[downstream_segment]) == 1 and downstream_segment != seg_num:
                potential_dead_ends += 1

        for upstream_segment in self.get_upstream_seg_nums(seg_num):
            if len(self.forward_links[upstream_segment]) == 1 and upstream_segment != seg_num:
                potential_dead_ends += 1

        return potential_dead_ends - self.dead_end_count(seg_num)

    def dead_end_change_if_path_deleted(self, path_segments):
        """
        Like the above function, but considers the whole path at once. It assumes that the path is
        simple and unbranching (i.e. could be merged into a single segment).
        This function does not check whether the path start and end both connect to the same
        segment. So if they form a hairpin loop, this function will return 0 even though the
        deletion of the path would create a dead end. This behaviour is intentionally left, as it
        helps to clean up such loops from the graph when they have been entirely used in bridges.
        """
        start = path_segments[0]
        end = path_segments[-1]

        potential_dead_ends = 0
        downstream_segments = self.get_downstream_seg_nums(end)
        for downstream_segment in downstream_segments:
            if len(self.reverse_links[downstream_segment]) == 1 and downstream_segment != end:
                potential_dead_ends += 1

        upstream_segments = self.get_upstream_seg_nums(start)
        for upstream_segment in upstream_segments:
            if len(self.forward_links[upstream_segment]) == 1 and upstream_segment != start:
                potential_dead_ends += 1

        dead_ends = 0
        if len(downstream_segments) == 0:
            dead_ends += 1
        if len(upstream_segments) == 0:
            dead_ends += 1
        return potential_dead_ends - dead_ends

    def clean(self, read_depth_filter, largest_component):
        """
        This function does various graph repairs, filters and normalisations to make it a bit
        nicer.
        """
        log.log('Repair multi way junctions  ' + get_dim_timestamp(), 3)
        self.repair_multi_way_junctions()
        log.log('Filter by read depth        ' + get_dim_timestamp(), 3)
        removed_count, removed_length = self.filter_by_read_depth(read_depth_filter)
        log.log('Filter homopolymer loops    ' + get_dim_timestamp(), 3)
        self.filter_homopolymer_loops()
        if largest_component:
            log.log('Keep largest component      ' + get_dim_timestamp(), 3)
            self.choose_largest_component()
        log.log('Merge all possible          ' + get_dim_timestamp(), 3)
        self.merge_all_possible(None, 2)
        log.log('Normalise read depths       ' + get_dim_timestamp(), 3)
        self.normalise_read_depths()
        log.log('Remove zero length segments ' + get_dim_timestamp(), 3)
        self.remove_zero_length_segs(suppress_log=True)
        log.log('Sort link order             ' + get_dim_timestamp(), 3)
        self.sort_link_order()
        log.log('Graph cleaning finished     ' + get_dim_timestamp(), 3)
        return removed_count, removed_length

    def final_clean(self):
        """
        This function cleans up the final assembled graph, in preparation for saving.
        """
        assert self.overlap == 0
        self.remove_zero_length_segs()
        self.merge_small_segments(5)
        self.reassign_read_depths()
        self.normalise_read_depths()
        self.renumber_segments()
        self.sort_link_order()
        self.paths = {}  # Don't need the paths anymore

    def repair_multi_way_junctions(self):
        """
        This function finds and fixes multi-way junctions in the graph, as these can mess up copy
        number determination. It fixes them by creating a new segment with no length (i.e with the
        overlap size) to bridge the connection.
        For example: A->B,C and D->B,C becomes A->E and D->E and E->B and E->C
        """
        seg_nums = list(self.segments) + [-x for x in self.segments]
        already_examined = set()
        for seg_num in seg_nums:
            if seg_num in already_examined:
                continue

            upstream_segs = {seg_num}
            downstream_segs = set()

            # Keep growing the upstream and downstream steps until they stop getting bigger.
            while True:
                upstream_size = len(upstream_segs)
                downstream_size = len(downstream_segs)
                for upstream_seg in upstream_segs:
                    downstream_segs.update(self.get_downstream_seg_nums(upstream_seg))
                for downstream_seg in downstream_segs:
                    upstream_segs.update(self.get_upstream_seg_nums(downstream_seg))
                if len(upstream_segs) == upstream_size and len(downstream_segs) == downstream_size:
                    break
            if len(upstream_segs) < 2 or len(downstream_segs) < 2:
                continue

            # We're looking at this group now so we can skip them in the future.
            already_examined.update(upstream_segs)
            already_examined.update([-s for s in downstream_segs])

            # Skip groups that are too large (otherwise checking each subset will take way too
            # long).
            if len(upstream_segs) > 8:
                continue

            # Examine every subset in the upstream segments of size 2 or larger. This includes the
            # full set.
            starting_seg_groups = []
            ending_seg_groups = []
            used_upstream_subsets = []
            for i in range(len(upstream_segs), 1, -1):
                upstream_subsets = set(itertools.combinations(upstream_segs, i))
                for upstream_subset in upstream_subsets:

                    # No need to examine this subset if we've already made use of a larger subset
                    # which contains it.
                    if any(set(upstream_subset).issubset(x) for x in used_upstream_subsets):
                        continue

                    # Now that we have some upstream segments to look at, see if there are more
                    # than one downstream segments which lead back to all of them.
                    downstream_subset = set()
                    for upstream_seg in upstream_subset:
                        downstream_subset.update(self.get_downstream_seg_nums(upstream_seg))
                    downstream_subset = [x for x in downstream_subset
                                         if all(y in self.get_upstream_seg_nums(x)
                                                for y in upstream_subset)]
                    if len(downstream_subset) < 2:
                        continue

                    # If we got here, then the combination of segments looks good!
                    used_upstream_subsets.append(set(upstream_subset))
                    starting_segs = sorted(list(upstream_subset))
                    ending_segs = sorted(list(downstream_subset))

                    # Sanity check: all upstream segments should lead to all downstream segments.
                    for start_seg in starting_segs:
                        for end_seg in ending_segs:
                            assert end_seg in self.forward_links[start_seg]
                            assert start_seg in self.reverse_links[end_seg]

                    # We'll set it aside to make a small bridge segment to replace the links
                    # between the start and end segments.
                    starting_seg_groups.append(starting_segs)
                    ending_seg_groups.append(ending_segs)

            assert len(starting_seg_groups) == len(ending_seg_groups)
            for starting_segs, ending_segs in zip(starting_seg_groups, ending_seg_groups):
                log.log('Multi-way junction:', 3)
                log.log('  start segs: ' + ', '.join([str(x) for x in starting_segs]), 3)
                log.log('  end segs:   ' + ', '.join([str(x) for x in ending_segs]), 3)

                # Create a new segment to bridge the starting and ending segments.
                bridge_num = self.get_next_available_seg_number()
                start_seg_depth_sum = sum(self.segments[abs(x)].depth for x in starting_segs)
                end_seg_depth_sum = sum(self.segments[abs(x)].depth for x in ending_segs)
                bridge_depth = (start_seg_depth_sum + end_seg_depth_sum) / 2.0
                bridge_seq = self.seq_from_signed_seg_num(ending_segs[0])[:self.overlap]
                bridge_seg = Segment(bridge_num, bridge_depth, bridge_seq, True)
                bridge_seg.build_other_sequence_if_necessary()
                self.segments[bridge_num] = bridge_seg
                log.log('   new seg:   ' + str(bridge_num), 3)

                # Now rebuild the links around the junction.
                for start_seg in starting_segs:
                    for end_seg in ending_segs:
                        self.remove_link(start_seg, end_seg)
                for start_seg in starting_segs:
                    self.add_link(start_seg, bridge_num)
                for end_seg in ending_segs:
                    self.add_link(bridge_num, end_seg)

                # Finally, we need to check to see if there were any paths through the junction.
                # If so, they need to be adjusted to contain the new segment.
                for name in self.paths:
                    for start_num in starting_segs:
                        for end_num in ending_segs:
                            self.paths[name] = insert_num_in_list(self.paths[name], start_num,
                                                                  end_num, bridge_num)
                            self.paths[name] = insert_num_in_list(self.paths[name], -end_num,
                                                                  -start_num, -bridge_num)
                log.log('', 3)

    def remove_unnecessary_links(self):
        """
        This function specifically looks for cases where two segments are directly linked and also
        linked via a 0 bp segment. If so, the direct connection is redundant and is deleted.
        """
        assert self.overlap == 0
        seg_nums = list(self.segments) + [-x for x in self.segments]
        for seg_num in seg_nums:
            down_segs = self.get_downstream_seg_nums(seg_num)
            zero_bp_segs = [x for x in down_segs if self.segments[abs(x)].get_length() == 0]
            for zero_bp_seg in zero_bp_segs:
                down_segs_2 = self.get_downstream_seg_nums(zero_bp_seg)
                common_down_segs = set(down_segs) & set(down_segs_2)
                for common_down_seg in common_down_segs:
                    try:
                        self.remove_link(seg_num, common_down_seg)
                    except ValueError:
                        pass

    def get_next_available_seg_number(self):
        """
        This function finds the largest used segment number and returns the next number.
        """
        return max(self.segments) + 1

    def get_depth_string(self, segment):
        """
        Given a particular segment, this function returns a string with the segment's copy depths
        (if it has any).
        """
        if segment.number not in self.copy_depths:
            return ''
        return '\n'.join(['%.3f' % x for x in self.copy_depths[segment.number]])

    def get_copy_number(self, segment):
        """
        Returns the segment's copy number (0 if copy number determination did not occur for this
        segment).
        """
        if segment.number not in self.copy_depths:
            return 0
        return len(self.copy_depths[segment.number])

    def get_copy_number_colour(self, segment):
        """
        Given a particular segment, this function returns a colour string based on the copy number.
        """
        copy_number = self.get_copy_number(segment)
        if copy_number == 0:
            return 'grey'
        elif copy_number == 1:
            return 'forestgreen'
        elif copy_number == 2:
            return 'gold'
        elif copy_number == 3:
            return 'darkorange'
        else:  # 4+
            return 'red'

    def is_seg_num_single_copy(self, seg_num):
        return seg_num in self.copy_depths and len(self.copy_depths[seg_num]) == 1

    def get_single_copy_segments(self):
        """
        Returns a list of the graph segments with a copy number of 1.
        """
        single_copy_segments = []
        for num, segment in self.segments.items():
            if self.is_seg_num_single_copy(num):
                single_copy_segments.append(segment)
        return single_copy_segments

    def get_no_copy_depth_segments(self):
        """
        Returns a list of the graph segments which failed to have a copy depth assigned.
        """
        no_copy_depth_segments = []
        for num, segment in self.segments.items():
            if num not in self.copy_depths or len(self.copy_depths[num]) == 0:
                no_copy_depth_segments.append(segment)
        return no_copy_depth_segments

    def get_path_sequence(self, path_segments):
        """
        Gets a linear (i.e. not circular) path sequence from the graph.
        """
        path_sequence = ''
        prev_segment_number = None
        for i, seg_num in enumerate(path_segments):
            segment = self.segments[abs(seg_num)]
            if seg_num > 0:
                seg_sequence = segment.forward_sequence
            else:
                seg_sequence = segment.reverse_sequence
            if i == 0:
                path_sequence = seg_sequence
            else:
                if seg_num not in self.forward_links[prev_segment_number]:
                    raise BadPath(str(path_segments) + ' is not a valid path')
                if self.overlap > 0 and \
                        path_sequence[-self.overlap:] != seg_sequence[:self.overlap]:
                    raise BadOverlaps('overlaps do not match when merging ' +
                                      str(prev_segment_number) + ' and ' + str(seg_num) +
                                      ' in path ' + str(path_segments))
                path_sequence += seg_sequence[self.overlap:]
            prev_segment_number = seg_num
        return path_sequence

    def apply_bridges(self, bridges, verbosity, min_bridge_qual):
        """
        Uses the supplied bridges to simplify the graph.
        """
        log.log_section_header('Applying bridges')
        log.log_explanation('Unicycler now applies to the graph in decreasing order of quality. '
                            'This ensures that when multiple, contradictory bridges exist, the '
                            'most supported option is used.')

        unbridged_graph = copy.deepcopy(self)

        # Each segment can have only one bridge per side, so we will track which segments have had
        # a bridge applied off one side or the other.
        right_bridged = set()
        left_bridged = set()
        seg_nums_used_in_bridges = []
        applied_bridges = []

        # Sort bridges first by type: LongReadBridge, SpadesContigBridge and then
        # LoopUnrollingBridge. Then sort by quality so within each type we apply the best bridges
        # first.
        sorted_bridges = sorted(bridges, key=lambda x: (x.get_type_score(), x.quality),
                                reverse=True)

        # The displayed table will only show applied bridges for verbosity 1. Higher verbosity
        # levels will show all bridges and whether or not they were applied.
        bridge_application_table = [['Bridge type', 'Start ' + get_right_arrow() + ' end', 'Path',
                                     'Quality']]
        if verbosity > 1:
            bridge_application_table[0].append('Result')
        table_row_colours = {}
        for bridge in sorted_bridges:

            # Is the bridge available to be applied? The first criterion is simple: neither the
            # start nor end can have already been bridged in the proposed direction.
            can_use_bridge = self.start_end_available_to_bridge(bridge.start_segment,
                                                                bridge.end_segment,
                                                                right_bridged, left_bridged)
            if can_use_bridge:
                # If a bridge has multiple equally-good graph paths, then we can choose which one
                # to use based on the availability of the path.
                if isinstance(bridge, LongReadBridge) and len(bridge.all_paths) > 1:
                    bridge.set_path_based_on_availability(self, unbridged_graph)

                # The second criterion for whether the bridge can be applied is a bit trickier. We
                # need to check whether either the start or end segment has already been used the
                # path of a different bridge. That alone isn't necessarily a problem, and since
                # single copy determination can make mistakes we want to allow for this sort of
                # thing. But it is a problem if the start or end segment has been used in a
                # bridge that happens to start or end in this bridge. That arrangement (two bridges,
                # each of which end inside the other's path) can break up the graph if they are
                # both applied, so don't apply this bridge if such a case exists.
                bridges_using_this_segment = []
                if abs(bridge.start_segment) in seg_nums_used_in_bridges:
                    for applied_bridge in applied_bridges:
                        segs_in_applied_bridge_path = set(abs(x) for x in applied_bridge.graph_path)
                        if abs(bridge.start_segment) in segs_in_applied_bridge_path:
                            bridges_using_this_segment.append(applied_bridge)
                if abs(bridge.end_segment) in seg_nums_used_in_bridges:
                    for applied_bridge in applied_bridges:
                        segs_in_applied_bridge_path = set(abs(x) for x in applied_bridge.graph_path)
                        if abs(bridge.end_segment) in segs_in_applied_bridge_path:
                            bridges_using_this_segment.append(applied_bridge)
                if bridges_using_this_segment:
                    segs_in_path = set(abs(x) for x in bridge.graph_path)
                    for bridge_using_this_segment in bridges_using_this_segment:
                        if abs(bridge_using_this_segment.start_segment) in segs_in_path or \
                                        abs(bridge_using_this_segment.end_segment) in segs_in_path:
                            can_use_bridge = False

            start_to_end = (str(bridge.start_segment) + ' ' + get_right_arrow()).rjust(7) + ' ' + \
                str(bridge.end_segment)
            bridge_application_table_row = [bridge.get_type_name(), start_to_end,
                                            ', '.join([str(x) for x in bridge.graph_path]),
                                            '%.3f' % bridge.quality]
            if can_use_bridge:
                # Even if there's no conflict with other bridges, the quality still needs to be
                # high enough for this bridge to be applicable.
                if bridge.quality >= min_bridge_qual:
                    self.apply_bridge(bridge, right_bridged, left_bridged, seg_nums_used_in_bridges)
                    seg_nums_used_in_bridges = remove_dupes_preserve_order(seg_nums_used_in_bridges)
                    applied_bridges.append(bridge)
                    if verbosity > 1:
                        bridge_application_table_row.append('applied')
                    bridge_application_table.append(bridge_application_table_row)
                elif verbosity > 1:
                    table_row_colours[len(bridge_application_table)] = 'dim'
                    bridge_application_table.append(bridge_application_table_row + ['rejected'])
            elif verbosity > 1:
                table_row_colours[len(bridge_application_table)] = 'dim'
                bridge_application_table.append(bridge_application_table_row + ['unused'])

        print_table(bridge_application_table, alignments='LLLRR', indent=0,
                    sub_colour={'applied': 'green', 'rejected': 'clear_red'},
                    row_colour=table_row_colours, max_col_width=40)
        return set(seg_nums_used_in_bridges)

    def apply_bridge(self, bridge, right_bridged, left_bridged, seg_nums_used_in_bridges):
        """
        Applies a whole bridge, start to end.
        """
        # Remove all existing links for the segments being bridged.
        start = bridge.start_segment
        end = bridge.end_segment
        if start in self.forward_links:
            for link in [x for x in self.forward_links[start]]:
                self.remove_link(start, link)
        if end in self.reverse_links:
            for link in [x for x in self.reverse_links[end]]:
                self.remove_link(link, end)

        # If the new bridge is a miniasm bridge, then we might need to trim a bit from the segments
        # being bridged.
        if isinstance(bridge, MiniasmBridge):
            start_seg = self.segments[abs(start)]
            if start > 0:
                start_seg.trim_from_end(bridge.start_overlap)
            else:
                start_seg.trim_from_start(bridge.start_overlap)
            end_seg = self.segments[abs(end)]
            if end > 0:
                end_seg.trim_from_start(bridge.end_overlap)
            else:
                end_seg.trim_from_end(bridge.end_overlap)

        # Create a new bridge segment.
        new_seg_num = self.get_next_available_seg_number()
        new_seg = Segment(new_seg_num, bridge.depth, bridge.bridge_sequence, True, bridge,
                          bridge.graph_path)
        new_seg.build_other_sequence_if_necessary()
        self.segments[new_seg_num] = new_seg

        # Link the bridge segment in to the start/end segments.
        self.add_link(start, new_seg_num)
        self.add_link(new_seg_num, end)

        # Add the bridge to the segment (which will reduce the segment's depth if it's a new bridge)
        for seg_num in list(set(bridge.graph_path)):
            self.add_bridge_to_segment(self.segments[abs(seg_num)], bridge)

        add_to_bridged_sets(bridge.start_segment, bridge.end_segment, right_bridged, left_bridged)
        seg_nums_used_in_bridges.extend([abs(x) for x in bridge.graph_path])

    def add_bridge_to_segment(self, segment, bridge):
        """
        Adds a bridge that uses the segment. This function checks whether this bridge is new,
        and if so, subtracts the appropriate depth from the segment.
        """
        full_bridge_path = [bridge.start_segment] + bridge.graph_path + [bridge.end_segment]
        bridge_str = '_' + '_'.join([str(x) for x in full_bridge_path]) + '_'

        # If this is the first used-in bridge, we don't need to check anything for redundancy.
        if not segment.used_in_bridges:
            segment.used_in_bridges.append(bridge_str)
            self.subtract_depth_from_segment(segment, bridge)

        # If there are already used-in bridges, then we need to check for redundancy.
        else:
            reverse_bridge_path = [-x for x in full_bridge_path[::-1]]
            reverse_bridge_str = '_' + '_'.join([str(x) for x in reverse_bridge_path]) + '_'
            new_used_in_bridges = []
            redundancy_found = False
            for used_in_bridge in segment.used_in_bridges:
                if bridge_str in used_in_bridge or reverse_bridge_str in used_in_bridge:
                    new_used_in_bridges.append(used_in_bridge)
                    redundancy_found = True
                elif used_in_bridge in bridge_str or used_in_bridge in reverse_bridge_str:
                    new_used_in_bridges.append(bridge_str)
                    redundancy_found = True
                else:
                    new_used_in_bridges.append(used_in_bridge)
            segment.used_in_bridges = new_used_in_bridges
            if not redundancy_found:
                segment.used_in_bridges.append(bridge_str)
                self.subtract_depth_from_segment(segment, bridge)

    def subtract_depth_from_segment(self, seg, bridge):
        """
        Removes the given depth from the segment. Allows depths to go into the negative.
        """
        seg_num = seg.number
        removed_depth = bridge.depth
        seg.depth -= removed_depth
        seg.original_depth = False
        if seg_num in self.copy_depths and self.copy_depths[seg_num]:
            removed_copy_depth = min(self.copy_depths[seg_num],
                                     key=lambda x: abs(x - removed_depth))
            del self.copy_depths[seg_num][self.copy_depths[seg_num].index(removed_copy_depth)]
        else:
            removed_copy_depth = None
        bridge.segments_reduced_depth.append((seg_num, removed_depth, removed_copy_depth))

    @staticmethod
    def start_end_available_to_bridge(start, end, right_bridged, left_bridged):
        """
        Checks whether the start and end segments can be bridged together (i.e. that they are both
        unbridged on the relevant sides and not yet used in a bridge).
        """
        if start > 0 and start in right_bridged:
            return False
        if start < 0 and -start in left_bridged:
            return False
        if end > 0 and end in left_bridged:
            return False
        if end < 0 and -end in right_bridged:
            return False
        return True

    def clean_up_after_bridging_1(self, anchor_segments, seg_nums_used_in_bridges):
        """
        This function is run after bridge application to clean up necessary segments. This is the
        first of two such functions, and this one takes care of the simpler aspects of cleaning.
        """
        log.log_section_header('Cleaning up leftover segments', 2)

        # For the purposes of cleaning up, a graph segment which is a bridge counts as a graph
        # segment used in a bridge.
        for seg_num, seg in self.segments.items():
            if seg.bridge is not None:
                seg_nums_used_in_bridges.add(seg_num)

        log.log('Segments eligible for deletion:', 2)
        log.log_number_list(sorted(list(seg_nums_used_in_bridges)), 2)

        anchor_seg_nums = set(x.number for x in anchor_segments)
        self.remove_unbridging_segments(anchor_seg_nums)
        self.remove_components_without_anchor_segments(anchor_seg_nums)
        self.remove_components_entirely_used_in_bridges(seg_nums_used_in_bridges)

    def clean_up_after_bridging_2(self, seg_nums_used_in_bridges, min_component_size,
                                  min_dead_end_size, unbridged_graph, anchor_segments):
        """
        This is the second of two post-bridging cleaning functions, and it takes care of the more
        complex aspects of cleaning: deleting segments that were used in bridges but are still
        part of important areas in the graph. This is critically important for areas in the graph
        which didn't get directly bridged. If all goes well, such areas will be cleaned up to the
        best possible paths - no more (which would be redundant) and no less (which could either
        break up the graph or lead to its oversimplification and a misassembly).
        Ideally, we should be deleting segments that are all 'used up' - i.e. they have been used
        bridges the same number of times as their copy number. This can be hard to get just
        right, however, because copy number determination and bridge paths aren't perfect.
        """
        removed_segments = []

        # Get all usedupness scores once, outside the loop, to save time in the loop.
        usedupness_scores = defaultdict(float)
        for seg_num in seg_nums_used_in_bridges:
            if seg_num in self.segments and seg_num in unbridged_graph.segments:
                usedupness_scores[seg_num] = self.get_usedupness_score(seg_num, unbridged_graph)

        # For the second pass, we also remove segments (or simple paths of segments) which can be
        # removed without creating any dead ends.
        while True:
            # First we remove as many segments as possible that are used in bridges and have dead
            # ends.
            while True:
                for seg_num in seg_nums_used_in_bridges:
                    if seg_num in self.segments and self.dead_end_count(seg_num) > 0:
                        self.remove_segments([seg_num])
                        removed_segments.append(seg_num)
                        break
                else:
                    break

            # When the code gets here, that means all possible used-in-bridge-dead-ends have been
            # removed. Now we want to remove segments (or groups of segments in simple paths) which
            # have been entirely used in bridges and can be removed without creating dead ends.

            # Group the segments based on simple paths. Segments not on a simple path will be in
            # their own group.
            path_groups = []
            segs_in_path_groups = set()
            for seg_num in seg_nums_used_in_bridges:
                if seg_num in self.segments and seg_num not in segs_in_path_groups:
                    # noinspection PyTypeChecker
                    path = self.get_simple_path(seg_num, None, 2)
                    if all(abs(x) in seg_nums_used_in_bridges for x in path):
                        path_groups.append(path)
                        segs_in_path_groups.update(path)

            # Sort the path groups by how likely it is that they are truly all 'used up'. These are
            # the ones we would like to remove first. Each segment in the path is scored on its
            # 'usedupness', and the path gets the lowest score of its constituent segments.
            scored_path_groups = []
            for path_group in path_groups:
                min_score = 100.0
                for path_seg in path_group:
                    min_score = min(min_score, usedupness_scores[abs(path_seg)])
                scored_path_groups.append((min_score, path_group))
            scored_path_groups = sorted(scored_path_groups, reverse=True, key=lambda x: x[0])

            for _, path in scored_path_groups:
                if self.dead_end_change_if_path_deleted(path) <= 0:
                    unsigned_path = [abs(x) for x in path]
                    self.remove_segments(unsigned_path)
                    removed_segments += unsigned_path
                    break
            else:
                break

        # It's possible at this point that there are bubbles remaining in the graph which are
        # mostly used up. If we can delete them without introducing dead ends, we do so.
        while True:
            potentially_deletable_paths = []
            for seg_num in self.segments:
                # noinspection PyTypeChecker
                path = self.get_simple_path(seg_num, None, 2)
                path_lengths = [max(1, self.segments[abs(x)].get_length() - self.overlap)
                                for x in path]
                path_usedupness = [usedupness_scores[abs(x)] for x in path]
                average_usedupness = weighted_average_list(path_usedupness, path_lengths)
                potentially_deletable_paths.append((average_usedupness, path))
            for usedupness, path in potentially_deletable_paths:
                if usedupness > settings.CLEANING_USEDUPNESS_THRESHOLD and \
                                self.dead_end_change_if_path_deleted(path) <= 0:
                    unsigned_path = [abs(x) for x in path]
                    self.remove_segments(unsigned_path)
                    removed_segments += unsigned_path
                    break
            else:
                break

        # It's also possible for entire graph components to be mostly used up, in which case we can
        # delete those as well.
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            component_lengths = [self.segments[abs(x)].get_length() for x in component_nums]
            component_usedupness = [usedupness_scores[abs(x)] for x in component_nums]
            average_usedupness = weighted_average_list(component_usedupness, component_lengths)
            if average_usedupness > settings.CLEANING_USEDUPNESS_THRESHOLD:
                self.remove_segments(component_nums)
                removed_segments += component_nums

        if removed_segments:
            removed_segments = sorted(list(set(removed_segments)))
            log.log('Removed segments used in bridges:', 2)
            log.log_number_list(removed_segments, 2)

        # Now that clean up is finished, we no longer want to allow depths below zero.
        for segment in self.segments.values():
            segment.depth = max(0.0, segment.depth)

        anchor_seg_nums = set(x.number for x in anchor_segments)
        self.remove_components_without_anchor_segments(anchor_seg_nums)
        self.remove_components_entirely_used_in_bridges(seg_nums_used_in_bridges)
        self.remove_unbridging_segments(anchor_seg_nums)
        self.remove_small_components(min_component_size)
        self.remove_small_dead_ends(min_dead_end_size)

    def remove_components_without_anchor_segments(self, anchor_seg_nums):
        """
        Deletes all graph components that contain no anchor segments.
        """
        segment_nums_to_remove = []
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            for seg_num in component_nums:
                if abs(seg_num) in anchor_seg_nums:
                    break
            else:
                segment_nums_to_remove += component_nums
        if segment_nums_to_remove:
            log.log('Removed components with no single copy segments:', 2)
            log.log_number_list(sorted(segment_nums_to_remove), 2)
        self.remove_segments(segment_nums_to_remove)

    def remove_components_entirely_used_in_bridges(self, seg_nums_used_in_bridges):
        """
        Deletes all graph components which have been entirely used in bridges.
        """
        segment_nums_to_remove = []
        connected_components = self.get_connected_components()
        for component_nums in connected_components:
            for seg_num in component_nums:
                if abs(seg_num) not in seg_nums_used_in_bridges:
                    break
            else:
                segment_nums_to_remove += component_nums
        if segment_nums_to_remove:
            log.log('Removed components used in bridges:', 2)
            log.log_number_list(sorted(segment_nums_to_remove), 2)
        self.remove_segments(segment_nums_to_remove)

    def remove_unbridging_segments(self, anchor_seg_nums):
        """
        Deletes any segments which cannot possibly connect two anchor segments.
        """
        segment_nums_to_remove = []
        for seg_num in self.segments:
            if seg_num in anchor_seg_nums:
                continue
            if not (self.search(seg_num, anchor_seg_nums) and
                    self.search(-seg_num, anchor_seg_nums)):
                segment_nums_to_remove.append(seg_num)
        if segment_nums_to_remove:
            log.log('Removed unbridging segments:', 2)
            log.log_number_list(segment_nums_to_remove, 2)
        self.remove_segments(segment_nums_to_remove)

    def get_usedupness_score(self, seg_num, unbridged_graph):
        """
        Returns a score for the segment which reflects how likely it is that the segment has been
        'used up' in bridges. E.g. a segment which originally had a depth of 2.0 and now has a
        depth of 0.03 is probably used up, but if its depth is now 0.94, it's less likely that
        it's used up.
        """
        original_depth = unbridged_graph.segments[seg_num].depth
        current_depth = self.segments[seg_num].depth
        depth_used = original_depth - current_depth

        # Since segment depths can get negative, depth_fraction_used can exceed 1.0.
        depth_fraction_used = depth_used / original_depth

        # A score penalty is applied based on the original depth. For example, a segment that
        # originally had 20x depth and is now down to 2x is less confidently used up than a segment
        # which originally had 2x depth and is now down to 0.2x.
        penalty = score_function(original_depth, 4.0)

        return depth_fraction_used - (penalty / 2.0)

    def find_all_simple_loops(self):
        """
        This function finds all cases of a simple loop in the graph: A->B->C->B->D.
        It returns them as a list of 4-tuples of segment numbers in this order:
        (start, end, middle, repeat).
        """
        simple_loops = []

        # We'll search specifically for the middle segments as they should be easy to spot.
        for middle in self.segments:
            if self.segments[middle].get_length() > settings.MAX_SIMPLE_LOOP_SIZE:
                continue

            # A middle segment will always have exactly one connection on each end which connect
            # to the same segment (the repeat segment).
            if middle not in self.forward_links or middle not in self.reverse_links:
                continue
            if len(self.forward_links[middle]) != 1 or len(self.reverse_links[middle]) != 1:
                continue
            if self.forward_links[middle][0] != self.reverse_links[middle][0]:
                continue
            repeat = self.forward_links[middle][0]

            # The repeat segment should have exactly two connections on each end. If less, then we
            # have a simple path which can be merged. If more, it's a more complex loop.
            if len(self.forward_links[repeat]) != 2 or len(self.reverse_links[repeat]) != 2:
                continue

            # Find the start and end segment numbers. It's okay if the start and the end are the
            # same, but we exclude any other screwy cases where the start or end is the middle or
            # repeat segment.
            start = self.reverse_links[repeat][0]
            if abs(start) == abs(middle):
                start = self.reverse_links[repeat][1]
            if abs(start) == abs(middle) or abs(start) == abs(repeat):
                continue

            end = self.forward_links[repeat][0]
            if abs(end) == abs(middle):
                end = self.forward_links[repeat][1]
            if abs(end) == abs(middle) or abs(end) == abs(repeat):
                continue

            simple_loops.append((start, end, middle, repeat))

        # Simple loops may just have the repeat node loop back to itself (i.e. no separate middle
        # segment). Look for these structures.
        for repeat in self.segments:
            if repeat not in self.forward_links or repeat not in self.reverse_links or \
                    len(self.forward_links[repeat]) != 2 or len(self.reverse_links[repeat]) != 2:
                continue

            # Make sure that the repeat loops back to itself.
            if repeat not in self.forward_links[repeat] or repeat not in self.reverse_links[repeat]:
                continue

            start_segs = list(self.reverse_links[repeat])
            start_segs.remove(repeat)
            end_segs = list(self.forward_links[repeat])
            end_segs.remove(repeat)
            if len(start_segs) != 1 or len(end_segs) != 1:
                continue
            start = start_segs[0]
            end = end_segs[0]
            if abs(start) == abs(repeat) or abs(end) == abs(repeat):
                continue

            simple_loops.append((start, end, None, repeat))

        return simple_loops

    def max_path_segment_count(self, seg_num, start_end_depth):
        """
        This function returns the maximum allowed number of times a segment can be in a bridge
        path. It uses both the segment's copy depth (if it has one) and the relative depth of the
        segment as compared to the start/end of the bridge.
        """
        if abs(seg_num) in self.copy_depths:
            count_by_copies = len(self.copy_depths[abs(seg_num)])
        else:
            count_by_copies = 1
        depth = self.segments[abs(seg_num)].depth
        count_by_depth = max(1, int(round(depth / start_end_depth)))
        return 2 * max(count_by_copies, count_by_depth)

    def get_path_length(self, path):
        """
        Returns the length of the given path.
        """
        if not path:
            return 0
        try:
            path_length = 0
            for seg in path:
                path_length += self.segments[abs(seg)].get_length()
            overlap_count = len(path) - 1
            path_length -= overlap_count * self.overlap
            return path_length
        except KeyError:
            return 0

    def get_bridge_path_length(self, path):
        """
        Like get_path_length, but if the path is empty it returns the graph overlap size (for a
        direct connection).
        """
        if not path:
            return self.overlap
        else:
            return self.get_path_length(path)

    def renumber_segments(self):
        """
        This function gives the longest segment the number 1, the second-longest the number 2, etc.
        """
        old_nums = [x.number for x in sorted(self.segments.values(), reverse=True,
                                             key=lambda x: x.get_length())]
        new_nums = list(range(1, len(old_nums) + 1))
        old_nums += [-x for x in old_nums]
        new_nums += [-x for x in new_nums]
        changes = dict(zip(old_nums, new_nums))

        new_segments = {}
        for seg_num, seg in self.segments.items():
            new_num = changes[seg_num]
            seg.number = new_num
            new_segments[new_num] = seg
        self.segments = new_segments

        new_forward_links = {}
        for seg_num, link_nums in self.forward_links.items():
            if link_nums:
                new_forward_links[changes[seg_num]] = [changes[x] for x in link_nums]
        self.forward_links = new_forward_links

        new_reverse_links = {}
        for seg_num, link_nums in self.reverse_links.items():
            if link_nums:
                new_reverse_links[changes[seg_num]] = [changes[x] for x in link_nums]
        self.reverse_links = new_reverse_links

        self.copy_depths = {changes[x]: y for x, y in self.copy_depths.items()}

        new_paths = {}
        for name, path_nums in self.paths.items():
            new_paths[name] = [changes[x] for x in path_nums]
        self.paths = new_paths

    def print_component_table(self):
        component_table = [['Component', 'Segments', 'Links', 'Length', 'N50',
                            'Longest segment', 'Status']]
        components = self.get_connected_components()
        if len(components) > 1:
            n50, _, _, _, _, longest = self.get_contig_stats()
            total_row = ['total',
                         int_to_str(len(self.segments)),
                         int_to_str(self.get_total_link_count()),
                         int_to_str(self.get_total_length()),
                         int_to_str(n50),
                         int_to_str(longest),
                         '']
            component_table.append(total_row)
        for i, component in enumerate(components):
            status = 'complete' if self.is_component_complete(component) else 'incomplete'
            component_len = sum(self.segments[x].get_length() for x in component)
            segment_count = len(component)
            link_count = self.get_component_link_count(component)
            n50, _, _, _, _, longest = self.get_contig_stats(seg_nums=component)
            component_row = [str(i+1),
                             int_to_str(segment_count),
                             int_to_str(link_count),
                             int_to_str(component_len),
                             int_to_str(n50),
                             int_to_str(longest),
                             status]
            component_table.append(component_row)
        print_table(component_table, alignments='RRRRRRR', indent=0,
                    sub_colour={' complete': 'green', ' incomplete': 'red'})

    def get_total_link_count(self):
        """
        Returns the total number of forward links in the graph, not counting rev comp duplicates.
        """
        links = set()
        for start, ends in self.forward_links.items():
            for end in ends:
                if (start, end) not in links and (-end, -start) not in links:
                    links.add((start, end))
        return len(links)

    def get_component_link_count(self, component_segs):
        """
        Returns the total number of forward links in the component, not counting rev comp
        duplicates. This function assumes the given segments make up a connected component - it
        doesn't check.
        """
        links = set()
        component_segs = set(component_segs)  # positive segment numbers
        for start, ends in self.forward_links.items():
            for end in ends:
                if abs(start) in component_segs and abs(end) in component_segs and \
                        (start, end) not in links and (-end, -start) not in links:
                    links.add((start, end))
        return len(links)

    def get_contig_stats(self, seg_nums=None):
        """
        Returns various contig length metrics.
        """
        if seg_nums is None:
            segs = self.segments.values()
        else:
            segs = [self.segments[x] for x in seg_nums]
        segment_lengths = sorted([x.get_length() for x in segs])
        if not segment_lengths:
            return 0, 0, 0, 0, 0, 0

        shortest = segment_lengths[0]
        longest = segment_lengths[-1]

        first_quartile_index = (len(segment_lengths) - 1) / 4
        median_index = (len(segment_lengths) - 1) / 2
        third_quartile_index = (len(segment_lengths) - 1) * 3 / 4

        first_quartile = int(round(value_from_fractional_index(segment_lengths,
                                                               first_quartile_index)))
        median = int(round(value_from_fractional_index(segment_lengths, median_index)))
        third_quartile = int(round(value_from_fractional_index(segment_lengths,
                                                               third_quartile_index)))

        half_total_length = sum(segment_lengths) / 2
        total_so_far = 0
        segment_lengths = segment_lengths[::-1]
        for length in segment_lengths:
            total_so_far += length
            if total_so_far >= half_total_length:
                n50 = length
                break
        else:
            n50 = 0

        return n50, shortest, first_quartile, median, third_quartile, longest

    def completed_circular_replicons(self):
        """
        Returns a list of graph components which are simple loops: one segment connected to itself
        to make a circular piece of DNA.
        """
        completed_components = []
        single_segment_components = [x for x in self.get_connected_components() if len(x) == 1]
        for component in single_segment_components:
            only_segment = component[0]
            if only_segment in self.forward_links and \
                    self.forward_links[only_segment] == [only_segment] and \
                    only_segment in self.reverse_links and \
                    self.reverse_links[only_segment] == [only_segment]:
                completed_components.append(only_segment)
        return completed_components

    def is_component_complete(self, component):
        """
        Given a list of unsigned segment numbers, this function returns whether the component is
        a completed circular replicon.
        """
        if len(component) != 1:
            return False
        seg = component[0]
        if self.get_downstream_seg_nums(seg) != [seg]:
            return False
        return self.get_upstream_seg_nums(seg) == [seg]

    def get_simple_path(self, starting_seg, single_copy_seg_nums, bridging_mode):
        """
        Starting with the given segment, this function tries to expand outward as far as possible
        while maintaining a simple (i.e. can be merged) path. If it can't expand at all, it will
        just return a list of the starting segment. At lower bridging modes, we only allow the
        merging of paths which are made up of single copy segments and bridges.
        """
        simple_path = [starting_seg]

        # Expand forward as much as possible.
        while True:
            if simple_path[-1] not in self.forward_links or \
                            len(self.forward_links[simple_path[-1]]) != 1:
                break
            potential = self.forward_links[simple_path[-1]][0]
            if potential in simple_path or -potential in simple_path:
                break
            abs_potential = abs(potential)
            if bridging_mode < 2 and not self.is_single_copy_or_bridge(abs_potential, bridging_mode,
                                                                       single_copy_seg_nums):
                break
            if len(self.reverse_links[potential]) == 1 and \
                    self.reverse_links[potential][0] == simple_path[-1]:
                simple_path.append(potential)
            else:
                break

        # Expand backward as much as possible.
        while True:
            if simple_path[0] not in self.reverse_links or \
                            len(self.reverse_links[simple_path[0]]) != 1:
                break
            potential = self.reverse_links[simple_path[0]][0]
            if potential in simple_path or -potential in simple_path:
                break
            abs_potential = abs(potential)
            if bridging_mode < 2 and not self.is_single_copy_or_bridge(abs_potential, bridging_mode,
                                                                       single_copy_seg_nums):
                break
            if len(self.forward_links[potential]) == 1 and \
                    self.forward_links[potential][0] == simple_path[0]:
                simple_path.insert(0, potential)
            else:
                break

        return simple_path

    def sort_link_order(self):
        """
        This function sorts the lists in links so path finding can be consistent from one run to
        the next.
        """
        for seg_num in self.forward_links:
            self.forward_links[seg_num].sort()
        for seg_num in self.reverse_links:
            self.reverse_links[seg_num].sort()

    def search(self, start, ends):
        """
        Conducts a DFS from the start segment to see if it leads to any of the end segments.
        The start segment is signed, i.e. positive start and negative start will conduct the
        search in different directions. The end segments are not signed, i.e. the search is
        successful if it reaches either orientation of an end.
        """
        end_set = set(ends)
        end_set.update(-x for x in ends)
        visited, stack = set(), [start]
        while stack:
            seg = stack.pop()
            if seg not in visited:
                visited.add(seg)
                if seg in self.forward_links:
                    for next_seg in self.forward_links[seg]:
                        if next_seg in end_set:
                            return True
                        if next_seg not in visited:
                            stack.append(next_seg)
        return False

    def get_path_availability(self, path):
        """
        Given a path, this function returns the fraction that is available. A segment is considered
        fully available it has a depth of above 0.5. Below that, availability drops towards 0. A
        single segment can have negative availability (if it has negative depth), but this function
        will only return a number from 0 to 1.
        """
        total_bases = 0
        available_bases = 0.0
        for seg_num in path:
            seg = self.segments[abs(seg_num)]
            if seg.depth >= 0.5:
                seg_availability = 1.0
            else:
                seg_availability = 2 * seg.depth
            seg_len = seg.get_length() - self.overlap
            total_bases += seg_len
            available_bases += seg_len * seg_availability
        if total_bases == 0:
            return 1.0
        else:
            return max(0.0, available_bases / total_bases)

    def get_estimated_sequence_len(self):
        """
        Returns an estimate sequence length, based on the depth. E.g. a segment of depth 1.0
        contributes it own length (minus overlaps), a segment of depth 2.5 contributes 2.5 times
        its own length, etc.
        """
        single_copy_depth = self.get_single_copy_depth()
        total_seq_len = 0.0
        for seg_num, seg in self.segments.items():
            seg_len = seg.get_length()
            if seg_num in self.forward_links:
                seg_len -= self.overlap / 2
            if seg_num in self.reverse_links:
                seg_len -= self.overlap / 2
            seg_len *= (seg.depth / single_copy_depth)
            total_seq_len += seg_len
        return total_seq_len

    def remove_all_overlaps(self):
        """
        This function removes all overlaps from the graph by shortening segments. It assumes that
        all overlaps in the graph are the same size.
        """
        if self.overlap == 0:
            log.log('Graph has no overlaps - overlap removal not needed')
            return

        # First we create a set of all graph edges, in both directions.
        all_edges = set()
        for start, ends in self.forward_links.items():
            for end in ends:
                all_edges.add((start, end))
                all_edges.add((-end, -start))

        # The overlap to be removed is an odd number, as SPAdes only uses odd k-mers. We'll split
        # this value approximately in half.
        large_half = int(math.ceil(self.overlap / 2))
        small_half = int(math.floor(self.overlap / 2))

        # This function will strive to put each edge into one of two groups:
        #   1) Trim more sequence from the end of the starting segment
        #   2) Trim more sequence from start of the ending segment

        # To do this, we determine which edges which must be in the same group as each other and
        # edges which must be in different groups.
        must_match = defaultdict(set)
        must_differ = defaultdict(set)

        # Firstly, each edge must be in the opposite group as its complement edge. E.g. if we
        # have an edge (5, -4) and trim more from the first segment, then we need to trim more from
        # the second segment of its complement edge (4, -5).
        for edge in all_edges:
            rev_edge = (-edge[1], -edge[0])
            must_differ[edge].add(rev_edge)
            must_differ[rev_edge].add(edge)

        # Edges which connect to the same side of a segment must be grouped together.
        pos_and_neg_seg_nums = list(self.segments) + [-x for x in self.segments]
        for seg in pos_and_neg_seg_nums:
            downstream_segs = self.get_downstream_seg_nums(seg)
            if len(downstream_segs) > 1:
                edge_1_for = (seg, downstream_segs[0])
                edge_1_rev = (-downstream_segs[0], -seg)
                for downstream_seg in downstream_segs[1:]:
                    edge_2_for = (seg, downstream_seg)
                    edge_2_rev = (-downstream_seg, -seg)
                    must_match[edge_1_for].add(edge_2_for)
                    must_match[edge_2_for].add(edge_1_for)
                    must_match[edge_1_rev].add(edge_2_rev)
                    must_match[edge_2_rev].add(edge_1_rev)
            upstream_segs = self.get_upstream_seg_nums(seg)
            if len(upstream_segs) > 1:
                edge_1_for = (upstream_segs[0], seg)
                edge_1_rev = (-seg, -upstream_segs[0])
                for upstream_seg in upstream_segs[1:]:
                    edge_2_for = (upstream_seg, seg)
                    edge_2_rev = (-seg, -upstream_seg)
                    must_match[edge_1_for].add(edge_2_for)
                    must_match[edge_2_for].add(edge_1_for)
                    must_match[edge_1_rev].add(edge_2_rev)
                    must_match[edge_2_rev].add(edge_1_rev)

        # Segments which are equal to the overlap size cannot have the larger trim applied to
        # both sides, so we require that edges on opposite sides of these segments to be grouped
        # together.
        small_seg_nums = [x for x in pos_and_neg_seg_nums
                          if self.segments[abs(x)].get_length() == self.overlap]
        for seg in small_seg_nums:
            downstream_segs = self.get_downstream_seg_nums(seg)
            upstream_segs = self.get_upstream_seg_nums(seg)
            if downstream_segs and upstream_segs:
                for downstream_seg in downstream_segs:
                    edge_1_for = (seg, downstream_seg)
                    edge_1_rev = (-downstream_seg, -seg)
                    for upstream_seg in upstream_segs:
                        edge_2_for = (upstream_seg, seg)
                        edge_2_rev = (-seg, -seg)
                        must_match[edge_1_for].add(edge_2_for)
                        must_match[edge_2_for].add(edge_1_for)
                        must_match[edge_1_rev].add(edge_2_rev)
                        must_match[edge_2_rev].add(edge_1_rev)

        # Now for the grouping! We order the edges for consistency from one run to the next.
        ordered_edges = list(all_edges)
        group_1 = set()
        group_2 = set()
        num_edges = len(ordered_edges)
        extra_verbose = log.logger.stdout_verbosity_level > 2
        if extra_verbose:
            log.log('Grouping graph edges based on overlap removal')
        edge_num = 0
        for edge in ordered_edges:
            if extra_verbose:
                log.log_progress_line(edge_num, num_edges)
                edge_num += 1
            if edge in group_1 or edge in group_2:
                continue

            new_group_1 = set()
            new_group_2 = set()

            # Put an edge in the first group (arbitrary decision).
            new_group_1.add(edge)

            while True:
                group_1_size_before = len(new_group_1)
                group_2_size_before = len(new_group_2)

                group_1_temp = set()
                group_2_temp = set()

                for group_1_edge in new_group_1:
                    for must_match_edge in must_match[group_1_edge]:
                        group_1_temp.add(must_match_edge)
                    for must_differ_edge in must_differ[group_1_edge]:
                        group_2_temp.add(must_differ_edge)

                for group_2_edge in new_group_2:
                    for must_match_edge in must_match[group_2_edge]:
                        group_2_temp.add(must_match_edge)
                    for must_differ_edge in must_differ[group_2_edge]:
                        group_1_temp.add(must_differ_edge)

                new_group_1.update(group_1_temp)
                new_group_2.update(group_2_temp)

                # We continue to loop until the groups stop growing.
                if len(new_group_1) == group_1_size_before and \
                        len(new_group_2) == group_2_size_before:
                    break

            group_1.update(new_group_1)
            group_2.update(new_group_2)

        if extra_verbose:
            log.log_progress_line(num_edges, num_edges, end_newline=True)

        # Make sure there are no conflicts.
        for group_1_edge in group_1:
            for must_match_edge in must_match[group_1_edge]:
                if must_match_edge in group_2:
                    raise CannotTrimOverlaps
            for must_differ_edge in must_differ[group_1_edge]:
                if must_differ_edge in group_1:
                    raise CannotTrimOverlaps
        for group_2_edge in group_2:
            for must_match_edge in must_match[group_2_edge]:
                if must_match_edge in group_1:
                    raise CannotTrimOverlaps
            for must_differ_edge in must_differ[group_2_edge]:
                if must_differ_edge in group_2:
                    raise CannotTrimOverlaps

        # If the code got here, that means that all edges have been grouped according to the
        # rules, so now we produce sets of what to do for each segment. Segments in the
        # large_trim_end set will have the larger overlap trimmed from their end. Segments not in
        # that set will have the shorter overlap trimmed from their end. And similarly for the
        # large_trim_start set. These two sets are not exclusive - a segment may be in neither,
        # just one or both.
        large_trim_end = set()
        large_trim_start = set()
        for edge in group_1:
            start_seg = edge[0]
            if start_seg > 0:
                large_trim_end.add(start_seg)
            else:
                large_trim_start.add(-start_seg)
        for edge in group_2:
            end_seg = edge[1]
            if end_seg > 0:
                large_trim_start.add(end_seg)
            else:
                large_trim_end.add(-end_seg)

        # Now we finally do the segment trimming!
        log.log('\nRemoving graph overlaps\n', 3)
        log.log('             Bases     Bases', 3)
        log.log('           trimmed   trimmed', 3)
        log.log(' Segment      from      from', 3)
        log.log('  number     start       end', 3)
        for seg_num, segment in self.segments.items():
            start_trim = large_half if seg_num in large_trim_start else small_half
            end_trim = large_half if seg_num in large_trim_end else small_half
            segment.trim_from_start(start_trim)
            segment.trim_from_end(end_trim)
            log.log(str(seg_num).rjust(8) + str(start_trim).rjust(10) + str(end_trim).rjust(10), 3)

        log.log('Graph overlaps removed')
        self.overlap = 0

    def get_downstream_seg_nums(self, seg_num):
        """
        Returns the list of downstream segments for the given segment (or an empty list if there
        are none).
        """
        if seg_num not in self.forward_links:
            return []
        else:
            return self.forward_links[seg_num]

    def get_upstream_seg_nums(self, seg_num):
        """
        Returns the list of upstream segments for the given segment (or an empty list if there
        are none).
        """
        if seg_num not in self.reverse_links:
            return []
        else:
            return self.reverse_links[seg_num]

    def remove_zero_length_segs(self, suppress_log=False):
        """
        This function removes zero-length segments from the graph (segments with a length equal
        to the graph overlap), but only if they they aren't serving a purpose (such as in a
        multi-way junction).
        """
        segs_to_remove = []
        for seg_num in sorted(self.segments):  # sort for consistency between runs
            seg = self.segments[seg_num]
            if seg.get_length() != self.overlap:
                continue
            if seg_num in self.forward_links:
                forward_connections = len(self.forward_links[seg_num])
            else:
                forward_connections = 0
            if seg_num in self.reverse_links:
                reverse_connections = len(self.reverse_links[seg_num])
            else:
                reverse_connections = 0

            # Don't delete junction points.
            if forward_connections > 1 and reverse_connections > 1:
                continue
            if forward_connections == 0 and reverse_connections > 1:
                continue
            if forward_connections > 1 and reverse_connections == 0:
                continue

            segs_to_remove.append(seg_num)

            if forward_connections == 1 and reverse_connections > 0:
                downstream_seg = self.forward_links[seg_num][0]
                for upstream_seg in self.reverse_links[seg_num]:
                    self.add_link(upstream_seg, downstream_seg)

            elif reverse_connections == 1 and forward_connections > 0:
                upstream_seg = self.reverse_links[seg_num][0]
                for downstream_seg in self.forward_links[seg_num]:
                    self.add_link(upstream_seg, downstream_seg)

        if segs_to_remove:
            self.remove_segments(segs_to_remove)
        if not suppress_log and segs_to_remove:
            log.log('')
            log.log('Removed zero-length segments:')
            log.log_number_list(segs_to_remove)
        if suppress_log and segs_to_remove:
            log.log('  Removed zero-length segments:', 2)
            log.log_number_list(segs_to_remove, 2)
        return len(segs_to_remove)

    def merge_small_segments(self, max_merge_size):
        """
        In some cases, small segments can be merged into neighbouring segments to simplify the
        graph somewhat. This function does just that!
        """
        assert self.overlap == 0
        merged_seg_nums = []
        while True:
            for seg_num in sorted(self.segments):  # sort for consistency between runs
                segment = self.segments[seg_num]
                if segment.get_length() > max_merge_size or segment.get_length() == 0:
                    continue
                downstream_segs = self.get_downstream_seg_nums(seg_num)
                upstream_segs = self.get_upstream_seg_nums(seg_num)

                # If the segment has one downstream segment and multiple upstream segments,
                # then we can merge it into the upstream segments.
                if len(downstream_segs) == 1 and len(upstream_segs) > 1 and \
                        all(self.lead_exclusively_to(x, seg_num) for x in upstream_segs):
                    for upstream_seg_num in upstream_segs:
                        upstream_seg = self.segments[abs(upstream_seg_num)]
                        if upstream_seg_num > 0:
                            upstream_seg.append_to_forward_sequence(segment.forward_sequence)
                        else:
                            upstream_seg.append_to_reverse_sequence(segment.forward_sequence)
                    segment.remove_sequence()
                    merged_seg_nums.append(seg_num)
                    break

                # If the segment has one upstream segment and multiple downstream segments,
                # then we can merge it into the downstream segments.
                if len(upstream_segs) == 1 and len(downstream_segs) > 1 and \
                        all(self.lead_exclusively_from(x, seg_num) for x in downstream_segs):
                    for downstream_seg_num in downstream_segs:
                        downstream_seg = self.segments[abs(downstream_seg_num)]
                        if downstream_seg_num > 0:
                            downstream_seg.prepend_to_forward_sequence(segment.forward_sequence)
                        else:
                            downstream_seg.prepend_to_reverse_sequence(segment.forward_sequence)
                    segment.remove_sequence()
                    merged_seg_nums.append(seg_num)
                    break
            else:
                break
            self.remove_zero_length_segs(suppress_log=True)

        if merged_seg_nums:
            log.log('\nMerged small segments:')
            log.log_number_list(merged_seg_nums)
            self.remove_zero_length_segs()
        return len(merged_seg_nums)

    def expand_repeats(self):
        """
        This function moves sequence into repeat segments, wherever possible.
        """
        for seg_num in sorted(self.segments):  # sort for consistency between runs

            def trim_amount_okay(seg_nums, trim_length):
                for num in seg_nums:
                    pos_seg_num = abs(num)
                    seg_count = [abs(x) for x in seg_nums].count(pos_seg_num)  # should be 1 or 2
                    if seg_count * trim_length > self.segments[pos_seg_num].get_length():
                        return False
                return True

            segment = self.segments[seg_num]
            inputs = sorted(self.get_upstream_seg_nums(seg_num))
            exclusive_inputs = sorted(self.get_exclusive_inputs_signed(seg_num))
            if len(inputs) > 1 and inputs == exclusive_inputs:
                common_end = os.path.commonprefix([self.seq_from_signed_seg_num(x)[::-1]
                                                   for x in inputs])[::-1]
                common_end_len = len(common_end)
                if common_end_len > 0 and trim_amount_okay(inputs, common_end_len):
                    segment.prepend_to_forward_sequence(common_end)
                    for in_seg in inputs:
                        if in_seg > 0:
                            self.segments[in_seg].trim_from_end(common_end_len)
                        else:
                            self.segments[-in_seg].trim_from_start(common_end_len)

            outputs = sorted(self.get_downstream_seg_nums(seg_num))
            exclusive_outputs = sorted(self.get_exclusive_outputs_signed(seg_num))
            if len(outputs) > 1 and outputs == exclusive_outputs:
                common_start = os.path.commonprefix([self.seq_from_signed_seg_num(x)
                                                     for x in outputs])
                common_start_len = len(common_start)
                if common_start_len > 0 and trim_amount_okay(outputs, common_start_len):
                    segment.append_to_forward_sequence(common_start)
                    for out_seg in outputs:
                        if out_seg > 0:
                            self.segments[out_seg].trim_from_start(common_start_len)
                        else:
                            self.segments[-out_seg].trim_from_end(common_start_len)

    def starts_with_dead_end(self, signed_seg_num):
        """
        Returns whether or not the given segment (in the direction implied by the sign) is not
        lead into by any other segments.
        """
        if signed_seg_num not in self.reverse_links:
            return True
        return not self.reverse_links[signed_seg_num]

    def ends_with_dead_end(self, signed_seg_num):
        """
        Returns whether or not the given segment (in the direction implied by the sign) does not
        lead to any other segments.
        """
        if signed_seg_num not in self.forward_links:
            return True
        return not self.forward_links[signed_seg_num]

    def is_single_copy_or_bridge(self, seg_num, bridging_mode, single_copy_seg_nums):
        """
        Returns True if the given segment number is a single copy segment or a bridge segment. For
        this function, what counts as a 'single copy segment' depends on the bridging mode. At
        conservative bridging mode, only segments which are bridges or original single copy
        segments are allowed. At normal bridging mode, we also allow segments which have become
        single copy due to all but one of their copy depths being used up.
        """
        # Bridging mode level 2 (bold) means merge everything.
        if bridging_mode == 2 or single_copy_seg_nums is None:
            return True

        # Bridges are always okay to merge.
        if self.segments[seg_num].bridge is not None:
            return True

        # Original single copy segments are always okay to merge
        if seg_num in single_copy_seg_nums:
            return True

        # If the code got here, then the segment isn't a bridge or an original single copy. For
        # bridging mode level 0 (conservative), this is unmergeable.
        if bridging_mode == 0:
            return False

        # If the code got here, then the bridging mode is level 1 (normal). If the segment has
        # become single copy, then it's okay to merge.
        return seg_num in self.copy_depths and len(self.copy_depths[seg_num]) == 1

    def find_simple_two_way_junctions(self, valid_segments):
        """
        This function returns a list of segment numbers that are simple two-way junctions.
        Simple two-way junctions are defined as cases where two single copy segments join and then
        split. The single copy segments must be in the valid_segments list. Also, the single
        copy segments do not have to be unique, e.g. a loop can count as a two-way junction (the
        same segment is in inputs and outputs.
        """
        valid_seg_nums = set([x.number for x in valid_segments])
        simple_two_way_junctions = []
        for segment in self.segments.values():
            if self.get_copy_number(segment) != 2:
                continue
            seg_num = segment.number
            input_count = (len(self.reverse_links[seg_num])
                           if seg_num in self.reverse_links else 0)
            output_count = (len(self.forward_links[seg_num])
                            if seg_num in self.forward_links else 0)
            if input_count != 2 or output_count != 2:
                continue
            exclusive_input_count = len(self.get_exclusive_inputs(seg_num))
            exclusive_output_count = len(self.get_exclusive_outputs(seg_num))
            if exclusive_input_count != 2 or exclusive_output_count != 2:
                continue
            neighbour_segments = self.get_connected_segments(seg_num)
            if any(x not in valid_seg_nums for x in neighbour_segments):
                continue
            simple_two_way_junctions.append(seg_num)
        return simple_two_way_junctions

    def rotate_circular_sequences(self, shift_fraction=0.70710678118655):
        """
        Rotates the sequence to a new starting point. It shifts by a non-rational (well, almost)
        fraction of the sequence length so repeated executions of this function don't result in
        repeated starting positions.
        """
        completed_replicons = self.completed_circular_replicons()
        for completed_replicon in completed_replicons:
            segment = self.segments[completed_replicon]
            shift = int(segment.get_length() * shift_fraction)
            segment.rotate_sequence(shift, False) def determine_copy_depth(graph):
    """
    Assigns a copy depth to each segment in the graph.
    """
    # Reset any existing copy depths.
    graph.copy_depths = {}

    log.log_section_header('Determining graph multiplicity')
    log.log_explanation('Multiplicity is the number of times a '
                        'sequence occurs in the underlying sequence. Single-copy contigs '
                        '(those with a multiplicity of one, occurring only once in the underlying '
                        'sequence) are particularly useful.',
                        verbosity=1, extra_empty_lines_after=0)

    log.log_explanation('Multiplicity determination begins by identifying contigs which are '
                        'clearly single-copy because they are of low depth and do not have more '
                        'than one link per side.',
                        verbosity=2)

    single_copy_depth = graph.get_single_copy_depth()

    # Assign single copy status to segments within the tolerance of the single copy depth.
    # Also, if the graph has manually set multiplicity values (using the ML tag), we also use them
    # to assign single copy status.
    max_depth = single_copy_depth + settings.INITIAL_SINGLE_COPY_TOLERANCE
    initial_single_copy_segments = []
    for segment in sorted([x for x in graph.segments.values()],
                          key=lambda x: x.get_length(), reverse=True):
        num = segment.number
        depth = segment.depth
        if (depth <= max_depth and okay_for_initial_single_copy(graph, segment)) or \
                (num in graph.manual_multiplicity and graph.manual_multiplicity[num] == 1):
            graph.copy_depths[segment.number] = [segment.depth]
            initial_single_copy_segments.append(segment.number)

    if initial_single_copy_segments:
        log.log('\nInitial single copy segments:', 2)
        log.log_number_list(initial_single_copy_segments, 2)
    else:
        log.log('Initial single copy segments: none', 2)

    log.log('', verbosity=2)
    log.log_explanation('Unicycler now uses a greedy algorithm to propagate multiplicity '
                        'through the graph. For example, if two single-copy contigs merge '
                        'together, the resulting contig will get a multiplicity of two. '
                        'When no more propagation is possible, additional single-copy contigs '
                        'are added and the process is repeated. This allows for multiplicity to be '
                        'correctly assigned to the chromosome (at the median depth) but also for '
                        'plasmids (which may be higher or lower in depth).',
                        verbosity=2)

    # Propagate copy depth as much as possible using those initial assignments.
    copy_depth_table = [['Input', '', 'Output']]
    determine_copy_depth_part_2(graph, settings.COPY_PROPAGATION_TOLERANCE, copy_depth_table)

    # Assign single copy to the largest available segment, propagate and repeat.
    while True:
        assignments = assign_single_copy_depth(graph, settings.MIN_SINGLE_COPY_LENGTH,
                                               copy_depth_table)
        determine_copy_depth_part_2(graph, settings.COPY_PROPAGATION_TOLERANCE, copy_depth_table)
        if not assignments:
            break

    # Now propagate with no tolerance threshold to complete the remaining segments.
    if log.logger.stdout_verbosity_level >= 3:
        copy_depth_table.append(['REMOVING PROPAGATION TOLERANCE', '', ''])
    determine_copy_depth_part_2(graph, 1.0, copy_depth_table)

    print_table(copy_depth_table, alignments='RLL', max_col_width=999, hide_header=True,
                indent=0, col_separation=1, verbosity=2) def create_simple_long_read_bridges(graph, out_dir, keep, threads, read_dict, long_read_filename,
                                    scoring_scheme, anchor_segments):
    """
    Create and return simple long read bridges.
    """
    log.log_section_header('Creating simple long read bridges')
    log.log_explanation('Unicycler uses long read alignments (from minimap) to resolve simple '
                        'repeat structures in the graph. This takes care of some "low-hanging '
                        'fruit" of the graph simplification.')

    bridging_dir = os.path.join(out_dir, 'simple_bridging')
    if not os.path.exists(bridging_dir):
        os.makedirs(bridging_dir)
    minimap_alignments = align_long_reads_to_assembly_graph(graph, long_read_filename,
                                                            bridging_dir, threads)
    start_overlap_reads, end_overlap_reads = build_start_end_overlap_sets(minimap_alignments)
    bridges = simple_bridge_two_way_junctions(graph, start_overlap_reads, end_overlap_reads,
                                              minimap_alignments, anchor_segments)
    bridges += simple_bridge_loops(graph, start_overlap_reads, end_overlap_reads,
                                   minimap_alignments, read_dict, scoring_scheme, threads,
                                   anchor_segments)
    if keep < 3:
        shutil.rmtree(bridging_dir, ignore_errors=True)
    return bridges def make_miniasm_string_graph(graph, read_dict, long_read_filename, scoring_scheme, read_nicknames,
                              counter, args, anchor_segments, existing_long_read_assembly):
    log.log_section_header('Assembling contigs and long reads with miniasm')
    if graph is not None:
        log.log_explanation('Unicycler uses miniasm to construct a string graph assembly using '
                            'both the short read contigs and the long reads. It will then use the '
                            'resulting string graph to produce bridges between contigs. This '
                            'method requires decent coverage of long reads and therefore may not '
                            'be fruitful if long reads are sparse. However, it does not rely on '
                            'the short read assembly graph having good connectivity and is able to '
                            'bridge an assembly graph even when it contains many dead ends.',
                            extra_empty_lines_after=0)
        log.log_explanation('Unicycler uses two types of "reads" as assembly input: anchor contigs '
                            'from the short-read assembly and actual long reads which '
                            'overlap two or more of these contigs. It then assembles them with '
                            'miniasm.')
        if existing_long_read_assembly:
            log.log_explanation('Since you directly provided a long read assembly (using the '
                                '--existing_long_read_assembly option), Unicycler will use that '
                                'instead of performing its own miniasm/Racon assembly. However, '
                                'Unicycler still runs miniasm as this can provide useful contig '
                                'trimming information.')

    miniasm_dir = os.path.join(args.out, 'miniasm_assembly')
    if not os.path.exists(miniasm_dir):
        os.makedirs(miniasm_dir)

    short_reads_available = graph is not None
    seg_nums_to_bridge = set(x.number for x in anchor_segments)

    assembly_reads_filename = os.path.join(miniasm_dir, '01_assembly_reads.fastq')
    mappings_filename = os.path.join(miniasm_dir, '02_mappings.paf')
    string_graph_filename = os.path.join(miniasm_dir, '10_final_string_graph.gfa')
    branching_paths_removed_filename = os.path.join(miniasm_dir, '11_branching_paths_removed.gfa')
    unitig_graph_filename = os.path.join(miniasm_dir, '12_unitig_graph.gfa')
    racon_polished_filename = os.path.join(miniasm_dir, '13_racon_polished.gfa')
    contigs_placed_filename = os.path.join(miniasm_dir, '14_contigs_placed.gfa')
    miniasm_read_list = os.path.join(miniasm_dir, 'all_reads.txt')

    assembly_read_names = get_miniasm_assembly_reads(graph, read_dict, long_read_filename,
                                                     miniasm_dir, args.threads)

    # TO DO: identify chimeric reads and throw them out. This was part of miniasm, but it was
    # removed due to 'not working as intended', so I pulled it out of my miniasm as well.

    save_assembly_reads_to_file(assembly_reads_filename, assembly_read_names, read_dict, graph,
                                seg_nums_to_bridge)

    # Do an all-vs-all alignment of the assembly FASTQ, for miniasm input. Contig-contig
    # alignments are excluded (because single-copy contigs, by definition, should not
    # significantly overlap each other).
    log.log('Finding overlaps with minimap... ', end='')
    minimap_alignments_str = minimap_align_reads(assembly_reads_filename,
                                                 assembly_reads_filename,
                                                 args.threads, 0, 'read vs read')
    overlap_count = 0
    with open(mappings_filename, 'wt') as mappings:
        for minimap_alignment_str in line_iterator(minimap_alignments_str):
            if minimap_alignment_str.count('CONTIG_') < 2:
                mappings.write(minimap_alignment_str)
                mappings.write('\n')
                overlap_count += 1

    if overlap_count == 0:
        log.log(red('failed'))
    else:
        log.log(green('success'))
        log.log('  ' + int_to_str(overlap_count) + ' overlaps\n')

    # TO DO: refine these overlaps? Perhaps using Unicycler-align? I suspect that the quality
    # of a miniasm assembly is highly dependent on the input overlaps.
    #
    # I could even try to do more sophisticated stuff, like using the alignments to identify
    # repeat regions, then finding these repeat regions in reads and tossing out alignments
    # which are contained only in repeat regions.

    # TO DO: Unicycler's mode (conservative, normal or bold) should appropriately affect the
    # miniasm settings.

    # Now actually do the miniasm assembly, which will create a GFA file of the string graph.
    log.log('Assembling reads with miniasm... ', end='')
    min_depth = 3
    miniasm_assembly(assembly_reads_filename, mappings_filename, miniasm_dir, min_depth)
    if not os.path.isfile(string_graph_filename):
        log.log(red('failed'))
        raise MiniasmFailure('miniasm failed to generate a string graph')
    string_graph = StringGraph(string_graph_filename)

    if len(string_graph.segments) == 0:
        log.log(red('empty result'))
        unitig_graph = None

    else:
        log.log(green('success'))
        log.log('  ' + int_to_str(len(string_graph.segments)) + ' segments, ' +
                int_to_str(len(string_graph.links) // 2) + ' links\n')

        if not short_reads_available and args.keep > 0:
            shutil.copyfile(string_graph_filename,
                            gfa_path(args.out, next(counter), 'string_graph'))

        string_graph.remove_branching_paths()
        string_graph.save_to_gfa(branching_paths_removed_filename, include_depth=False)

        log.log('Merging segments into unitigs:')
        unitig_graph = merge_string_graph_segments_into_unitig_graph(string_graph,
                                                                     read_nicknames)
        circular_count = unitig_graph.get_circular_segment_count()
        linear_count = unitig_graph.get_linear_segment_count()
        unitig_graph_size = unitig_graph.get_total_segment_length()
        if circular_count > 0:
            log.log('  ' + int_to_str(circular_count) + ' circular unitig' +
                    ('' if circular_count == 1 else 's'))
        if linear_count > 0:
            log.log('  ' + int_to_str(linear_count) + ' linear unitig' +
                    ('' if linear_count == 1 else 's'))
        log.log('  total size = ' + int_to_str(unitig_graph_size) + ' bp')

        unitig_graph.save_to_gfa(unitig_graph_filename, include_depth=False)
        if not short_reads_available and args.keep > 0:
            unitig_graph.save_to_gfa(gfa_path(args.out, next(counter), 'unitig_graph'),
                                     include_depth=False)

        # If the miniasm assembly looks too small, then we don't bother polishing it or using
        # it for bridging.
        if short_reads_available:
            estimated_genome_size = graph.get_estimated_sequence_len()
            target_size = estimated_genome_size * \
                settings.REQUIRED_MINIASM_ASSEMBLY_SIZE_FOR_BRIDGING
            if unitig_graph_size < target_size:
                log.log('')
                log.log(red('miniasm assembly too small for bridging'))
                unitig_graph = None

        if unitig_graph is not None:
            # If the user supplied a long read assembly, then we use that instead of doing a
            # Racon polish.
            if existing_long_read_assembly:
                log.log('')
                log.log('Using provided long read assembly instead of running Racon: ' +
                        existing_long_read_assembly)
                log.log('')
                unitig_graph = StringGraph(existing_long_read_assembly)
            else:
                polish_unitigs_with_racon(unitig_graph, miniasm_dir, read_dict, graph,
                                          args.racon_path, args.threads, scoring_scheme,
                                          seg_nums_to_bridge)
                unitig_graph.save_to_gfa(racon_polished_filename)
                if not short_reads_available and args.keep > 0:
                    unitig_graph.save_to_gfa(gfa_path(args.out, next(counter),
                                                      'racon_polished'))
            if short_reads_available and args.keep > 0:
                unitig_graph.save_to_gfa(gfa_path(args.out, next(counter),
                                                  'long_read_assembly'))

    if unitig_graph is not None and short_reads_available:
        log.log('')
        trim_dead_ends_based_on_miniasm_trimming(graph, miniasm_read_list)
        unitig_graph = place_contigs(miniasm_dir, graph, unitig_graph, args.threads,
                                     scoring_scheme, seg_nums_to_bridge)
        unitig_graph.save_to_gfa(contigs_placed_filename, include_depth=False)

    if args.keep < 3:
        shutil.rmtree(miniasm_dir, ignore_errors=True)
    return unitig_graph def create_miniasm_bridges(graph, string_graph, anchor_segments, scoring_scheme, verbosity,
                           min_bridge_qual):
    """
    Makes bridges between single copy segments using the miniasm string graph.
    """
    log.log_section_header('Creating miniasm/Racon bridges')
    log.log_explanation('Now that the miniasm/Racon string graph is complete, Unicycler will '
                        'use it to build bridges between anchor segments.', verbosity=1)
    bridges = []
    anchor_seg_nums = set(x.number for x in anchor_segments)

    string_graph_bridge_segments = sorted([x for x in string_graph.segments
                                           if x.startswith('BRIDGE_') or
                                           x.startswith('OVERLAPPING_BRIDGE_')])

    filtered_string_graph_bridge_segments = []
    for bridge_seg_name in string_graph_bridge_segments:
        pos_seg_name = bridge_seg_name + '+'
        preceding_segments = string_graph.get_preceding_segments(pos_seg_name)
        following_segments = string_graph.get_following_segments(pos_seg_name)
        if len(preceding_segments) != 1:
            continue
        if len(following_segments) != 1:
            continue
        preceding_seg_name = preceding_segments[0]
        following_seg_name = following_segments[0]
        if not preceding_seg_name.startswith('CONTIG_') or \
                not following_seg_name.startswith('CONTIG_'):
            continue
        filtered_string_graph_bridge_segments.append(bridge_seg_name)

    # We want to display this table one row at a time, so we have to fix all of the column widths
    # at the start.
    bridge_count = len(filtered_string_graph_bridge_segments)
    alignments, col_widths = get_bridge_table_parameters(graph, bridge_count, verbosity,
                                                         'MiniasmBridge')
    print_bridge_table_header(alignments, col_widths, verbosity, 'MiniasmBridge')
    completed_count = 0

    for bridge_seg_name in filtered_string_graph_bridge_segments:
        bridge_seg = string_graph.segments[bridge_seg_name]
        pos_seg_name = bridge_seg_name + '+'
        preceding_seg_name = string_graph.get_preceding_segments(pos_seg_name)[0]
        following_seg_name = string_graph.get_following_segments(pos_seg_name)[0]

        first_link = string_graph.links[(preceding_seg_name, pos_seg_name)]
        second_link = string_graph.links[(pos_seg_name, following_seg_name)]
        preceding_seg_name = preceding_seg_name[7:]
        following_seg_name = following_seg_name[7:]
        preceding_segment_number = int(preceding_seg_name[:-1]) * \
            (1 if preceding_seg_name[-1] == '+' else -1)
        following_segment_number = int(following_seg_name[:-1]) * \
            (1 if following_seg_name[-1] == '+' else -1)
        assert abs(preceding_segment_number) in anchor_seg_nums
        assert abs(following_segment_number) in anchor_seg_nums

        start_overlap = first_link.seg_1_overlap
        end_overlap = second_link.seg_2_overlap
        output = []
        bridge = MiniasmBridge(graph, preceding_segment_number, following_segment_number,
                               bridge_seg.forward_sequence, start_overlap, end_overlap,
                               scoring_scheme, output)
        bridges.append(bridge)

        completed_count += 1
        print_bridge_table_row(alignments, col_widths, output, completed_count,
                               bridge_count, min_bridge_qual, verbosity, 'MiniasmBridge')

    # Now that the bridges are finalised, we split bridges that contain anchor segments in their
    # path such that all bridges start and end on an anchor segment but contain no anchor segments
    # in their path.
    split_bridges = []
    for bridge in bridges:
        if not bridge.graph_path:
            split_bridges.append(bridge)
        else:  # bridge has a path
            if not any(abs(x) in anchor_seg_nums for x in bridge.graph_path):  # already good
                split_bridges.append(bridge)

            # If the bridge path contains one or more anchor segments, it must be split!
            else:

                # A bridge with a path shouldn't have overlaps.
                assert bridge.start_overlap == 0
                assert bridge.end_overlap == 0

                full_path = [bridge.start_segment] + bridge.graph_path + [bridge.end_segment]
                anchor_indices = []
                for i, seg_num in enumerate(full_path):
                    if abs(seg_num) in anchor_seg_nums:
                        anchor_indices.append(i)
                anchor_indices = sorted(anchor_indices)
                for i in range(len(anchor_indices) - 1):
                    start_i = anchor_indices[i]
                    end_i = anchor_indices[i+1]
                    start_seg_num = full_path[start_i]
                    end_seg_num = full_path[end_i]
                    new_path = full_path[start_i+1:end_i]
                    bridge_sequence = graph.get_path_sequence(new_path)
                    split_bridge = MiniasmBridge(graph, start_seg_num, end_seg_num, bridge_sequence,
                                                 0, 0, scoring_scheme, [])
                    split_bridge.graph_path = new_path
                    split_bridge.all_paths = [new_path]
                    split_bridge.quality = bridge.quality
                    split_bridges.append(split_bridge)

    return bridges def create_long_read_bridges(graph, read_dict, read_names, anchor_segments, verbosity,
                             min_scaled_score, threads, scoring_scheme, min_alignment_length,
                             expected_linear_seqs, min_bridge_qual):
    """
    Makes bridges between single copy segments using the alignments in the long reads.
    """
    log.log_section_header('Building long read bridges')
    log.log_explanation('Unicycler uses the long read alignments to produce bridges between '
                        'anchor segments. These bridges can be formed using as few as one long '
                        'read, giving Unicycler the ability to bridge the graph even when '
                        'long-read depth is low.')

    anchor_seg_nums = set(x.number for x in anchor_segments)

    # This dictionary will collect the read sequences which span between two single copy segments.
    # Key = tuple of signed segment numbers (the segments being bridged)
    # Value = list of tuples containing the bridging sequence and the single copy segment
    #         alignments.
    spanning_read_seqs = defaultdict(list)

    for read_name in read_names:
        read = read_dict[read_name]
        alignments = get_single_copy_alignments(read, anchor_seg_nums, min_scaled_score)
        if len(alignments) < 2:
            continue

        # If the code got here, then we have some alignments to anchor segments. We grab
        # neighbouring pairs of alignments, starting with the highest scoring ones and work our
        # way down. This means that we should have a pair for each neighbouring alignment, but
        # potentially also more distant pairs if the alignments are strong.
        already_added = set()
        sorted_alignments = sorted(alignments, key=lambda x: x.raw_score, reverse=True)
        available_alignments = []
        for alignment in sorted_alignments:

            # If the alignment being added is to a reference that has already been added but in the
            # opposite direction, then we don't include it. E.g. we don't add an alignment for 10
            # if we already have an alignment for -10. This is because there's no legitimate way
            # for a single copy segment to appear in the same read in two different directions. The
            # same direction is okay, as that can happen with a circular piece of DNA, but opposite
            # directions implies multi-copy.
            opposite_num = -alignment.get_signed_ref_num()
            if opposite_num in set(x.get_signed_ref_num() for x in available_alignments):
                continue

            available_alignments.append(alignment)
            available_alignments = sorted(available_alignments,
                                          key=lambda x: x.read_start_positive_strand())
            if len(available_alignments) < 2:
                continue

            for i in range(len(available_alignments)):
                if i < len(available_alignments) - 1:
                    alignment_1 = available_alignments[i]
                    alignment_2 = available_alignments[i + 1]

                # Special case: when the first and last alignments are to the same graph segment,
                # make a bridge for them, even if they aren't a particularly high scoring pair of
                # alignments. This can help to circularise plasmids which are very tied up with
                # other, similar plasmids.
                elif available_alignments[0].ref.name == available_alignments[-1].ref.name:
                    alignment_1 = available_alignments[0]
                    alignment_2 = available_alignments[-1]
                else:
                    continue

                # Standardise the order so we don't end up with both directions (e.g. 5 to -6 and
                # 6 to -5) in spanning_read_seqs.
                seg_nums, flipped = flip_number_order(alignment_1.get_signed_ref_num(),
                                                      alignment_2.get_signed_ref_num())
                if seg_nums not in already_added:
                    bridge_start = alignment_1.read_end_positive_strand()
                    bridge_end = alignment_2.read_start_positive_strand()

                    if bridge_end > bridge_start:
                        bridge_seq = read.sequence[bridge_start:bridge_end]
                        bridge_qual = read.qualities[bridge_start:bridge_end]
                        if flipped:
                            bridge_seq = reverse_complement(bridge_seq)
                            bridge_qual = bridge_qual[::-1]
                    else:
                        bridge_seq = bridge_end - bridge_start  # 0 or a negative number
                        bridge_qual = ''

                    spanning_read_seqs[seg_nums].append((bridge_seq, bridge_qual, alignment_1,
                                                         alignment_2))
                    already_added.add(seg_nums)

    # If a bridge already exists for a spanning sequence, we add the sequence to the bridge. If
    # not, we create a new bridge and add it.
    new_bridges = []
    for seg_nums, span in spanning_read_seqs.items():
        start, end = seg_nums

        # If the start and end are the same and already exclusively connect (i.e. if this segment
        # is already circular), then skip - there's no need to bridge.
        if start == end and graph.get_downstream_seg_nums(start) == [start] and \
                graph.get_upstream_seg_nums(start) == [start]:
            continue

        new_bridge = LongReadBridge(graph, start, end)
        new_bridge.reads += span
        new_bridges.append(new_bridge)

    new_bridges = sorted(new_bridges, key=lambda x: (x.start_segment, x.end_segment))

    # During finalisation, we will compare the expected read count to the actual read count for
    # each bridge. To do this, we'll need the lengths of all reads (excluding those with no
    # alignments). We also need an estimate of the genome size.
    read_lengths = defaultdict(int)
    for read_name in read_names:
        read = read_dict[read_name]
        if read.alignments:
            read_lengths[read.get_length()] += 1
    estimated_genome_size = graph.get_estimated_sequence_len()

    # Now we need to finalise the bridges. This is the intensive step, as it involves creating a
    # consensus sequence, finding graph paths and doing alignments between the consensus and the
    # graph paths. We can therefore use threads to make this faster.
    num_long_read_bridges = len(new_bridges)

    # We want to display this table one row at a time, so we have to fix all of the column widths
    # at the start.
    alignments, col_widths = get_bridge_table_parameters(graph, num_long_read_bridges, verbosity,
                                                         'LongReadBridge')
    print_bridge_table_header(alignments, col_widths, verbosity, 'LongReadBridge')
    completed_count = 0

    # Use a simple loop if we only have one thread.
    if threads == 1:
        for bridge in new_bridges:
            output = bridge.finalise(scoring_scheme, min_alignment_length, read_lengths,
                                     estimated_genome_size, expected_linear_seqs)
            completed_count += 1
            print_bridge_table_row(alignments, col_widths, output, completed_count,
                                   num_long_read_bridges, min_bridge_qual, verbosity,
                                   'LongReadBridge')

    # Use a thread pool if we have more than one thread.
    else:
        pool = ThreadPool(threads)
        arg_list = []

        # Sort the bridges based on how long they're predicted to take to finalise. This will make
        # the big ones runs first which helps to more efficiently use the CPU cores.
        # E.g. if the biggest bridge was at the end, we'd be left waiting for it to finish with
        # only one core (bad), but if it was at the start, other work could be done in parallel.
        long_read_bridges = sorted(new_bridges, reverse=True,
                                   key=lambda x: x.predicted_time_to_finalise())
        for bridge in long_read_bridges:
            arg_list.append((bridge, scoring_scheme, min_alignment_length, read_lengths,
                             estimated_genome_size, expected_linear_seqs))
        for output in pool.imap_unordered(finalise_bridge, arg_list):
            completed_count += 1
            print_bridge_table_row(alignments, col_widths, output, completed_count,
                                   num_long_read_bridges, min_bridge_qual, verbosity,
                                   'LongReadBridge')

    # Now that the bridges are finalised, we split bridges that contain anchor segments in their
    # path such that all bridges start and end on an anchor segment but contain no anchor segments
    # in their path.
    split_bridges = []
    for bridge in new_bridges:
        if not bridge.graph_path:
            split_bridges.append(bridge)
        else:  # bridge has a path
            if not any(abs(x) in anchor_seg_nums for x in bridge.graph_path):  # already good
                split_bridges.append(bridge)

            # If the bridge path contains one or more anchor segments, it must be split!
            else:
                full_path = [bridge.start_segment] + bridge.graph_path + [bridge.end_segment]
                anchor_indices = []
                for i, seg_num in enumerate(full_path):
                    if abs(seg_num) in anchor_seg_nums:
                        anchor_indices.append(i)
                anchor_indices = sorted(anchor_indices)
                for i in range(len(anchor_indices) - 1):
                    start_i = anchor_indices[i]
                    end_i = anchor_indices[i+1]
                    start_seg_num = full_path[start_i]
                    end_seg_num = full_path[end_i]
                    new_path = full_path[start_i+1:end_i]
                    split_bridge = LongReadBridge(graph, start_seg_num, end_seg_num)
                    split_bridge.graph_path = new_path
                    split_bridge.all_paths = [new_path]
                    split_bridge.bridge_sequence = graph.get_path_sequence(new_path)
                    split_bridge.quality = bridge.quality
                    split_bridges.append(split_bridge)

    return split_bridges def create_spades_contig_bridges(graph, anchor_segments):
    """
    Builds graph bridges using the SPAdes contig paths.
    """
    log.log_section_header('Creating SPAdes contig bridges')
    log.log_explanation('SPAdes uses paired-end information to perform repeat resolution (RR) and '
                        'produce contigs from the assembly graph. SPAdes saves the graph paths '
                        'corresponding to these contigs in the contigs.paths file. When one of '
                        'these paths contains two or more anchor contigs, Unicycler can '
                        'create a bridge from the path.', verbosity=1)

    bridge_path_set = set()
    single_copy_numbers = [x.number for x in anchor_segments]
    for segment in anchor_segments:
        for path in graph.paths.values():
            flipped_path = [-x for x in reversed(path)]
            contig_bridges = find_contig_bridges(segment.number, path, single_copy_numbers)
            contig_bridges += find_contig_bridges(segment.number, flipped_path, single_copy_numbers)
            for contig_bridge in contig_bridges:
                flipped_contig_bridge = [-x for x in reversed(contig_bridge)]
                contig_bridge_str = ','.join([str(x) for x in contig_bridge])
                flipped_contig_bridge_str = ','.join([str(x) for x in flipped_contig_bridge])
                if contig_bridge_str not in bridge_path_set and \
                        flipped_contig_bridge_str not in bridge_path_set:
                    if contig_bridge[0] < 0 and contig_bridge[-1] < 0:
                        bridge_path_set.add(flipped_contig_bridge_str)
                    else:
                        bridge_path_set.add(contig_bridge_str)

    bridge_path_list = sorted(list([[int(y) for y in x.split(',')] for x in bridge_path_set]))

    # If multiple bridge paths start with or end with the same segment, that implies a conflict
    # between SPADes' paths and our single copy determination. Throw these bridges out.
    bridge_paths_by_start = {}
    bridge_paths_by_end = {}
    for path in bridge_path_list:
        start = path[0]
        end = path[-1]
        if start not in bridge_paths_by_start:
            bridge_paths_by_start[start] = []
        if end not in bridge_paths_by_end:
            bridge_paths_by_end[end] = []
        if -end not in bridge_paths_by_start:
            bridge_paths_by_start[-end] = []
        if -start not in bridge_paths_by_end:
            bridge_paths_by_end[-start] = []
        bridge_paths_by_start[start].append(path)
        bridge_paths_by_end[end].append(path)
        bridge_paths_by_start[-end].append(path)
        bridge_paths_by_end[-start].append(path)
    conflicting_paths = []
    for grouped_paths in bridge_paths_by_start.values():
        if len(grouped_paths) > 1:
            conflicting_paths += grouped_paths
    for grouped_paths in bridge_paths_by_end.values():
        if len(grouped_paths) > 1:
            conflicting_paths += grouped_paths
    conflicting_paths_no_dupes = []
    for path in conflicting_paths:
        if path not in conflicting_paths_no_dupes:
            conflicting_paths_no_dupes.append(path)
    conflicting_paths = conflicting_paths_no_dupes
    final_bridge_paths = [x for x in bridge_path_list if x not in conflicting_paths]

    bridges = [SpadesContigBridge(spades_contig_path=x, graph=graph) for x in final_bridge_paths]

    if bridges:
        bridge_table = [['Start', 'Path', 'End', 'Bridge quality']]
        spaded_arrow = ' ' + get_right_arrow() + ' '
        max_path_str_len = 5
        for bridge in bridges:
            path_str = spaded_arrow.join(str(seg) for seg in bridge.graph_path)
            max_path_str_len = max(max_path_str_len, len(path_str))
            bridge_table.append([str(bridge.start_segment), path_str, str(bridge.end_segment),
                                 float_to_str(bridge.quality, 1)])
        print_table(bridge_table, alignments='RCLR', left_align_header=False, indent=0,
                    fixed_col_widths=[5, max_path_str_len, 5, 7])
    else:
        log.log('No SPAdes contig bridges')

    return bridges def create_loop_unrolling_bridges(graph, anchor_segments):
    """
    This function creates loop unrolling bridges using the information in SPAdes paths.
    """
    log.log_section_header('Creating loop unrolling bridges')
    log.log_explanation('When a SPAdes contig path connects an anchor contig with the '
                        'middle contig of a simple loop, Unicycler concludes that the sequences '
                        'are contiguous (i.e. the loop is not a separate piece of DNA). It then '
                        'uses the read depth of the middle and repeat contigs to guess the number '
                        'of times to traverse the loop and makes a bridge.', verbosity=1)

    bridges = []
    simple_loops = graph.find_all_simple_loops()
    simple_loops = [x for x in simple_loops if x[2] is not None]
    seg_nums_to_bridge = set(x.number for x in anchor_segments)

    # A simple loop can either be caused by a repeat in one sequence (probably more typical) or by
    # a separate circular sequence which has some common sequence (less typical, but still very
    # possible: plasmids). We only want to unroll the former group, so we look for cases where the
    # loop's start or end is in a SPAdes contig path along with the middle. That implies that they
    # are on the same piece of DNA and can be unrolled.
    for start, end, middle, repeat in simple_loops:

        # We only want where the start and end are to-be-bridged segments but the middle is not.
        if abs(start) not in seg_nums_to_bridge:
            continue
        if abs(end) not in seg_nums_to_bridge:
            continue
        if abs(repeat) in seg_nums_to_bridge:
            continue

        joined = False
        for path in graph.paths.values():
            flipped_path = [-x for x in reversed(path)]
            if (start in path and middle in path) or \
                    (end in path and middle in path) or \
                    (start in flipped_path and middle in flipped_path) or \
                    (end in flipped_path and middle in flipped_path):
                joined = True
                break

        # If we've found evidence the simply loop is a single piece of DNA, then we'll make a loop
        # unrolling bridge!
        if joined:
            bridges.append(LoopUnrollingBridge(graph, start, end, middle, repeat))

    if bridges:
        bridge_table = [['Start', 'Repeat', 'Middle', 'End', 'Loop count by repeat',
                         'Loop count by middle', 'Loop count', 'Bridge quality']]
        for bridge in bridges:
            bridge_table.append([str(bridge.start_segment), str(bridge.repeat_segment),
                                 str(bridge.middle_segment), str(bridge.end_segment),
                                 float_to_str(bridge.loop_count_by_repeat, 2),
                                 float_to_str(bridge.loop_count_by_middle, 2),
                                 str(bridge.loop_count), float_to_str(bridge.quality, 1)])
        print_table(bridge_table, alignments='RRRRRRRR', left_align_header=False, indent=0,
                    fixed_col_widths=[5, 6, 6, 5, 10, 10, 5, 7])
    else:
        log.log('No loop unrolling bridges made')

    return bridges def int_to_str(num, max_num=0):
    """
    Converts a number to a string. Will add left padding based on the max value to ensure numbers
    align well.
    """
    if num is None:
        num_str = 'n/a'
    else:
        num_str = '{:,}'.format(num)
    max_str = '{:,}'.format(int(max_num))
    return num_str.rjust(len(max_str)) def float_to_str(num, decimals, max_num=0):
    """
    Converts a number to a string. Will add left padding based on the max value to ensure numbers
    align well. Will also add commas to numbers over 1000.
    """
    if decimals == 0:
        return int_to_str(int(round(num)), max_num=max_num)
    if num is None:
        num_str = 'n/a'
    else:
        num_str = '%.' + str(decimals) + 'f'
        num_str = num_str % num
        parts = num_str.split('.')
        before_decimal = parts[0]
        after_decimal = parts[1]
        num_str = int_to_str(int(before_decimal)) + '.' + after_decimal
    if max_num > 0:
        max_str = float_to_str(max_num, decimals)
        num_str = num_str.rjust(len(max_str))
    return num_str def quit_with_error(message):
    """
    Displays the given message and ends the program's execution.
    """
    log.log('Error: ' + message, 0, stderr=True)
    sys.exit(1) def get_percentile(unsorted_list, percentile):
    """
    Returns a percentile of a list of numbers. Doesn't assume the list has already been sorted.
    Implements the nearest rank method:
    https://en.wikipedia.org/wiki/Percentile#The_Nearest_Rank_method
    """
    return get_percentile_sorted(sorted(unsorted_list), percentile) def bold(text):
    return BOLD + text + END_FORMATTING def check_input_files(args):
    """
    Checks to make sure all files in the list are present and either program, as needed.
    """
    if args.short1:
        check_file_exists(args.short1)
    if args.short2:
        check_file_exists(args.short2)
    if args.unpaired:
        check_file_exists(args.unpaired)
    if args.long:
        check_file_exists(args.long)
    if args.existing_long_read_assembly:
        check_file_exists(args.existing_long_read_assembly)
    if args.short1 and args.short2 and args.short1 == args.short2:
        quit_with_error('first and second read pair files cannot be the same file') class MyHelpFormatter(argparse.HelpFormatter):
    """
    This is a custom formatter class for argparse. It allows for some custom formatting,
    in particular for the help texts with multiple options (like bridging mode and verbosity level).
    http://stackoverflow.com/questions/3853722
    """
    def __init__(self, prog):
        terminal_width = shutil.get_terminal_size().columns
        os.environ['COLUMNS'] = str(terminal_width)
        max_help_position = min(max(24, terminal_width // 3), 40)
        try:
            self.colours = int(subprocess.check_output(['tput', 'colors']).decode().strip())
        except (ValueError, subprocess.CalledProcessError, FileNotFoundError, AttributeError):
            self.colours = 1
        super().__init__(prog, max_help_position=max_help_position)

    def _get_help_string(self, action):
        """
        Override this function to add default values, but only when 'default' is not already in the
        help text.
        """
        help_text = action.help
        if action.default != argparse.SUPPRESS and 'default' not in help_text.lower() and \
                action.default is not None:
            help_text += ' (default: ' + str(action.default) + ')'
        return help_text

    def start_section(self, heading):
        """
        Override this method to add bold underlining to section headers.
        """
        if self.colours > 1:
            heading = BOLD + heading + END_FORMATTING
        super().start_section(heading)

    def _split_lines(self, text, width):
        """
        Override this method to add special behaviour for help texts that start with:
          'B|' - loop text to the column of the equals sign
          'R|' - loop text one option per line
        """
        if text.startswith('B|') or text.startswith('R|'):
            text_lines = text[2:].splitlines()
            wrapped_text_lines = []
            for line in text_lines:
                if len(line) <= width:
                    wrapped_text_lines.append(line)
                else:
                    wrap_column = 2
                    if text.startswith('B|'):
                        line_parts = line.split()
                        wrap_column += line.find('=')
                        join = ''
                        current_line = '  ' + line_parts[0]
                    else:  # text.startswith('R|')
                        line_parts = line.split(', ')
                        join = ','
                        current_line = line_parts[0]
                    for part in line_parts[1:]:
                        if len(current_line) + len(join) + 1 + len(part) <= width:
                            current_line += join + ' ' + part
                        else:
                            wrapped_text_lines.append(current_line + join)
                            current_line = ' ' * wrap_column + part
                    wrapped_text_lines.append(current_line)
            return wrapped_text_lines
        else:
            return argparse.HelpFormatter._split_lines(self, text, width)

    def _fill_text(self, text, width, indent):
        if text.startswith('R|'):
            return ''.join(indent + line for line in text[2:].splitlines(keepends=True))
        else:
            return argparse.HelpFormatter._fill_text(self, text, width, indent)

    def _format_action(self, action):
        """
        Override this method to make help descriptions dim.
        """
        # determine the required width and the entry label
        help_position = min(self._action_max_length + 2,
                            self._max_help_position)
        help_width = self._width - help_position
        action_width = help_position - self._current_indent - 2
        action_header = self._format_action_invocation(action)

        # ho nelp; start on same line and add a final newline
        if not action.help:
            tup = self._current_indent, '', action_header
            action_header = '%*s%s\n' % tup
            indent_first = 0

        # short action name; start on the same line and pad two spaces
        elif len(action_header) <= action_width:
            tup = self._current_indent, '', action_width, action_header
            action_header = '%*s%-*s  ' % tup
            indent_first = 0

        # long action name; start on the next line
        else:
            tup = self._current_indent, '', action_header
            action_header = '%*s%s\n' % tup
            indent_first = help_position

        # collect the pieces of the action help
        parts = [action_header]

        # if there was help for the action, add lines of help text
        if action.help:
            help_text = self._expand_help(action)
            help_lines = self._split_lines(help_text, help_width)
            first_line = help_lines[0]
            if self.colours > 8:
                first_line = DIM + first_line + END_FORMATTING
            parts.append('%*s%s\n' % (indent_first, '', first_line))
            for line in help_lines[1:]:
                if self.colours > 8:
                    line = DIM + line + END_FORMATTING
                parts.append('%*s%s\n' % (help_position, '', line))

        # or add a newline if the description doesn't end with one
        elif not action_header.endswith('\n'):
            parts.append('\n')

        # if there are any sub-actions, add their help as well
        for subaction in self._iter_indented_subactions(action):
            parts.append(self._format_action(subaction))

        # return a single string
        return self._join_parts(parts) def print_table(table, alignments='', max_col_width=30, col_separation=3, indent=2,
                row_colour=None, sub_colour=None, row_extra_text=None, leading_newline=False,
                subsequent_indent='', return_str=False, header_format='underline',
                hide_header=False, fixed_col_widths=None, left_align_header=True,
                bottom_align_header=True, verbosity=1):
    """
    Args:
        table: a list of lists of strings (one row is one list, all rows should be the same length)
        alignments: a string of L and R, indicating the alignment for each row
        max_col_width: values longer than this will be wrapped
        col_separation: the number of spaces between columns
        indent: the number of spaces between the table and the left side of the terminal
        row_colour: a dictionary of row indices and their colour names
        sub_colour: a dictionary of values to colour names for which the text colour will be set
        row_extra_text: a dictionary of row indices and extra text to display after the row
        leading_newline: if True, the function will print a blank line above the table
        subsequent_indent: this string will be added to the start of wrapped text lines
        return_str: if True, this function will return a string of the table instead of printing it
        header_format: the formatting (colour, underline, etc) of the header line
        hide_header: if True, the header is not printed
        fixed_col_widths: a list to specify exact column widths (automatic if not used)
        left_align_header: if False, the header will follow the column alignments
        bottom_align_header: if False, the header will align to the top, like other rows
        verbosity: the table will only be logged if the logger verbosity is >= this value
    """
    column_count = len(table[0])
    table = [x[:column_count] for x in table]
    table = [x + [''] * (column_count - len(x)) for x in table]
    if row_colour is None:
        row_colour = {}
    if sub_colour is None:
        sub_colour = {}
    if row_extra_text is None:
        row_extra_text = {}
    if leading_newline:
        log.log('', verbosity)

    # Ensure the alignments string is the same length as the column count
    alignments += 'L' * (column_count - len(alignments))
    alignments = alignments[:column_count]

    if fixed_col_widths is not None:
        col_widths = fixed_col_widths
    else:
        col_widths = [0] * column_count
        for row in table:
            col_widths = [min(max(col_widths[i], len_without_format(x)), max_col_width)
                          for i, x in enumerate(row)]
    separator = ' ' * col_separation
    indenter = ' ' * indent
    full_table_str = ''
    for i, row in enumerate(table):
        row = [str(x) for x in row]
        if hide_header and i == 0:
            continue

        if fixed_col_widths is not None:
            wrapped_row = []
            for col, fixed_width in zip(row, fixed_col_widths):
                wrapper = textwrap.TextWrapper(subsequent_indent=subsequent_indent,
                                               width=fixed_width)
                wrapped_row.append(wrapper.wrap(col))
        else:
            wrapper = textwrap.TextWrapper(subsequent_indent=subsequent_indent, width=max_col_width)
            wrapped_row = [wrapper.wrap(x) for x in row]

        row_rows = max(len(x) for x in wrapped_row)
        if i == 0 and bottom_align_header:
            wrapped_row = [[''] * (row_rows - len(x)) + x for x in wrapped_row]

        for j in range(row_rows):
            row_line = [x[j] if j < len(x) else '' for x in wrapped_row]
            aligned_row = []
            for value, col_width, alignment in zip(row_line, col_widths, alignments):
                if alignment == 'L' or (i == 0 and left_align_header):
                    aligned_row.append(value.ljust(col_width))
                elif alignment == 'C':
                    aligned_row.append(value.center(col_width))
                else:
                    aligned_row.append(value.rjust(col_width))
            row_str = separator.join(aligned_row)
            if i in row_extra_text:
                row_str += row_extra_text[i]
            if i == 0 and header_format:
                row_str = colour(row_str, header_format)
            if i in row_colour:
                row_str = colour(row_str, row_colour[i])
            for text, colour_name in sub_colour.items():
                row_str = row_str.replace(text, colour(text, colour_name))
            if j < row_rows - 1 and UNDERLINE in row_str:
                row_str = re.sub('\033\[4m', '', row_str)
            if return_str:
                full_table_str += indenter + row_str + '\n'
            else:
                log.log(indenter + row_str, verbosity)
    if return_str:
        return full_table_str def get_ascii_art():
    ascii_art = (bold_red("       __\n") +
                 bold_red("       \ \___\n") +
                 bold_red("        \ ___\\\n") +
                 bold_red("        //\n") +
                 bold_red("   ____//      ") +
                 bold_yellow("_    _         _                     _\n") +
                 bold_red(" //_  //\\\\    ") +
                 bold_yellow("| |  | |       |_|                   | |\n") +
                 bold_red("//  \\//  \\\\   ") +
                 bold_yellow("| |  | | _ __   _   ___  _   _   ___ | |  ___  _ __\n") +
                 bold_red("||  (O)  ||   ") +
                 bold_yellow("| |  | || '_ \ | | / __|| | | | / __|| | / _ \| '__|\n") +
                 bold_red("\\\\    \_ //   ") +
                 bold_yellow("| |__| || | | || || (__ | |_| || (__ | ||  __/| |\n") +
                 bold_red(" \\\\_____//     ") +
                 bold_yellow("\____/ |_| |_||_| \___| \__, | \___||_| \___||_|\n") +
                 bold_yellow("                                        __/ |\n") +
                 bold_yellow("                                       |___/"))
    return ascii_art def get_default_thread_count():
    return min(multiprocessing.cpu_count(), settings.MAX_AUTO_THREAD_COUNT) def spades_path_and_version(spades_path):
    found_spades_path = shutil.which(spades_path)
    if found_spades_path is None:
        return spades_path, '', 'not found'
    command = [found_spades_path, '-v']
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
    out, _ = process.communicate()
    out = out.decode()

    if not out or 'Verification of expression' in out:
        return found_spades_path, '', 'bad'
    version = spades_version_from_spades_output(out)

    if 'python version' in out and 'is not supported' in out:
        return found_spades_path, '', 'Python problem'

    # Make sure SPAdes is 3.13.1 or later
    try:
        status = spades_status_from_version(version)
    except (ValueError, IndexError):
        version, status = '?', 'too old'

    return found_spades_path, version, status def makeblastdb_path_and_version(makeblastdb_path):
    found_makeblastdb_path = shutil.which(makeblastdb_path)
    if found_makeblastdb_path is None:
        return makeblastdb_path, '', 'not found'
    command = [found_makeblastdb_path, '-version']
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
    out, _ = process.communicate()
    version = out.decode().split('makeblastdb: ')[-1].split()[0]
    try:
        int(version.split('.')[0]), int(version.split('.')[1])
    except (ValueError, IndexError):
        version, status = '?', 'too old'
    return found_makeblastdb_path, version, 'good' def tblastn_path_and_version(tblastn_path):
    found_tblastn_path = shutil.which(tblastn_path)
    if found_tblastn_path is None:
        return tblastn_path, '', 'not found'
    command = [found_tblastn_path, '-version']
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
    out, _ = process.communicate()
    version = out.decode().split('tblastn: ')[-1].split()[0]
    try:
        int(version.split('.')[0]), int(version.split('.')[1])
    except (ValueError, IndexError):
        version, status = '?', 'too old'
    return found_tblastn_path, version, 'good' def racon_path_and_version(racon_path):
    found_racon_path = shutil.which(racon_path)
    if found_racon_path is None:
        return racon_path, '', 'not found'
    command = [found_racon_path]
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
    out, _ = process.communicate()
    out = out.decode().lower()
    if 'racon' not in out or 'options' not in out:
        return found_racon_path, '-', 'bad'

    return found_racon_path, racon_version(found_racon_path), 'good' def gfa_path(out_dir, file_num, name):
    return os.path.join(out_dir, str(file_num).zfill(3) + '_' + name + '.gfa') def red(text):
    return RED + text + END_FORMATTING def get_best_spades_graph(short1, short2, short_unpaired, out_dir, read_depth_filter, verbosity,
                          spades_path, threads, keep, kmer_count, min_k_frac, max_k_frac, kmers,
                          expected_linear_seqs, largest_component, spades_graph_prefix,
                          spades_options):
    """
    This function tries a SPAdes assembly at different k-mers and returns the best one.
    """
    spades_dir = os.path.join(out_dir, 'spades_assembly')
    if not os.path.exists(spades_dir):
        os.makedirs(spades_dir)

    threads = min(threads, 32)  # SPAdes can possibly crash if given too many threads.
    check_fastqs(short1, short2, short_unpaired)
    reads = (short1, short2, short_unpaired)

    kmer_range = get_kmer_range(kmers, short1, short2, short_unpaired, spades_dir, kmer_count,
                                min_k_frac, max_k_frac, spades_path)

    log.log_section_header('SPAdes assemblies')
    log.log_explanation('Unicycler now uses SPAdes to assemble the short reads. It scores the '
                        'assembly graph for each k-mer using the number of contigs (fewer is '
                        'better) and the number of dead ends (fewer is better). The score '
                        'function is 1/(c*(d+2)), where c is the contig count and d is the '
                        'dead end count.')

    if verbosity > 1:
        spades_results_table = [['K-mer', 'Contigs', 'Links', 'Total length', 'N50',
                                 'Longest contig', 'Dead ends', 'Score']]
    else:
        spades_results_table = [['K-mer', 'Contigs', 'Dead ends', 'Score']]

    graph_files, insert_size_mean, insert_size_deviation = \
        run_spades_all_kmers(reads, spades_dir, kmer_range, threads, spades_path,
                             spades_graph_prefix, spades_options)

    existing_graph_files = [x for x in graph_files if x is not None]
    if not existing_graph_files:
        quit_with_error('SPAdes failed to produce assemblies. '
                        'See spades_assembly/spades.log for more info.')
    median_segment_count = statistics.median(count_segments_in_gfa(x)
                                             for x in existing_graph_files)

    best_score, best_kmer, best_graph_filename = 0.0, 0, ''
    for graph_file, kmer in zip(graph_files, kmer_range):
        table_line = [int_to_str(kmer)]

        if graph_file is None:
            table_line += [''] * (7 if verbosity > 1 else 2)
            table_line.append('failed')
            spades_results_table.append(table_line)
            continue

        assembly_graph = AssemblyGraph(graph_file, kmer, insert_size_mean=insert_size_mean,
                                       insert_size_deviation=insert_size_deviation)

        # If this graph has way too many segments, then we will just skip it because very complex
        # graphs take forever to clean up.
        # TO DO: I can remove this awkward hack if I make the graph cleaning more efficient.
        if len(assembly_graph.segments) > 4 * median_segment_count:
            table_line += [''] * (6 if verbosity > 1 else 2)
            table_line.append('too complex')
            spades_results_table.append(table_line)
            continue

        log.log('\nCleaning k{} graph'.format(kmer), 2)
        assembly_graph.clean(read_depth_filter, largest_component)
        clean_graph_filename = os.path.join(spades_dir, ('k%03d' % kmer) + '_assembly_graph.gfa')
        assembly_graph.save_to_gfa(clean_graph_filename, verbosity=2)

        segment_count = len(assembly_graph.segments)
        dead_ends = assembly_graph.total_dead_end_count()

        # If the user is expecting some linear sequences, then the dead end count can be adjusted
        # down so expected dead ends don't penalise this k-mer.
        adjusted_dead_ends = max(0, dead_ends - (2 * expected_linear_seqs))
        if segment_count == 0:
            score = 0.0
        else:
            score = 1.0 / (segment_count * (adjusted_dead_ends + 2))

        # Prepare the table line for this k-mer graph.
        table_line += [int_to_str(segment_count)]
        if verbosity > 1:
            n50, shortest, _, median, _, longest = assembly_graph.get_contig_stats()
            table_line += [int_to_str(assembly_graph.get_total_link_count()),
                           int_to_str(assembly_graph.get_total_length()),
                           int_to_str(n50), int_to_str(longest)]
        table_line += [int_to_str(dead_ends), '{:.2e}'.format(score)]
        spades_results_table.append(table_line)

        if score > best_score:
            best_kmer, best_score, best_graph_filename = kmer, score, graph_file

    log.log('', 2)

    if not best_kmer:
        quit_with_error('none of the SPAdes graphs were suitable for scaffolding in Unicycler')

    assembly_graph = AssemblyGraph(best_graph_filename, best_kmer,
                                   insert_size_mean=insert_size_mean,
                                   insert_size_deviation=insert_size_deviation)
    removed_count, removed_length = assembly_graph.clean(read_depth_filter, largest_component)
    clean_graph_filename = os.path.join(spades_dir, 'k' + str(best_kmer) + '_assembly_graph.gfa')
    assembly_graph.save_to_gfa(clean_graph_filename, verbosity=2)

    if best_score == 0.0:
        quit_with_error('none of the SPAdes assemblies produced assembled sequence')

    # Print the SPAdes result table, highlighting the best k-mer in green.
    log.log_section_header('SPAdes assembly graph summary', 2)
    best_kmer_row = [x[0] for x in spades_results_table].index(int_to_str(best_kmer))
    print_table(spades_results_table, alignments='RRRRRRRR', indent=0,
                row_colour={best_kmer_row: 'green'},
                row_extra_text={best_kmer_row: ' ' + get_left_arrow() + 'best'})

    # Report on the results of the read depth filter (can help with identifying levels of
    # contamination).
    log.log('\nRead depth filter: removed {} contigs totalling {} bp'.format(removed_count,
                                                                             removed_length))
    # Clean up.
    if keep < 1:
        for g in graph_files:
            if os.path.isfile(g):
                log.log('Deleting ' + g)
                os.remove(g)
    if keep < 3 and os.path.isdir(spades_dir):
        log.log('Deleting ' + spades_dir + '/')
        shutil.rmtree(spades_dir, ignore_errors=True)

    return assembly_graph def find_start_gene(sequence, start_genes_fasta, identity_threshold, coverage_threshold, blast_dir,
                    makeblastdb_path, tblastn_path):
    """
    This function uses tblastn to look for start genes in the sequence. It returns the first gene
    (using the order in the file) which meets the identity and coverage thresholds, as well as
    the position of that gene (including which strand it is on).
    This function assumes that the sequence is circular with no overlap.
    """
    # Prepare the replicon sequence. In order to get a solid, single BLAST hit in cases where the
    # gene overlaps from the end to the start, we have to duplicate some of the replicon sequence
    # for the BLAST database.
    seq_len = len(sequence)
    start_genes_fasta = os.path.abspath(start_genes_fasta)
    queries = load_fasta(start_genes_fasta)
    if not queries:
        raise CannotFindStart
    longest_query = max(len(x[1]) for x in queries)
    longest_query *= 3  # amino acids to nucleotides
    dup_length = min(seq_len, longest_query)
    sequence = sequence + sequence[:dup_length]

    # BLAST has serious issues with paths that contain spaces. This page explains some of it:
    #   https://www.ncbi.nlm.nih.gov/books/NBK279669/
    # But I couldn't make it all work for makeblastdb (spaces made it require -out, and it never
    # accepted spaces in the -out path, no matter how I used quotes). So we will just move into the
    # temporary directory to run the BLAST commands.
    starting_dir = os.getcwd()
    os.chdir(blast_dir)

    # Create a FASTA file of the replicon sequence.
    replicon_fasta_filename = 'replicon.fasta'
    replicon_fasta = open(replicon_fasta_filename, 'w')
    replicon_fasta.write('>replicon\n')
    replicon_fasta.write(sequence)
    replicon_fasta.write('\n')
    replicon_fasta.close()

    # Build the BLAST database.
    command = [makeblastdb_path, '-dbtype', 'nucl', '-in', replicon_fasta_filename]
    log.log('  ' + ' '.join(command), 2)
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    _, err = process.communicate()
    if err:
        log.log('\nmakeblastdb encountered an error:\n' + err.decode())
        os.chdir(starting_dir)
        raise CannotFindStart

    # Run the tblastn search.
    command = [tblastn_path, '-db', replicon_fasta_filename, '-query', start_genes_fasta, '-outfmt',
               '6 qseqid sstart send pident qlen qseq qstart bitscore', '-num_threads', '1']
    log.log('  ' + ' '.join(command), 2)
    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    blast_out, blast_err = process.communicate()
    process.wait()
    if blast_err:
        log.log('\nBLAST encountered an error:\n' + blast_err.decode())

    # Find the best hit in the results.
    best_hit, best_bitscore = None, 0
    for line in blast_out.decode().splitlines():
        hit = BlastHit(line, seq_len)
        if hit.pident >= identity_threshold and hit.query_cov >= coverage_threshold and \
                hit.qstart == 0 and hit.bitscore > best_bitscore:
            best_hit = hit
            best_bitscore = hit.bitscore

    os.chdir(starting_dir)

    if best_bitscore:
        return best_hit
    else:
        raise CannotFindStart class CannotFindStart(Exception):
    pass SAM_WRITE_LOCK = threading.Lock()
VERBOSITY = 0
        VERBOSITY = verbosity
def fix_up_arguments(args):
def semi_global_align_long_reads(references, ref_fasta, read_dict, read_names, reads_fastq,
                                 threads, scoring_scheme, low_score_threshold_list, keep_bad,
                                 min_align_length, sam_filename, full_command, allowed_overlap,
                                 sensitivity_level, contamination_fasta, verbosity=None,
                                 stdout_header='Aligning reads', display_low_score=True,
                                 single_copy_segment_names=None):
def get_percent_contamination(read_dict):
def print_alignment_summary_table(read_dict, verbosity, using_contamination):
def load_sam_alignments(sam_filename, read_dict, reference_dict, scoring_scheme):
def seqan_alignment_one_arg(all_args):
def seqan_alignment(read, reference_dict, scoring_scheme, ref_seqs_ptr, low_score_threshold,
                    keep_bad, min_align_length, sam_filename, allowed_overlap,
                    minimap_alignments, sensitivity_level, single_copy_segment_names):
def group_reads_by_fraction_aligned(read_dict):
def get_auto_score_threshold(scoring_scheme, std_devs_over_mean): def get_read_nickname_dict(read_names):
    """
    Read names can be quite long, so for the sake of output brevity, this function tries to come
    up with some shorter nicknames for the reads.
    """
    max_read_name_len = max(len(name) for name in read_names)
    for nickname_length in range(1, max_read_name_len):
        nicknames = set()
        for name in read_names:
            nickname = name[:nickname_length]
            if nickname in nicknames:
                break
            nicknames.add(nickname)
        else:
            return {name: name[:nickname_length] for name in read_names}

    # If we couldn't find a length for shorter nicknames, then the nicknames are just the full
    # names. Oh well.
    return {name: name for name in read_names} def load_long_reads(filename, silent=False, section_header='Loading reads', output_dir=None):
    """
    This function loads in long reads from a FASTQ file and returns a dictionary where key = read
    name and value = Read object. It also returns a list of read names, in the order they are in
    the file.
    """
    # Read files can be either FASTA or FASTQ and optionally gzipped.
    try:
        file_type = get_sequence_file_type(filename)
    except ValueError:
        file_type = ''
        quit_with_error(filename + ' is not in either FASTA or FASTQ format')
    if get_compression_type(filename) == 'gz':
        open_func = gzip.open
    else:  # plain text
        open_func = open

    if not silent:
        log.log_section_header(section_header)

    read_dict = {}
    read_names = []
    total_bases = 0
    last_progress = 0.0
    step = settings.LOADING_READS_PROGRESS_STEP
    duplicate_read_names_found = False

    if file_type == 'FASTQ':
        with open_func(filename, 'rt') as fastq:
            num_reads = sum(1 for _ in fastq) // 4
    else:  # file_type == 'FASTA'
        with open_func(filename, 'rt') as fasta:
            num_reads = sum(1 for line in fasta if line.startswith('>'))
    if not num_reads:
        quit_with_error('There are no read sequences in ' + filename)
    if not silent:
        log.log_progress_line(0, num_reads)

    if file_type == 'FASTQ':
        with open_func(filename, 'rt') as fastq:
            for line in fastq:
                stripped_line = line.strip()
                if len(stripped_line) == 0:
                    continue
                if not stripped_line.startswith('@'):
                    continue
                original_name = stripped_line[1:].split()[0]
                sequence = next(fastq).strip()
                _ = next(fastq)
                qualities = next(fastq).strip()

                # Don't allow duplicate read names, so add a trailing number when they occur.
                name = original_name
                duplicate_name_number = 1
                while name in read_dict:
                    duplicate_read_names_found = True
                    duplicate_name_number += 1
                    name = original_name + '_' + str(duplicate_name_number)

                read_dict[name] = Read(name, sequence, qualities)
                read_names.append(name)
                total_bases += len(sequence)
                progress = 100.0 * len(read_dict) / num_reads
                progress_rounded_down = math.floor(progress / step) * step
                if progress == 100.0 or progress_rounded_down > last_progress:
                    if not silent:
                        log.log_progress_line(len(read_dict), num_reads, total_bases)
                    last_progress = progress_rounded_down

    else:  # file_type == 'FASTA'
        with open_func(filename, 'rt') as fasta:
            name = ''
            sequence = ''
            last_progress = 0.0
            for line in fasta:
                line = line.strip()
                if not line:
                    continue
                if line.startswith('>'):  # Header line = start of new contig
                    if name:
                        read_dict[name] = Read(name, sequence, None)
                        read_names.append(name)
                        total_bases += len(sequence)
                        progress = 100.0 * len(read_dict) / num_reads
                        progress_rounded_down = math.floor(progress / step) * step
                        if progress == 100.0 or progress_rounded_down > last_progress:
                            if not silent:
                                log.log_progress_line(len(read_dict), num_reads, total_bases)
                            last_progress = progress_rounded_down
                        sequence = ''

                    # Don't allow duplicate read names, so add a trailing number when they occur.
                    original_name = get_nice_header(line[1:])
                    name = original_name
                    duplicate_name_number = 1
                    while name in read_dict:
                        duplicate_read_names_found = True
                        duplicate_name_number += 1
                        name = original_name + '_' + str(duplicate_name_number)

                else:
                    sequence += line
            if name:
                read_dict[name] = Read(name, sequence, None)
                read_names.append(name)
                total_bases += len(sequence)
                if not silent:
                    log.log_progress_line(len(read_dict), num_reads, total_bases)

    if not silent:
        log.log_progress_line(len(read_dict), len(read_dict), total_bases, end_newline=True)

    # If there were duplicate read names, then we save the reads back out to file with their fixed
    # names.
    if duplicate_read_names_found:
        no_dup_filename = strip_read_extensions(filename) + '_no_duplicates'
        if file_type == 'FASTQ':
            no_dup_filename += '.fastq.gz'
        else:  # file_type == 'FASTA'
            no_dup_filename += '.fasta.gz'

        # If an output directory was provided, we put the no duplicate read file there.
        if output_dir is not None:
            no_dup_filename = os.path.join(output_dir, no_dup_filename)

        # If an output directory wasn't provided, we put the no duplicate read file in the same
        # directory as the normal read file.
        else:
            no_dup_filename = os.path.join(os.path.dirname(os.path.abspath(filename)),
                                           no_dup_filename)

        if not silent:
            log.log('\nDuplicate read names found. Saving duplicate-free file:')
            log.log(no_dup_filename)
        with gzip.open(no_dup_filename, 'wb') as f:
            for read_name in read_names:
                read = read_dict[read_name]
                if file_type == 'FASTQ':
                    f.write(read.get_fastq().encode())
                else:  # file_type == 'FASTA'
                    f.write(read.get_fasta().encode())

    else:
        no_dup_filename = filename

    return read_dict, read_names, no_dup_filename
import argparse
import os
import sys
import shutil
import random
import itertools
import multiprocessing
from .alignment import AlignmentScoringScheme
from .assembly_graph import AssemblyGraph
from .assembly_graph_copy_depth import determine_copy_depth
from .bridge_long_read_simple import create_simple_long_read_bridges
from .miniasm_assembly import make_miniasm_string_graph
from .bridge_miniasm import create_miniasm_bridges
from .bridge_long_read import create_long_read_bridges
from .bridge_spades_contig import create_spades_contig_bridges
from .bridge_loop_unroll import create_loop_unrolling_bridges
from .misc import int_to_str, float_to_str, quit_with_error, get_percentile, bold, \
    check_input_files, MyHelpFormatter, print_table, get_ascii_art, \
    get_default_thread_count, spades_path_and_version, makeblastdb_path_and_version, \
    tblastn_path_and_version, racon_path_and_version, gfa_path, red
from .spades_func import get_best_spades_graph
from .blast_func import find_start_gene, CannotFindStart
from .unicycler_align import fix_up_arguments, semi_global_align_long_reads, load_references, \
    load_sam_alignments, print_alignment_summary_table
from .read_ref import get_read_nickname_dict, load_long_reads
from . import log
from . import settings
from .version import __version__
    else:
        racon_path, racon_version, racon_status = racon_path_and_version(args.racon_path)
    racon_row = ['racon', racon_version, racon_status]
    if args.verbosity > 1:
        racon_row.append(racon_path)
    program_table.append(racon_row)

    # Rotation dependencies
    if args.no_rotate:
        makeblastdb_path, makeblastdb_version, makeblastdb_status = '', '', 'not used'
        tblastn_path, tblastn_version, tblastn_status = '', '', 'not used'
    else:
        makeblastdb_path, makeblastdb_version, makeblastdb_status = \
            makeblastdb_path_and_version(args.makeblastdb_path)
        tblastn_path, tblastn_version, tblastn_status = tblastn_path_and_version(args.tblastn_path)
    makeblastdb_row = ['makeblastdb', makeblastdb_version, makeblastdb_status]
    tblastn_row = ['tblastn', tblastn_version, tblastn_status]
    if args.verbosity > 1:
        makeblastdb_row.append(makeblastdb_path)
        tblastn_row.append(tblastn_path)
    program_table.append(makeblastdb_row)
    program_table.append(tblastn_row)

    row_colours = {}
    for i, row in enumerate(program_table):
        if 'not used' in row:
            row_colours[i] = 'dim'
        elif ('too old' in row or 'too new' in row or 'not found' in row or 'bad' in row or
              'Python problem' in row):
            row_colours[i] = 'red'

    print_table(program_table, alignments='LLLL', row_colour=row_colours, max_col_width=60,
                sub_colour={'good': 'green'})

    quit_if_dependency_problem(spades_status, racon_status, makeblastdb_status, tblastn_status,
                               args)


def quit_if_dependency_problem(spades_status, racon_status, makeblastdb_status, tblastn_status,
                               args):
    if all(x == 'good' or x == 'not used'
           for x in [spades_status, racon_status, makeblastdb_status, tblastn_status]):
        return

    log.log('')
    if spades_status == 'not found':
        quit_with_error('could not find SPAdes at ' + args.spades_path)
    if spades_status == 'too old' or spades_status == 'too new':
        quit_with_error('Unicycler requires SPAdes v3.13.1 or later')
    if spades_status == 'Python problem':
        quit_with_error('SPAdes cannot run due to an incompatible Python version')
    if spades_status == 'bad':
        quit_with_error('SPAdes was found but does not produce output (make sure to use '
                        '"spades.py" location, not "spades")')
    if makeblastdb_status == 'not found':
        quit_with_error('could not find makeblastdb - either specify its location using '
                        '--makeblastdb_path or use --no_rotate to remove BLAST dependency')
    if tblastn_status == 'not found':
        quit_with_error('could not find tblastn - either specify its location using '
                        '--tblastn_path or use --no_rotate to remove BLAST dependency')
    if racon_status == 'not found':
        quit_with_error('could not find racon - either specify its location using --racon_path '
                        'or use --no_miniasm to remove Racon dependency')
    if racon_status == 'bad':
        quit_with_error('Racon was found but does not produce output')

    # Code should never get here!
    quit_with_error('Unspecified error with Unicycler dependencies')


def rotate_completed_replicons(graph, args, counter):
    completed_replicons = graph.completed_circular_replicons()
    if len(completed_replicons) > 0:
        log.log_section_header('Rotating completed replicons')
        log.log_explanation('Any completed circular contigs (i.e. single contigs which have one '
                            'link connecting end to start) can have their start position changed '
                            'without altering the sequence. For consistency, Unicycler now '
                            'searches for a starting gene (dnaA or repA) in each such contig, and '
                            'if one is found, the contig is rotated to start with that gene on '
                            'the forward strand.')

        rotation_result_table = [['Segment', 'Length', 'Depth', 'Starting gene', 'Position',
                                  'Strand', 'Identity', 'Coverage']]
        blast_dir = os.path.join(args.out, 'blast')
        if not os.path.exists(blast_dir):
            os.makedirs(blast_dir)
        completed_replicons = sorted(completed_replicons, reverse=True,
                                     key=lambda x: graph.segments[x].get_length())
        rotation_count = 0
        for completed_replicon in completed_replicons:
            segment = graph.segments[completed_replicon]
            sequence = segment.forward_sequence

            try:
                seg_name = str(segment.number)
            except AttributeError:
                seg_name = segment.full_name

            log.log('Segment ' + seg_name + ':', 2)
            rotation_result_row = [seg_name, int_to_str(len(sequence)),
                                   float_to_str(segment.depth, 2) + 'x']
            try:
                blast_hit = find_start_gene(sequence, args.start_genes, args.start_gene_id,
                                            args.start_gene_cov, blast_dir, args.makeblastdb_path,
                                            args.tblastn_path)
            except CannotFindStart:
                rotation_result_row += ['none found', '', '', '', '']
            else:
                rotation_result_row += [blast_hit.qseqid, int_to_str(blast_hit.start_pos),
                                        'reverse' if blast_hit.flip else 'forward',
                                        '%.1f' % blast_hit.pident + '%',
                                        '%.1f' % blast_hit.query_cov + '%']
                segment.rotate_sequence(blast_hit.start_pos, blast_hit.flip)
                rotation_count += 1
            rotation_result_table.append(rotation_result_row)

        log.log('', 2)
        print_table(rotation_result_table, alignments='RRRLRLRR', indent=0,
                    sub_colour={'none found': 'red'})
        if rotation_count and args.keep > 0:

            graph.save_to_gfa(gfa_path(args.out, next(counter), 'rotated'), newline=True)