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