parameters: (graph, read_dict, read_names, anchor_segments, verbosity,
                               min_scaled_score, threads, scoring_scheme, min_alignment_length,
                               expected_linear_seqs, min_bridge_qual)
return_statement: return split_bridges
docstring: Makes bridges between single copy segments using the alignments in the long reads.
docstring_summary: Makes bridges between single copy segments using the alignments in the long reads.
docstring_tokens: Makes bridges between single copy segments using the alignments in the long reads .
function: 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
function_tokens: 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 ) : 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