]> git.donarmstrong.com Git - biopieces.git/blob - code_ruby/test/maasha/test_seq.rb
rewrite of FASTQ internals
[biopieces.git] / code_ruby / test / maasha / test_seq.rb
1 #!/usr/bin/env ruby
2 $:.unshift File.join(File.dirname(__FILE__), '..', '..')
3
4 # Copyright (C) 2011 Martin A. Hansen.
5
6 # This program is free software; you can redistribute it and/or
7 # modify it under the terms of the GNU General Public License
8 # as published by the Free Software Foundation; either version 2
9 # of the License, or (at your option) any later version.
10
11 # This program is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 # GNU General Public License for more details.
15
16 # You should have received a copy of the GNU General Public License
17 # along with this program; if not, write to the Free Software
18 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19
20 # http://www.gnu.org/copyleft/gpl.html
21
22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
23
24 # This software is part of the Biopieces framework (www.biopieces.org).
25
26 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
27
28 require 'maasha/seq'
29 require 'test/unit'
30 require 'test/helper'
31
32 class TestSeq < Test::Unit::TestCase 
33   def setup
34     @entry = Seq.new
35   end
36
37   test "Seq.new_bp returns correctly" do
38     record = {:SEQ_NAME => "test", :SEQ => "ATCG", :SEQ_TYPE => "dna", :SCORES => "hhhh"}
39     seq    = Seq.new_bp(record)
40     assert_equal("test", seq.seq_name)
41     assert_equal("ATCG", seq.seq)
42     assert_equal("dna",  seq.type)
43     assert_equal("hhhh", seq.qual)
44   end
45
46   test "#is_dna? with no sequence type returns false" do
47     assert(@entry.is_dna? == false)
48   end
49
50   test "#is_dna? with dna sequence type returns true" do
51     @entry.type = 'dna'
52     assert(@entry.is_dna? == true)
53   end
54
55   test "#is_rna? with no sequence type returns false" do
56     assert(@entry.is_rna? == false)
57   end
58
59   test "#is_rna? with rna sequence type returns true" do
60     @entry.type = 'rna'
61     assert(@entry.is_rna? == true)
62   end
63
64   test "#is_protein? with no sequence type returns false" do
65     assert(@entry.is_protein? == false)
66   end
67
68   test "#is_protein? with protein sequence type returns true" do
69     @entry.type = 'protein'
70     assert_equal(true, @entry.is_protein?)
71   end
72
73   test "#type_guess without sequence raises" do
74     assert_raise(SeqError) { @entry.type_guess }
75   end
76
77   test "#type_guess with protein returns protein" do
78     @entry.seq = 'atcatcrFgatcg'
79     assert_equal('protein', @entry.type_guess)
80   end
81
82   test "#type_guess with rna returns rna" do
83     @entry.seq = 'atcatcrUgatcg'
84     assert_equal('rna', @entry.type_guess)
85   end
86
87   test "#type_guess with dna returns dna" do
88     @entry.seq = 'atcatcgatcg'
89     assert_equal('dna', @entry.type_guess)
90   end
91
92   test "#type_guess! without sequence raises" do
93     assert_raise(SeqError) { @entry.type_guess! }
94   end
95
96   test "#type_guess! with protein returns protein" do
97     @entry.seq = 'atcatcrFgatcg'
98     @entry.type_guess!
99     assert_equal('protein', @entry.type)
100   end
101
102   test "#type_guess! with rna returns rna" do
103     @entry.seq = 'atcatcrUgatcg'
104     @entry.type_guess!
105     assert_equal('rna', @entry.type)
106   end
107
108   test "#type_guess! with dna returns dna" do
109     @entry.seq = 'atcatcgatcg'
110     @entry.type_guess!
111     assert_equal('dna', @entry.type)
112   end
113
114   test "#length returns corretly" do
115     @entry.seq = 'ATCG'
116     assert_equal(4, @entry.length)
117   end
118
119   test "#indels returns correctly" do
120     @entry.seq = 'ATCG.-~_'
121     assert_equal(4, @entry.indels)
122   end
123
124   test "#to_rna with no sequence raises" do
125     @entry.type = 'dna'
126     assert_raise(SeqError) { @entry.to_rna }
127   end
128
129   test "#to_rna with bad type raises" do
130     @entry.seq  = 'ATCG'
131     @entry.type = 'rna'
132     assert_raise(SeqError) { @entry.to_rna }
133   end
134
135   test "#to_rna transcribes correctly" do
136     @entry.seq  = 'ATCGatcg'
137     @entry.type = 'dna'
138     assert_equal("AUCGaucg", @entry.to_rna)
139   end
140
141   test "#to_rna changes entry type to rna" do
142     @entry.seq  = 'ATCGatcg'
143     @entry.type = 'dna'
144     @entry.to_rna
145     assert_equal("rna", @entry.type)
146   end
147
148   test "#to_dna with no sequence raises" do
149     @entry.type = 'rna'
150     assert_raise(SeqError) { @entry.to_dna }
151   end
152
153   test "#to_dna with bad type raises" do
154     @entry.seq  = 'AUCG'
155     @entry.type = 'dna'
156     assert_raise(SeqError) { @entry.to_dna }
157   end
158
159   test "#to_dna transcribes correctly" do
160     @entry.seq  = 'AUCGaucg'
161     @entry.type = 'rna'
162     assert_equal("ATCGatcg", @entry.to_dna)
163   end
164
165   test "#to_dna changes entry type to dna" do
166     @entry.seq  = 'AUCGaucg'
167     @entry.type = 'rna'
168     @entry.to_dna
169     assert_equal("dna", @entry.type)
170   end
171
172   test "#to_bp returns correct record" do
173     @entry.seq_name = 'test'
174     @entry.seq      = 'ATCG'
175     assert_equal({:SEQ_NAME=>"test", :SEQ=>"ATCG", :SEQ_LEN=>4}, @entry.to_bp)
176   end
177
178   test "#to_bp with missing seq_name raises" do
179     @entry.seq = 'ATCG'
180     assert_raise(SeqError) { @entry.to_bp }
181   end
182
183   test "#to_bp with missing sequence raises" do
184     @entry.seq_name = 'test'
185     assert_raise(SeqError) { @entry.to_bp }
186   end
187
188   test "#to_fasta with missing seq_name raises" do
189     @entry.seq = 'ATCG'
190     assert_raise(SeqError) { @entry.to_fasta }
191   end
192
193   test "#to_fasta with empty seq_name raises" do
194     @entry.seq_name = ''
195     @entry.seq      = 'ATCG'
196     assert_raise(SeqError) { @entry.to_fasta }
197   end
198
199   test "#to_fasta with missing seq raises" do
200     @entry.seq_name = 'test'
201     assert_raise(SeqError) { @entry.to_fasta }
202   end
203
204   test "#to_fasta with empty seq raises" do
205     @entry.seq_name = 'test'
206     @entry.seq      = ''
207     assert_raise(SeqError) { @entry.to_fasta }
208   end
209
210   test "#to_fasta returns correct entry" do
211     @entry.seq_name = 'test'
212     @entry.seq      = 'ATCG'
213     assert_equal(">test\nATCG\n", @entry.to_fasta)
214   end
215
216   test "#to_fasta wraps correctly" do
217     entry = Seq.new("test", "ATCG")
218     assert_equal(">test\nAT\nCG\n", entry.to_fasta(2))
219   end
220
221   test "#to_fastq returns correct entry" do
222     @entry.seq_name = 'test'
223     @entry.seq      = 'ATCG'
224     @entry.qual     = 'hhhh'
225     assert_equal("@test\nATCG\n+\nhhhh\n", @entry.to_fastq)
226   end
227
228   test "#to_key with bad residue raises" do
229     entry = Seq.new("test", "AUCG")
230     assert_raise(SeqError) { entry.to_key }
231   end
232
233   test "#to_key returns correctly" do
234     entry = Seq.new("test", "ATCG")
235     assert_equal(54, entry.to_key)
236   end
237
238   test "#reverse returns correctly" do
239     @entry.seq = "ATCG"
240     new_entry  = @entry.reverse
241     assert_equal("GCTA", new_entry.seq)
242     assert_equal("ATCG", @entry.seq)
243   end
244
245   test "#reverse! returns correctly" do
246     @entry.seq = "ATCG"
247     @entry.reverse!
248     assert_equal("GCTA", @entry.seq)
249   end
250
251   test "#complement with no sequence raises" do
252     @entry.type = 'dna'
253     assert_raise(SeqError) { @entry.complement }
254   end
255
256   test "#complement with bad type raises" do
257     @entry.seq  = 'ATCG'
258     @entry.type = 'protein'
259     assert_raise(SeqError) { @entry.complement }
260   end
261
262   test "#complement for DNA is correct" do
263     @entry.seq  = 'ATCGatcg'
264     @entry.type = 'dna'
265     comp        = @entry.complement
266     assert_equal("TAGCtagc", comp.seq)
267     assert_equal("ATCGatcg", @entry.seq)
268   end
269
270   test "#complement for RNA is correct" do
271     @entry.seq  = 'AUCGaucg'
272     @entry.type = 'rna'
273     comp        = @entry.complement
274     assert_equal("UAGCuagc", comp.seq)
275     assert_equal("AUCGaucg", @entry.seq)
276   end
277
278   test "#complement! with no sequence raises" do
279     @entry.type = 'dna'
280     assert_raise(SeqError) { @entry.complement! }
281   end
282
283   test "#complement! with bad type raises" do
284     @entry.seq  = 'ATCG'
285     @entry.type = 'protein'
286     assert_raise(SeqError) { @entry.complement! }
287   end
288
289   test "#complement! for DNA is correct" do
290     @entry.seq  = 'ATCGatcg'
291     @entry.type = 'dna'
292     assert_equal("TAGCtagc", @entry.complement!.seq)
293   end
294
295   test "#complement! for RNA is correct" do
296     @entry.seq  = 'AUCGaucg'
297     @entry.type = 'rna'
298     assert_equal("UAGCuagc", @entry.complement!.seq)
299   end
300
301
302   test "#hamming distance returns correctly" do
303     seq1 = Seq.new("test1", "ATCG")
304     seq2 = Seq.new("test2", "atgg")
305     assert_equal(1, seq1.hamming_distance(seq2))
306   end
307
308   test "#generate with length < 1 raises" do
309     assert_raise(SeqError) { @entry.generate(-10, "dna") }
310     assert_raise(SeqError) { @entry.generate(0, "dna") }
311   end
312
313   test "#generate with bad type raises" do
314     assert_raise(SeqError) { @entry.generate(10, "foo") }
315   end
316
317   test "#generate with ok type dont raise" do
318     %w[dna DNA rna RNA protein Protein].each do |type|
319       assert_nothing_raised { @entry.generate(10, type) }
320     end
321   end
322
323   test "#shuffle returns correctly" do
324     orig       = "actgactgactgatcgatcgatcgatcgtactg" 
325     @entry.seq = "actgactgactgatcgatcgatcgatcgtactg"
326     entry_shuf = @entry.shuffle
327     assert_equal(orig, @entry.seq)
328     assert_not_equal(@entry.seq, entry_shuf.seq)
329   end
330
331   test "#shuffle! returns correctly" do
332     @entry.seq = "actgactgactgatcgatcgatcgatcgtactg"
333     assert_not_equal(@entry.seq, @entry.shuffle!.seq)
334   end
335
336   test "#subseq with start < 0 raises" do
337     @entry.seq = "ATCG"
338     assert_raise(SeqError) { @entry.subseq(-1, 1) }
339   end
340
341   test "#subseq with start plus length gt seq raises" do
342     @entry.seq = "ATCG"
343     assert_raise(SeqError) { @entry.subseq(0, 5) }
344   end
345
346   test "#subseq returns correct sequence" do
347     @entry.seq  = "ATCG"
348     assert_equal("AT", @entry.subseq(0, 2).seq)
349     assert_equal("CG", @entry.subseq(2, 2).seq)
350   end
351
352   test "#subseq without length returns correct sequence" do
353     @entry.seq  = "ATCG"
354     assert_equal("ATCG", @entry.subseq(0).seq)
355     assert_equal("CG",   @entry.subseq(2).seq)
356   end
357
358   test "#subseq returns correct qual" do
359     @entry.seq  = "ATCG"
360     @entry.qual = "abcd"
361     assert_equal("ab", @entry.subseq(0, 2).qual)
362     assert_equal("cd", @entry.subseq(2, 2).qual)
363   end
364
365   test "#subseq without length returns correct qual" do
366     @entry.seq  = "ATCG"
367     @entry.qual = "abcd"
368     assert_equal("abcd", @entry.subseq(0).qual)
369     assert_equal("cd",   @entry.subseq(2).qual)
370   end
371
372   test "#subseq! with start < 0 raises" do
373     @entry.seq = "ATCG"
374     assert_raise(SeqError) { @entry.subseq!(-1, 1) }
375   end
376
377   test "#subseq! with start plus length > seq.length raises" do
378     @entry.seq = "ATCG"
379     assert_raise(SeqError) { @entry.subseq!(0, 5) }
380   end
381
382   test "#subseq! returns correct sequence" do
383     @entry.seq  = "ATCG"
384     @entry.subseq!(0, 2)
385     assert_equal("AT", @entry.seq)
386     @entry.seq  = "ATCG"
387     @entry.subseq!(2, 2)
388     assert_equal("CG", @entry.seq)
389   end
390
391   test "#subseq! without length returns correct sequence" do
392     @entry.seq  = "ATCG"
393     @entry.subseq!(0)
394     assert_equal("ATCG", @entry.seq)
395     @entry.seq  = "ATCG"
396     @entry.subseq!(2)
397     assert_equal("CG", @entry.seq)
398   end
399
400   test "#subseq! with pos and length returns correct qual" do
401     @entry.seq  = "ATCG"
402     @entry.qual = "abcd"
403     @entry.subseq!(0, 2)
404     assert_equal("ab", @entry.qual)
405     @entry.seq  = "ATCG"
406     @entry.qual = "abcd"
407     @entry.subseq!(2, 2)
408     assert_equal("cd", @entry.qual)
409   end
410
411   test "#subseq! with pos returns correct qual" do
412     @entry.seq  = "ATCG"
413     @entry.qual = "abcd"
414     @entry.subseq!(0)
415     assert_equal("abcd", @entry.qual)
416     @entry.seq  = "ATCG"
417     @entry.qual = "abcd"
418     @entry.subseq!(2)
419     assert_equal("cd", @entry.qual)
420   end
421
422   test "#subseq_rand returns correct sequence" do
423     @entry.seq  = "ATCG"
424     assert_equal("ATCG", @entry.subseq_rand(4).seq)
425   end
426
427   test "#indels_remove without qual returns correctly" do
428     @entry.seq  = "A-T.CG~CG"
429     @entry.qual = nil
430     assert_equal("ATCGCG", @entry.indels_remove.seq)
431   end
432
433   test "#indels_remove with qual returns correctly" do
434     @entry.seq  = "A-T.CG~CG"
435     @entry.qual = "a@b@cd@fg"
436     assert_equal("ATCGCG", @entry.indels_remove.seq)
437     assert_equal("abcdfg", @entry.indels_remove.qual)
438   end
439
440   test "#composition returns correctly" do
441     @entry.seq = "AAAATTTCCG"
442     assert_equal(4, @entry.composition["A"])
443     assert_equal(3, @entry.composition["T"])
444     assert_equal(2, @entry.composition["C"])
445     assert_equal(1, @entry.composition["G"])
446     assert_equal(0, @entry.composition["X"])
447   end
448
449   test "#homopol_max returns 0 with empty sequence" do
450     @entry.seq = ""
451     assert_equal(0, @entry.homopol_max)
452   end
453
454   test "#homopol_max returns 0 with nil sequence" do
455     @entry.seq = nil
456     assert_equal(0, @entry.homopol_max)
457   end
458
459   test "#homopol_max returns 0 when not found" do
460     @entry.seq = "AtTcCcGggGnnNnn"
461     assert_equal(0, @entry.homopol_max(6))
462   end
463
464   test "#homopol_max returns correctly" do
465     @entry.seq = "AtTcCcGggGnnNnn"
466     assert_equal(5, @entry.homopol_max(3))
467   end
468
469   test "#hard_mask returns correctly" do
470     @entry.seq = "--AAAANn"
471     assert_equal(33.33, @entry.hard_mask)
472   end
473
474   test "#soft_mask returns correctly" do
475     @entry.seq = "--AAAa"
476     assert_equal(25.00, @entry.soft_mask)
477   end
478
479   test "#mask_seq_hard! with nil seq raises" do
480     @entry.seq  = nil
481     @entry.qual = ""
482
483     assert_raise(SeqError) { @entry.mask_seq_hard!(20) }
484   end
485
486   test "#mask_seq_hard! with nil qual raises" do
487     @entry.seq  = ""
488     @entry.qual = nil
489
490     assert_raise(SeqError) { @entry.mask_seq_hard!(20) }
491   end
492
493   test "#mask_seq_hard! with bad cutoff raises" do
494     assert_raise(SeqError) { @entry.mask_seq_hard!(-1) }
495     assert_raise(SeqError) { @entry.mask_seq_hard!(41) }
496   end
497
498   test "#mask_seq_hard! with OK cutoff dont raise" do
499     @entry.seq  = "ATCG"
500     @entry.qual = "RSTU"
501
502     assert_nothing_raised { @entry.mask_seq_hard!(0) }
503     assert_nothing_raised { @entry.mask_seq_hard!(40) }
504   end
505
506   test "#mask_seq_hard! returns correctly" do
507     @entry.seq  = "-ATCG"
508     @entry.qual = "33456"
509
510     assert_equal("-NNCG", @entry.mask_seq_hard!(20).seq)
511   end
512
513   test "#mask_seq_soft! with nil seq raises" do
514     @entry.seq  = nil
515     @entry.qual = ""
516
517     assert_raise(SeqError) { @entry.mask_seq_soft!(20) }
518   end
519
520   test "#mask_seq_soft! with nil qual raises" do
521     @entry.seq  = ""
522     @entry.qual = nil
523
524     assert_raise(SeqError) { @entry.mask_seq_soft!(20) }
525   end
526
527   test "#mask_seq_soft! with bad cutoff raises" do
528     assert_raise(SeqError) { @entry.mask_seq_soft!(-1) }
529     assert_raise(SeqError) { @entry.mask_seq_soft!(41) }
530   end
531
532   test "#mask_seq_soft! with OK cutoff dont raise" do
533     @entry.seq  = "ATCG"
534     @entry.qual = "RSTU"
535
536     assert_nothing_raised { @entry.mask_seq_soft!(0) }
537     assert_nothing_raised { @entry.mask_seq_soft!(40) }
538   end
539
540   test "#mask_seq_soft! returns correctly" do
541     @entry.seq  = "-ATCG"
542     @entry.qual = "33456"
543
544     assert_equal("-atCG", @entry.mask_seq_soft!(20).seq)
545   end
546
547   # qual score detection
548
549   test "#qual_base33? returns correctly" do
550     # self.qual.match(/[!-:]/)
551     @entry.qual = '!"#$%&\'()*+,-./0123456789:'
552     assert_equal(true,  @entry.qual_base33? )
553     @entry.qual = 32.chr
554     assert_equal(false, @entry.qual_base33? )
555     @entry.qual = 59.chr
556     assert_equal(false, @entry.qual_base33? )
557   end
558
559   test "#qual_base64? returns correctly" do
560     # self.qual.match(/[K-h]/)
561     @entry.qual = 'KLMNOPQRSTUVWXYZ[\]^_`abcdefgh'
562     assert_equal(true,  @entry.qual_base64? )
563     @entry.qual = 74.chr
564     assert_equal(false, @entry.qual_base64? )
565     @entry.qual = 105.chr
566     assert_equal(false, @entry.qual_base64? )
567   end
568
569   test "#qual_valid? with nil qual raises" do
570     assert_raise(SeqError) { @entry.qual_valid?(:base_33) }
571     assert_raise(SeqError) { @entry.qual_valid?(:base_64) }
572   end
573
574   test "#qual_valid? with bad encoding raises" do
575     @entry.qual = "abc"
576     assert_raise(SeqError) { @entry.qual_valid?("foobar") }
577   end
578
579   test "#qual_valid? with OK range returns correctly" do
580     @entry.qual = ((Seq::SCORE_MIN + 33).chr .. (Seq::SCORE_MAX + 33).chr).to_a.join
581     assert_equal(true,  @entry.qual_valid?(:base_33))
582     @entry.qual = ((Seq::SCORE_MIN + 64).chr .. (Seq::SCORE_MAX + 64).chr).to_a.join
583     assert_equal(true,  @entry.qual_valid?(:base_64))
584   end
585
586   test "#qual_valid? with bad range returns correctly" do
587     @entry.qual = ((Seq::SCORE_MIN + 33 - 1).chr .. (Seq::SCORE_MAX + 33).chr).to_a.join
588     assert_equal(false,  @entry.qual_valid?(:base_33))
589     @entry.qual = ((Seq::SCORE_MIN + 33).chr .. (Seq::SCORE_MAX + 33 + 1).chr).to_a.join
590     assert_equal(false,  @entry.qual_valid?(:base_33))
591
592     @entry.qual = ((Seq::SCORE_MIN + 64 - 1).chr .. (Seq::SCORE_MAX + 64).chr).to_a.join
593     assert_equal(false,  @entry.qual_valid?(:base_64))
594     @entry.qual = ((Seq::SCORE_MIN + 64).chr .. (Seq::SCORE_MAX + 64 + 1).chr).to_a.join
595     assert_equal(false,  @entry.qual_valid?(:base_64))
596   end
597
598   # convert sanger to ...
599
600   test "#qual_convert! from base33 to base33 returns OK" do
601     @entry.qual = 'BCDEFGHI'
602     assert_equal('BCDEFGHI', @entry.qual_convert!(:base_33, :base_33).qual)
603   end
604
605   test "#qual_convert! from base33 to base64 returns OK" do
606     @entry.qual = 'BCDEFGHI'
607     assert_equal('abcdefgh', @entry.qual_convert!(:base_33, :base_64).qual)
608   end
609
610   test "#qual_convert! from base64 to base64 returns OK" do
611     @entry.qual = 'BCDEFGHI'
612     assert_equal('BCDEFGHI', @entry.qual_convert!(:base_64, :base_64).qual)
613   end
614
615   test "#qual_convert! from base64 to base33 returns OK" do
616     @entry.qual = 'abcdefgh'
617     assert_equal('BCDEFGHI', @entry.qual_convert!(:base_64, :base_33).qual)
618   end
619
620   test "#qual_coerce! returns correctly" do
621     @entry.qual = ('!' .. '~').to_a.join
622     assert_equal("!\"\#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII", @entry.qual_coerce!(:base_33).qual)
623     @entry.qual = ('!' .. '~').to_a.join
624     assert_equal("!\"\#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZh\\h^_`abcdefghhhhhhhhhhhhhhhhhhhhhhh", @entry.qual_coerce!(:base_64).qual)
625   end
626
627   test "#scores_mean without qual raises" do
628     @entry.qual = nil
629     assert_raise(SeqError) { @entry.scores_mean }
630   end
631
632   test "#scores_mean returns correctly" do
633     @entry.qual = '!!II'
634     assert_equal(20.0, @entry.scores_mean)
635   end
636 end