-
Notifications
You must be signed in to change notification settings - Fork 1
/
multiT_runLastal.pl
70 lines (62 loc) · 1.87 KB
/
multiT_runLastal.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
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/perl -w
use threads;
use threads::shared;
die "Usage: perl $0 query.fasta database.fasta threads identity out.aln\n"
if(!defined($ARGV[0]) || !defined($ARGV[1]) || !defined($ARGV[2]) || !defined($ARGV[4]));
my $query = $ARGV[0];
my $dataB = $ARGV[1];
my $maxThreads = int $ARGV[2]/2;
my $identity = (defined $ARGV[3])?$ARGV[3]:60;
my $outaln = $ARGV[4];
system("rm -rf faByChr");
system("mkdir faByChr");
system("rm -rf tmp");
system("mkdir tmp");
#system("rm dblast*");
#system("lastdb -uNEAR -R01 dblast $dataB");
my %infordb;
my %namedb;
my $count = 0;
open(IN, $query) or die"";
$/='>';
<IN>;
while(<IN>){
chomp;
my ($name,$seq) = split(/\n/,$_,2);
$name =~ s/\s+.*//g;
$seq =~ s/\s+//g;
$count++;
$infordb{$count}->{'name'} = $name;
$infordb{$count}->{'seq'} = $seq;
my $rname = $count;
$namedb{$rname} = $name;
$infordb{$count}->{'rname'} = $rname;
my $outfile = $rname.".fasta";
open(my $out, ">faByChr/$outfile") or die"";
print $out ">$rname\n$seq\n";
close $out;
}
close IN;
system("rm command.txt");
system("ls faByChr/*fasta |xargs -I\{\} basename \{\} .fasta|xargs -I\{\} echo \"lastal -P 1 dblast faByChr/{}.fasta -a 180 -b 1 -q 80 -r 30|maf-convert blasttab >tmp/\{\}.out\" >>command.txt");
system("ParaFly -c command.txt -CPU $maxThreads");
open(OUT, "> $outaln") or die"";
open(IN, "cat tmp/*.out |awk '\$3>88' |") or die"";
my $content = <IN>;
@linedb = split(/\n/,$content);
foreach my $line (@linedb){
my @data = split(/\s+/,$line);
$data[0] = $namedb{$data[0]} if(exists($namedb{$data[0]}));
foreach my $i(0..$#data){
print OUT "$data[$i] ";
}
print OUT "\n";
}
close IN;
close OUT;
system("rm -rf tmp/");
#sub run_lastal{
# my $scaf = shift;
# my $cmd = "lastal -P 2 dblast faByChr/".$scaf.".fasta -a 180 -b 1 -q 80 -r 30 |maf-convert blasttab |awk '\$3>$identity' > tmp/".$scaf.".out";
# system($cmd);
# }