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