parameters: (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)
return_statement: return assembly_graph
docstring: This function tries a SPAdes assembly at different k-mers and returns the best one.
docstring_summary: This function tries a SPAdes assembly at different k-mers and returns the best one.
docstring_tokens: This function tries a SPAdes assembly at different k - mers and returns the best one .
function: 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
function_tokens: 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 ) : 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