-
Notifications
You must be signed in to change notification settings - Fork 30
/
Copy pathdbm_read_fa.pl
executable file
·59 lines (50 loc) · 1.14 KB
/
dbm_read_fa.pl
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
#!/usr/bin/perl -w
#
#Author: Ruan Jue
#
use strict;
use DB_File;
my $dbf = shift or die("Usage: $0 <dbm_file> [read_name1 ...]\n");
if($dbf!~/\.dbm$/){
$dbf .= ".dbm" if(-e "$dbf.dbm");
}
my @names = @ARGV;
if(@names == 0){
while(<>){
chomp;
push(@names, $_);
}
}
my @tags = ();
foreach my $tag (@names){
if($tag=~/^(.+?)(\[([+-])(:(-?\d+),(-?\d+))?\])$/){
push(@tags, [$1, $3 eq '+'? 1:2, (defined $5)? $5:1, (defined $6)? $6:-1, 1]);
} else {
push(@tags, [$tag, 1, 1, -1, 0]);
}
}
my %seqs;
tie %seqs, 'DB_File', $dbf, O_RDONLY or die "Cannot open $dbf: $!";
foreach my $tag (@tags){
if(exists $seqs{$tag->[0]}){
my $seq = $seqs{$tag->[0]};
$tag->[3] = length($seq) if($tag->[3] < 1);
if($tag->[4]){
print ">", join("_", $tag->[0], $tag->[1] == 1? "F":"R", $tag->[2], $tag->[3]), "\n";
} else {
print ">$tag->[0]\n";
}
if($tag->[2] < $tag->[3]){
my $ss = substr($seq, $tag->[2] - 1, $tag->[3] - $tag->[2] + 1);
if($tag->[1] == 2){
$ss =~tr/ACGTacgt/TGCAtgca/;
$ss = reverse $ss;
}
while($ss=~/(.{1,100})/g){ print "$1\n"; }
}
} else {
warn("'$tag->[0]' was not found\n");
}
}
untie %seqs;
1;