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
|
# Mart ja Mikk may 2017
# Database building script for PlasmidSeeker
# Takes in multifasta file with all plasmid sequences as input
# Database in the format of list files, each plasmid separate, with names.txt for keeping real names
use warnings;
use strict;
use Getopt::Long;
# GLOBALS
my $fastafile;
my $glistmaker = "GenomeTester4/glistmaker";
my $dir_location = "plasmid_db";
my $help;
my $version;
my $time = localtime;
# PARAMETERS
my $word = 20;
my $threads = 32;
# OPTIONS
GetOptions(
'i=s' => \$fastafile,
'd=s' => \$dir_location,
'w=i' => \$word,
'h' => \$help,
'v' => \$version,
't=s' => \$threads,
) or die printHelp()."\n";
if ($version){ die "PlasmidSeeker builder v1.0 (16 May 2017)\n"; }
# Check presence of options and files
if(!$dir_location || !$fastafile) { die printHelp()."\n"; }
##########
### SUBS
##########
#General help
sub printHelp {
print "Usage: $0 -i <PLASMIDS FASTA FILE>\n";
print "Options:
-i\t Input fasta file with all plasmid sequences
-d\t Database directory (default „plasmid_db“, will be created if does not exist)
-t\t Number of threads used (default 32)
-w\t K-mer length used (default 20)\n
-h\t Print this help
-v\t Print version of the program\n";
return "";
}
##########
### MAIN
##########
###### Step 1 - extracts all plasmids as single fna files to a tmp directory
print STDERR "Getting all individual plasmid sequences and putting them to $dir_location\_fna\n";
open(FNA,'<',$fastafile) or die("Please specify plasmid fasta file!");
system "mkdir $dir_location\_fna";
my $count = 1;
while(<FNA>) {
if($_ =~ /^>/) {
close FILE;
open(FILE,'>',"$dir_location\_fna\/plasmid_$count\.fna");
$count++;
}
print FILE $_;
}
###### Step 2 - converts each fna file to a list file (multithreaded)
print STDERR "Converting plasmid fna sequences to list files in $dir_location\n";
opendir my $dir, "$dir_location\_fna" or die "Cannot open temporary fna dir $!";
my @fna_files = grep { $_ ne '.' && $_ ne '..' } readdir $dir;
closedir $dir;
system "mkdir $dir_location";
$count = 0;
my $cmd;
my $total = 0;
foreach (@fna_files) {
$count++;
$cmd .= "$glistmaker $dir_location\_fna\/$_ -w $word -o $dir_location\/$_ & ";
#Multithread execution if thread limit reached
if ($count == $threads) {
$total = $total + $threads;
print STDERR "Done: $total\r";
system "$cmd wait";
$count = 0;
$cmd = "";
}
}
# Check for unexecuted commands (last nodes if number of threads is not a multiple of node count)
if ($cmd ne "") {
system "$cmd wait";
}
###### Step 3 - build names.txt file to give all plasmids real names - default taken from FASTA headers
print STDERR "Creating names.txt file in $dir_location\n";
open(FILE,'>',"$dir_location/names.txt");
$count = scalar(@fna_files);
print FILE "# Database: $dir_location\tPlasmids total: $count\tBuilt on: $time\tK-mer length: $word\n#\n"; # Print all info to header
foreach(@fna_files) {
chomp;
my $cmd = "head -n 1 $dir_location\_fna\/$_";
my $head = qx/$cmd/;
#$head = (split(/\|\s+/,$head))[1]; Split only in RefSeq format
chomp $head;
print FILE "$_\t$head\n";
}
close FILE;
## FINAL - remove fna files and dir, leave only db lists
system "rm $dir_location\_fna/*.fna";
system "rmdir $dir_location\_fna";
|
