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
|
package Fasta_retriever;
use strict;
use warnings;
use Carp;
use threads;
use threads::shared;
my $LOCKVAR :shared;
our $DEBUG = 0;
sub new {
my ($packagename) = shift;
my $filename = shift;
unless ($filename) {
confess "Error, need filename as param";
}
my $self = { filename => $filename,
acc_to_pos_index => undef,
fh => undef,
};
my %acc_to_pos_index :shared;
$self->{acc_to_pos_index} = \%acc_to_pos_index;
bless ($self, $packagename);
$self->_init();
return($self);
}
sub _init {
my $self = shift;
my $filename = $self->{filename};
# use a samtools faidx index if available
my $index_file = "$filename.fai";
if (-s $index_file) {
open(my $fh, $index_file) or die "Error, cannot open file: $index_file";
while(<$fh>) {
chomp;
my @x = split(/\t/);
my $acc = $x[0];
my $file_pos = $x[2];
$self->{acc_to_pos_index}->{$acc} = $file_pos;
}
close $fh;
}
else {
print STDERR "-missing faidx file: $index_file, extracting positions directly.\n";
print STDERR "-Fasta_retriever:: begin initializing for $filename\n";
open (my $fh, $filename) or die $!;
$self->{fh} = $fh;
while (<$fh>) {
if (/>(\S+)/) {
my $acc = $1;
my $file_pos = tell($fh);
$self->{acc_to_pos_index}->{$acc} = $file_pos;
}
}
print STDERR "-Fasta_retriever:: done initializing for $filename\n";
}
return;
}
sub refresh_fh {
my $self = shift;
open (my $fh, $self->{filename}) or die "Error, cannot open file : " . $self->{filename};
$self->{fh} = $fh;
return $fh;
}
sub get_seq {
my $self = shift;
my $acc = shift;
unless (defined $acc) {
confess "Error, need acc as param";
}
my $seq = "";
{
lock $LOCKVAR;
my $file_pos = $self->{acc_to_pos_index}->{$acc} or confess "Error, no seek pos for acc: $acc";
my $fh = $self->refresh_fh();
seek($fh, $file_pos, 0);
print STDERR "seeking $acc -> $file_pos\n" if $DEBUG;
while (<$fh>) {
if (/^>/) {
print STDERR " reached $_, stopping\n" if $DEBUG;
last;
}
$seq .= $_;
}
print STDERR "-done seeking $acc\n\n" if $DEBUG;
}
$seq =~ s/\s+//g;
return($seq);
}
1; #EOM
|
