t/SeqIO/phd.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 => 21);
  11
  12         use_ok('Bio::SeqIO::phd');
  13 }
  14
  15 my $DEBUG = test_debug();
  16
  17 print("Checking to see if Bio::Seq::Quality objects can be created from a file...\n") if ($DEBUG);
  18 my $in_phd  = Bio::SeqIO->new('-file' => test_input_file('test.phd'),
  19                               '-format'  => 'phd',
  20                               '-verbose' => $DEBUG);
  21 isa_ok($in_phd,'Bio::SeqIO::phd');
  22
  23
  24
  25 my $phd = $in_phd->next_seq();
  26 is($phd->quality_levels,'99',"Did you get the 'QUALITY_LEVELS' comment?");
  27 isa_ok($phd,"Bio::Seq::Quality");
  28
  29
  30 if( $DEBUG ) {
  31     my $position = 6;
  32     print("I saw these in phredfile.phd:\n\n");
  33     print $_->tagname,": ",$_->display_text || 0," \n"
  34         for ($phd->annotation->get_Annotations('header'));
  35
  36     print("What is the base at position $position (using subseq)?\n");
  37     print($phd->subseq($position,$position)."\n");
  38     print("What is the base at position $position (using baseat)?\n");
  39     print($phd->baseat($position)."\n");
  40     print("What is the quality at $position? (using subqual)\n");
  41
  42 my @qualsretr = @{$phd->subqual($position,$position)};
  43     print($qualsretr[0]."\n");
  44     print("What is the quality at $position? (using qualat)\n");
  45     print($phd->qualat($position)."\n");
  46     print("What is the trace at $position? (using trace_index_at)\n");
  47     print($phd->trace_index_at($position)."\n");
  48     print("What is the trace at $position? (using subtrace)\n");
  49     my @tracesretr = @{$phd->subtrace($position,$position)};
  50     print($tracesretr[0]."\n");
  51 }
  52
  53 print("OK. Now testing write_phd...\n") if($DEBUG);
  54
  55 my $outfile = test_output_file();
  56 my $out_phd = Bio::SeqIO->new(-file => ">$outfile",
  57                               '-format' => 'phd');
  58 isa_ok($out_phd,"Bio::SeqIO::phd");
  59
  60 $out_phd->write_seq($phd);
  61
  62 ok( -s $outfile);
  63
  64 # Bug 2120
  65
  66 my @qual = q(9 9 12 12 8 8 9 8 8 8 9);
  67 my @trace = q(113 121 130 145 153 169 177 203 210 218 234);
  68
  69 $in_phd  = Bio::SeqIO->new('-file' => test_input_file('bug2120.phd'),
  70                               '-format'  => 'phd',
  71                               '-verbose' => $DEBUG);
  72
  73 my $seq = $in_phd->next_seq;
  74 is($seq->subseq(10,20),'gggggccttat','$seq->subseq()');
  75 my @seq_qual =$seq->subqual_text(10,20);
  76 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
  77 my @seq_trace = $seq->subtrace_text(10,20);
  78 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
  79
  80 if($DEBUG) {
  81     print "\nDefault header ... \n\n";
  82     use Bio::Seq::Quality;
  83     my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
  84     $out_phd->_fh(\*STDOUT);
  85     $out_phd->write_header($seq);
  86     print "Complete output\n\n";
  87     $out_phd->write_seq($seq);
  88 }
  89
  90 print("Testing the header manipulation\n") if($DEBUG);
  91 is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
  92 $phd->chromat_file('ML4924R.esd');
  93 is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
  94 $phd->touch();
  95
  96 # Commented out 1/17/09.
  97 # This isn't exactly a stable regression test as the comparison tests
  98 # localtime() called from two different timepoints. They can differ if the calls
  99 # occurred before and after a change in seconds, for example.
 100
 101 #my $localtime = localtime();
 102 #is($phd->time, "$localtime", $phd->time.':'.$localtime);
 103
 104 if ($DEBUG){
 105     print "Testing the sequence ...\n";
 106     print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
 107     my $revcom = $phd->revcom;
 108     print ">revcom\n",$revcom->seq,"\n";
 109     print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
 110     print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
 111     my $trunc = $phd->trunc(10,20);
 112     print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
 113 }
 114
 115 # Multiple seqs in one file
 116
 117 $in_phd  = Bio::SeqIO->new('-file' => test_input_file('multi.phd'),
 118                               '-format'  => 'phd',
 119                               '-verbose' => $DEBUG);
 120
 121 @qual = qq(9 9 15 17 17 22 22 25 25 22 22);
 122 @trace = qq(98 105 119 128 143 148 162 173 185 197 202);
 123
 124 $seq = $in_phd->next_seq;
 125 is($seq->id, 'ML4924F');
 126 is($seq->subseq(10,20),'tctcgagggta','$seq->subseq()');
 127 @seq_qual =$seq->subqual_text(10,20);
 128 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
 129 @seq_trace = $seq->subtrace_text(10,20);
 130 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
 131
 132 @qual = qq(11 9 6 6 9 19 20 32 34 34 39);
 133 @trace = qq(98 104 122 128 140 147 159 167 178 190 200);
 134
 135 $seq = $in_phd->next_seq;
 136 is($seq->id, 'ML4924R');
 137 is($seq->subseq(10,20),'gcctgcaggta','$seq->subseq()');
 138 @seq_qual =$seq->subqual_text(10,20);
 139 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
 140 @seq_trace = $seq->subtrace_text(10,20);
 141 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
 142
 143 #if($DEBUG) {
 144 #    print "\nDefault header ... \n\n";
 145 #    use Bio::Seq::Quality;
 146 #    my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
 147 #    $out_phd->_fh(\*STDOUT);
 148 #    $out_phd->write_header($seq);
 149 #    print "Complete output\n\n";
 150 #    $out_phd->write_seq($seq);
 151 #}
 152
 153 ##print("Testing the header manipulation\n") if($DEBUG);
 154 #is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
 155 #$phd->chromat_file('ML4924R.esd');
 156 #is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
 157 #$phd->touch();
 158 #my $localtime = localtime();
 159 #is($phd->time, "$localtime");
 160 #if ($DEBUG){
 161 #    print "Testing the sequence ...\n";
 162 #    print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
 163 #    my $revcom = $phd->revcom;
 164 #    print ">revcom\n",$revcom->seq,"\n";
 165 #    print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
 166 #    print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
 167 #    my $trunc = $phd->trunc(10,20);
 168 #    print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
 169 #}
 170 #
 171
 172 # Whole-read tags in the file
 173 $in_phd  = Bio::SeqIO->new('-file' => test_input_file('multiseq_tags.phd'),
 174                               '-format'  => 'phd',
 175                               '-verbose' => $DEBUG);
 176 isa_ok($in_phd,'Bio::SeqIO::phd');
 177 my @seqs = ();
 178 while (my $seq = $in_phd->next_seq){
 179   push @seqs, $seq;
 180 }
 181 is( scalar @seqs, 2 );