Always filter :largecontigs by length

Before, datasets initiated from assembly bypassed this filter. Now all loaded assemblies are processed to remove contigs shorter than 1 kbp

Author: lmrodriguezr

lib/miga/common/format.rb

@@ -23,8 +23,9 @@ def tabulate(header, values, tabular = false)
   end
 
   ##
-  # Cleans a FastA file in place.
-  def clean_fasta_file(file)
+  # Cleans a FastA file in place, removing all sequences shorter than
+  # +min_len+
+  def clean_fasta_file(file, min_len = 1)
     tmp_fh = nil
     tmp_path = nil
     begin
@@ -39,19 +40,24 @@ def clean_fasta_file(file)
         tmp_path = tmp_fh.path
         fh = File.open(file, 'r')
       end
-      buffer = ''.dup
+      next_seq = ['', '']
       fh.each_line do |ln|
         ln.chomp!
         if ln =~ /^>\s*(\S+)(.*)/
           id, df = $1, $2
-          tmp_fh.print buffer.wrap_width(80)
-          buffer = ''.dup
-          tmp_fh.puts ">#{id.gsub(/[^A-Za-z0-9_\|\.]/, '_')}#{df}"
+          if next_seq[1].length >= min_len
+            tmp_fh.puts next_seq[0]
+            tmp_fh.print next_seq[1].wrap_width(80)
+          end
+          next_seq = [">#{id.gsub(/[^A-Za-z0-9_\|\.]/, '_')}#{df}", '']
         else
-          buffer << ln.gsub(/[^A-Za-z\.\-]/, '')
+          next_seq[1] += ln.gsub(/[^A-Za-z\.\-]/, '')
         end
       end
-      tmp_fh.print buffer.wrap_width(80)
+      if next_seq[1].length >= min_len
+        tmp_fh.puts next_seq[0]
+        tmp_fh.print next_seq[1].wrap_width(80)
+      end
       tmp_fh.close
       fh.close
       FileUtils.mv(tmp_path, file)

lib/miga/dataset/result/add.rb

@@ -84,7 +84,7 @@ def add_result_assembly(base, opts)
     opts[:is_clean] ||= false
     r.clean! if opts[:is_clean]
     unless r.clean?
-      MiGA::MiGA.clean_fasta_file(r.file_path(:largecontigs))
+      MiGA::MiGA.clean_fasta_file(r.file_path(:largecontigs), 1000)
       r.clean!
     end
     r

lib/miga/version.rb

@@ -12,15 +12,15 @@ module MiGA
   # - String indicating release status:
   #   - rc* release candidate, not released as gem
   #   - [0-9]+ stable release, released as gem
-  VERSION = [1.3, 18, 0].freeze
+  VERSION = [1.3, 19, 0].freeze
 
   ##
   # Nickname for the current major.minor version.
   VERSION_NAME = 'mezzotint'
 
   ##
   # Date of the current gem relese.
-  VERSION_DATE = Date.new(2024, 7, 12)
+  VERSION_DATE = Date.new(2024, 7, 17)
 
   ##
   # References of MiGA

scripts/assembly.bash

@@ -58,7 +58,7 @@ if [[ -s "$DATASET/scaffold.fa" ]] ; then
 else
   ln -s "$DATASET/contig.fa" "$DATASET.AllContigs.fna"
 fi
-FastA.length.pl "$DATASET.AllContigs.fna" | awk '$2>=1000{print $1}' \
+FastA.length.pl "$DATASET.AllContigs.fna" | awk '$2 >= 1000 { print $1 }' \
   | FastA.filter.pl /dev/stdin "$DATASET.AllContigs.fna" \
   > "$DATASET.LargeContigs.fna"
 

test/result_stats_test.rb

@@ -83,23 +83,37 @@ def test_assembly
     # Prepare result
     dir = 'data/05.assembly'
     fa = file_path(dir, '.LargeContigs.fna')
-    File.open(fa, 'w') { |fh| fh.puts '>1', 'ACTAC' }
+    File.open(fa, 'w') { |fh| fh.puts('>1', 'ACTAC' * 500) }
     touch_done(dir)
     r = dataset.add_result(:assembly)
 
     # Test assertions
     assert_equal({}, r[:stats])
     r.compute_stats
     assert_equal(1, r[:stats][:contigs])
-    assert_equal([5, 'bp'], r[:stats][:total_length])
+    assert_equal([2500, 'bp'], r[:stats][:total_length])
     assert_equal([40.0, '%'], r[:stats][:g_c_content])
   end
 
+  def test_large_contigs
+    # Prepare result
+    dir = 'data/05.assembly'
+    fa = file_path(dir, '.LargeContigs.fna')
+    File.open(fa, 'w') { |fh| fh.puts('>1', 'ACTAC' * 50) }
+    touch_done(dir)
+    r = dataset.add_result(:assembly)
+
+    # Test assertions
+    assert_equal({}, r[:stats])
+    r.compute_stats
+    assert_equal(0, r[:stats][:contigs])
+  end
+
   def test_cds
     # Prepare result
     dir = 'data/06.cds'
     fa = file_path(dir, '.faa')
-    File.open(fa, 'w') { |fh| fh.puts '>1', 'M' }
+    File.open(fa, 'w') { |fh| fh.puts('>1', 'M' * 500) }
     gff = file_path(dir, '.gff3.gz')
     Zlib::GzipWriter.open(gff) do |fh|
       fh.puts '# Model Data: a=b;transl_table=11;'
@@ -111,7 +125,7 @@ def test_cds
     assert_equal({}, r[:stats])
     r.compute_stats
     assert_equal(1, r[:stats][:predicted_proteins])
-    assert_equal([1.0, 'aa'], r[:stats][:average_length])
+    assert_equal([500.0, 'aa'], r[:stats][:average_length])
     assert_nil(r[:stats][:coding_density])
     test_assembly
     r.compute_stats