t/Align/SimpleAlign.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use Bio::Root::Test;
   8
   9     test_begin( -tests => 206 );
  10
  11     use_ok('Bio::SimpleAlign');
  12     use_ok('Bio::AlignIO');
  13     use_ok('Bio::SeqFeature::Generic');
  14     use_ok('Bio::Location::Simple');
  15     use_ok('Bio::Location::Split');
  16 }
  17
  18 my $DEBUG = test_debug();
  19
  20 my ( $str, $aln, @seqs, $seq );
  21
  22 $str = Bio::AlignIO->new( -file => test_input_file('testaln.pfam') );
  23 isa_ok( $str, 'Bio::AlignIO' );
  24 $aln = $str->next_aln();
  25 is $aln->get_seq_by_pos(1)->get_nse, '1433_LYCES/9-246', "pfam input test";
  26
  27 my $aln1 = $aln->remove_columns( ['mismatch'] );
  28 is(
  29     $aln1->match_line,
  30     '::*::::*:**:*:*:***:**.***::*.*::**::**:***..**:'
  31       . '*:*.::::*:.:*.*.**:***.**:*.:.**::**.*:***********:::*:.:*:**.*::*:'
  32       . '.*.:*:**:****************::',
  33     'match_line'
  34 );
  35
  36 my $aln2 = $aln->select( 1, 3 );
  37 isa_ok( $aln2, 'Bio::Align::AlignI' );
  38 is( $aln2->num_sequences, 3, 'num_sequences' );
  39
  40 # test select non contiguous-sorted by default
  41 $aln2 = $aln->select_noncont( 6, 7, 8, 9, 10, 1, 2, 3, 4, 5 );
  42 is( $aln2->num_sequences, 10, 'num_sequences' );
  43 is(
  44     $aln2->get_seq_by_pos(2)->id,
  45     $aln->get_seq_by_pos(2)->id,
  46     'select_noncont'
  47 );
  48 is(
  49     $aln2->get_seq_by_pos(8)->id,
  50     $aln->get_seq_by_pos(8)->id,
  51     'select_noncont'
  52 );
  53
  54 # test select non contiguous-nosort option
  55 $aln2 = $aln->select_noncont( 'nosort', 6, 7, 8, 9, 10, 1, 2, 3, 4, 5 );
  56 is( $aln2->num_sequences, 10, 'num_sequences' );
  57 is(
  58     $aln2->get_seq_by_pos(2)->id,
  59     $aln->get_seq_by_pos(7)->id,
  60     'select_noncont'
  61 );
  62 is(
  63     $aln2->get_seq_by_pos(8)->id,
  64     $aln->get_seq_by_pos(3)->id,
  65     'select_noncont'
  66 );
  67
  68 # test select non contiguous by name
  69 my $aln3 = $aln->select_noncont_by_name('1433_LYCES','BMH1_YEAST','143T_HUMAN');
  70 is( $aln3->num_sequences, 3, 'select_noncont_by_name' );
  71 my @seqs3 = $aln3->each_seq();
  72 is $seqs3[0]->id, '1433_LYCES', 'select_noncont_by_name';
  73 is $seqs3[1]->id, 'BMH1_YEAST', 'select_noncont_by_name';
  74 is $seqs3[2]->id, '143T_HUMAN', 'select_noncont_by_name';
  75
  76
  77 @seqs = $aln->each_seq();
  78 is scalar @seqs, 16, 'each_seq';
  79 is $seqs[0]->get_nse, '1433_LYCES/9-246', 'get_nse';
  80 is $seqs[0]->id,      '1433_LYCES',       'id';
  81 is $seqs[0]->num_gaps, 3,                  'num_gaps';
  82 @seqs = $aln->each_alphabetically();
  83 is scalar @seqs, 16, 'each_alphabetically';
  84
  85 is $aln->column_from_residue_number( '1433_LYCES', 10 ), 2,
  86   'column_from_residue_number';
  87 is $aln->displayname( '1433_LYCES/9-246', 'my_seq' ), 'my_seq',
  88   'display_name get/set';
  89 is $aln->displayname('1433_LYCES/9-246'), 'my_seq', 'display_name get';
  90 is substr( $aln->consensus_string(50), 0, 60 ),
  91   "RE??VY?AKLAEQAERYEEMV??MK?VAE??????ELSVEERNLLSVAYKNVIGARRASW",
  92   'consensus_string';
  93 is substr( $aln->consensus_string(100), 0, 60 ),
  94   "?????????L????E????M???M????????????L??E?RNL?SV?YKN??G??R??W",
  95   'consensus_string';
  96 is substr( $aln->consensus_string(0), 0, 60 ),
  97   "REDLVYLAKLAEQAERYEEMVEFMKKVAELGAPAEELSVEERNLLSVAYKNVIGARRASW",
  98   'consensus_string';
  99
 100 ok( @seqs = $aln->each_seq_with_id('143T_HUMAN') );
 101 is scalar @seqs, 1, 'each_seq_with_id';
 102
 103 is $aln->is_flush, 1, 'is_flush';
 104 ok( $aln->id('x') && $aln->id eq 'x', 'id get/set' );
 105
 106 is $aln->length,        242,  'length';
 107 is $aln->num_residues,  3769, 'num_residues';
 108 is $aln->num_sequences, 16,   'num_sequences';
 109 is( sprintf( "%.2f", $aln->overall_percentage_identity() ),
 110     33.06, 'overall_percentage_identity' );
 111 is( sprintf( "%.2f", $aln->overall_percentage_identity('align') ),
 112     33.06, 'overall_percentage_identity (align)' );
 113 is( sprintf( "%.2f", $aln->overall_percentage_identity('short') ),
 114     35.24, 'overall_percentage_identity (short)' );
 115 is( sprintf( "%.2f", $aln->overall_percentage_identity('long') ),
 116     33.47, 'overall_percentage_identity (long)' );
 117 is( sprintf( "%.2f", $aln->average_percentage_identity() ),
 118     66.91, 'average_percentage_identity' );
 119
 120 ok $aln->set_displayname_count;
 121 is $aln->displayname('1433_LYCES/9-246'), '1433_LYCES_1',
 122   'set_displayname_count';
 123 ok $aln->set_displayname_flat;
 124 is $aln->displayname('1433_LYCES/9-246'), '1433_LYCES', 'set_displayname_flat';
 125 ok $aln->set_displayname_normal;
 126 is $aln->displayname('1433_LYCES/9-246'), '1433_LYCES/9-246',
 127   'set_displayname_normal';
 128 ok $aln->uppercase;
 129 ok $aln->map_chars( '\.', '-' );
 130 @seqs = $aln->each_seq_with_id('143T_HUMAN');
 131 is substr( $seqs[0]->seq, 0, 60 ),
 132   'KTELIQKAKLAEQAERYDDMATCMKAVTEQGA---ELSNEERNLLSVAYKNVVGGRRSAW',
 133   'uppercase, map_chars';
 134
 135 is(
 136     $aln->match_line,
 137     '       ::*::::*  : *   *:           *: *:***:**.***::*.'
 138       . ' *::**::**:***      .  .      **  :* :*   .  :: ::   *:  .     :* .*. **:'
 139       . '***.** :*.            :  .*  *   :   : **.*:***********:::* : .: *  :** .'
 140       . '*::*: .*. : *: **:****************::     ',
 141     'match_line'
 142 );
 143 ok $aln->remove_seq( $seqs[0] ), 'remove_seqs';
 144 is $aln->num_sequences, 15, 'remove_seqs';
 145 ok $aln->add_seq( $seqs[0] ), 'add_seq';
 146 is $aln->num_sequences, 16, 'add_seq';
 147 ok $seq = $aln->get_seq_by_pos(1), 'get_seq_by_pos';
 148 is( $seq->id, '1433_LYCES', 'get_seq_by_pos' );
 149 ok( ( $aln->missing_char(), 'P' ) and ( $aln->missing_char('X'), 'X' ) );
 150 ok( ( $aln->match_char(),   '.' ) and ( $aln->match_char('-'),   '-' ) );
 151 ok( ( $aln->gap_char(),     '-' ) and ( $aln->gap_char('.'),     '.' ) );
 152
 153 is $aln->purge(0.7), 12, 'purge';
 154 is $aln->num_sequences, 4, 'purge';
 155
 156 SKIP: {
 157     test_skip( -tests => 24, -requires_module => 'IO::String' );
 158
 159     my $string;
 160     my $out = IO::String->new($string);
 161
 162     my $s1 = Bio::LocatableSeq->new(
 163         -id       => 'AAA',
 164         -seq      => 'aawtat-tn-',
 165         -start    => 1,
 166         -end      => 8,
 167         -alphabet => 'dna'
 168     );
 169     my $s2 = Bio::LocatableSeq->new(
 170         -id       => 'BBB',
 171         -seq      => '-aaaat-tt-',
 172         -start    => 1,
 173         -end      => 7,
 174         -alphabet => 'dna'
 175     );
 176     $a = Bio::SimpleAlign->new();
 177     $a->add_seq($s1);
 178     $a->add_seq($s2);
 179
 180     is( $a->consensus_iupac, "aAWWAT-TN-", 'IO::String consensus_iupac' );
 181     $s1->seq('aaaaattttt');
 182     $s1->alphabet('dna');
 183     $s1->end(10);
 184     $s2->seq('-aaaatttt-');
 185     $s2->end(8);
 186     $a = Bio::SimpleAlign->new();
 187     $a->add_seq($s1);
 188     $a->add_seq($s2);
 189
 190     my $strout = Bio::AlignIO->new( -fh => $out, -format => 'pfam' );
 191     $strout->write_aln($a);
 192     is(
 193         $string,
 194         "AAA/1-10    aaaaattttt\n" . "BBB/1-8     -aaaatttt-\n",
 195         'IO::String write_aln normal'
 196     );
 197
 198     $out->setpos(0);
 199     $string = '';
 200     my $b = $a->slice( 2, 9 );
 201     $strout->write_aln($b);
 202     is $string,
 203       "AAA/2-9    aaaatttt\n" . "BBB/1-8    aaaatttt\n",
 204       'IO::String write_aln slice';
 205
 206     $out->setpos(0);
 207     $string = '';
 208     $b = $a->slice( 9, 10 );
 209     $strout->write_aln($b);
 210     is $string,
 211       "AAA/9-10    tt\n" . "BBB/8-8     t-\n",
 212       'IO::String write_aln slice';
 213
 214     $a->verbose(-1);
 215     $out->setpos(0);
 216     $string = '';
 217     $b = $a->slice( 1, 2 );
 218     $strout->write_aln($b);
 219     is $string,
 220       "AAA/1-2    aa\n" . "BBB/1-1    -a\n",
 221       'IO::String write_aln slice';
 222
 223     # not sure what coordinates this should return...
 224     $a->verbose(-1);
 225     $out->setpos(0);
 226     $string = '';
 227     $b = $a->slice( 1, 1, 1 );
 228     $strout->write_aln($b);
 229     is $string,
 230       "AAA/1-1    a\n" . "BBB/1-0    -\n",
 231       'IO::String write_aln slice';
 232
 233     $a->verbose(-1);
 234     $out->setpos(0);
 235     $string = '';
 236     $b = $a->slice( 2, 2 );
 237     $strout->write_aln($b);
 238     is $string,
 239       "AAA/2-2    a\n" . "BBB/1-1    a\n",
 240       'IO::String write_aln slice';
 241
 242     eval { $b = $a->slice( 11, 13 ); };
 243
 244     like( $@, qr/EX/ );
 245
 246     # remove_columns by position
 247     $out->setpos(0);
 248     $string = '';
 249     $str    = Bio::AlignIO->new( -file => test_input_file('mini-align.aln') );
 250     $aln1   = $str->next_aln;
 251     $aln2   = $aln1->remove_columns( [ 0, 0 ] );
 252     $strout->write_aln($aln2);
 253     is $string,
 254         "P84139/2-33              NEGEHQIKLDELFEKLLRARLIFKNKDVLRRC\n"
 255       . "P814153/2-33             NEGMHQIKLDVLFEKLLRARLIFKNKDVLRRC\n"
 256       . "BAB68554/1-14            ------------------AMLIFKDKQLLQQC\n"
 257       . "gb|443893|124775/1-32    MRFRFQIKVPPAVEGARPALLIFKSRPELGGC\n",
 258       'remove_columns by position';
 259
 260     # and when arguments are entered in "wrong order"?
 261     $out->setpos(0);
 262     $string = '';
 263     my $aln3 = $aln1->remove_columns( [ 1, 1 ], [ 30, 30 ], [ 5, 6 ] );
 264     $strout->write_aln($aln3);
 265     is $string,
 266         "P84139/1-33              MEGEIKLDELFEKLLRARLIFKNKDVLRC\n"
 267       . "P814153/1-33             MEGMIKLDVLFEKLLRARLIFKNKDVLRC\n"
 268       . "BAB68554/1-14            ----------------AMLIFKDKQLLQC\n"
 269       . "gb|443893|124775/2-32    -RFRIKVPPAVEGARPALLIFKSRPELGC\n",
 270       'remove_columns by position (wrong order)';
 271
 272     my %cigars = $aln1->cigar_line;
 273     is $cigars{'gb|443893|124775/1-32'}, '19,19:21,24:29,29:32,32',
 274       'cigar_line';
 275     is $cigars{'P814153/1-33'},  '20,20:22,25:30,30:33,33', 'cigar_line';
 276     is $cigars{'BAB68554/1-14'}, '1,1:3,6:11,11:14,14',     'cigar_line';
 277     is $cigars{'P84139/1-33'},   '20,20:22,25:30,30:33,33', 'cigar_line';
 278
 279     # sort_alphabetically
 280     my $s3 = Bio::LocatableSeq->new(
 281         -id       => 'ABB',
 282         -seq      => '-attat-tt-',
 283         -start    => 1,
 284         -end      => 7,
 285         -alphabet => 'dna'
 286     );
 287     $a->add_seq($s3);
 288
 289     is $a->get_seq_by_pos(2)->id, "BBB", 'sort_alphabetically - before';
 290     ok $a->sort_alphabetically;
 291     is $a->get_seq_by_pos(2)->id, "ABB", 'sort_alphabetically - after';
 292
 293     $b = $a->remove_gaps();
 294     is $b->consensus_string, "aaaattt", 'remove_gaps';
 295
 296     $s1->seq('aaaaattt--');
 297
 298     $b = $a->remove_gaps( undef, 'all_gaps_only' );
 299     is $b->consensus_string, "aaaaatttt", 'remove_gaps all_gaps_only';
 300
 301     # test set_new_reference:
 302     $str = Bio::AlignIO->new( -file => test_input_file('testaln.clustalw') );
 303     $aln = $str->next_aln();
 304     my $new_aln = $aln->set_new_reference(3);
 305     $a       = $new_aln->get_seq_by_pos(1)->display_id;
 306     $new_aln = $aln->set_new_reference('P851414');
 307     $b       = $new_aln->get_seq_by_pos(1)->display_id;
 308     is $a, 'P851414', 'set_new_reference';
 309     is $b, 'P851414', 'set_new_reference';
 310
 311     # test uniq_seq:
 312     $str = Bio::AlignIO->new(
 313         -verbose => $DEBUG,
 314         -file    => test_input_file('testaln2.fasta')
 315     );
 316     $aln     = $str->next_aln();
 317     $new_aln = $aln->uniq_seq();
 318     $a       = $new_aln->num_sequences;
 319     is $a, 11, 'uniq_seq';
 320
 321     # check if slice works well with a LocateableSeq in its negative strand
 322     my $seq1 = Bio::LocatableSeq->new(
 323         -SEQ    => "ATGCTG-ATG",
 324         -START  => 1,
 325         -END    => 9,
 326         -ID     => "test1",
 327         -STRAND => -1
 328     );
 329
 330     my $seq2 = Bio::LocatableSeq->new(
 331         -SEQ    => "A-GCTGCATG",
 332         -START  => 1,
 333         -END    => 9,
 334         -ID     => "test2",
 335         -STRAND => 1
 336     );
 337
 338     $string = '';
 339     my $aln_negative = Bio::SimpleAlign->new();
 340     $aln_negative->add_seq($seq1);
 341     $aln_negative->add_seq($seq2);
 342     my $start_column =
 343       $aln_negative->column_from_residue_number(
 344         $aln_negative->get_seq_by_pos(1)->display_id, 2 );
 345     my $end_column =
 346       $aln_negative->column_from_residue_number(
 347         $aln_negative->get_seq_by_pos(1)->display_id, 5 );
 348     $aln_negative = $aln_negative->slice( $end_column, $start_column );
 349     my $seq_negative = $aln_negative->get_seq_by_pos(1);
 350     is( $seq_negative->start, 2, "bug 2099" );
 351     is( $seq_negative->end,   5, "bug 2099" );
 352
 353     # bug 2793
 354     my $s11 = Bio::LocatableSeq->new( -id => 'testseq1', -seq => 'AAA' );
 355     my $s21 = Bio::LocatableSeq->new( -id => 'testseq2', -seq => 'CCC' );
 356     $a = Bio::SimpleAlign->new();
 357     ok( $a->add_seq( $s11, 1 ), "bug 2793" );
 358     is( $a->get_seq_by_pos(1)->seq, 'AAA', "bug 2793" );
 359     ok( $a->add_seq( $s21, 2 ), "bug 2793" );
 360     is( $a->get_seq_by_pos(2)->seq, 'CCC', "bug 2793" );
 361     throws_ok { $a->add_seq( $s21, 0 ) } qr/must be >= 1/, 'Bad sequence, bad!';
 362 }
 363
 364 # test for Bio::SimpleAlign annotation method and
 365 # Bio::FeatureHolder stuff
 366
 367 $aln = Bio::SimpleAlign->new;
 368 isa_ok($aln,"Bio::AnnotatableI");
 369
 370 for my $seqset ( [qw(one AGAGGAT)], [qw(two AGACGAT)], [qw(three AGAGGTT)] ) {
 371     $aln->add_seq(
 372         Bio::LocatableSeq->new(
 373             -id  => $seqset->[0],
 374             -seq => $seqset->[1]
 375         )
 376     );
 377 }
 378
 379 is $aln->num_sequences, 3, 'added 3 seqs';
 380
 381 $aln->add_SeqFeature(
 382     Bio::SeqFeature::Generic->new(
 383         -start       => 1,
 384         -end         => 1,
 385         -primary_tag => 'charLabel',
 386     )
 387 );
 388 $aln->add_SeqFeature(
 389     Bio::SeqFeature::Generic->new(
 390         -start       => 3,
 391         -end         => 3,
 392         -primary_tag => 'charLabel',
 393
 394     )
 395 );
 396 is( $aln->feature_count, 2, 'first 2 features added' );
 397
 398 my $splitloc = Bio::Location::Split->new;
 399 $splitloc->add_sub_Location(
 400     Bio::Location::Simple->new(
 401         -start => 2,
 402         -end   => 3
 403     )
 404 );
 405
 406 $splitloc->add_sub_Location(
 407     Bio::Location::Simple->new(
 408         -start => 5,
 409         -end   => 6
 410     )
 411 );
 412
 413 $aln->add_SeqFeature(
 414     Bio::SeqFeature::Generic->new(
 415         -location    => $splitloc,
 416         -primary_tag => 'charLabel',
 417     )
 418 );
 419
 420 is( $aln->feature_count, 3, '3rd feature added' );
 421
 422 #do slices and masks as defined by the feature
 423 my $i          = 0;
 424 my @slice_lens = qw(1 1 2 2);
 425 for my $feature ( $aln->get_SeqFeatures ) {
 426     for my $loc ( $feature->location->each_Location ) {
 427         my $masked = $aln->mask_columns( $loc->start, $loc->end, '?');
 428         $masked->verbose(2);
 429         lives_ok {my $fslice = $masked->slice( $loc->start, $loc->end )};
 430
 431         $masked->verbose(-1);
 432         my $fslice = $masked->slice( $loc->start, $loc->end );
 433         is( $fslice->length, $slice_lens[ $i++ ], "slice $i len" );
 434         for my $s ( $fslice->each_seq ) {
 435             like( $s->seq, qr/^\?+$/, 'correct masked seq' );
 436         }
 437     }
 438 }
 439
 440 # test set_displayname_safe & restore_displayname:
 441 $str = Bio::AlignIO->new( -file => test_input_file('pep-266.aln') );
 442 $aln = $str->next_aln();
 443 is $aln->get_seq_by_pos(3)->display_id, 'Smik_Contig1103.1',
 444   'initial display id ok';
 445 my ( $new_aln, $ref ) = $aln->set_displayname_safe();
 446 is $new_aln->get_seq_by_pos(3)->display_id, 'S000000003', 'safe display id ok';
 447 my $restored_aln = $new_aln->restore_displayname($ref);
 448 is $restored_aln->get_seq_by_pos(3)->display_id, 'Smik_Contig1103.1',
 449   'restored display id ok';
 450
 451 # test sort_by_list:
 452 $str = Bio::AlignIO->new( -file => test_input_file('testaln.clustalw') );
 453 my $list_file = test_input_file('testaln.list');
 454 $aln     = $str->next_aln();
 455 $new_aln = $aln->sort_by_list($list_file);
 456 $a       = $new_aln->get_seq_by_pos(1)->display_id;
 457 is $a, 'BAB68554', 'sort by list ok';
 458
 459 # test for Binary/Morphological/Mixed data
 460
 461 # sort_by_start
 462
 463 # test sort_by_list:
 464
 465 my $s1 = Bio::LocatableSeq->new(
 466     -id       => 'AAA',
 467     -seq      => 'aawtat-tn-',
 468     -start    => 12,
 469     -end      => 19,
 470     -alphabet => 'dna'
 471 );
 472 my $s2 = Bio::LocatableSeq->new(
 473     -id       => 'BBB',
 474     -seq      => '-aaaat-tt-',
 475     -start    => 1,
 476     -end      => 7,
 477     -alphabet => 'dna'
 478 );
 479 my $s3 = Bio::LocatableSeq->new(
 480     -id       => 'BBB',
 481     -seq      => '-aaaat-tt-',
 482     -start    => 31,
 483     -end      => 37,
 484     -alphabet => 'dna'
 485 );
 486 $a = Bio::SimpleAlign->new();
 487 $a->add_seq($s1);
 488 $a->add_seq($s2);
 489 $a->add_seq($s3);
 490
 491 @seqs = $a->each_seq;
 492 is( $seqs[0]->start, 12 );
 493 is( $seqs[1]->start, 1 );
 494 is( $seqs[2]->start, 31 );
 495
 496 $a->sort_by_start;
 497 @seqs = $a->each_seq;
 498
 499 is( $seqs[0]->start, 1 );
 500 is( $seqs[1]->start, 12 );
 501 is( $seqs[2]->start, 31 );
 502
 503 my %testdata = (
 504     'allele1' => 'GGATCCATT[C/C]CTACT',
 505     'allele2' => 'GGAT[C/-][C/-]ATT[C/C]CT[A/C]CT',
 506     'allele3' => 'G[G/C]ATCCATT[C/G]CTACT',
 507     'allele4' => 'GGATCCATT[C/G]CTACT',
 508     'allele5' => 'GGATCCATT[C/G]CTAC[T/A]',
 509     'allele6' => 'GGATCCATT[C/G]CTA[C/G][T/A]',
 510     'testseq' => 'GGATCCATT[C/G]CTACT'
 511 );
 512
 513 my $alnin = Bio::AlignIO->new(
 514     -format => 'fasta',
 515     -file   => test_input_file('alleles.fas')
 516 );
 517
 518 $aln = $alnin->next_aln;
 519
 520 my $ct = 0;
 521
 522 # compare all to test seq
 523
 524 for my $ls ( sort keys %testdata ) {
 525     $ct++;
 526     my $str = $aln->bracket_string(
 527         -refseq  => 'testseq',
 528         -allele1 => 'allele1',
 529         -allele2 => $ls,
 530     );
 531     is( $str, $testdata{$ls}, "BIC:$str" );
 532 }
 533
 534 %testdata = (
 535     'allele1' => 'GGATCCATT{C.C}CTACT',
 536     'allele2' => 'GGAT{C.-}{C.-}ATT{C.C}CT{A.C}CT',
 537     'allele3' => 'G{G.C}ATCCATT{C.G}CTACT',
 538     'allele4' => 'GGATCCATT{C.G}CTACT',
 539     'allele5' => 'GGATCCATT{C.G}CTAC{T.A}',
 540     'allele6' => 'GGATCCATT{C.G}CTA{C.G}{T.A}',
 541     'testseq' => 'GGATCCATT{C.G}CTACT'
 542 );
 543
 544 for my $ls ( sort keys %testdata ) {
 545     $ct++;
 546     my $str = $aln->bracket_string(
 547         -refseq     => 'testseq',
 548         -allele1    => 'allele1',
 549         -allele2    => $ls,
 550         -delimiters => '{}',
 551         -separator  => '.'
 552     );
 553     is( $str, $testdata{$ls}, "BIC:$str" );
 554 }
 555
 556 # is _remove_col really working correctly?
 557 my $a = Bio::LocatableSeq->new(
 558     -id    => 'a',
 559     -strand => 1,
 560     -seq   => 'atcgatcgatcgatcg',
 561     -start => 5,
 562     -end   => 20
 563 );
 564 my $b = Bio::LocatableSeq->new(
 565     -id    => 'b',
 566     -strand => 1,
 567     -seq   => '-tcgatc-atcgatcg',
 568     -start => 30,
 569     -end   => 43
 570 );
 571 my $c = Bio::LocatableSeq->new(
 572     -id    => 'c',
 573     -strand => -1,
 574     -seq   => 'atcgatcgatc-atc-',
 575     -start => 50,
 576     -end   => 63
 577 );
 578 my $d = Bio::LocatableSeq->new(
 579     -id    => 'd',
 580     -strand => -1,
 581     -seq   => '--cgatcgatcgat--',
 582     -start => 80,
 583     -end   => 91
 584 );
 585 my $e = Bio::LocatableSeq->new(
 586     -id    => 'e',
 587     -strand => 1,
 588     -seq   => '-t-gatcgatcga-c-',
 589     -start => 100,
 590     -end   => 111
 591 );
 592 $aln = Bio::SimpleAlign->new();
 593 $aln->add_seq($a);
 594 $aln->add_seq($b);
 595 $aln->add_seq($c);
 596
 597 my $gapless = $aln->remove_gaps();
 598 foreach my $seq ( $gapless->each_seq ) {
 599     if ( $seq->id eq 'a' ) {
 600         is $seq->start, 6;
 601         is $seq->end,   19;
 602         is $seq->seq,   'tcgatcatcatc';
 603     }
 604     elsif ( $seq->id eq 'b' ) {
 605         is $seq->start, 30;
 606         is $seq->end,   42;
 607         is $seq->seq,   'tcgatcatcatc';
 608     }
 609     elsif ( $seq->id eq 'c' ) {
 610         is $seq->start, 51;
 611         is $seq->end,   63;
 612         is $seq->seq,   'tcgatcatcatc';
 613     }
 614 }
 615
 616 $aln->add_seq($d);
 617 $aln->add_seq($e);
 618 $gapless = $aln->remove_gaps();
 619 foreach my $seq ( $gapless->each_seq ) {
 620     if ( $seq->id eq 'a' ) {
 621         is $seq->start, 8;
 622         is $seq->end,   17;
 623         is $seq->seq,   'gatcatca';
 624     }
 625     elsif ( $seq->id eq 'b' ) {
 626         is $seq->start, 32;
 627         is $seq->end,   40;
 628         is $seq->seq,   'gatcatca';
 629     }
 630     elsif ( $seq->id eq 'c' ) {
 631         is $seq->start, 53;
 632         is $seq->end,   61;
 633         is $seq->seq,   'gatcatca';
 634     }
 635     elsif ( $seq->id eq 'd' ) {
 636         is $seq->start, 81;
 637         is $seq->end,   90;
 638         is $seq->seq,   'gatcatca';
 639     }
 640     elsif ( $seq->id eq 'e' ) {
 641         is $seq->start, 101;
 642         is $seq->end,   110;
 643         is $seq->seq,   'gatcatca';
 644     }
 645 }
 646
 647 # bug 3077
 648
 649 my $slice = $aln->slice(1,4);
 650
 651 for my $seq ($slice->each_seq) {
 652     if ( $seq->id eq 'a' ) {
 653         is $seq->start, 5;
 654         is $seq->end,   8;
 655         is $seq->strand, 1;
 656         is $seq->seq,   'atcg';
 657     }
 658     elsif ( $seq->id eq 'b' ) {
 659         is $seq->start, 30;
 660         is $seq->end,   32;
 661         is $seq->strand, 1;
 662         is $seq->seq,   '-tcg';
 663     }
 664     elsif ( $seq->id eq 'c' ) {
 665         is $seq->start, 60;
 666         is $seq->end,   63;
 667         is $seq->strand, -1;
 668         is $seq->seq,   'atcg';
 669     }
 670     elsif ( $seq->id eq 'd' ) {
 671         is $seq->start, 90;
 672         is $seq->end,   91;
 673         is $seq->strand, -1;
 674         is $seq->seq,   '--cg';
 675     }
 676     elsif ( $seq->id eq 'e' ) {
 677         is $seq->start, 100;
 678         is $seq->end,   101;
 679         is $seq->strand, 1;
 680         is $seq->seq,   '-t-g';
 681     }
 682 }
 683
 684 my $f = Bio::LocatableSeq->new(
 685     -id    => 'f',
 686     -seq   => 'a-cgatcgatcgat-g',
 687     -start => 30,
 688     -end   => 43
 689 );
 690 $aln = Bio::SimpleAlign->new();
 691 $aln->add_seq($a);
 692 $aln->add_seq($f);
 693
 694 $gapless = $aln->remove_gaps();
 695 foreach my $seq ( $gapless->each_seq ) {
 696     if ( $seq->id eq 'a' ) {
 697         is $seq->start, 5;
 698         is $seq->end,   20;
 699         is $seq->seq,   'acgatcgatcgatg';
 700     }
 701     elsif ( $seq->id eq 'f' ) {
 702         is $seq->start, 30;
 703         is $seq->end,   43;
 704         is $seq->seq,   'acgatcgatcgatg';
 705     }
 706 }
 707
 708 my $g =
 709   Bio::LocatableSeq->new( -id => 'g', -seq => 'atgc', -start => 5, -end => 8 );
 710 my $h = Bio::LocatableSeq->new(
 711     -id    => 'h',
 712     -seq   => '-tcg',
 713     -start => 30,
 714     -end   => 32
 715 );
 716 $aln = Bio::SimpleAlign->new();
 717 $aln->add_seq($g);
 718 $aln->add_seq($h);
 719
 720 # test for new method in API get_seq_by_id
 721 my $retrieved = $aln->get_seq_by_id('g');
 722 is( defined $retrieved, 1 );
 723 my $removed = $aln->remove_columns( [ 1, 3 ] );
 724 foreach my $seq ( $removed->each_seq ) {
 725     if ( $seq->id eq 'g' ) {
 726         is $seq->start, 5;
 727         is $seq->end,   5;
 728         is $seq->seq,   'a';
 729     }
 730     elsif ( $seq->id eq 'h' ) {
 731         is $seq->start, 0;
 732         is $seq->end,   0;
 733         is $seq->seq,   '-';
 734     }
 735 }
 736
 737 # work out mask_columns(), see bug 2842
 738 SKIP: {
 739     test_skip(-tests => 6, -requires_module => 'IO::String');
 740     my $io = Bio::AlignIO->new( -file => test_input_file("testaln.clustalw") );
 741     my $aln = $io->next_aln();
 742     isa_ok( $aln, 'Bio::SimpleAlign' );
 743     my $consensus = <<EOU;
 744 MNEGEHQIKLDELFEKLLRARKIFKNKDVLRHSWEPKDLPHRHEQIEALAQILV
 745 PVLRGETMKIIFCGHHACELGEDRGTKGFVIDELKDVDEDRNGKVDVIEINCEH
 746 MDTHYRVLPNIAKLFDDCTGIGVPMHGGPTDEVTAKLKQVIDMKERFVIIVLDE
 747 IDKLVKKSGDEVLYSLTRINTELKRAKVSVIGISNDLKFKEYLDPRVLSSLSEE
 748 EVVFPPYDANQLRDILTQRAEEAFYPGVLDEGVIPLCAALAAREHGDARKALDL
 749 LRVAGEIAEREGASKVTEKHVWKAQEKIEQDMMEEVIKTLPLQSKVLLYAIVLL
 750 DENGDLPANTGDVYAVYRELCEYIDLEPLTQRRISDLINELDMLGIINAKVVSK
 751 GRYGRTKEIRLMVTSYKIRNVLRYDYSIQPLLTISLKSEQRRLI
 752 EOU
 753     $consensus =~ s/\n//g;
 754
 755     is( $aln->consensus_string, $consensus, 'consensus string looks ok' );
 756
 757
 758     my @cons_got = $aln->consensus_conservation;
 759     # 422 positions, mostly two of six sequences conserved, set as default
 760     my @cons_expect = (100 * 2/6) x 422;
 761     # Exceptionally columns as a mask, manually determined (1-based columns)
 762     $cons_expect[$_-1] = 100 * 1/6 for (5,12,41,70,82,310,390);
 763     $cons_expect[$_-1] = 100 * 3/6 for (27,30,32,36,47,49,61,66,69,71,77,79,
 764         81,91,96,97,105,114,115,117,118,121,122,129,140,146,156,159,160,162,
 765         183,197,217,221,229,242,247,248,261,266,282,287,295,316,323,329,335,337,344,);
 766     $cons_expect[$_-1] = 100 * 4/6 for (84,93,99,100,102,107,108,112,113,119,150,);
 767     $cons_expect[$_-1] = 100 * 5/6 for (81,110);
 768     # Format for string comparison
 769     @cons_expect = map { sprintf "%4.1f", $_ } @cons_expect;
 770     @cons_got = map { sprintf "%4.1f", $_ } @cons_got;
 771     is(length($aln->consensus_string), scalar(@cons_got),"conservation length");
 772     is_deeply(\@cons_got, \@cons_expect, "conservation scores");
 773
 774
 775     is( aln2str( $aln => 'pfam' ), <<EOA, 'looks like correct unmasked alignment (from clustalw)' );
 776 P84139/1-420              MNEGEHQIKLDELFEKLLRARKIFKNKDVLRHSYTPKDLPLRHEQIETLAQILVPVLRGETPSNIFVYG-KTGTGKTVTVK-FVTEELKRISEKYNIPVDVIYINCEIVDTHYRVLANIVNYFKDETGIGVPMVGWPTDEVYAKLKQVIDMKERFVIIVLDEIDKLVKKSGDEVLYSLTRINTELKRAKVSVIGISNDLKFKEYLDPRVLSSLSEEEVVFPPYDANQLRDILTQRAEEAFYPGVLDEGVIPLCAALAAREHGDARKALDLLRVAGEIAEREGASKVTEKHVWKAQEKIEQDMMEEVIKTRPLQSKVLLYAIVLLDENGDLPANTGDVYAVYRELCEYIDLEPLTQRRISDLINELDMLGIINAKVVSKGRYGRTKEIRLNVTSYKIRNVLRYDYSIQPLLTISLKSEQRRLI
 777 P814153/1-420             MNEGMHQIKLDVLFEKLLRARKIFKNKDVLRHSYTPKDLPHRHEQIETLAQILVPVLRGETPSNIFVYG-KTGTGKTVTVK-FVTEELKRISEKYNIPVDVIYINCEIVDTHYRVLANIVNYFKDETGIEVPMVGWPTDEVYAKLKQVIDMKERFVIIVLDEIDKLVKKSGDEVLYSLTRINTELKRAKVSVIGISNDLKFKEYLDPRVLSSLSEEEVVFPPYDANQLRDILTQRAEEAFYPGVLDEGVIPLCAALAAREHGDARKALDLLRVAGEIAEREGASKVTEKHVWKAQEKIEQDMMEEVIKTLPLQSKVLLYAIVLLDENGDLPANTGDVYAVYRELCEYIDLEPLTQRRISDLINELDMLGIINAKVVSKGRYGRTKEIRLMVTSYKIRNVLRYDYSIQPLLTISLKSEQRRLI
 778 P851414/1-60              -------------------------------------------------------------MKIVWCGH-ACFLVEDRGTK-ILIDPYPDVDEDRIGKVDYILQTHEHMD-HYGKTPLIAKLSD----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 779 P841414/1-60              -------------------------------------------------------------MKIVWCGH-ACFLVEDRGTK-ILIDPYPDVDEDRIGKVDYILVTHEHMD-HYGKTPLIAKLSD----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 780 BAB68554/1-141            --------------------MLTEDDKQLIQHVWEKVLEHQEDFGAEALERMFIVYPSTKTYFPHFDLHHDSEQIRHHGKK-VVGALGDAVKHIDNLSATLSELSNLHCY-NLRVDPVNFKLLSHCFQVVLGAHLG--REYTPQVQVAYDKFLAAVSAVLAEKYR-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 781 gb|443893|124775/1-331    -MRFRFGVVVPPAVAGARPELLVVGSRPELG-RWEPRGAVRLRPAGTAAGDGALALQEPGLWLGEVELA-AEEAAQDGAEPGRVDTFWYKFLKREPGGELSWEGNGPHHDRCCTYNENNLVDGVYCLPIG---HWGEATGHTNEMKHTTDFYFNIAGHQAMHYSRILPNIWLGSCPRQVEHVTIKLKHELGITAVMN-FQTEWDIVQNSSGCNRYPEPMTPDTMIKLYREEGLAYIWMP-TPDMSTEGRVQMLPQAVCLLHALLEKGHIVY-----VHCNAGVGRSTAAVCGWLQYVMGWNLRKVQYFLMAKRPAVYIDEEALARAQEDFFQKFGKVRSSVCSL------------------------------------------------------------------------------
 782 EOA
 783
 784     my $newaln = $aln->mask_columns(12,20,'?');
 785     is( aln2str( $newaln, 'pfam' ), <<EOA, 'looks like correct masked alignment (from clustalw)' );
 786 P84139/1-420              MNEGEHQIKLD?????????RKIFKNKDVLRHSYTPKDLPLRHEQIETLAQILVPVLRGETPSNIFVYG-KTGTGKTVTVK-FVTEELKRISEKYNIPVDVIYINCEIVDTHYRVLANIVNYFKDETGIGVPMVGWPTDEVYAKLKQVIDMKERFVIIVLDEIDKLVKKSGDEVLYSLTRINTELKRAKVSVIGISNDLKFKEYLDPRVLSSLSEEEVVFPPYDANQLRDILTQRAEEAFYPGVLDEGVIPLCAALAAREHGDARKALDLLRVAGEIAEREGASKVTEKHVWKAQEKIEQDMMEEVIKTRPLQSKVLLYAIVLLDENGDLPANTGDVYAVYRELCEYIDLEPLTQRRISDLINELDMLGIINAKVVSKGRYGRTKEIRLNVTSYKIRNVLRYDYSIQPLLTISLKSEQRRLI
 787 P814153/1-420             MNEGMHQIKLD?????????RKIFKNKDVLRHSYTPKDLPHRHEQIETLAQILVPVLRGETPSNIFVYG-KTGTGKTVTVK-FVTEELKRISEKYNIPVDVIYINCEIVDTHYRVLANIVNYFKDETGIEVPMVGWPTDEVYAKLKQVIDMKERFVIIVLDEIDKLVKKSGDEVLYSLTRINTELKRAKVSVIGISNDLKFKEYLDPRVLSSLSEEEVVFPPYDANQLRDILTQRAEEAFYPGVLDEGVIPLCAALAAREHGDARKALDLLRVAGEIAEREGASKVTEKHVWKAQEKIEQDMMEEVIKTLPLQSKVLLYAIVLLDENGDLPANTGDVYAVYRELCEYIDLEPLTQRRISDLINELDMLGIINAKVVSKGRYGRTKEIRLMVTSYKIRNVLRYDYSIQPLLTISLKSEQRRLI
 788 P851414/1-60              -------------------------------------------------------------MKIVWCGH-ACFLVEDRGTK-ILIDPYPDVDEDRIGKVDYILQTHEHMD-HYGKTPLIAKLSD----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 789 P841414/1-60              -------------------------------------------------------------MKIVWCGH-ACFLVEDRGTK-ILIDPYPDVDEDRIGKVDYILVTHEHMD-HYGKTPLIAKLSD----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 790 BAB68554/1-141            --------------------MLTEDDKQLIQHVWEKVLEHQEDFGAEALERMFIVYPSTKTYFPHFDLHHDSEQIRHHGKK-VVGALGDAVKHIDNLSATLSELSNLHCY-NLRVDPVNFKLLSHCFQVVLGAHLG--REYTPQVQVAYDKFLAAVSAVLAEKYR-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 791 gb|443893|124775/1-331    -MRFRFGVVVP?????????LLVVGSRPELG-RWEPRGAVRLRPAGTAAGDGALALQEPGLWLGEVELA-AEEAAQDGAEPGRVDTFWYKFLKREPGGELSWEGNGPHHDRCCTYNENNLVDGVYCLPIG---HWGEATGHTNEMKHTTDFYFNIAGHQAMHYSRILPNIWLGSCPRQVEHVTIKLKHELGITAVMN-FQTEWDIVQNSSGCNRYPEPMTPDTMIKLYREEGLAYIWMP-TPDMSTEGRVQMLPQAVCLLHALLEKGHIVY-----VHCNAGVGRSTAAVCGWLQYVMGWNLRKVQYFLMAKRPAVYIDEEALARAQEDFFQKFGKVRSSVCSL------------------------------------------------------------------------------
 792 EOA
 793
 794 ###### test with phylip
 795
 796     my $phylip_str = <<EOF;
 797  3 37
 798 seq1         AAAATGGGGG TGGT------ GGTACCT--- -------
 799 seq2         -----GGCGG TGGTGNNNNG GGTTCCCTNN NNNNNNN
 800 new          AAAATGGNGG TGGTN----N GGTNCCNTNN NNNNNNN
 801
 802 EOF
 803
 804     my $phylip_masked = <<EOF;
 805  3 37
 806 seq1         AAAATGGGGG TGGT------ GGTACCT--- -------
 807 seq2         -----GGCGG TGGT?????? GGTTCCCTNN NNNNNNN
 808 new          AAAATGGNGG TGGT?----? GGTNCCNTNN NNNNNNN
 809
 810 EOF
 811
 812     my $phy_fh = IO::String->new( $phylip_str );
 813
 814     my $in = Bio::AlignIO->new( -fh => $phy_fh, -format => 'phylip' );
 815     unified_diff;
 816
 817     $aln = $in->next_aln();
 818     eq_or_diff( aln2str( $aln, 'phylip' ), $phylip_str );
 819
 820     $newaln = $aln->mask_columns(15,20,'?');
 821     eq_or_diff( aln2str( $newaln,'phylip' ), $phylip_masked, 'align after looks ok' );
 822 }
 823
 824 ######## SUBROUTINES
 825
 826 sub aln2str {
 827     my ( $aln, $fmt ) = @_;
 828     my $out;
 829     my $out_fh = IO::String->new( $out );
 830     my $alignio_out = Bio::AlignIO->new(-fh => $out_fh, -format => $fmt);
 831     $alignio_out->write_aln( $aln );
 832     return $out;
 833 }