ABI
view release on metacpan or search on metacpan
$abi = ABI->new("filename"); # same thing
=cut
sub new {
my $class = shift;
my $self = {};
bless $self, ref($class) || $class;
$self->_init(@_);
#print "****", $self->{_mac_header}, "\n";
return $self;
}
sub _init {
my ( $self, @args ) = @_;
my ($file) = $self->_rearrange( ["FILE"], @args );
if ( !defined($file) ) {
croak "Can't open the input file\n";
} else {
$self->set_file_handle($file);
$self->{DATA11} = undef;
$self->{DATA12} = undef;
$self->{PBAS1} = undef;
$self->{PBAS2} = undef;
$self->{FWO} = undef;
$self->{_mac_header} = 0;
$self->{_maximum_trace} = 0;
if ( $self->_is_abi() ) {
#print "ABI FILE\n";
$self->_set_index();
$self->_set_base_calls();
$self->_set_corrected_base_calls();
$self->_set_seq();
$self->_set_corrected_seq();
$self->_set_traces();
$self->_set_max_trace();
$self->_set_sample_name();
close( $self->{_fh} );
}
# What we intend to do is loop through the @{$order} variable,
# and for each value, we use that as a key into our associative
# array, pushing the value at that key onto our return array.
my ($key);
foreach $key ( @{$order} ) {
my ($value) = $param{$key};
delete $param{$key};
push( @return_array, $value );
}
# print "\n_rearrange() after processing:\n";
# my $i; for ($i=0;$i<@return_array;$i++) { printf "%20s => %s\n", ${$order}[$i], $return_array[$i]; } <STDIN>;
return (@return_array);
}
sub _is_abi {
my $self = shift;
my $fh = $self->{"_fh"};
my $buf;
seek( $fh, 0, 0 );
read( $fh, $buf, 3 );
$self->{_mac_header} = 128;
}
sub _set_index {
my $self = shift;
my $data_counter = 0;
my $pbas_counter = 0;
my $ploc_counter = 0;
my ( $num_records, $buf );
#print $self->{_fh}, "\n";
#print $self->{_mac_header}, "\n";
seek( $self->{_fh}, $self->{_abs_index} + $self->{_mac_header}, 0 );
read( $self->{_fh}, $buf, 4 );
$self->{_index} = unpack( "N", $buf );
#print $self->{_index};
seek( $self->{_fh}, $self->{_abs_index} - 8 + $self->{_mac_header}, 0 );
read( $self->{_fh}, $buf, 4 );
$num_records = unpack( "N", $buf );
for ( my $i = 0 ; $i <= $num_records - 1 ; $i++ ) {
seek( $self->{_fh}, $self->{_index} + ( $i * 28 ), 0 );
read( $self->{_fh}, $buf, 4 );
if ( $buf eq "FWO_" ) {
$self->{FWO} = $self->{_index} + ( $i * 28 ) + 20;
}
if ( $buf eq "DATA" ) {
sub _set_base_calls {
my $self = shift;
my $buf;
my $length = $self->{_seq_length} * 2;
my $fh = $self->{_fh};
seek( $fh, $self->{PLOC}, 0 );
read( $fh, $buf, $length );
@{ $self->{_basecalls} } = unpack( "n" x $length, $buf );
# print "@{$self->{_basecalls}}" , "\n";
}
sub _set_corrected_base_calls {
my $self = shift;
my $buf;
my $length = $self->{_seq_length_corrected} * 2;
my $fh = $self->{_fh};
seek( $fh, $self->{PLOC1}, 0 );
read( $fh, $buf, $length );
@{ $self->{_basecalls_corrected} } = unpack( "n" x $length, $buf );
sub _set_seq {
my $self = shift;
my $buf;
my $length = $self->{_seq_length};
my $fh = $self->{_fh};
seek( $fh, $self->{PBAS2}, 0 );
read( $fh, $buf, $length );
$self->{_sequence} = $buf;
#my @seq = unpack( "C" x $length, $buf);
#print $buf, "\n";
}
sub _set_corrected_seq {
my $self = shift;
my $buf;
my $length = $self->{_seq_length_corrected};
my $fh = $self->{_fh};
seek( $fh, $self->{PBAS1}, 0 );
read( $fh, $buf, $length );
$self->{_sequence_corrected} = $buf;
my $self = shift;
my $buf;
my ( @pointers, @A, @G, @C, @T );
my (@datas) =
( $self->{DATA9}, $self->{DATA10}, $self->{DATA11}, $self->{DATA12} );
my $fh = $self->{_fh};
seek( $fh, $self->{FWO}, 0 );
read( $fh, $buf, 4 );
my @order = split( //, $buf );
#print "@order", "\n";
for ( my $i = 0 ; $i < 4 ; $i++ ) {
if ( $order[$i] =~ /A/i ) {
$pointers[0] = $datas[$i];
} elsif ( $order[$i] =~ /C/i ) {
$pointers[1] = $datas[$i];
} elsif ( $order[$i] =~ /G/i ) {
$pointers[2] = $datas[$i];
} elsif ( $order[$i] =~ /T/i ) {
$pointers[3] = $datas[$i];
} else {
Usage : my $seq_length = $abi->get_corrected_sequence_length();
Function : Returns the length of the edited sequence.
Args : Nothing
Returns : A scalar
=cut
sub get_corrected_sequence_length {
my $self = shift;
#print STDERR "**ABI**",$self->{_seq_length_corrected},"\n";
return $self->{_seq_length_corrected};
}
=head2 get_trace_length()
Title : get_trace_length()
Usage : my $trace_length = $abi->get_trace_length();
Function : Returns the trace length
Args : Nothing
Returns : A scalar
# Before `make install' is performed this script should be runnable with
# `make test'. After `make install' it should work as `perl test.pl'
#########################
# change 'tests => 1' to 'tests => last_test_to_print'
#use FindBin;
#use lib "$FindBin::Bin/..";
use Test::More qw(no_plan);
#BEGIN { plan tests => 1 };
require_ok("ABI");
my $abi = ABI->new("t/TEST_modified.ab1");
isa_ok( $abi, "ABI" );
my @base_calls = $abi->get_base_calls;
#print STDERR "@base_calls";
my @b_calls =
qw(20 33 46 56 67 82 89 102 111 124 136 151 163 174 185 192 204 215
225 233 246 258 274 285 297 308 324 340 349 356 369 383 393 403 413
422 434 447 458 467 477 487 499 513 523 529 544 552 562 570 581 593
603 615 624 634 644 657 668 680 690 702 714 726 736 746 756 766 776
784 796 808 820 830 842 854 863 874 884 897 906 918 930 940 949 960
972 983 996 1005 1017 1027 1039 1049 1063 1074 1086 1097 1109 1117 1131
1143 1153 1163 1174 1186 1196 1208 1218 1232 1243 1255 1266 1278 1288
1298 1310 1323 1333 1345 1357 1368 1380 1392 1402 1413 1425 1436 1448
1459 1469 1481 1493 1504 1513 1524 1537 1547 1559 1570 1583 1595 1606
GTCGAGGGGCTCAAGCGTTTGCGCGAGATCCTTGAGCGCAACGGTTTGTCCAGTGAGTCGCTG
TTTCAGTAACAGGCATCCTGCTCGCTAAAAAGCCCCGAAATATTCGGGGCTTTTTTGTGCCCG
CAGAATCTGGACCGCTGCTGCCAAGGGGTTTTTTTGAGTGCGTGCGGGTGACCGGTCAGTCTC
AAAAGTGCAGTCAGGCAGGGGTTGGAACTTTATCTGTCATGGGCTGTAAGCCTTTGCTTACCT
TTNATGTAAGCCAAGGGCGAAAACAGGCTTGCGGATAGNTTCGCTTCTGACTTTTCATAGGTT
GNAACTGATTGAAATTTAAACATTNTNATTGTTNTGNTAAGAN";
$seq =~ s/\n//g;
$seq =~ s/\s//g;
is( $seq, $abi->get_sequence() );
#print STDERR $abi->get_sequence_length();
is( 733, $abi->get_sequence_length() );
my @trace = qw(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 3 6 8 11 12 10 7 4 1 0 12 42 92 161 251 341 416 460 462 423
349 254 163 90 38 8 0 0 0 0 0 0 0 0 1 4 8 13 20 26 34 44 56 72 90 110 128
146 163 183 207 242 287 342 405 475 548 622 698 775 852 929 1008 1093 1187
1294 1417 1556 1600 1600 1600 1600 1600 1600 1600 1600 1600 1565 1395 1265
1174 1101 1020 914 775 611 439 281 159 75 25 2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 13 73 182 339 533 739 898 976 956 843 663 455 276 138 48
4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 23 101 236 420 642 857
0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 9 19 33 51 70 88 105 119 129 135 137 137 134 130 125 120 115 111 108 104 101 97 94 91 88 86 83 80 76 71 65 58 51 44 38 34 30 29 29 32 36 40 45 50 56
60 64 66 65 62 57 50 41 34 27 23 22 24 28 31 34 36 36 36 35 35 35 36 36 36 34 29 24 19 17 18 25 38 54 73 91 108 121 131 137 140 139 136 131 124 116 109 104 100 99 101 103 105 106 103 98 90 79 65 50
34 21 10 3 0 0 0 4 9 15 21 26 29 27 24 21 18 18 21 27 34 43 49 54 55 56 55 56 57 59 62 63 64 65 65 65 65 65 63 62 59 56 52 48 42 36 29 22 18 16 17 23 33 44 56 66 71 73 71 65
57 46 35 24 14 7 2 0 0 0 0 1 4 7 11 14 16 15 13 12 13 14 19 26 34 43 53 64 74 82 90 94 96 95 92 85 76 65 54 43 35 28 25 22 19 17 13 10 6 3 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4 10 19 28 37 43 45 40 30 22 12 4 5 22 35
);
is_deeply(\@trace,\@array);
is (8853, $abi->get_trace_length());
#print STDERR $abi->get_trace_length();
#my $abi1 = ABI->new("t/TEST_modified.abi");
is(734, $abi->get_corrected_sequence_length());
@base_calls = $abi->get_corrected_base_calls();
@b_calls =
qw(20 33 46 52 56 67 82 89 102 111 124 136 151 163 174 185 192 204 215
225 233 246 258 274 285 297 308 324 340 349 356 369 383 393 403 413
422 434 447 458 467 477 487 499 513 523 529 544 552 562 570 581 593
603 615 624 634 644 657 668 680 690 702 714 726 736 746 756 766 776
TTTCAGTAACAGGCATCCTGCTCGCTAAAAAGCCCCGAAATATTCGGGGCTTTTTTGTGCCCG
CAGAATCTGGACCGCTGCTGCCAAGGGGTTTTTTTGAGTGCGTGCGGGTGACCGGTCAGTCTC
AAAAGTGCAGTCAGGCAGGGGTTGGAACTTTATCTGTCATGGGCTGTAAGCCTTTGCTTACCT
TTNATGTAAGCCAAGGGCGAAAACAGGCTTGCGGATAGNTTCGCTTCTGACTTTTCATAGGTT
GNAACTGATTGAAATTTAAACATTNTNATTGTTNTGNTAAGAN";
$seq =~ s/\n//g;
$seq =~ s/\s//g;
is ($seq,$abi->get_corrected_sequence);
#print STDERR "\n", $abi->{_index},"\n";
#print STDERR "\n", $abi->{_seq_length_corrected},"\n";
#print "@array", "\n";
#ok(1); # If we made it this far, we're ok.
#########################
# Insert your test code below, the Test module is use()ed here so read
# its man page ( perldoc Test ) for help writing this test script.
( run in 1.924 second using v1.01-cache-2.11-cpan-de7293f3b23 )