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