9 test_begin( -tests => 890,
10 -requires_module => 'DB_File' );
13 use_ok('Bio::LocatableSeq');
14 use_ok('Bio::Seq::Quality');
15 use_ok('Bio::Assembly::IO');
16 use_ok('Bio::Assembly::Singlet');
23 my ($seq_id, $seq_str, $qual, $start, $end, $strand) = ('seq1', 'CAGT-GGT',
24 '0 1 2 3 4 5 6 7', 1, 7, -1, );
26 my $seq = Bio::PrimarySeq->new(-seq => $seq_str, -id => $seq_id);
27 my $singlet_id = 'singlet1';
28 ok my $singlet = Bio::Assembly::Singlet->new( -id => $singlet_id, -seqref =>
29 $seq), 'singlet from Bio::PrimarySeq';
30 isa_ok $singlet, 'Bio::Assembly::Contig';
31 isa_ok $singlet, 'Bio::Assembly::Singlet';
32 is $singlet->id, $singlet_id;
33 ok my $consensus = $singlet->get_consensus_sequence;
34 isa_ok $consensus, 'Bio::LocatableSeq';
35 is $consensus->seq, $seq_str;
36 is $consensus->start, 1;
37 is $consensus->end, 7;
38 is $consensus->strand, 1;
39 is $singlet->get_consensus_length, 8;
40 is $singlet->get_consensus_quality, undef;
41 ok my $seqref = $singlet->seqref;
42 is $seqref->id, $seq_id;
43 is $seqref->seq, $seq_str;
46 is $seqref->length, 8;
48 $seq = Bio::Seq::Quality->new(-seq => $seq_str, -id => $seq_id, -qual => $qual);
49 ok $singlet = Bio::Assembly::Singlet->new( -id => $singlet_id, -seqref => $seq),
50 'singlet from Bio::Seq::Quality';
51 is $singlet->get_consensus_length, 8;
52 ok $consensus = $singlet->get_consensus_quality;
53 isa_ok $consensus, 'Bio::Seq::QualI';
54 is join(' ', @{$consensus->qual}), $qual;
56 $seq = Bio::LocatableSeq->new( -seq => $seq_str, -id => $seq_id, -start =>
57 $start, -end => $end, -strand => $strand );
58 ok $singlet = Bio::Assembly::Singlet->new( -id => $singlet_id, -seqref => $seq),
59 'singlet from LocatableSeq';
60 ok $consensus = $singlet->get_consensus_sequence;
61 is $consensus->start, 1;
62 is $consensus->end, 7;
63 is $consensus->strand, -1;
64 ok $seqref = $singlet->seqref;
67 is $seqref->strand, -1;
69 ($start, $end) = (20, 26);
70 $seq = Bio::LocatableSeq->new( -seq => $seq_str, -id => $seq_id, -start =>
71 $start, -end => $end, -strand => $strand );
72 ok $singlet = Bio::Assembly::Singlet->new( -id => $singlet_id, -seqref => $seq),
73 'singlet from LocatableSeq with set coordinates';
74 ok $consensus = $singlet->get_consensus_sequence;
75 is $consensus->start, 1;
76 is $consensus->end, 7;
77 is $consensus->strand, -1;
78 ok $seqref = $singlet->seqref;
81 is $seqref->strand, -1;
92 ok my $aio = Bio::Assembly::IO->new(
93 -file => test_input_file('assembly_with_singlets.ace'),
94 -format => 'ace-consed',
96 is $aio->variant, 'consed', 'consed';
97 ok $aio = Bio::Assembly::IO->new(
98 -file => test_input_file('assembly_with_singlets.ace'),
101 is $aio->variant, 'consed';
102 ok $aio->variant('454');
103 is $aio->variant, '454';
109 my $in = Bio::Assembly::IO->new
110 (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
112 isa_ok $in, 'Bio::Assembly::IO';
113 while (my $contig = $in->next_contig) {
114 isa_ok($contig, 'Bio::Assembly::Contig');
117 $in = Bio::Assembly::IO->new
118 (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
120 isa_ok $in, 'Bio::Assembly::IO';
123 local $TODO = "phrap parser doesn't include the sequence string in the sequence objects.";
125 eval {$sc = $in->next_assembly};
130 $in = Bio::Assembly::IO->new
131 (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
133 ok $sc = $in->next_assembly;
134 isa_ok $sc, 'Bio::Assembly::Scaffold';
141 is $sc->id, "NoName";
142 is $sc->id('test'), "test";
144 isa_ok($sc->annotation, 'Bio::AnnotationCollectionI');
145 is $sc->annotation->get_all_annotation_keys, 0,"no annotations in Annotation collection?";
146 is $sc->get_nof_contigs, 1;
147 is $sc->get_nof_sequences_in_contigs, 2;
148 is $sc->get_nof_singlets, 2, "get_nof_singlets";
149 is $sc->get_contig_seq_ids, 2, "get_contig_seq_ids";
150 is $sc->get_contig_ids, 1, "get_contig_ids";
151 is $sc->get_singlet_ids, 2, "get_singlet_ids";
153 my @phrap_contigs = $sc->all_contigs();
154 isa_ok $phrap_contigs[0], "Bio::Assembly::Contig",'the contig is a Bio::Assembly::Contig';
155 my @singlets = $sc->all_singlets();
156 isa_ok $singlets[0], "Bio::Assembly::Contig", 'the singlet is a Bio::Assembly::Contig';
157 isa_ok $singlets[0], "Bio::Assembly::Singlet", 'the singlet is a Bio::Assembly::Singlet';
160 ok @contig_seq_ids = $sc->get_contig_seq_ids, "get_contig_seq_ids";
161 is @contig_seq_ids, 2;
162 for my $contig_seq_id (@contig_seq_ids) {
163 ok not $contig_seq_id =~ m/contig/i;
166 ok @contig_ids = $sc->get_contig_ids, "get_contig_ids";
168 for my $contig_id (@contig_ids) {
169 ok not $contig_id =~ m/contig/i;
172 ok @singlet_ids = $sc->get_singlet_ids, "get_singlet_ids";
174 for my $singlet_id (@singlet_ids) {
175 ok not $singlet_id =~ m/contig/i;
178 ok @all_seq_ids = $sc->get_all_seq_ids, "get_all_seq_ids";
179 for my $seq_id (@all_seq_ids) {
180 ok not $seq_id =~ m/contig/i;
185 # Testing ContigAnalysis
192 # ACE Consed variant (default)
193 $aio = Bio::Assembly::IO->new(
194 -file => test_input_file('consed_project','edit_dir','test_project.fasta.screen.ace.2'),
198 my $assembly = $aio->next_assembly();
200 my @contigs = $assembly->all_contigs();
202 my $direction = $contigs[0]->strand;
205 my $features = $contigs[0]->get_features_collection;
207 my @contig_features = $features->features;
208 is @contig_features, 61, 'contig features'; # 59 contig features + 2 seqfeatures
210 my @annotations = $features->get_features_by_type('Annotation');
214 foreach my $an (@annotations) {
215 if ($an->has_tag('extra_info')) {
217 is (($an->get_tag_values('extra_info'))[0], "contig extra\ninfo\n");
219 elsif ($an->has_tag('comment')){
221 is (($an->get_tag_values('comment'))[0], "contig tag\ncomment\n");
226 is $assembly->get_nof_contigs, 1;
227 is $assembly->get_nof_sequences_in_contigs, 2;
228 is $assembly->get_nof_singlets, 0, "get_nof_singlets";
229 is $assembly->get_contig_seq_ids, 2, "get_contig_seq_ids";
230 is $assembly->get_contig_ids, 1, "get_contig_ids";
231 is $assembly->get_singlet_ids, 0, "get_singlet_ids";
233 $aio = Bio::Assembly::IO->new(
234 -file => test_input_file('assembly_with_singlets.ace'),
238 while (my $obj = $aio->next_contig) {
239 isa_ok $obj, 'Bio::Assembly::Contig'; # Singlets are contigs too
242 $aio = Bio::Assembly::IO->new(
243 -file => test_input_file('assembly_with_singlets.ace'),
247 $assembly = $aio->next_assembly();
248 is $assembly->get_nof_contigs, 3;
249 my @ace_contigs = $assembly->all_contigs();
250 my $ace_contig = $ace_contigs[0];
251 isa_ok $ace_contig, "Bio::Assembly::Contig",'the contig is a Bio::Assembly::Contig';
253 ok my @test_reads = $ace_contig->get_seq_ids();
254 is scalar @test_reads, 2;
255 is $test_reads[0], '5704073';
256 is $test_reads[1], '5762101';
258 is $assembly->get_nof_sequences_in_contigs, 6;
259 is $assembly->get_nof_singlets, 33, "get_nof_singlets";
260 @singlets = $assembly->all_singlets();
261 isa_ok $singlets[0], "Bio::Assembly::Contig", 'the singlet is a Bio::Assembly::Contig';
262 isa_ok $singlets[0], "Bio::Assembly::Singlet", 'the singlet is a Bio::Assembly::Singlet';
263 ok @contig_seq_ids = $assembly->get_contig_seq_ids, "get_contig_seq_ids";
264 is @contig_seq_ids, 6;
265 for my $contig_seq_id (@contig_seq_ids) {
266 ok not $contig_seq_id =~ m/contig/i;
268 ok @contig_ids = $assembly->get_contig_ids, "get_contig_ids";
270 for my $contig_id (@contig_ids) {
271 ok $contig_id =~ m/contig/i;
273 ok @singlet_ids = $assembly->get_singlet_ids, "get_singlet_ids";
275 for my $singlet_id (@singlet_ids) {
276 ok $singlet_id =~ m/contig/i;
278 ok @all_seq_ids = $assembly->get_all_seq_ids, "get_all_seq_ids";
279 for my $seq_id (@all_seq_ids) {
280 ok not $seq_id =~ m/contig/i;
282 is(@all_seq_ids, 39);
285 ok $aio = Bio::Assembly::IO->new(
286 -file=>test_input_file('singlet_w_CT.ace'),
291 $aio = Bio::Assembly::IO->new(
292 -file=>test_input_file('27-contig_Newbler.ace'),
295 $assembly = $aio->next_assembly();
296 @contigs = $assembly->all_contigs();
297 # All read positions should be >0
298 my $contig = $contigs[0];
299 my $min_aln_coord = undef;
300 for my $read ($contig->each_seq) {
301 my ($feat) = $contig->get_features_collection->get_features_by_type("_aligned_coord:".$read->id);
302 my $aln_coord_start = $feat->location->start;
303 if ( (not defined $min_aln_coord) or ($aln_coord_start < $min_aln_coord) ) {
304 $min_aln_coord = $aln_coord_start;
307 is ($min_aln_coord, 1, '454 ACE variant coordinates check');
308 # The ends of the consensus should be padded
309 my $left_pad_length = 29;
310 my $sequence_length = 203;
311 my $right_pad_length = 81;
312 my $consensus_length = $left_pad_length + $sequence_length + $right_pad_length;
313 my $cons_seq = $contig->get_consensus_sequence->seq;
314 is length $cons_seq, $consensus_length;
315 $cons_seq =~ m/^(-*).*?(-*)$/;
316 is length $1, $left_pad_length, '454 ACE variant consensus check';
317 is length $2, $right_pad_length;
318 my $cons_qual = $contig->get_consensus_quality->qual;
319 is scalar @$cons_qual, $consensus_length;
320 $cons_qual = join ' ', @{$contig->get_consensus_quality->qual};
321 my $lpad = $left_pad_length x '0 ';
322 my $rpad = $right_pad_length x '0 ';
323 $cons_qual =~ m/^($lpad).*($rpad)$/;
328 my $asm_infile = '27-contig_Newbler.ace';
329 my $asm_outfile = test_output_file().'.ace';
330 my $asm_out = Bio::Assembly::IO->new(
331 -file => ">$asm_outfile",
335 ok $asm_in = Bio::Assembly::IO->new(
336 -file => test_input_file($asm_infile),
339 )->next_assembly, 'writing in the ACE format';
340 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1 );
342 $asm_infile = 'assembly_with_singlets.ace';
343 ok $asm_in = Bio::Assembly::IO->new(
344 -file => test_input_file($asm_infile),
347 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1 );
349 $asm_infile = 'reference_ace.ace';
350 ok $asm_in = Bio::Assembly::IO->new(
351 -file => test_input_file($asm_infile),
354 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1 );
358 # Testing TIGR format
361 # Importing an assembly
363 $asm_in = Bio::Assembly::IO->new(
364 -file => test_input_file("sample_dataset.tigr "),
367 while (my $obj = $asm_in->next_contig) {
368 isa_ok $obj, 'Bio::Assembly::Contig'; # Singlets are contigs too
371 $asm_in = Bio::Assembly::IO->new(
372 -file => test_input_file("sample_dataset.tigr "),
376 my $scaf_in = $asm_in->next_assembly;
377 isa_ok $scaf_in, 'Bio::Assembly::Scaffold';
378 is $scaf_in->id, 'NoName';
379 is $scaf_in->get_nof_contigs, 13;
380 is $scaf_in->get_nof_sequences_in_contigs, 36;
381 is $scaf_in->get_nof_singlets, 1;
382 my @contigseqids = sort qw(sdsu|SDSU1_RFPERU_001_A09.x01.phd.1
383 sdsu|SDSU1_RFPERU_001_B03.x01.phd.1 sdsu|SDSU1_RFPERU_001_B04.x01.phd.1
384 sdsu|SDSU1_RFPERU_001_E04.x01.phd.1 sdsu|SDSU_RFPERU_002_A01.x01.phd.1
385 sdsu|SDSU_RFPERU_002_B07.x01.phd.1 sdsu|SDSU_RFPERU_002_C12.x01.phd.1
386 sdsu|SDSU_RFPERU_002_D08.x01.phd.1 sdsu|SDSU_RFPERU_002_H12.x01.phd.1
387 sdsu|SDSU_RFPERU_003_G09.x01.phd.1 sdsu|SDSU_RFPERU_004_H12.x01.phd.1
388 sdsu|SDSU_RFPERU_005_F02.x01.phd.1 sdsu|SDSU_RFPERU_006_D03.x01.phd.1
389 sdsu|SDSU_RFPERU_006_E04.x01.phd.1 sdsu|SDSU_RFPERU_006_E05.x01.phd.1
390 sdsu|SDSU_RFPERU_006_H08.x01.phd.1 sdsu|SDSU_RFPERU_007_E09.x01.phd.1
391 sdsu|SDSU_RFPERU_007_F06.x01.phd.1 sdsu|SDSU_RFPERU_008_B02.x01.phd.1
392 sdsu|SDSU_RFPERU_009_E07.x01.phd.1 sdsu|SDSU_RFPERU_010_B05.x01.phd.1
393 sdsu|SDSU_RFPERU_010_B06.x01.phd.1 sdsu|SDSU_RFPERU_010_C09.x01.phd.1
394 sdsu|SDSU_RFPERU_010_D10.x01.phd.1 sdsu|SDSU_RFPERU_012_H02.x01.phd.1
395 sdsu|SDSU_RFPERU_013_B05.x01.phd.1 sdsu|SDSU_RFPERU_013_C07.x01.phd.1
396 sdsu|SDSU_RFPERU_013_C08.x01.phd.1 sdsu|SDSU_RFPERU_013_G10.x01.phd.1
397 sdsu|SDSU_RFPERU_013_H05.x01.phd.1 sdsu|SDSU_RFPERU_014_H06.x01.phd.1
398 sdsu|SDSU_RFPERU_015_A05.x01.phd.1 sdsu|SDSU_RFPERU_015_C06.x01.phd.1
399 sdsu|SDSU_RFPERU_015_E04.x01.phd.1 sdsu|SDSU_RFPERU_015_G04.x01.phd.1
400 sdsu|SDSU_RFPERU_015_H03.x01.phd.1);
401 my @contigids = sort qw(106 144 148 17 185 2 210 36 453 500 613 668 93);
402 my @singletids = sort qw(123);
403 my @singletseqids = sort qw(asdf);
404 is_deeply [sort $scaf_in->get_contig_seq_ids], \@contigseqids;
405 is_deeply [sort $scaf_in->get_contig_ids], \@contigids ;
406 is_deeply [sort $scaf_in->get_singlet_ids], \@singletids ;
407 isa_ok $scaf_in->get_seq_by_id('sdsu|SDSU1_RFPERU_001_A09.x01.phd.1'),'Bio::LocatableSeq';
408 $contig = $scaf_in->get_contig_by_id('106');
409 isa_ok $contig,'Bio::Assembly::Contig';
411 # check Contig object SeqFeature::Collection
412 # should add more specific Contig tests...
413 my @sfs = $contig->get_features_collection->features; # 5 contig features + 2 seqfeatures
415 is $sfs[1]->seq_id(), undef; # should this be undef?
416 ok $contig->get_features_collection->get_features_by_type('_aligned_coord:sdsu|SDSU_RFPERU_006_E04.x01.phd.1');
417 isa_ok $scaf_in->annotation, 'Bio::AnnotationCollectionI';
418 is $scaf_in->annotation->get_all_annotation_keys, 0, "no annotations in Annotation collection?";
421 # Exporting an assembly
422 $asm_outfile = test_output_file();
423 $asm_out = Bio::Assembly::IO->new(
424 -file=> ">$asm_outfile",
427 ok $asm_out->write_assembly( -scaffold => $scaf_in), 'writing in the TIGR format';
434 my $file = 'test.maq';
435 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
436 -format => 'maq' ), "init maq IO object";
437 ok $assembly = $aio->next_assembly, "get maq assy";
438 is $assembly->get_nof_contigs, 11, "got all contigs";
439 ok open(my $tf, test_input_file($file)), "read test file as text";
441 is $assembly->get_nof_contig_seqs, scalar @lines, "recorded all maq reads";
442 ok !$assembly->get_nof_singlets, "no singlets";
444 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
446 isa_ok $aio, 'Bio::Assembly::IO';
447 while (my $contig = $aio->next_contig) {
448 isa_ok $contig, 'Bio::Assembly::Contig';
452 # Testing maq with singlets
454 $file = 'test_singlets.maq';
455 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
457 ok $assembly = $aio->next_assembly, "get maq assy";
458 isa_ok $aio, 'Bio::Assembly::IO';
461 ok @contig_seq_ids = $assembly->get_contig_seq_ids, "get_contig_seq_ids";
462 is @contig_seq_ids, 246;
463 for my $contig_seq_id (@contig_seq_ids) {
464 ok not $contig_seq_id =~ m/maq_assy/i;
467 ok @contig_ids = $assembly->get_contig_ids, "get_contig_ids";
469 for my $contig_id (@contig_ids) {
470 ok $contig_id =~ m/maq_assy/i;
473 ok @singlet_ids = $assembly->get_singlet_ids, "get_singlet_ids";
475 for my $singlet_id (@singlet_ids) {
476 ok $singlet_id =~ m/maq_assy/i;
479 ok @all_seq_ids = $assembly->get_all_seq_ids, "get_all_seq_ids";
480 for my $seq_id (@all_seq_ids) {
481 ok not $seq_id =~ m/maq_assy/i;
483 is @all_seq_ids, 250;
485 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
487 while (my $contig = $aio->next_contig) {
488 isa_ok $contig, 'Bio::Assembly::Contig';