t/scf.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib 't/lib';
   8     use BioperlTest;
   9
  10     test_begin(-tests => 35);
  11
  12         use_ok('Bio::SeqIO');
  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 is (scalar(@indices), 761);
  45
  46 warn("Now checking version3...\n") if $verbose;
  47 my $in_scf_v3 = Bio::SeqIO->new(-file => test_input_file('version3.scf'),
  48                                                                                   -format => 'scf',
  49                                                                                   -verbose => $verbose);
  50
  51 my $v3 = $in_scf_v3->next_seq();
  52 isa_ok($v3, 'Bio::Seq::SequenceTrace');
  53 my $ind = $v3->peak_indices();
  54 my @ff = @$ind;
  55
  56 @indices = @{$v3->peak_indices()};
  57 is (scalar(@indices), 1106);
  58
  59 my %header = %{$in_scf_v3->get_header()};
  60 is $header{bases}, 1106;
  61 is $header{samples},  14107;
  62
  63 # is the Bio::Seq::SequenceTrace AnnotatableI?
  64 my $ac = $v3->annotation();
  65
  66 isa_ok($ac,"Bio::Annotation::Collection");
  67
  68 my @name_comments = grep {$_->tagname() eq 'NAME'}
  69   $ac->get_Annotations('comment');
  70
  71 is $name_comments[0]->as_text(), 'Comment: IIABP1D4373';
  72
  73 # also get comments this way...
  74 $ac = $in_scf_v3->get_comments();
  75
  76 isa_ok($ac,"Bio::Annotation::Collection");
  77
  78 @name_comments = grep {$_->tagname() eq 'NAME'}
  79   $ac->get_Annotations('comment');
  80
  81 is $name_comments[0]->as_text(), 'Comment: IIABP1D4373';
  82
  83 my @conv_comments = grep {$_->tagname() eq 'CONV'}
  84   $ac->get_Annotations('comment');
  85
  86 is $conv_comments[0]->as_text(), 'Comment: phred version=0.990722.h';
  87
  88 # is the SequenceTrace object annotated?
  89 my $st_ac = $swq->annotation();
  90
  91 isa_ok ($st_ac, "Bio::Annotation::Collection");
  92
  93 my @ann =   $st_ac->get_Annotations();
  94
  95 is $ann[0]->tagname, 'SIGN';
  96 is $ann[2]->text, 'SRC3700';
  97 is $ann[5]->tagname, 'LANE';
  98 is $ann[5]->text, 89;
  99 is $ann[6]->text, 'phred version=0.980904.e';
 100 is $ann[8]->text, 'ABI 373A or 377';
 101
 102 my $outfile = test_output_file();
 103 my $out_scf = Bio::SeqIO->new(-file => ">$outfile",
 104                                                                                 -format => 'scf',
 105                                                                                 -verbose => $verbose);
 106
 107 # Bug 2196 - commentless scf
 108
 109 my $in = Bio::SeqIO->new(-file => test_input_file('13-pilE-F.scf'),
 110                                                           -format => 'scf',
 111                                                           -verbose => $verbose);
 112
 113 my $seq = $in->next_seq;
 114
 115 ok ($seq);
 116
 117 isa_ok($seq, 'Bio::Seq::SequenceTrace');
 118
 119 $ac = $seq->annotation;
 120
 121 isa_ok($ac, 'Bio::Annotation::Collection');
 122
 123 @name_comments = grep {$_->tagname() eq 'NAME'}
 124   $ac->get_Annotations('comment');
 125
 126 is $name_comments[0], undef;
 127
 128 @conv_comments = grep {$_->tagname() eq 'CONV'}
 129   $ac->get_Annotations('comment');
 130
 131 is $conv_comments[0], undef;
 132
 133 # the new way
 134
 135 warn("Now testing the _writing_ of scfs\n") if $verbose;
 136
 137 $out_scf->write_seq(-target     =>      $v3,
 138                                                   -MACH         =>      'CSM sequence-o-matic 5000',
 139                                                   -TPSW         =>      'trace processing software',
 140                                                   -BCSW         =>      'basecalling software',
 141                                                   -DATF         =>      'AM_Version=2.00',
 142                                                   -DATN         =>      'a22c.alf',
 143                                                   -CONV         =>      'Bioperl-scf.pm');
 144
 145 ok( -s $outfile && ! -z "$outfile" );
 146
 147 # TODO? tests below don't do much
 148
 149 $out_scf = Bio::SeqIO->new(-verbose => 1,
 150                                                         -file => ">$outfile",
 151                                                         -format => 'scf');
 152
 153 $swq = Bio::Seq::Quality->new(-seq =>'ATCGATCGAA',
 154                                                                                 -qual =>"10 20 30 40 50 20 10 30 40 50",
 155                                                                                 -alphabet =>'dna');
 156
 157 my $trace = Bio::Seq::SequenceTrace->new(-swq => $swq);
 158
 159 $out_scf->write_seq(    -target =>      $trace,
 160                                                         -MACH           =>      'CSM sequence-o-matic 5000',
 161                                                         -TPSW           =>      'trace processing software',
 162                                                         -BCSW           =>      'basecalling software',
 163                                                         -DATF           =>      'AM_Version=2.00',
 164                                                         -DATN           =>      'a22c.alf',
 165                                                         -CONV           =>      'Bioperl-scf.pm' );
 166
 167 warn("Trying to write an scf with a subset of a real scf...\n") if $verbose;
 168 $out_scf = Bio::SeqIO->new(-verbose => 1,
 169                                                                         -file => ">$outfile",
 170                                                                         -format => 'scf');
 171
 172 $in_scf_v3 = Bio::SeqIO->new(-file => test_input_file('version3.scf'),
 173                                                                           -format => 'scf',
 174                                                                           -verbose => $verbose);
 175 $v3 = $in_scf_v3->next_seq();
 176
 177 my $sub_v3 = $v3->sub_trace_object(5,50);
 178
 179 #warn("The subtrace object is this:\n") if $DEBUG;
 180
 181 $out_scf->write_seq(-target => $sub_v3 );
 182
 183 my $in_scf_v2 = Bio::SeqIO->new(-file => test_input_file('version2.scf'),
 184                                                                                   -format => 'scf',
 185                                                                                   -verbose => $verbose);
 186 $v3 = $in_scf_v2->next_seq();
 187 ok($v3);
 188
 189 $out_scf = Bio::SeqIO->new(-file   => ">$outfile",
 190                            -format => "scf");
 191 $out_scf->write_seq( -target  => $v3,
 192                      -version => 2 );
 193
 194 # now some version 2 things...