t/Assembly/core.t

   1 use strict;
   2 use warnings;
   3 my %ASSEMBLY_TESTS;
   4
   5 BEGIN {
   6     use lib '.';
   7     use Bio::Root::Test;
   8
   9     test_begin( -tests => 890,
  10                 -requires_module => 'DB_File' );
  11
  12     use_ok('Bio::Seq');
  13     use_ok('Bio::LocatableSeq');
  14     use_ok('Bio::Seq::Quality');
  15     use_ok('Bio::Assembly::IO');
  16     use_ok('Bio::Assembly::Singlet');
  17 }
  18 use Bio::Root::IO;
  19
  20 #
  21 # Testing Singlet
  22 #
  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, );
  25
  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;
  44 is $seqref->start, 1;
  45 is $seqref->end, 7;
  46 is $seqref->length, 8;
  47
  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;
  55
  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;
  65 is $seqref->start, 1;
  66 is $seqref->end, 7;
  67 is $seqref->strand, -1;
  68
  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;
  79 is $seqref->start, 1;
  80 is $seqref->end, 7;
  81 is $seqref->strand, -1;
  82
  83 #
  84 # Testing Contig
  85 #
  86
  87 #
  88 # Testing IO
  89 #
  90
  91 # ACE variants
  92 ok my $aio = Bio::Assembly::IO->new(
  93     -file   => test_input_file('assembly_with_singlets.ace'),
  94     -format => 'ace-consed',
  95 );
  96 is $aio->variant, 'consed', 'consed';
  97 ok $aio = Bio::Assembly::IO->new(
  98     -file   => test_input_file('assembly_with_singlets.ace'),
  99     -format => 'ace',
 100 );
 101 is $aio->variant, 'consed';
 102 ok $aio->variant('454');
 103 is $aio->variant, '454';
 104
 105 #
 106 # Some PHRAP input
 107 #
 108
 109 my $in = Bio::Assembly::IO->new
 110     (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
 111      -verbose => -1);
 112 isa_ok $in, 'Bio::Assembly::IO';
 113 while (my $contig = $in->next_contig) {
 114     isa_ok($contig, 'Bio::Assembly::Contig');
 115 }
 116
 117 $in = Bio::Assembly::IO->new
 118     (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
 119      -verbose => -1);
 120 isa_ok $in, 'Bio::Assembly::IO';
 121 my $sc;
 122 TODO: {
 123     local $TODO = "phrap parser doesn't include the sequence string in the sequence objects.";
 124     $in->verbose(2);
 125     eval {$sc = $in->next_assembly};
 126     ok !$@;
 127 }
 128
 129 $in->verbose(-1);
 130 $in = Bio::Assembly::IO->new
 131     (-file => test_input_file('consed_project','edit_dir','test_project.phrap.out'),
 132      -verbose => -1);
 133 ok $sc = $in->next_assembly;
 134 isa_ok $sc, 'Bio::Assembly::Scaffold';
 135
 136
 137 #
 138 # Testing Scaffold
 139 #
 140
 141 is $sc->id, "NoName";
 142 is $sc->id('test'), "test";
 143
 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";
 152
 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';
 158
 159 my @contig_seq_ids;
 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;
 164 }
 165 my @contig_ids;
 166 ok @contig_ids = $sc->get_contig_ids, "get_contig_ids";
 167 is @contig_ids, 1;
 168 for my $contig_id (@contig_ids) {
 169   ok not $contig_id =~ m/contig/i;
 170 }
 171 my @singlet_ids;
 172 ok @singlet_ids = $sc->get_singlet_ids, "get_singlet_ids";
 173 is @singlet_ids, 2;
 174 for my $singlet_id (@singlet_ids) {
 175   ok not $singlet_id =~ m/contig/i;
 176 }
 177 my @all_seq_ids;
 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;
 181 }
 182 is @all_seq_ids, 4;
 183
 184 #
 185 # Testing ContigAnalysis
 186 #
 187
 188 #
 189 # Testing ACE
 190 #
 191
 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'),
 195     -format => 'ace'
 196 );
 197
 198 my $assembly = $aio->next_assembly();
 199
 200 my @contigs = $assembly->all_contigs();
 201
 202 my $direction = $contigs[0]->strand;
 203 is $direction, 1;
 204
 205 my $features =  $contigs[0]->get_features_collection;
 206
 207 my @contig_features = $features->features;
 208 is @contig_features, 61, 'contig features'; # 59 contig features + 2 seqfeatures
 209
 210 my @annotations = $features->get_features_by_type('Annotation');
 211 is @annotations, 2;
 212
 213 my $had_tag = 0;
 214 foreach my $an (@annotations) {
 215     if ($an->has_tag('extra_info')) {
 216         $had_tag++;
 217         is (($an->get_tag_values('extra_info'))[0], "contig extra\ninfo\n");
 218     }
 219     elsif ($an->has_tag('comment')){
 220         $had_tag++;
 221         is (($an->get_tag_values('comment'))[0], "contig tag\ncomment\n");
 222     }
 223 }
 224 is $had_tag, 2;
 225
 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";
 232
 233 $aio = Bio::Assembly::IO->new(
 234     -file   => test_input_file('assembly_with_singlets.ace'),
 235     -format => 'ace',
 236 );
 237
 238 while (my $obj = $aio->next_contig) {
 239     isa_ok $obj, 'Bio::Assembly::Contig'; # Singlets are contigs too
 240 }
 241
 242 $aio = Bio::Assembly::IO->new(
 243     -file   => test_input_file('assembly_with_singlets.ace'),
 244     -format => 'ace',
 245 );
 246
 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';
 252
 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';
 257
 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;
 267 }
 268 ok @contig_ids = $assembly->get_contig_ids, "get_contig_ids";
 269 is @contig_ids, 3;
 270 for my $contig_id (@contig_ids) {
 271   ok $contig_id =~ m/contig/i;
 272 }
 273 ok @singlet_ids = $assembly->get_singlet_ids, "get_singlet_ids";
 274 is @singlet_ids, 33;
 275 for my $singlet_id (@singlet_ids) {
 276   ok $singlet_id =~ m/contig/i;
 277 }
 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;
 281 }
 282 is(@all_seq_ids, 39);
 283
 284 # bug 2758
 285 ok $aio = Bio::Assembly::IO->new(
 286     -file=>test_input_file('singlet_w_CT.ace'),
 287     -format=>'ace'
 288 );
 289
 290 # ACE 454 variant
 291 $aio = Bio::Assembly::IO->new(
 292     -file=>test_input_file('27-contig_Newbler.ace'),
 293     -format=>'ace-454'
 294 );
 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;
 305    }
 306 }
 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)$/;
 324 ok defined $1;
 325 ok defined $2;
 326
 327 # Writing ACE files
 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",
 332     -format  =>'ace',
 333 );
 334 my $asm_in;
 335 ok $asm_in = Bio::Assembly::IO->new(
 336     -file    => test_input_file($asm_infile),
 337     -format  => 'ace',
 338     -variant => '454',
 339 )->next_assembly, 'writing in the ACE format';
 340 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1 );
 341
 342 $asm_infile = 'assembly_with_singlets.ace';
 343 ok $asm_in = Bio::Assembly::IO->new(
 344     -file   => test_input_file($asm_infile),
 345     -format => 'ace',
 346 )->next_assembly;
 347 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1  );
 348
 349 $asm_infile = 'reference_ace.ace';
 350 ok $asm_in = Bio::Assembly::IO->new(
 351     -file   => test_input_file($asm_infile),
 352     -format => 'ace',
 353 )->next_assembly;
 354 ok $asm_out->write_assembly( -scaffold => $asm_in, -singlets => 1  );
 355
 356
 357 #
 358 # Testing TIGR format
 359 #
 360
 361 # Importing an assembly
 362
 363 $asm_in = Bio::Assembly::IO->new(
 364     -file => test_input_file("sample_dataset.tigr "),
 365     -format=>'tigr'
 366 );
 367 while (my $obj = $asm_in->next_contig) {
 368     isa_ok $obj, 'Bio::Assembly::Contig'; # Singlets are contigs too
 369 }
 370
 371 $asm_in = Bio::Assembly::IO->new(
 372     -file => test_input_file("sample_dataset.tigr "),
 373     -format=>'tigr'
 374 );
 375
 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';
 410
 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
 414 is scalar @sfs, 7;
 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?";
 419
 420
 421 # Exporting an assembly
 422 $asm_outfile = test_output_file();
 423 $asm_out = Bio::Assembly::IO->new(
 424     -file=> ">$asm_outfile",
 425     -format=>'tigr'
 426 );
 427 ok $asm_out->write_assembly( -scaffold => $scaf_in), 'writing in the TIGR format';
 428
 429
 430 #
 431 # Testing maq
 432 # /maj
 433 #
 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";
 440 my @lines = <$tf>;
 441 is $assembly->get_nof_contig_seqs, scalar @lines, "recorded all maq reads";
 442 ok !$assembly->get_nof_singlets, "no singlets";
 443
 444 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
 445                                   -format => 'maq' );
 446 isa_ok $aio, 'Bio::Assembly::IO';
 447 while (my $contig = $aio->next_contig) {
 448     isa_ok $contig, 'Bio::Assembly::Contig';
 449 }
 450
 451 #
 452 # Testing maq with singlets
 453 #
 454 $file = 'test_singlets.maq';
 455 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
 456                                   -format => 'maq' );
 457 ok $assembly = $aio->next_assembly, "get maq assy";
 458 isa_ok $aio, 'Bio::Assembly::IO';
 459
 460
 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;
 465 }
 466
 467 ok @contig_ids = $assembly->get_contig_ids, "get_contig_ids";
 468 is @contig_ids, 37;
 469 for my $contig_id (@contig_ids) {
 470   ok $contig_id =~ m/maq_assy/i;
 471 }
 472
 473 ok @singlet_ids = $assembly->get_singlet_ids, "get_singlet_ids";
 474 is @singlet_ids, 4;
 475 for my $singlet_id (@singlet_ids) {
 476   ok $singlet_id =~ m/maq_assy/i;
 477 }
 478
 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;
 482 }
 483 is @all_seq_ids, 250;
 484
 485 ok $aio = Bio::Assembly::IO->new( -file => test_input_file($file),
 486                                   -format => 'maq' );
 487 while (my $contig = $aio->next_contig) {
 488     isa_ok $contig, 'Bio::Assembly::Contig';
 489 }
 490
 491 exit;