t/Seq/LocatableSeq.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 => 119);
  10
  11     use_ok('Bio::LocatableSeq');
  12     use_ok('Bio::AlignIO');
  13 }
  14
  15 my ($str, $aln, $seq, $loc);
  16
  17 # basic tests
  18
  19 ok $seq = Bio::LocatableSeq->new(
  20                  -seq => '--atg---gta--',
  21                  -strand => 1,
  22                  -alphabet => 'dna'
  23                  );
  24 is $seq->alphabet, 'dna';
  25 is $seq->start, 1;
  26 is $seq->end, 6;
  27 is $seq->strand, 1;
  28 is $seq->num_gaps, 1;
  29 is $seq->column_from_residue_number(4), 9;
  30 is $seq->column_from_residue_number(3), 5;
  31
  32 ok $loc = $seq->location_from_column(4);
  33 isa_ok $loc,'Bio::Location::Simple';
  34 is $loc->to_FTstring, 2;
  35
  36 ok $loc = $seq->location_from_column(6);
  37 isa_ok $loc,'Bio::Location::Simple';
  38 is $loc->start, 3;
  39 is $loc->location_type, 'IN-BETWEEN';
  40 is $loc->to_FTstring, '3^4';
  41
  42 is $loc = $seq->location_from_column(2), undef;
  43 TODO: {
  44   local $TODO = "Need to fix columns before start of seq w/ start > 1";
  45   $seq->start(90);
  46   is $loc = $seq->location_from_column(2), undef;
  47 }
  48
  49 $str = Bio::AlignIO->new(-file=> test_input_file('testaln.pfam'));
  50 ok defined($str);
  51 isa_ok $str,'Bio::AlignIO';
  52 $aln = $str->next_aln();
  53 ok $seq = $aln->get_seq_by_pos(1);
  54 is ref($seq), 'Bio::LocatableSeq';
  55
  56 is $seq->get_nse, '1433_LYCES/9-246';
  57 is $seq->id, '1433_LYCES';
  58
  59 # test invalid sequence
  60
  61 throws_ok{ $seq = Bio::LocatableSeq->new( -seq => '//!\\' ) } qr/.+/;
  62
  63 # test revcom and trunc
  64
  65 $seq = Bio::LocatableSeq->new(
  66                  -seq => '--atg---gta--',
  67                  -strand => 1,
  68                  -alphabet => 'dna'
  69                  );
  70
  71 my $seq2 = $seq->trunc(1,9);
  72 is $seq2->seq, '--atg---g';
  73 is $seq2->start, 1;
  74 is $seq2->end, 4;
  75 is $seq2->strand, $seq->strand;
  76
  77 $seq2 = $seq->trunc(3,8);
  78 is $seq2->seq, 'atg---';
  79 is $seq2->start, 1;
  80 is $seq2->end, 3;
  81
  82 is $seq->strand(-1), -1;
  83 is $seq->start, 1;
  84 is $seq->end, 6;
  85 $seq2 = $seq->trunc(3,8);
  86 is $seq2->seq, 'atg---';
  87 is $seq2->start, 4;
  88 is $seq2->end, 6;
  89 $seq2 = $seq->revcom();
  90 is $seq2->seq, '--tac---cat--';
  91 is $seq2->start, $seq->start;
  92 is $seq2->end, $seq->end;
  93 is $seq2->strand, $seq->strand * -1;
  94 is $seq2->column_from_residue_number(4), 9;
  95 is $seq2->column_from_residue_number(3), 5;
  96
  97 # test column-mapping for -1 strand sequence
  98 $seq = Bio::LocatableSeq->new(
  99                  -seq => '--atg---gtaa-',
 100                  -strand => -1,
 101                  -alphabet => 'dna'
 102                  );
 103 is $seq->column_from_residue_number(5),5;
 104 is $seq->column_from_residue_number(4),9;
 105 ok $loc = $seq->location_from_column(4);
 106 isa_ok $loc,'Bio::Location::Simple';
 107 is $loc->to_FTstring, 6;
 108 ok $loc = $seq->location_from_column(6);
 109 isa_ok $loc,'Bio::Location::Simple';
 110 is $loc->start, 4;
 111 is $loc->location_type, 'IN-BETWEEN';
 112 is $loc->to_FTstring, '4^5';
 113
 114
 115 # more tests for trunc() with strand -1
 116
 117
 118 ok $seq = Bio::LocatableSeq->new(
 119                  -seq => '--atg---gta--',
 120                  -strand => -1,
 121                  -alphabet => 'dna'
 122                  );
 123 is $seq->alphabet, 'dna';
 124 is $seq->start, 1;
 125 is $seq->end, 6;
 126 is $seq->strand, -1;
 127 is $seq->num_gaps, 1;
 128 is $seq->column_from_residue_number(4), 5;
 129
 130
 131 ok $seq2 = $seq->trunc(1,9);
 132 is $seq2->seq, '--atg---g';
 133 is $seq2->start, 3;
 134 is $seq2->end, 6;
 135 is $seq2->strand, $seq->strand;
 136
 137 is $seq->location_from_column(3)->start, 6;
 138 is $seq->location_from_column(11)->start, 1;
 139 is $seq->location_from_column(9)->start, 3;
 140
 141
 142
 143 ok $seq2 = $seq->trunc(7,12);
 144 is $seq2->seq, '--gta-';
 145 is $seq2->start, 1;
 146 is $seq2->end, 3;
 147
 148
 149 ok $seq2 = $seq->trunc(2,6);
 150 is $seq2->seq, '-atg-';
 151 is $seq2->start, 4;
 152 is $seq2->end, 6;
 153
 154 ok $seq2 = $seq->trunc(4,7);
 155 is $seq2->seq, 'tg--';
 156 is $seq2->start, 4;
 157 is $seq2->end, 5;
 158
 159 ok $seq = Bio::LocatableSeq->new();
 160 is $seq->seq, undef;
 161 is $seq->start, undef;
 162 is $seq->end, undef;
 163 my $nse;
 164 eval{$nse = $seq->get_nse};
 165 ok($@);
 166 is ($nse, undef);
 167 $seq->force_nse(1);
 168 eval{$nse = $seq->get_nse};
 169 ok(!$@);
 170 is ($nse, '/0-0');
 171
 172 # test mapping
 173
 174 # mapping only supported for 1 => 1, 3 => 1, or 1 => 3 mapping relationships
 175
 176 eval{$seq = Bio::LocatableSeq->new(
 177                  -mapping => [40 => 2],
 178                  );};
 179
 180 ok($@);
 181 like($@, qr/Mapping values other than 1 or 3 are not currently supported/);
 182
 183 eval{$seq = Bio::LocatableSeq->new(
 184                  -mapping => [3 => 3],
 185                  );};
 186
 187 ok($@);
 188
 189 # sequence is translated to protein, retains original DNA coordinates
 190 # mapping is 1 residue for every 3 coordinate positions
 191 $seq = Bio::LocatableSeq->new(
 192                  -seq => 'KKKAIDLVGVDKARENRQAIYLGASAIAEF',
 193                  -strand => -1,
 194                  -mapping => [1 => 3],
 195                  -start => 1,
 196                  -end => 90,
 197                  -alphabet => 'dna'
 198                  );
 199
 200 is $seq->seq, 'KKKAIDLVGVDKARENRQAIYLGASAIAEF';
 201 is $seq->start, 1;
 202 is $seq->end, 90;
 203
 204 # sequence is reverse-translated to DNA, retains original protein coordinates
 205 # mapping is 3 residues for every 1 coordinate positions
 206 $seq = Bio::LocatableSeq->new(
 207                  -seq => 'aaraaraargcnathgayytngtnggngtngayaargcnmgngaraaymgncargcnathtayytnggngcnwsngcnathgcngartty',
 208                  -strand => -1,
 209                  -mapping => [3 => 1],
 210                  -start => 1,
 211                  -end => 30,
 212                  -alphabet => 'protein'
 213                  );
 214
 215 is $seq->seq, 'aaraaraargcnathgayytngtnggngtngayaargcnmgngaraaymgncargcnathtayytnggngcnwsngcnathgcngartty';
 216 is $seq->start, 1;
 217 is $seq->end, 30;
 218
 219 # frameshifts (FASTA-like)
 220 # support for this is preliminary
 221 # this is a real example from a TFASTY report
 222
 223 $seq = Bio::LocatableSeq->new(
 224                  -seq => 'MGSSSTDRELLSAADVGRTVSRIAHQIIEKTALDDPAERTRVVLLGIPTRGVILATRLAAKIKEFAGEDVPHGALDITLYRDDLNFKPPRPLEATSIPAF\GGVDDAIVILVDDVLYSGRSVRSALDALRDIGRPRIVQLAVLVDRGHRELPI--/DYVGKNVPTSRSESVHVLLSEHDDRDGVVISK',
 225                  -strand => 1,
 226                  -mapping => [1 => 3],
 227                  -start => 1,
 228                  -end => 552,
 229                  -frameshifts => { # position, frameshift
 230                     298 => -1,
 231                     455 => 1
 232                     },
 233                  -alphabet => 'dna'
 234                  );
 235
 236 is $seq->seq, 'MGSSSTDRELLSAADVGRTVSRIAHQIIEKTALDDPAERTRVVLLGIPTRGVILATRLAAKIKEFAGEDVPHGALDITLYRDDLNFKPPRPLEATSIPAF\GGVDDAIVILVDDVLYSGRSVRSALDALRDIGRPRIVQLAVLVDRGHRELPI--/DYVGKNVPTSRSESVHVLLSEHDDRDGVVISK';
 237 is $seq->start, 1;
 238 is $seq->end, 552;
 239 $seq->verbose(2);
 240 eval { $seq->end(554);};
 241 ok $@;
 242 like $@, qr/Overriding value \[554\] with value 552/;
 243
 244 lives_ok { $seq = Bio::LocatableSeq->new(
 245                  -seq => 'LSYC*',
 246                  -strand => 0,
 247                  -start => 1,
 248                  -end => 5,
 249                  -verbose => 2
 250                  );} '* is counted in length';
 251
 252 throws_ok { $seq = Bio::LocatableSeq->new(
 253                  -seq => 'LSYC*',
 254                  -strand => 0,
 255                  -start => 1,
 256                  -end => 6,
 257                  -verbose => 2
 258                  );} qr/Overriding value \[6\] with value 5/, '* is counted in length, but end is wrong';
 259
 260 # setting symbols (class variables) - demonstrate scoping issues when using
 261 # globals with and w/o localization.  To be fixed in a future BioPerl version
 262
 263 # see bug 2715
 264 my $temp;
 265
 266 {
 267     $temp = $Bio::LocatableSeq::GAP_SYMBOLS;
 268     $Bio::LocatableSeq::GAP_SYMBOLS = '-\?';
 269     $seq = Bio::LocatableSeq->new(
 270                      -seq => '??atg-?-gta-?',
 271                      -strand => 1,
 272                      -start => 10,
 273                      -end => 15,
 274                      -alphabet => 'dna',
 275                      );
 276     is $Bio::LocatableSeq::GAP_SYMBOLS, '-\?';
 277     is $seq->start, 10;
 278     is $seq->end, 15;
 279 }
 280
 281 is $Bio::LocatableSeq::GAP_SYMBOLS, '-\?';
 282 is $seq->end(15), 15;
 283 $Bio::LocatableSeq::GAP_SYMBOLS = $temp;
 284 is $Bio::LocatableSeq::GAP_SYMBOLS, '\-\.=~';
 285
 286 {
 287     local $Bio::LocatableSeq::GAP_SYMBOLS = '-\?';
 288     $seq = Bio::LocatableSeq->new(
 289                      -seq => '??atg-?-gta-?',
 290                      -strand => 1,
 291                      -start => 10,
 292                      -end => 15,
 293                      -alphabet => 'dna',
 294                      );
 295     is $Bio::LocatableSeq::GAP_SYMBOLS, '-\?';
 296     is $seq->start, 10;
 297     is $seq->end, 15;
 298 }
 299
 300 is $seq->end, 15;
 301
 302 # note, recalling the end() method uses old $GAP_SYMBOLS, which
 303 # no longer are set (this argues for locally set symbols)
 304 TODO: {
 305     local $TODO = 'Bio::LocatableSeq global variables have scoping issues';
 306     is $Bio::LocatableSeq::GAP_SYMBOLS, '-\?';
 307     # this should be 15
 308     isnt $seq->end(19), 19;
 309 }
 310