class Segment(object):
    """
    This hold a graph segment with a number, depth, direction and sequence.
    """
    def __init__(self, number, depth, sequence, positive, bridge=None, graph_path=None,
                 original_depth=True):
        self.number = number
        self.depth = depth
        self.original_depth = original_depth
        self.forward_sequence = ''
        self.reverse_sequence = ''
        self.bridge = bridge
        self.graph_path = graph_path
        if positive:
            self.forward_sequence = sequence
        else:
            self.reverse_sequence = sequence
        self.used_in_bridges = []

    def __repr__(self):
        if len(self.forward_sequence) > 6:
            seq_string = self.forward_sequence[:3] + '...' + self.forward_sequence[-3:]
        else:
            seq_string = self.forward_sequence
        return str(self.number) + ' (' + seq_string + ')'

    def add_sequence(self, sequence, positive):
        if positive:
            self.forward_sequence = sequence
        else:
            self.reverse_sequence = sequence

    def build_other_sequence_if_necessary(self):
        if not self.forward_sequence:
            self.forward_sequence = reverse_complement(self.reverse_sequence)
        if not self.reverse_sequence:
            self.reverse_sequence = reverse_complement(self.forward_sequence)

    def get_length(self):
        return len(self.forward_sequence)

    def get_length_no_overlap(self, overlap):
        return len(self.forward_sequence) - overlap

    def is_homopolymer(self):
        """
        Returns True if the segment's sequence is made up of only one base.
        """
        if len(self.forward_sequence) == 0:
            return False
        first_base = self.forward_sequence[0].lower()
        for base in self.forward_sequence[1:]:
            if base.lower() != first_base:
                return False
        return True

    def gfa_segment_line(self):
        """
        Returns an entire S line for GFA output, including the newline.
        """
        s_line = 'S\t'
        s_line += str(self.number) + '\t'
        s_line += self.forward_sequence + '\t'
        s_line += 'LN:i:' + str(self.get_length()) + '\t'
        s_line += 'dp:f:' + str(self.depth) + '\n'
        return s_line

    def get_fasta_name_and_description_line(self, circular_seg_nums=None):
        """
        Returns the segment's fasta line, including the '>' and the newline.
        """
        name_and_description = ''.join(['>', str(self.number), ' length=', str(self.get_length()),
                                        ' depth=', '%.2f' % self.depth, 'x'])
        if circular_seg_nums and self.number in circular_seg_nums:
            name_and_description += ' circular=true'
        return name_and_description + '\n'

    def save_to_fasta(self, fasta_filename):
        """
        Saves the segment's sequence to FASTA file.
        """
        fasta = open(fasta_filename, 'w')
        fasta.write(self.get_fasta_name_and_description_line())
        fasta.write(add_line_breaks_to_sequence(self.forward_sequence))
        fasta.close()

    def get_seg_type_label(self):
        """
        Given a particular segment, this function returns a label string based its type.
        """
        if self.bridge is None:
            return ''
        if isinstance(self.bridge, SpadesContigBridge):
            label = 'SPAdes contig bridge'
        elif isinstance(self.bridge, LoopUnrollingBridge):
            label = 'Loop unrolling bridge'
        elif isinstance(self.bridge, LongReadBridge):
            label = 'Long read bridge'
        elif isinstance(self.bridge, SimpleLongReadBridge):
            label = 'Simple long read bridge'
        elif isinstance(self.bridge, MiniasmBridge):
            label = 'Miniasm bridge'
        else:
            raise TypeError("unknown bridge type")
        if self.graph_path:
            graph_path_str = ', '.join([str(x) for x in self.graph_path])
            graph_path_str = '\\n'.join(textwrap.wrap(graph_path_str, 40))
            label += ':\\n' + graph_path_str
        return label

    def trim_from_end(self, amount):
        """
        Removes the specified number of bases from the end of the segment sequence.
        """
        assert self.get_length() >= amount
        if amount == 0:
            return
        self.forward_sequence = self.forward_sequence[:-amount]
        self.reverse_sequence = self.reverse_sequence[amount:]

    def trim_from_start(self, amount):
        """
        Removes the specified number of bases from the end of the segment sequence.
        """
        assert self.get_length() >= amount
        if amount == 0:
            return
        self.forward_sequence = self.forward_sequence[amount:]
        self.reverse_sequence = self.reverse_sequence[:-amount]

    def append_to_forward_sequence(self, additional_seq):
        """
        Adds the given sequence to the end of the forward sequence (and updates the reverse
        sequence accordingly).
        """
        self.forward_sequence = self.forward_sequence + additional_seq
        self.reverse_sequence = reverse_complement(self.forward_sequence)

    def append_to_reverse_sequence(self, additional_seq):
        """
        Adds the given sequence to the end of the reverse sequence (and updates the forward
        sequence accordingly).
        """
        self.reverse_sequence = self.reverse_sequence + additional_seq
        self.forward_sequence = reverse_complement(self.reverse_sequence)

    def prepend_to_forward_sequence(self, additional_seq):
        """
        Adds the given sequence to the end of the forward sequence (and updates the reverse
        sequence accordingly).
        """
        self.forward_sequence = additional_seq + self.forward_sequence
        self.reverse_sequence = reverse_complement(self.forward_sequence)

    def prepend_to_reverse_sequence(self, additional_seq):
        """
        Adds the given sequence to the end of the reverse sequence (and updates the forward
        sequence accordingly).
        """
        self.reverse_sequence = additional_seq + self.reverse_sequence
        self.forward_sequence = reverse_complement(self.reverse_sequence)

    def remove_sequence(self):
        """
        Gets rid of the segment sequence entirely, turning it into a zero-length segment.
        """
        self.forward_sequence = ''
        self.reverse_sequence = ''

    def rotate_sequence(self, start_pos, flip):
        """
        Rotates the sequence so it begins at start_pos. If flip is True, it also switches the
        forward and reverse strands. This function assumes that the segment is a circular
        completed replicon with no overlap.
        """
        unrotated_seq = self.forward_sequence
        rotated_seq = unrotated_seq[start_pos:] + unrotated_seq[:start_pos]
        rev_comp_rotated_seq = reverse_complement(rotated_seq)

        if flip:
            self.forward_sequence = rev_comp_rotated_seq
            self.reverse_sequence = rotated_seq
        else:
            self.forward_sequence = rotated_seq
            self.reverse_sequence = rev_comp_rotated_seq 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 weighted_average_list(nums, weights):
    """
    A simple weighted mean of a list of numbers.
    """
    weight_sum = sum(weights)
    if weight_sum == 0.0:
        weights = [1.0] * len(nums)
    return sum(num * (weights[i] / weight_sum) for i, num in enumerate(nums)) def score_function(val, half_score_val):
    """
    For inputs of 0.0 and greater, this function returns a value between 0.0 and 1.0, approaching
    1.0 with large values. The half_score_val argument is the point at which the function returns
    0.5. If it's large the function approaches 1.0 more slowly, if it's small the function
    approaches 1.0 more quickly.
    """
    return 1.0 - (half_score_val / (half_score_val + val)) def add_line_breaks_to_sequence(sequence, line_length=0):
    """
    Wraps sequences to the defined length.  All resulting sequences end in a line break.
    """
    if not sequence:
        return '\n'
    if line_length <= 0:
        line_length = settings.BASES_PER_FASTA_LINE
    seq_with_breaks = ''
    pos = 0
    while pos < len(sequence):
        seq_with_breaks += sequence[pos:pos+line_length] + '\n'
        pos += line_length
    return seq_with_breaks 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_dim_timestamp():
    return dim('(' + get_timestamp() + ')') def get_right_arrow():
    """
    This function returns either a Unicode right arrow or '->', depending on the system encoding.
    """
    try:
        '\u2192'.encode(sys.stdout.encoding)
    except (AttributeError, UnicodeEncodeError):
        return '->'
    else:
        return '\u2192' def remove_dupes_preserve_order(lst):
    seen = set()
    return [x for x in lst if not (x in seen or seen.add(x))] class LongReadBridge(object):
    """
    This class describes a bridge created from long read alignments.
    """
    def __init__(self, graph, start, end):

        # The numbers of the two single copy segments which are being bridged.
        self.start_segment = start
        self.end_segment = end

        # The individual reads contributing to the bridge. The sequences/qualities are not for the
        # entire read, just the part in the bridge. In the case of overlapping alignments, this has
        # the overlap size instead of the sequence.
        self.reads = []

        # The consensus of all read sequences. If there is only one read, then this is the same as
        # that read's sequence. If there are multiple reads, this is hopefully a more accurate
        # sequence. If this is a case where the start and end segments overlap, then there will not
        # be a consensus sequence and consensus_length will be the mean of read_lengths (a negative
        # value).
        self.consensus_sequence = ''

        # The path through the unbridged graph, if one was found.
        self.graph_path = []
        self.all_paths = []

        # The bridge sequence, gotten from the graph path if a good path was found. Otherwise it's
        # from the consensus sequence.
        self.bridge_sequence = ''

        # The bridge depth, a weighted mean of the start and end depths.
        self.depth = get_mean_depth(graph.segments[abs(self.start_segment)],
                                    graph.segments[abs(self.end_segment)], graph)

        # A score used to determine the order of bridge application.
        self.quality = 1.0

        # When a bridge is applied, the segments in the bridge may have their depth reduced
        # accordingly. This member stores which segments have had their depth reduced and by how
        # much due to this bridge's application. It is stored so if this bridge is later deleted,
        # we can restore the depth to the segments.
        self.segments_reduced_depth = []

        self.graph = graph

    def __repr__(self):
        return 'long read bridge: ' + get_bridge_str(self) + \
               ' (quality = ' + float_to_str(self.quality, 2) + ')'

    def predicted_time_to_finalise(self):
        """
        This function very roughly predicts how long the bridge will take to finalise. It's not
        meant to be particularly accurate, but can hopefully be used to roughly order the bridges
        from slow to fast.
        """
        total_seq_length = 0
        seq_count = 0
        for read in self.reads:
            if not isinstance(read[0], int):
                total_seq_length += len(read[0])
                seq_count += 1
        if not seq_count:
            mean_seq_length = 0.0
        else:
            mean_seq_length = total_seq_length / seq_count

        if seq_count > 1:
            predicted_consensus_time = ((1.34e-9 * (total_seq_length ** 2)) +
                                        (2.76e-5 * total_seq_length))
        else:
            predicted_consensus_time = 0.0
        predicted_path_time = ((1.78e-7 * (mean_seq_length ** 2)) + (3.75e-3 * mean_seq_length))

        return predicted_consensus_time + predicted_path_time

    def finalise(self, scoring_scheme, min_alignment_length, read_lengths, estimated_genome_size,
                 expected_linear_seqs):
        """
        Determines the consensus sequence for the bridge, attempts to find it in the graph and
        assigns a quality score to the bridge. This is the big performance-intensive step of long
        read bridging!
        """
        start_seg = self.graph.segments[abs(self.start_segment)]
        end_seg = self.graph.segments[abs(self.end_segment)]

        output = [str(self.start_segment), str(self.end_segment), str(len(self.reads))]

        start_alignment_scaled_scores = [x[2].scaled_score for x in self.reads]
        end_alignment_scaled_scores = [x[3].scaled_score for x in self.reads]
        best_overall_scaled_score = min(max(start_alignment_scaled_scores),
                                        max(end_alignment_scaled_scores))
        alignment_scaled_scores = start_alignment_scaled_scores + end_alignment_scaled_scores
        mean_alignment_scaled_score = statistics.mean(alignment_scaled_scores)
        read_to_ref_ratios = [x[2].get_read_to_ref_ratio() for x in self.reads] + \
                             [x[3].get_read_to_ref_ratio() for x in self.reads]
        mean_read_to_ref_ratio = statistics.mean(read_to_ref_ratios)

        # Partition the full span reads into two groups: those with negative numbers (implying that
        # the two segments overlap) and those with actual sequences.
        reads_without_seq = []
        reads_with_seq = []
        for read in self.reads:
            if isinstance(read[0], int):
                reads_without_seq.append(read)
            else:
                reads_with_seq.append(read)

        # There shouldn't usually be both full spans with sequence and without. If there are some
        # of each, we'll throw out the minority group.
        if reads_with_seq and reads_without_seq:
            if len(reads_without_seq) > len(reads_with_seq):
                reads_with_seq = []
            else:
                reads_without_seq = []

        # For reads with sequence, we perform a MSA and get a consensus sequence.
        if reads_with_seq:

            self.consensus_sequence = get_consensus_sequence(reads_with_seq, scoring_scheme,
                                                             output)

            # We now make an expected scaled score for an alignment between the consensus and a
            # graph path. I.e. when we find a path in the graph for this consensus, this is about
            # how well it should align. This goes up with alignment scores (better alignment scores
            # imply better read sequences) and read count (more reads can give a better consensus).
            #     https://www.desmos.com/calculator
            #     y=100\cdot \left(\left(1-\frac{a}{100}\right)
            #       \left(1-\frac{3}{2+x}\right)+\frac{a}{100}\right)
            #     y = expected_scaled_score, x = num_span_reads, a = mean_alignment_scaled_score
            num_span_reads = len(self.reads)
            expected_scaled_score = 100.0 * ((1.0 - mean_alignment_scaled_score / 100.0) *
                                             (1.0 - (3.0 / (2.0 + num_span_reads))) +
                                             mean_alignment_scaled_score / 100.0)
            expected_scaled_score = max(expected_scaled_score, best_overall_scaled_score)

            # We can also predict the ratio in length between the consensus sequence and the
            # graph path. For a low number of input reads, this is close to the mean ratio for
            # the read alignments, but it approaches 1 as we have more reads and expect our
            # consensus to be better.
            #     https://www.desmos.com/calculator
            #     y=\left(a-1\right)\left(\frac{b}{x+b-1}\right)+1
            expected_consensus_to_ref_ratio = 1.0 + (mean_read_to_ref_ratio - 1.0) * \
                                                    (4 / (4 + num_span_reads - 1))
            target_path_length = int(round((len(self.consensus_sequence) /
                                            expected_consensus_to_ref_ratio)))

        # For reads without sequence, we simply need a mean distance.
        else:
            self.consensus_sequence = ''
            mean_overlap = int(round(sum(x[0] for x in reads_without_seq) / len(reads_without_seq)))
            output.append(str(mean_overlap))
            output.append('')
            target_path_length = 0
            expected_scaled_score = 100.0

        output.append(str(target_path_length))

        path_start_time = time.time()
        self.all_paths, progressive_path_search = \
            get_best_paths_for_seq(self.graph, self.start_segment, self.end_segment,
                                   target_path_length, self.consensus_sequence, scoring_scheme,
                                   expected_scaled_score)
        path_time = time.time() - path_start_time

        output.append(str(len(self.all_paths)))
        output.append('progressive' if progressive_path_search else 'exhaustive')
        output.append(float_to_str(path_time, 1))

        # If paths were found, use a path sequence for the bridge.
        if self.all_paths:
            if self.all_paths[0][0]:
                output.append(', '.join(str(x) for x in self.all_paths[0][0]))
            else:
                output.append('direct connection')
            best_path_len = self.graph.get_bridge_path_length(self.all_paths[0][0])
            output.append(str(best_path_len))
            output.append(float_to_str(self.all_paths[0][1], 1))
            output.append(float_to_str(self.all_paths[0][3], 2))
            output.append(str(self.all_paths[0][2]))

            self.graph_path = self.all_paths[0][0]
            self.bridge_sequence = self.graph.get_path_sequence(self.graph_path)

            # We start this bridge's quality using a function that takes into account the
            # actual, expected and minimum acceptable scores. If the actual scaled score is 100,
            # this function gives 1. If it is the expected value, this function gives a number
            # around 0.7. Then as the actual score approaches the minimum acceptable score, this
            # function approaches 0.
            #     https://www.desmos.com/calculator
            #     y=\left(\frac{1}{1+2^{a-x}}\right)^{0.5}
            #     y = self.quality, x = actual_scaled_score, a = expected_scaled_score,
            #     b = min_acceptable_scaled_score
            actual_scaled_score = self.all_paths[0][3]
            self.quality = math.sqrt(1.0 /
                                     (1.0 + 2.0 ** (expected_scaled_score - actual_scaled_score)))

        # If a path wasn't found, the consensus sequence is the bridge.
        else:
            self.graph_path = []
            output += ['', '', '', '', '']

            if self.consensus_sequence:
                self.bridge_sequence = self.consensus_sequence
            else:
                self.bridge_sequence = ''

            # The quality of non-graph-path-supported bridges depends on the number of dead ends
            # and whether or not linear sequences are expected.
            dead_end_count = 0
            if self.graph.ends_with_dead_end(self.start_segment):
                dead_end_count += 1
            if self.graph.starts_with_dead_end(self.end_segment):
                dead_end_count += 1

            if expected_linear_seqs:
                if dead_end_count == 2:
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_TWO_DEAD_ENDS_WITH_LINEAR_SEQS
                elif dead_end_count == 1:
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_ONE_DEAD_END_WITH_LINEAR_SEQS
                else:  # dead_end_count == 0
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_NO_DEAD_ENDS_WITH_LINEAR_SEQS
            else:
                if dead_end_count == 2:
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_TWO_DEAD_ENDS
                elif dead_end_count == 1:
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_ONE_DEAD_END
                else:  # dead_end_count == 0
                    self.quality = settings.PATHLESS_BRIDGE_QUAL_NO_DEAD_ENDS

            # Bridge length affects quality too: short bridges are better.
            bridge_len = max(0, len(self.bridge_sequence))
            half_qual_len = settings.LONG_READ_BRIDGE_HALF_QUAL_LENGTH
            self.quality *= half_qual_len / (bridge_len + half_qual_len)

        # Expected read count is determined using the read lengths and bridge size. For a given
        # read length and bridge, there are an estimable number of positions where a read of that
        # length would be able to contribute to the bridge. This is used to get the probability
        # that any read would create a bridge, and totalling those up gives us our estimated count.
        min_read_len = (2 * min_alignment_length) + len(self.bridge_sequence)
        total_possible_placements = 0
        for read_len, count in read_lengths.items():
            if read_len < min_read_len:
                continue
            possible_read_placements = read_len - min_read_len + 1
            possible_read_placements *= count
            possible_read_placements *= max(self.depth, 1)
            total_possible_placements += possible_read_placements
        expected_read_count = total_possible_placements / estimated_genome_size
        actual_read_count = len(self.reads)

        # Adjust the expected read count down, especially for higher values.
        # TO DO: reevaluate this step - is it necessary?
        expected_read_count = reduce_expected_count(expected_read_count, 30, 0.5)

        # The start segment and end segment should agree in depth. If they don't, that's very bad,
        # as it implies that they aren't actually single copy or on the same piece of DNA.
        depth_agreement_factor = get_depth_agreement_factor(start_seg.depth, end_seg.depth)
        self.quality *= depth_agreement_factor

        # The number of reads which contribute to a bridge is a big deal, so the read count factor
        # scales linearly. This value is capped at 1, which means that bridges with too few reads
        # are punished but bridges with excess reads are not rewarded.
        try:
            read_count_factor = min(1.0, actual_read_count / expected_read_count)
            self.quality *= read_count_factor
        except ZeroDivisionError:
            pass

        # The length of alignments to the start/end segments is positively correlated with quality
        # to reward bridges with long alignments. Specifically, we want there to be at least one
        # spanning read with a long alignment to the start segment and at least one spanning read
        # with a long alignment to the end segment.
        longest_start_alignment = max(x[2].get_aligned_ref_length() for x in self.reads)
        longest_end_alignment = max(x[3].get_aligned_ref_length() for x in self.reads)
        alignment_length = min(longest_start_alignment, longest_end_alignment)
        align_length_factor = score_function(alignment_length, min_alignment_length * 4)
        self.quality *= align_length_factor

        # The mean alignment score to the start/end segments is positively correlated with quality,
        # so bridges with high quality alignments are rewarded. Specifically, we want there to be at
        # least one spanning read with a high quality alignment to the start segment and at least
        # one spanning read with a high quality alignment to the end segment.
        best_start_alignment = max(x[2].scaled_score for x in self.reads)
        best_end_alignment = max(x[3].scaled_score for x in self.reads)
        alignment_quality = min(best_start_alignment, best_end_alignment)
        align_score_factor = alignment_quality / 100.0
        self.quality *= align_score_factor

        # Bridges between long start/end segments are rewarded, as they are more likely to actually
        # be single copy. We apply a length factor for both the start and the end segments,
        # and then apply the smaller of two again. This is to punish cases where both segments
        # are not long.
        start_length_factor = score_function(start_seg.get_length(), min_alignment_length * 4)
        self.quality *= start_length_factor
        end_length_factor = score_function(end_seg.get_length(), min_alignment_length * 4)
        self.quality *= end_length_factor
        smaller_length_factor = min(start_length_factor, end_length_factor)
        self.quality *= smaller_length_factor

        # We finalise the quality to a range of 0 to 100. We also use the sqrt function to pull
        # the scores up a bit (otherwise they tend to hang near the bottom of the range).
        self.quality = 100.0 * math.sqrt(self.quality)

        # noinspection PyTypeChecker
        output.append(self.quality)

        return output

    def set_path_based_on_availability(self, graph, unbridged_graph):
        """
        This function will change a bridge's graph path based on what's currently available. This
        is to handle the case where a bridge has multiple possible graph paths, but its first
        choice isn't available anymore (because other bridges used up those segments).
        It has to balance the path quality with the path availability to make a choice.
        """
        best_path = self.all_paths[0][0]
        best_sequence = unbridged_graph.get_path_sequence(best_path)
        best_scaled_score = self.all_paths[0][3]
        best_availability = graph.get_path_availability(best_path)
        for i in range(1, len(self.all_paths)):
            potential_path = self.all_paths[i][0]
            potential_scaled_score = self.all_paths[i][3]
            potential_availability = graph.get_path_availability(potential_path)

            # relative_score measures how much worse this path aligned than the current best.
            # Differences matter more close to scores of 100. E.g. 99 to 97 is a big drop in
            # score, but 79 to 77 is not as large of a drop.
            if potential_scaled_score == 100.0:
                relative_score = 1.0
            else:
                relative_score = (100.0 - best_scaled_score) / (100.0 - potential_scaled_score)
                relative_score = min(1.0, relative_score)

            # relative_availability measures how much more available this path is than the current
            # best. We use 1.1 (instead of 1.0) in the equation to attenuate the affect of
            # availability a bit (because score is more important).
            relative_availability = (1.1 - best_availability) / (1.1 - potential_availability)
            relative_availability = min(2.0, relative_availability)

            # If this path looks better than our current best (considering both score and
            # availability), then it becomes the new best.
            if relative_score * relative_availability > 1.0:
                best_path = potential_path
                best_sequence = unbridged_graph.get_path_sequence(potential_path)
                best_scaled_score = potential_scaled_score
                best_availability = potential_availability

        self.graph_path = best_path
        self.bridge_sequence = best_sequence

    @staticmethod
    def get_type_score():
        """
        Returns a score indicating the relative importance of the bridge types:
        LongReadBridge = 2, SpadesContigBridge = 1, LoopUnrollingBridge = 0
        """
        return 2

    @staticmethod
    def get_type_name():
        """
        Returns the name of the bridge type.
        """
        return 'long read' class MiniasmBridge(object):
    """
    This class describes a bridge created from long read alignments.
    """
    def __init__(self, graph, start, end, bridge_sequence, start_overlap, end_overlap,
                 scoring_scheme, output, do_path_search=True):

        # The numbers of the two single copy segments which are being bridged.
        self.start_segment = start
        self.end_segment = end
        output += [str(start), str(end), '']

        # In some miniasm bridges where the contigs are very close, there will be overlap between
        # the contigs and the bridge. I.e. when the bridge is applied, contig sequence will be
        # trimmed off.
        self.start_overlap = start_overlap
        self.end_overlap = end_overlap

        # The bridge depth, a weighted mean of the start and end depths.
        self.depth = get_mean_depth(graph.segments[abs(self.start_segment)],
                                    graph.segments[abs(self.end_segment)], graph)

        self.segments_reduced_depth = []

        # When a MiniasmBridge is being made by splitting an already-made bridge, then we don't
        # do the path search.
        if not do_path_search:
            self.bridge_sequence = bridge_sequence
            self.all_paths = []
            self.graph_path = []
            self.quality = 1.0

        else:
            # Look for a graph path corresponding to the bridge sequence.
            target_path_length = len(bridge_sequence)
            output += [str(target_path_length), '', str(target_path_length)]
            path_start_time = time.time()
            self.all_paths, progressive_path_search = \
                get_best_paths_for_seq(graph, self.start_segment, self.end_segment,
                                       target_path_length, bridge_sequence, scoring_scheme, 90.0)
            path_time = time.time() - path_start_time

            output.append(str(len(self.all_paths)))
            output.append('progressive' if progressive_path_search else 'exhaustive')
            output.append(float_to_str(path_time, 1))

            if self.all_paths:
                self.graph_path = self.all_paths[0][0]
                raw_score = self.all_paths[0][1]
                scaled_score = self.all_paths[0][3]
                len_discrepancy = self.all_paths[0][2]
                if self.graph_path:
                    output.append(', '.join(str(x) for x in self.graph_path))
                else:
                    output.append('direct connection')
                best_path_len = graph.get_bridge_path_length(self.graph_path)
                output.append(str(best_path_len))
                output.append(float_to_str(raw_score, 1))
                output.append(float_to_str(scaled_score, 2))
                output.append(str(len_discrepancy))

            else:
                self.graph_path = []
                output += ['', '', '', '', '']
                scaled_score = 0.0

            # If a very good match was found, use the graph path and give the bridge a very high
            # quality score.
            if scaled_score > settings.MINIASM_BRIDGE_SCALED_SCORE_TO_USE_GRAPH_PATH:
                self.bridge_sequence = graph.get_path_sequence(self.graph_path)
                self.quality = settings.MINIASM_BRIDGE_QUAL_WITH_GRAPH_PATH

            # Otherwise, just use the miniasm bridge sequence for the bridge.
            else:
                self.bridge_sequence = bridge_sequence
                if graph.ends_with_dead_end(self.start_segment) or \
                        graph.starts_with_dead_end(self.end_segment):
                    self.quality = settings.MINIASM_BRIDGE_QUAL_WITH_DEAD_END
                else:
                    self.quality = settings.MINIASM_BRIDGE_QUAL_WITHOUT_PATH_OR_DEAD_END

            # Depth agreement affects bridge quality.
            start_seg = graph.segments[abs(self.start_segment)]
            end_seg = graph.segments[abs(self.end_segment)]
            self.quality *= get_depth_agreement_factor(start_seg.depth, end_seg.depth)

            # Bridge length affects quality too: short bridges are better.
            bridge_len = max(0, len(self.bridge_sequence))
            half_qual_len = settings.MINIASM_BRIDGE_HALF_QUAL_LENGTH
            self.quality *= half_qual_len / (bridge_len + half_qual_len)

            self.quality = 100.0 * math.sqrt(self.quality)
            output.append(self.quality)

    def __repr__(self):
        return 'miniasm bridge: ' + get_bridge_str(self) + \
               ' (quality = ' + float_to_str(self.quality, 2) + ')'

    @staticmethod
    def get_type_score():
        """
        Returns a score indicating the relative importance of the bridge types:
        LongReadBridge = 2, SpadesContigBridge = 1, LoopUnrollingBridge = 0
        """
        return 2

    @staticmethod
    def get_type_name():
        """
        Returns the name of the bridge type.
        """
        return 'miniasm'
import math
import copy
import os
import itertools
from collections import deque, defaultdict
from .assembly_graph_segment import Segment
from .misc import int_to_str, float_to_str, weighted_average_list, score_function, \
    add_line_breaks_to_sequence, print_table, get_dim_timestamp, get_right_arrow, \
    remove_dupes_preserve_order
from .bridge_long_read import LongReadBridge
from .bridge_miniasm import MiniasmBridge
from . import settings
from . import log
        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))