t/Seq/LocatableSeq.t

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