t/SeqIO/scf.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 => 80);
  10
  11     use_ok('Bio::SeqIO::scf');
  12     use_ok('Bio::Seq::SequenceTrace');
  13 }
  14
  15 my $verbose = test_debug();
  16
  17 ok my $in_scf = Bio::SeqIO->new(-file => test_input_file('chad100.scf'),
  18                                 -format => 'scf',
  19                                 -verbose => $verbose);
  20
  21 my $swq = $in_scf->next_seq();
  22
  23 isa_ok($swq,"Bio::Seq::SequenceTrace");
  24
  25 cmp_ok (length($swq->seq()), '>', 10);
  26 my $qualities = join(' ',@{$swq->qual()});
  27
  28 cmp_ok (length($qualities), '>', 10);
  29 my $id = $swq->id();
  30 is ($swq->id(), "ML4942R");
  31
  32 my $a_channel = $swq->trace("a");
  33 cmp_ok (scalar(@$a_channel), '>', 10);
  34 my $c_channel = $swq->trace("c");
  35 cmp_ok (scalar(@$c_channel), '>', 10);
  36 my $g_channel = $swq->trace("g");
  37 cmp_ok (scalar(@$g_channel), '>', 10);
  38 my $t_channel = $swq->trace("t");
  39 cmp_ok (scalar(@$t_channel), '>', 10);
  40
  41 my $ref = $swq->peak_indices();
  42 my @indices = @$ref;
  43 my $indexcount = 761;
  44 is (scalar(@indices), $indexcount);
  45
  46 #use Data::Dumper;
  47 #----------------------------------------
  48 isa_ok $swq->seq_obj, 'Bio::Seq::Quality';
  49 isa_ok $swq->qual_obj, 'Bio::Seq::Quality';
  50 is $swq->alphabet, 'dna', 'alphabet';
  51
  52 is $swq->display_id, 'ML4942R', 'display_id';
  53 like $swq->primary_id, qr/HASH/, 'primary_id is the stringified memory position';
  54 is $swq->primary_id('ABC'), 'ABC', 'set primary_id';
  55
  56 is $swq->accession_number, 'unknown', 'accession_number';
  57 is $swq->desc, undef, 'desc';
  58 is $swq->desc('test'), 'test', 'desc';
  59 is $swq->id, 'ML4942R', 'id';
  60 is $swq->id('test'), 'test', 'id';
  61 is length($swq->seq), $indexcount, 'seq';
  62
  63
  64
  65 my $len = 7;
  66 my $start = $swq->length-$len+1;
  67 my $end = $swq->length;
  68
  69 is $swq->subseq($start,$end), 'cctcaag', 'subseq';
  70 is $swq->baseat($start), 'c', 'baseat';
  71 is $swq->qualat($start), '18', 'qualat';
  72
  73 is $swq->trace_value_at('a',$start), '482', 'trace_value_at';
  74
  75 TODO: {
  76     local $TODO = 'documentation and code for accuracies() do not match' if 1;
  77     is $swq->accuracies('a',$start), '482', 'accuracies';
  78 }
  79 my $qualstring = join(' ',@{$swq->subqual($start,$end)});
  80 is ($qualstring, '18 18 21 15 8 8 8');
  81
  82 my $refs = $swq->sub_peak_index($start,$end);
  83 is @$refs, $len, 'sub_peak_index';
  84 is $swq->peak_index_at($start), 8819, 'peak_index_at';
  85
  86 my $indices_at_end = join(' ',@{$swq->sub_peak_index($start,$end)});
  87 is($indices_at_end, '8819 8831 8843 8853 8862 8873 8891');
  88
  89 my $swq_end = $swq->trace_length();
  90 my $swq_start = $swq_end - $len +1;
  91 my $subtrace_a = join(' ',@{$swq->sub_trace('a',$swq_start,$swq_end)});
  92 is $subtrace_a, '13 3 0 0 75 274 431';
  93
  94 my $swq2 = $swq->sub_trace_object(1,5);
  95 #$traces2->verbose(-1);
  96
  97 isa_ok($swq2, 'Bio::Seq::SequenceTrace');
  98
  99 $swq2->_synthesize_traces(), 1; # this should not be a private method! Heikki
 100 $swq2->set_accuracies(), 1;
 101
 102 is $swq->accuracy_at('a',1), '755', 'accuracy_at';
 103
 104 #----------------------------------------
 105
 106
 107 warn("Now checking version3...\n") if $verbose;
 108 my $in_scf_v3 = Bio::SeqIO->new(-file => test_input_file('version3.scf'),
 109                                 -format => 'scf',
 110                                 -verbose => $verbose);
 111
 112 my $v3 = $in_scf_v3->next_seq();
 113 isa_ok($v3, 'Bio::Seq::SequenceTrace');
 114 my $ind = $v3->peak_indices();
 115 my @ff = @$ind;
 116
 117 @indices = @{$v3->peak_indices()};
 118 is (scalar(@indices), 1106);
 119
 120 my %header = %{$in_scf_v3->get_header()};
 121 is $header{bases}, 1106;
 122 is $header{samples},  14107;
 123
 124 # is the Bio::Seq::SequenceTrace AnnotatableI?
 125 my $ac = $v3->annotation();
 126
 127 isa_ok($ac,"Bio::Annotation::Collection");
 128
 129 my @name_comments = grep {$_->tagname() eq 'NAME'}
 130 $ac->get_Annotations('comment');
 131
 132 is $name_comments[0]->as_text(), 'Comment: IIABP1D4373';
 133
 134 # also get comments this way...
 135 $ac = $in_scf_v3->get_comments();
 136
 137 isa_ok($ac,"Bio::Annotation::Collection");
 138
 139 @name_comments = grep {$_->tagname() eq 'NAME'}
 140 $ac->get_Annotations('comment');
 141
 142 is $name_comments[0]->as_text(), 'Comment: IIABP1D4373';
 143
 144 my @conv_comments = grep {$_->tagname() eq 'CONV'}
 145 $ac->get_Annotations('comment');
 146
 147 is $conv_comments[0]->as_text(), 'Comment: phred version=0.990722.h';
 148
 149 # is the SequenceTrace object annotated?
 150 my $st_ac = $swq->annotation();
 151
 152 isa_ok ($st_ac, "Bio::Annotation::Collection");
 153
 154 my @ann =   $st_ac->get_Annotations();
 155
 156 is $ann[0]->tagname, 'SIGN';
 157 is $ann[2]->text, 'SRC3700';
 158 is $ann[5]->tagname, 'LANE';
 159 is $ann[5]->text, 89;
 160 is $ann[6]->text, 'phred version=0.980904.e';
 161 is $ann[8]->text, 'ABI 373A or 377';
 162
 163 my $outfile = test_output_file();
 164 my $out_scf = Bio::SeqIO->new(-file => ">$outfile",
 165                               -format => 'scf',
 166                               -verbose => $verbose);
 167
 168 # Bug 2196 - commentless scf
 169
 170 my $in = Bio::SeqIO->new(-file => test_input_file('13-pilE-F.scf'),
 171                          -format => 'scf',
 172                          -verbose => $verbose);
 173
 174 my $seq = $in->next_seq;
 175
 176 ok ($seq);
 177
 178 isa_ok($seq, 'Bio::Seq::SequenceTrace');
 179
 180 $ac = $seq->annotation;
 181
 182 isa_ok($ac, 'Bio::Annotation::Collection');
 183
 184 @name_comments = grep {$_->tagname() eq 'NAME'}
 185 $ac->get_Annotations('comment');
 186
 187 is $name_comments[0], undef;
 188
 189 @conv_comments = grep {$_->tagname() eq 'CONV'}
 190 $ac->get_Annotations('comment');
 191
 192 is $conv_comments[0], undef;
 193
 194 # the new way
 195
 196 warn("Now testing the _writing_ of scfs\n") if $verbose;
 197
 198 $out_scf->write_seq(-target     =>      $v3,
 199                     -MACH               =>      'CSM sequence-o-matic 5000',
 200                     -TPSW               =>      'trace processing software',
 201                     -BCSW               =>      'basecalling software',
 202                     -DATF               =>      'AM_Version=2.00',
 203                     -DATN               =>      'a22c.alf',
 204                     -CONV               =>      'Bioperl-scf.pm');
 205
 206 ok( -s $outfile && ! -z "$outfile" );
 207
 208 # TODO? tests below don't do much
 209
 210 $out_scf = Bio::SeqIO->new(-verbose => 1,
 211                            -file => ">$outfile",
 212                            -format => 'scf');
 213
 214 $swq = Bio::Seq::Quality->new(-seq =>'ATCGATCGAA',
 215                               -qual =>"10 20 30 40 50 20 10 30 40 50",
 216                               -alphabet =>'dna');
 217
 218 my $trace = Bio::Seq::SequenceTrace->new(-swq => $swq);
 219
 220 $out_scf->write_seq(    -target         =>      $trace,
 221                         -MACH           =>      'CSM sequence-o-matic 5000',
 222                         -TPSW           =>      'trace processing software',
 223                         -BCSW           =>      'basecalling software',
 224                         -DATF           =>      'AM_Version=2.00',
 225                         -DATN           =>      'a22c.alf',
 226                         -CONV           =>      'Bioperl-scf.pm' );
 227
 228 warn("Trying to write an scf with a subset of a real scf...\n") if $verbose;
 229 $out_scf = Bio::SeqIO->new(-verbose => 1,
 230                            -file => ">$outfile",
 231                            -format => 'scf');
 232
 233 $in_scf_v3 = Bio::SeqIO->new(-file => test_input_file('version3.scf'),
 234                              -format => 'scf',
 235                              -verbose => $verbose);
 236 $v3 = $in_scf_v3->next_seq();
 237
 238 my $sub_v3 = $v3->sub_trace_object(5,50);
 239
 240 #warn("The subtrace object is this:\n") if $DEBUG;
 241
 242 $out_scf->write_seq(-target => $sub_v3 );
 243
 244 my $in_scf_v2 = Bio::SeqIO->new(-file => test_input_file('version2.scf'),
 245                                 -format => 'scf',
 246                                 -verbose => $verbose);
 247 $v3 = $in_scf_v2->next_seq();
 248 ok($v3);
 249
 250 $out_scf = Bio::SeqIO->new(-file   => ">$outfile",
 251                            -format => "scf");
 252 $out_scf->write_seq( -target  => $v3,
 253                      -version => 2 );
 254
 255 # simple round trip tests (bug 2881)
 256
 257 my %file_map = (
 258         # filename         # write_seq args
 259         'chad100.scf'   => 1,
 260         '13-pilE-F.scf' => 1,
 261         'version2.scf'  => 1,
 262         'version3.scf'  => 1
 263         );
 264
 265 for my $f (sort keys %file_map) {
 266         my $outfile = test_output_file();
 267         my $in = Bio::SeqIO->new(-file => test_input_file($f),
 268                                                          -format => 'scf');
 269         my $out = Bio::SeqIO->new(-file => ">$outfile",
 270                                                          -format => 'scf');
 271
 272         my $seq1 = $in->next_seq();
 273         isa_ok($seq1, 'Bio::Seq::SequenceTrace');
 274
 275         ok($out->write_seq(-target => $seq1));
 276
 277         my $in2 = Bio::SeqIO->new(-file => "<$outfile",
 278                                                           -format => 'scf');
 279         my $seq2 = $in2->next_seq();
 280         isa_ok($seq2, 'Bio::Seq::SequenceTrace');
 281         if ($seq1->display_id) {
 282                 TODO: {
 283                         local $TODO = "display_id doesn't round trip yet";
 284                         is($seq1->display_id, $seq2->display_id, 'display_id matches');
 285                 }
 286         }
 287         is_deeply($seq1->qual, $seq2->qual, 'qual scores match');
 288 }
 289
 290 # synthesizing traces roundtrip (bug 2881):
 291
 292 my @sequences=('A',
 293                'ATGGAGCTCATCAAAGAATCGACTCATATATCCATCCCTGAACGGCTGACTCACATTAATGGTTGA');
 294 foreach my $sequence (@sequences) {
 295     my $qualstr=join ' ', map { 65 } (1 .. length($sequence));
 296     my $seq_qual=Bio::Seq::Quality->new(-seq=>$sequence, -qual=>$qualstr);
 297
 298     my $outfile=test_output_file();
 299     my $out=Bio::SeqIO->new(-file=>">$outfile", -format=>'scf');
 300     $out->write_seq(-target=>$seq_qual);
 301
 302     my $in=Bio::SeqIO->new(-file=>$outfile, -format=>'scf');
 303     my $in_seq=$in->next_seq();
 304
 305     is_deeply($seq_qual, $in_seq->{swq}, 'Bio::Sequence::Quality matches');
 306 }
 307