ABI
view release on metacpan or search on metacpan
}
$self->{_sequence} = "";
$self->{_sequence_corrected} = "";
$self->{_sample} = "";
$self->{A} = [];
$self->{T} = [];
$self->{G} = [];
$self->{C} = [];
$self->{_basecalls} = [];
$self->{_basecalls_corrected} = [];
$self->{_trace_length} = 0;
$self->{_seq_length} = 0;
$self->{_seq_length_corrected} = 0;
$self->{_abs_index} = 26;
$self->{_index} = undef;
$self->{PLOC1} = undef;
$self->{PLOC} = undef;
$self->{_a_start} = undef;
$self->{_g_start} = undef;
$self->{_c_start} = undef;
$self->{_t_start} = undef;
$self->{DATA9} = undef;
$self->{DATA10} = undef;
if ( $ploc_counter == 2 ) {
$self->{PLOC} = $self->{_index} + ( $i * 28 ) + 20;
}
}
if ( $buf eq "SMPL" ) {
$self->{SMPL} = $self->{_index} + ( $i * 28 ) + 20;
}
}
seek( $self->{_fh}, $self->{DATA12} - 8, 0 );
read( $self->{_fh}, $buf, 4 );
$self->{_trace_length} = unpack( "N", $buf );
seek( $self->{_fh}, $self->{PBAS2} - 4, 0 );
read( $self->{_fh}, $buf, 4 );
$self->{_seq_length} = unpack( "N", $buf );
seek( $self->{_fh}, $self->{PBAS1} - 4, 0 );
read( $self->{_fh}, $buf, 4 );
$self->{_seq_length_corrected} = unpack( "N", $buf );
$self->{PLOC} = $self->_get_int( $self->{PLOC} ) + $self->{_mac_header};
$self->{PLOC1} = $self->_get_int( $self->{PLOC1} ) + $self->{_mac_header};
$self->{DATA9} = $self->_get_int( $self->{DATA9} ) + $self->{_mac_header};
$self->{DATA10} = $self->_get_int( $self->{DATA10} ) + $self->{_mac_header};
$self->{DATA11} = $self->_get_int( $self->{DATA11} ) + $self->{_mac_header};
$self->{DATA12} = $self->_get_int( $self->{DATA12} ) + $self->{_mac_header};
$self->{PBAS1} = $self->_get_int( $self->{PBAS1} ) + $self->{_mac_header};
$self->{PBAS2} = $self->_get_int( $self->{PBAS2} ) + $self->{_mac_header};
$self->{SMPL} += $self->{_mac_header};
}
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;
}
sub _set_traces {
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};
} elsif ( $order[$i] =~ /G/i ) {
$pointers[2] = $datas[$i];
} elsif ( $order[$i] =~ /T/i ) {
$pointers[3] = $datas[$i];
} else {
croak "Wrong traces\n";
}
}
for ( my $i = 0 ; $i < 4 ; $i++ ) {
seek( $fh, $pointers[$i], 0 );
read( $fh, $buf, $self->{_trace_length} * 2 );
if ( $i == 0 ) {
@A = unpack( "n" x $self->{_trace_length}, $buf );
}
if ( $i == 1 ) {
@C = unpack( "n" x $self->{_trace_length}, $buf );
}
if ( $i == 2 ) {
@G = unpack( "n" x $self->{_trace_length}, $buf );
}
if ( $i == 3 ) {
@T = unpack( "n" x $self->{_trace_length}, $buf );
}
}
@{ $self->{A} } = @A;
@{ $self->{G} } = @G;
@{ $self->{T} } = @T;
@{ $self->{C} } = @C;
}
sub _get_int {
my $self = shift;
}
$self->{_maximum_trace} = $max;
}
sub _set_sample_name {
my $self = shift;
my $buf;
my $fh = $self->{_fh};
seek( $fh, $self->{SMPL}, 0 );
read( $fh, $buf, 1 );
my $length = unpack( "C", $buf );
read( $fh, $buf, $length );
$self->{_sample} = $buf;
}
=head1 METHODS
=head2 get_max_trace()
Title : get_max_trace()
Usage : $max = $abi->get_max_trace();
Function : Returns the maximum trace value of all the traces.
Args : Nothing
Returns : A scalar
=cut
sub get_corrected_sequence {
my $self = shift;
return $self->{_sequence_corrected};
}
=head2 get_sequence_length()
Title : get_sequence_length()
Usage : my $seq_length = $abi->get_sequence_length();
Function : Returns the sequence length of the orginal unedited sequence.
Args : Nothing
Returns : A scalar
=cut
sub get_sequence_length {
my $self = shift;
return $self->{_seq_length};
}
=head2 get_corrected_sequence_length()
Title : get_corrected_sequence_length()
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
=cut
sub get_trace_length {
my $self = shift;
return $self->{_trace_length};
}
=head2 get_base_calls()
Title : get_base_calls()
Usage : my @base_calls = $abi->get_base_calls();
Function : Returns the called bases by the base caller. This method will return the unedited
original basecalls created by the basecaller.
Args : Nothing
Returns : An array
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
1007 1063 1022 913 784 681 633 641 678 706 694 625 507 368 248 183 197
82 79 75 71 66 60 55 51 50 50 52 54 56 55 52 47 40 32 24 16 9 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 2 2 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 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
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
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 0.488 second using v1.01-cache-2.11-cpan-65fba6d93b7 )