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