|
0 |
--- a/bin/barrnap
|
|
1 |
+++ b/bin/barrnap
|
|
2 |
@@ -4,8 +4,9 @@
|
|
3 |
use Time::Piece;
|
|
4 |
use List::Util qw(max);
|
|
5 |
use FindBin;
|
|
6 |
+use File::Temp qw/ tempfile tempdir /;
|
|
7 |
|
|
8 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
9 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
10 |
# global variables
|
|
11 |
|
|
12 |
my $EXE = $FindBin::RealScript;
|
|
13 |
@@ -29,13 +30,13 @@
|
|
14 |
);
|
|
15 |
my $MAXLEN = int( 1.2 * max(values %LENG) );
|
|
16 |
|
|
17 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
18 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
19 |
# command line options
|
|
20 |
|
|
21 |
my(@Options, $quiet, $kingdom, $threads, $evalue, $lencutoff, $reject, $incseq);
|
|
22 |
setOptions();
|
|
23 |
|
|
24 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
25 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
26 |
# check all is well
|
|
27 |
|
|
28 |
msg("This is $EXE $VERSION");
|
|
29 |
@@ -64,23 +65,38 @@
|
|
30 |
my $hmmdb = "$DBDIR/$kdom.hmm";
|
|
31 |
err("Can't find database: $hmmdb") unless -r $hmmdb;
|
|
32 |
msg("Using database: $hmmdb");
|
|
33 |
+my $nonfreehmmdb = "$DBDIR/nonfree/$kdom.hmm";
|
|
34 |
|
|
35 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
36 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
37 |
# run the external command
|
|
38 |
|
|
39 |
my $fasta = shift @ARGV;
|
|
40 |
$fasta && -r $fasta or err("Usage: $EXE <file.fasta>");
|
|
41 |
|
|
42 |
+if (-e $nonfreehmmdb) {
|
|
43 |
+ # prepare full pHMM file (with nonfree if it exists)
|
|
44 |
+ open ( HMMDB, "<", $hmmdb ) or die "Could not open file $hmmdb: $!";
|
|
45 |
+ open ( NONFREE, "<", $nonfreehmmdb );
|
|
46 |
+ my ($tmpfile, $tmpfilename) = tempfile();
|
|
47 |
+ while ( my $line = <HMMDB> ) {
|
|
48 |
+ print $tmpfile $line;
|
|
49 |
+ }
|
|
50 |
+ while ( my $line = <NONFREE> ) {
|
|
51 |
+ print $tmpfile $line;
|
|
52 |
+ }
|
|
53 |
+ $hmmdb = $tmpfilename
|
|
54 |
+}
|
|
55 |
+
|
|
56 |
msg("Scanning $fasta for $kdom rRNA genes... please wait");
|
|
57 |
my $cmd = "$NHMMER --cpu $threads -E $evalue --w_length $MAXLEN -o /dev/null --tblout /dev/stdout \Q$hmmdb\E \Q$fasta\E";
|
|
58 |
msg("Command: $cmd");
|
|
59 |
my @hits = qx($cmd 2>&1);
|
|
60 |
|
|
61 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
62 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
63 |
# process the output
|
|
64 |
|
|
65 |
my @feat;
|
|
66 |
-HIT:
|
|
67 |
+HIT:
|
|
68 |
foreach (@hits) {
|
|
69 |
chomp;
|
|
70 |
err("nhmmer failed to run - $_") if m/fail|error|core dump|bus error/i;
|
|
71 |
@@ -102,7 +118,7 @@
|
|
72 |
# check for incomplete alignments
|
|
73 |
my $note = '';
|
|
74 |
my $len = $end-$begin+1;
|
|
75 |
-
|
|
76 |
+
|
|
77 |
if ( $len < int($reject * $LENG{$gene}) ) {
|
|
78 |
msg("Rejecting short $len nt predicted $gene. Adjust via --reject option.");
|
|
79 |
next HIT;
|
|
80 |
@@ -120,7 +136,7 @@
|
|
81 |
];
|
|
82 |
}
|
|
83 |
|
|
84 |
-# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
85 |
+# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
86 |
# output a sorted GFF3
|
|
87 |
|
|
88 |
sub gff_sort {
|
|
89 |
@@ -168,13 +184,13 @@
|
|
90 |
|
|
91 |
sub show_citation {
|
|
92 |
print STDERR << "EOCITE";
|
|
93 |
-
|
|
94 |
+
|
|
95 |
If you use Barrnap in your work, please cite:
|
|
96 |
|
|
97 |
Seemann T (2013)
|
|
98 |
$EXE $VERSION : $DESC
|
|
99 |
$URL
|
|
100 |
-
|
|
101 |
+
|
|
102 |
Thank you.
|
|
103 |
|
|
104 |
EOCITE
|
|
105 |
@@ -193,7 +209,7 @@
|
|
106 |
{OPT=>"help", VAR=>\&usage, DESC=>"This help"},
|
|
107 |
{OPT=>"version", VAR=>\&version, DESC=>"Print version and exit"},
|
|
108 |
{OPT=>"citation",VAR=>\&show_citation, DESC=>"Print citation for referencing $EXE"},
|
|
109 |
- {OPT=>"kingdom=s", VAR=>\$kingdom, DEFAULT=>'bac',
|
|
110 |
+ {OPT=>"kingdom=s", VAR=>\$kingdom, DEFAULT=>'bac',
|
|
111 |
DESC=>"Kingdom: ".join(' ', values %KINGDOM) },
|
|
112 |
{OPT=>"quiet!", VAR=>\$quiet, DEFAULT=>0, DESC=>"No screen output"},
|
|
113 |
{OPT=>"threads=i", VAR=>\$threads, DEFAULT=>8, DESC=>"Number of threads/cores/CPUs to use"},
|
|
114 |
@@ -226,15 +242,15 @@
|
|
115 |
my $def = defined($_->{DEFAULT}) ? " (default '$_->{DEFAULT}')" : "";
|
|
116 |
$def = ($def ? ' (default OFF)' : '(default ON)') if $_->{OPT} =~ m/!$/;
|
|
117 |
my $opt = $_->{OPT};
|
|
118 |
- $opt =~ s/!$//;
|
|
119 |
- $opt =~ s/=s$/ [X]/;
|
|
120 |
+ $opt =~ s/!$//;
|
|
121 |
+ $opt =~ s/=s$/ [X]/;
|
|
122 |
$opt =~ s/=i$/ [N]/;
|
|
123 |
$opt =~ s/=f$/ [n.n]/;
|
|
124 |
printf STDERR " --%-15s %s%s\n", $opt, $_->{DESC}, $def;
|
|
125 |
}
|
|
126 |
else {
|
|
127 |
print STDERR "$_\n";
|
|
128 |
- }
|
|
129 |
+ }
|
|
130 |
}
|
|
131 |
exit(1);
|
|
132 |
}
|