Class: Bio::FinishM::RoundUp
- Inherits:
-
Object
- Object
- Bio::FinishM::RoundUp
- Includes:
- Logging
- Defined in:
- lib/finishm/roundup.rb
Constant Summary collapse
- DEFAULT_OPTIONS =
{ :contig_end_length => 200, :graph_search_leash_length => 20000, :unscaffold_first => false, :recoherence_kmer => 1, :debug => false, :gapfill_only => false, :max_explore_nodes => 10000, :max_gapfill_paths => 10, :gapfill_with_max_coverage => false, }
Instance Method Summary collapse
- #add_options(optparse_object, options) ⇒ Object
- #gapfill_a_scaffold(gapfiller, printer, master_graph, genome, scaffold_index, scaffold_direction, superscaffold_name, report, variants_file, options) ⇒ Object
- #piece_together_gapfill(printer, master_graph, first_sequence, aconn, second_sequence, gap_length, max_gapfill_paths, options = {}) ⇒ Object
- #revcom(seq) ⇒ Object
- #run(options, argv = []) ⇒ Object
- #setup_output_directory(given_directory) ⇒ Object
- #validate_options(options, argv) ⇒ Object
- #wander_a_genome(wanderer, genome, master_probed_graph, options, report) ⇒ Object
Methods included from Logging
Instance Method Details
#add_options(optparse_object, options) ⇒ Object
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 |
# File 'lib/finishm/roundup.rb', line 16 def (optparse_object, ) optparse_object. = "\nUsage: finishm roundup --genomes <genome1.fasta>[,<genome2.fasta>,...] --fastq-gz <reads..> --output-directory <directory> Takes one or more genomes and tries to improve their quality by reducing the number of scaffolds and N characters they contain. Example: finishm roundup --genomes genome1.fasta,genome2.fasta --fastq-gz reads.1.fq.gz,reads.2.fq.gz --output-directory finishm_roundup_results That will create a collapsed de-Bruijn graph from reads.1.fq.gz and reads.2.fq.gz, then try to find connections between the starts and the ends of the contigs in genome1.fasta. If any connections between contigs are mutually exclusive, then they are incorporated into scaffolds together, and gapfilling is attempted. The final sequences are output in the finishm_roundup_results directory in FASTA format. The procedure is then repeated for genome2.fasta. \n\n" .merge!(DEFAULT_OPTIONS) optparse_object.separator "\nRequired arguments:\n\n" optparse_object.on("--genomes FASTA_1[,FASTA_2...]", Array, "fasta files of genomes to be improved [required]") do |arg| [:assembly_files] = arg end optparse_object.on("--output-directory PATH", "Output results to this directory [required]") do |arg| [:output_directory] = arg end optparse_object.separator "\nThere must be some definition of reads too:\n\n" #TODO improve this help Bio::FinishM::ReadInput.new.(optparse_object, ) optparse_object.separator "\nOptional arguments:\n\n" optparse_object.on("--overhang NUM", Integer, "Start assembling this far from the ends of the contigs [default: #{[:contig_end_length] }]") do |arg| [:contig_end_length] = arg.to_i end optparse_object.on("--recoherence-kmer NUM", Integer, "Use a kmer longer than the original velvet one, to help remove bubbles and circular paths [default: none]") do |arg| [:recoherence_kmer] = arg end optparse_object.on("--leash-length NUM", Integer, "Don't explore too far in the graph, only this far and not much more [default: #{[:graph_search_leash_length] }]") do |arg| [:graph_search_leash_length] = arg end optparse_object.on("--unscaffold-first", "Break the scaffolds in the contigs file apart, and then wander between the resultant contigs. This option is only relevant to the wander step; gapfilling is attempted on all gaps by default. [default: #{[:unscaffold_first] }]") do [:unscaffold_first] = true end optparse_object.on("--gapfill-only", "Don't wander, just gapfill [default: #{[:gapfill_only] }]") do [:gapfill_only] = true end optparse_object.on("--max-gapfill-paths NUM", Integer, "When this number of paths is exceeded, don't gapfill, print as Ns [default: #{[:max_gapfill_paths] }]") do |arg| [:max_gapfill_paths] = arg end optparse_object.on("--max-explore-nodes NUM", Integer, "Only explore this many nodes. If max is reached, do not make connections. [default: #{[:max_explore_nodes] }]") do |arg| [:max_explore_nodes] = arg end optparse_object.on("--gapfill-with-max-coverage", "When gapfilling, take the path with maximal coverage and do not print variants [default: #{[:gapfill_with_max_coverage] }]") do [:gapfill_with_max_coverage] = true end optparse_object.on("--debug", "Build the graph, then drop to a pry console. [default: #{[:debug] }]") do [:debug] = true end optparse_object.on("--probe NUM",Integer,"debug mode: explore from this probe number (1-based index), gapfill only, no wander. [default: explore from all probes}]") do |arg| [:interesting_probes] = [arg-1] #convert to 0-based index [:gapfill_only] = true end #optparse_object.on("--proceed-on-short-contigs", "By default, when overly short contigs are encountered, finishm croaks. This option stops the croaking [default: #{options[:proceed_on_short_contigs] }]") do # options[:proceed_on_short_contigs] = true #end Bio::FinishM::GraphGenerator.new. optparse_object, end |
#gapfill_a_scaffold(gapfiller, printer, master_graph, genome, scaffold_index, scaffold_direction, superscaffold_name, report, variants_file, options) ⇒ Object
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 |
# File 'lib/finishm/roundup.rb', line 295 def gapfill_a_scaffold(gapfiller, printer, master_graph, genome, scaffold_index, scaffold_direction, superscaffold_name, report, variants_file, ) connections = [] genome.each_gap_probe_pair(scaffold_index) do |probe1, probe2| log.info "Gapfilling between probes #{probe1.number+1} and #{probe2.number+1}.." next unless [:interesting_probes].nil? or [:interesting_probes].include?(probe1.number) or [:interesting_probes].include?(probe2.number) connections.push gapfiller.gapfill(master_graph, probe1.index, probe2.index, ) end log.debug "Found #{connections.length} connections gapfilling in scaffold #{scaffold_index}" if log.debug? all_variants = [] num_gapfills = 0 scaffold = genome.scaffolds[scaffold_index] gapfilled_sequence = genome.scaffolds[scaffold_index].contigs[0].sequence connections.each_with_index do |aconn, i| rhs_sequence = scaffold.contigs[i+1].sequence gapfilled_sequence, variants, gapfilled = piece_together_gapfill( printer, master_graph, gapfilled_sequence, aconn, rhs_sequence, genome.gap_length(scaffold_index, i), [:max_gapfill_paths] ) if gapfilled num_gapfills += 1 variants.each{|v| all_variants << v} to_log = "Filled a gap on genome #{genome.filename}: scaffold #{scaffold.name}: #{scaffold.contigs[i].scaffold_position_end+1}-#{scaffold.contigs[i+1].scaffold_position_start-1}" report.puts to_log log.info to_log end end if scaffold_direction == false gapfilled_sequence = revcom(gapfilled_sequence) all_variants.each do |variant| variant.position = gapfilled_sequence.length - variant.position variant.reverse! end end all_variants.each do |variant| variant.reference_name = superscaffold_name variants_file.puts variant.vcf(gapfilled_sequence) end return gapfilled_sequence, num_gapfills, all_variants end |
#piece_together_gapfill(printer, master_graph, first_sequence, aconn, second_sequence, gap_length, max_gapfill_paths, options = {}) ⇒ Object
338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 |
# File 'lib/finishm/roundup.rb', line 338 def piece_together_gapfill(printer, master_graph, first_sequence, aconn, second_sequence, gap_length, max_gapfill_paths, ={}) scaffold_sequence = nil gapfilled = -1 if aconn.paths.length != 0 and aconn.paths.length <= max_gapfill_paths aconn.collapse_paths_to_maximal_coverage_path! if [:gapfill_with_max_coverage] scaffold_sequence, variants = printer.ready_two_contigs_and_connections( master_graph.graph, first_sequence, aconn, second_sequence, master_graph.velvet_sequences ) if !scaffold_sequence.nil? #sometimes it is just impossible even if there is paths scaffold_sequence.gsub!('-','') #remove gaps i.e. where the consensus is a gap gapfilled = true end end if gapfilled != true # No paths found, or gapfilling failed at the final step Just fill with Ns like it was before scaffold_sequence = first_sequence + 'N'*gap_length + second_sequence gapfilled = false end return scaffold_sequence, variants, gapfilled end |
#revcom(seq) ⇒ Object
360 361 362 |
# File 'lib/finishm/roundup.rb', line 360 def revcom(seq) Bio::Sequence::NA.new(seq).reverse_complement.to_s.upcase end |
#run(options, argv = []) ⇒ Object
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 |
# File 'lib/finishm/roundup.rb', line 105 def run(, argv=[]) # Make sure output directory is writeable to avoid late croaking output_directory = setup_output_directory [:output_directory] # Gather the probes from each genome supplied genomes = Bio::FinishM::InputGenome.parse_genome_fasta_files( [:assembly_files], [:contig_end_length], ) # Generate one velvet assembly to rule them all (forging the assembly is hard work..) probe_sequences = genomes.collect{|genome| genome.probe_sequences}.flatten # Generate the graph with the probe sequences in it. read_input = Bio::FinishM::ReadInput.new read_input. master_graph = Bio::FinishM::GraphGenerator.new.generate_graph(probe_sequences, read_input, ) binding.pry if [:debug] # For each genome, wander, gapfill, then output wanderer = Bio::FinishM::Wanderer.new gapfiller = Bio::FinishM::GapFiller.new printer = Bio::AssemblyGraphAlgorithms::ContigPrinter.new genomes.each do |genome| # wander using just the probes on the ends of the scaffolds connected_scaffolds = nil all_connections = [] gaps_filled_in_genome = 0 wandered_probe_indices = nil File.open(File.join(output_directory, File.basename(genome.filename)+".report.txt"),'w') do |report| report.puts "#{Time.now} FinishM report for roundup run with: #{.inspect}" if [:gapfill_only] log.info "Skipping wander, gapfilling only" connected_scaffolds = Bio::FinishM::ConnectionInterpreter.new([], (0...genome.scaffolds.length)).scaffolds([]) else log.debug "Wandering.." connected_scaffolds, all_connections, wandered_probe_indices = wander_a_genome(wanderer, genome, master_graph, , report) end # Write out all the connections File.open(File.join(output_directory, File.basename(genome.filename)+".connections.csv"),'w') do |con_file| all_connections.each do |connection| con_file.puts connection end end output_path = File.join(output_directory, File.basename(genome.filename)+".scaffolds.fasta") variants_path = File.join(output_directory, File.basename(genome.filename)+".at_least_half_completely_wrong.vcf") num_circular_scaffolds = 0 File.open(output_path, 'w') do |output_file| File.open(variants_path,'w') do |variants_file| variants_file.puts %w(#CHROM POS ID REF ALT QUAL FILTER INFO).join("\t") # gapfill between # (1) interpreted_connections # (2) gaps that were present before above wander connected_scaffolds.each_with_index do |cross_scaffold_connection, connected_scaffold_index| superscaffold_name = "scaffold#{connected_scaffold_index+1}" pretend_contig = cross_scaffold_connection.contigs[0] first_scaffold_index = pretend_contig.sequence_index # Gapfill contigs within the scaffold on the extreme LHS scaffold_sequence, num_gaps, variants = gapfill_a_scaffold(gapfiller, printer, master_graph, genome, first_scaffold_index, pretend_contig.direction, superscaffold_name, report, variants_file, ) gaps_filled_in_genome += num_gaps last_contig = nil cross_scaffold_connection.contigs.each_with_index do |contig, superscaffold_index| unless last_contig.nil? #skip the first contig - it be done last_name = genome.scaffolds[last_contig.sequence_index].name current_name = genome.scaffolds[contig.sequence_index].name log.debug "Connecting #{last_name} and #{current_name}" if log.debug? # Ready the contig on the RHS of this join # Gapfill within the scaffold on the RHS of the new gap rhs_sequence, num_gaps, variants = gapfill_a_scaffold(gapfiller, printer, master_graph, genome, contig.sequence_index, contig.direction, superscaffold_name, report, variants_file, ) gaps_filled_in_genome += num_gaps # Gapfill across the new gap between scaffolds aconn = gapfiller.gapfill(master_graph, last_contig.direction == true ? genome.last_probe(last_contig.sequence_index).index : genome.first_probe(last_contig.sequence_index).index, contig.direction == true ? genome.first_probe(contig.sequence_index).index : genome.last_probe(contig.sequence_index).index, ) second_sequence = genome.scaffolds[contig.sequence_index].contigs[0].sequence log.debug "Found #{aconn.paths.length} connections between #{last_name} and #{current_name}" if log.debug? connected = false if aconn.paths.length > 0 aconn.collapse_paths_to_maximal_coverage_path! if [:gapfill_with_max_coverage] scaffold_sequence2, variants = printer.ready_two_contigs_and_connections( master_graph.graph, scaffold_sequence, aconn, rhs_sequence, master_graph.velvet_sequences ) if !scaffold_sequence2.nil? scaffold_sequence = scaffold_sequence2 connected = true # Print variants # TODO: need to change coordinates of variants, particularly when >2 contigs are joined? variants.each do |variant| variant.reference_name = superscaffold_name variants_file.puts variant.vcf(scaffold_sequence) end end end if !connected # when this occurs, it is due to there being a circuit in the path, so no paths are printed. # (at least for now) TODO: this could be improved. # Just arbitrarily put in 100 N characters, to denote a join, but no gapfill # (or, it could be impossible to join because of low coverage resulting in inability to get sequence from the node trail) scaffold_sequence = scaffold_sequence+('N'*100)+rhs_sequence end end last_contig = contig end #Output the scaffold to the output directory descriptor = nil if cross_scaffold_connection.circular? descriptor = 'circular' num_circular_scaffolds += 1 else descriptor = 'scaffold' end scaffold_names = cross_scaffold_connection.contigs.collect do |contig| genome.scaffolds[contig.sequence_index].name end output_file.puts ">#{superscaffold_name} #{descriptor} #{scaffold_names.join(':') }" scaffold_sequence.gsub! '-', '' #remove dashes since these make things fail downstream output_file.puts scaffold_sequence end num_connected_scaffolds = genome.scaffolds.length - connected_scaffolds.length log.info "Wrote #{connected_scaffolds.length} scaffolds to #{output_path}, after scaffolding #{num_connected_scaffolds} scaffolds together (#{num_circular_scaffolds} were circular) and filling #{gaps_filled_in_genome} gaps." end end end end end |
#setup_output_directory(given_directory) ⇒ Object
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 |
# File 'lib/finishm/roundup.rb', line 251 def setup_output_directory(given_directory) output_directory = File.absolute_path(given_directory) log.debug "Using output directory: #{output_directory}" if log.debug? if File.exist?(output_directory) if !File.directory?(output_directory) log.error "Specified --output-directory #{output_directory} exists but is a file and not a directory. Cannot continue." exit 1 elsif !File.writable?(output_directory) log.error "Specified --output-directory #{output_directory} is not writeable. Cannot continue." exit 1 else log.debug "Already existing output directory #{output_directory} seems usable" end else # Creating a new output directory Dir.mkdir(output_directory) end return output_directory end |
#validate_options(options, argv) ⇒ Object
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 |
# File 'lib/finishm/roundup.rb', line 85 def (, argv) #TODO: give a better description of the error that has occurred #TODO: require reads options if argv.length != 0 return "Dangling argument(s) found e.g. #{argv[0] }" else [ :assembly_files, :output_directory, ].each do |sym| if [sym].nil? return "No option found to specify #{sym}." end end #if return nil from here, options all were parsed successfully return Bio::FinishM::ReadInput.new.(, []) end end |
#wander_a_genome(wanderer, genome, master_probed_graph, options, report) ⇒ Object
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 |
# File 'lib/finishm/roundup.rb', line 273 def wander_a_genome(wanderer, genome, master_probed_graph, , report) # Create new finishm_graph with only probes from the ends of the scaffolds of this genome probe_indices = [] genome.each_scaffold_end_numbered_probe{|probe| probe_indices.push(probe.number)} genome_graph = master_probed_graph.subgraph(probe_indices) num_scaffolds = genome.scaffolds.length all_connections = wanderer.probed_graph_to_connections(genome_graph, ) interpreter = Bio::FinishM::ConnectionInterpreter.new(all_connections, (0...num_scaffolds)) connections = interpreter.doubly_single_contig_connections report.puts "Found #{connections.length} connections between contigs that can be used for scaffolding" unconnected_probes = interpreter.unconnected_probes report.puts "Found #{unconnected_probes.length} contig ends that did not connect to any others" unconnected_probes.each do |probe| report.puts "Did not connect to any other probes: #{probe.inspect}" end return interpreter.scaffolds(connections), all_connections, probe_indices end |