Commit 705a4a19 authored by Sebastien's avatar Sebastien

Paralellisation parseBlast via un outil companion et des objets storage

Utilise le NB CPUS blast pour nombre de tranches
parent e83287a5
#!/usr/bin/perl
use strict;
use Bio::SearchIO;
use Storable;
use FindBin;
MAIN:
{
my $blastfile = $ARGV[0];
my $seqlen_storable = $ARGV[1];
my %blast_flag = %{retrieve ($ARGV[2])}; ;
my $pv_cutoff = $ARGV[3];
my %seq_len = %{retrieve ($seqlen_storable)};
my $parseoutfile = "$blastfile.parse";
open (PARSEOUT,">$parseoutfile");
my $searchio = Bio::SearchIO->new(-file => $blastfile,
-format => $blast_flag{'format'}) or die("Blast parsing failed!");
while (my $result = $searchio->next_result()) {
my $queryid=$result->query_name;
my $querylen;
if (defined $seq_len{$queryid}) {
$querylen=$seq_len{$queryid};
} else {
$querylen=$result->query_length; # query length is not stored in BLAST m8 format, so querylen will be 0
}
while( my $hit = $result->next_hit ) {
next unless numeric_pvalue($hit->significance) <= $pv_cutoff;
my $subjectid=$hit->name;
my $subjectlen;
if (defined $seq_len{$subjectid}) {
$subjectlen=$seq_len{$subjectid};
} else {
$subjectlen=$hit->length; # subject length is not stored in BLAST m8 format, so querylen will be 0
}
my $pvalue=numeric_pvalue($hit->significance);
if ($blast_flag{'hsp'}) {
my $simspanid=1;
my $simspan='';
my (@percentidentity,@hsplength);
while( my $hsp = $hit->next_hsp ) {
my $querystart=$hsp->start('query');
my $queryend=$hsp->end('query');
my $subjectstart=$hsp->start('sbjct');
my $subjectend=$hsp->end('sbjct');
$percentidentity[$simspanid]=$hsp->percent_identity;
$hsplength[$simspanid]=$hsp->length('hit');
$simspan.="$simspanid:$querystart-$queryend:$subjectstart-$subjectend.";
$simspanid++;
}
my $sum_identical=0;
my $sum_length=0;
for (my $i=1;$i<$simspanid;$i++) {
$sum_identical+=$percentidentity[$i]*$hsplength[$i];
$sum_length+=$hsplength[$i];
}
my $percentIdent=int($sum_identical/$sum_length);
print PARSEOUT "$queryid;$querylen;$subjectid;$subjectlen;$pvalue;$percentIdent;$simspan\n";
} else {
print PARSEOUT "$queryid;$querylen;$subjectid;$subjectlen;$pvalue;0;NULL\n";
}
}
}
close(PARSEOUT);
}
# Make pvalue numeric, used by subroutine blast_parse
# One Arguments:
# 1. String Variable: pvalue
# Last modified: 07/19/04
sub numeric_pvalue {
my $p=$_[0];
if ($p=~/^e-(\d+)/) {return "1e-".$1;}
else {return $p}
} # numeric_pvalue
......@@ -19,6 +19,7 @@ use lib "$FindBin::RealBin";
use Getopt::Long;
use File::Basename;
use Storable;
use orthomcl_module;
......
......@@ -503,62 +503,46 @@ sub blast_parse {
my %seq_len = %{$_[3]};
my %blast_flag = %{$_[4]};
open (PARSEOUT,">$parseoutfile");
write_log("\nParsing blast result!\n");
my $searchio = Bio::SearchIO->new(-file => $blastfile,
-format => $blast_flag{'format'}) or dieWithUnexpectedError("Blast parsing failed!");
store (\%seq_len, "$parseoutfile.seqlen");
store (\%blast_flag, "$parseoutfile.blast_flag");
my $split_cmd = "split -d -n l/$BLAST_NOCPU $blastfile $blastfile.split.";
system ($split_cmd);
my @a_split_files = `ls $blastfile.split.*`;
chomp @a_split_files;
#PARALELLISER:
# splitter le $blastfile en NB CPUS
# splitter le $blastfile en NB CPUS $BLAST_NOCPU
# ecrire un utils.pl pour faire ce filtre en ligne de commande sur chaque tranche
# merger les fichiers filtrés en ajoutant le $similarityid
my $similarityid=1;
while (my $result = $searchio->next_result()) {
my $queryid=$result->query_name;
my $querylen;
if (defined $seq_len{$queryid}) {
$querylen=$seq_len{$queryid};
} else {
$querylen=$result->query_length; # query length is not stored in BLAST m8 format, so querylen will be 0
}
while( my $hit = $result->next_hit ) {
next unless numeric_pvalue($hit->significance) <= $pv_cutoff;
my $subjectid=$hit->name;
my $subjectlen;
if (defined $seq_len{$subjectid}) {
$subjectlen=$seq_len{$subjectid};
} else {
$subjectlen=$hit->length; # subject length is not stored in BLAST m8 format, so querylen will be 0
}
my $pvalue=numeric_pvalue($hit->significance);
if ($blast_flag{'hsp'}) {
my $simspanid=1;
my $simspan='';
my (@percentidentity,@hsplength);
while( my $hsp = $hit->next_hsp ) {
my $querystart=$hsp->start('query');
my $queryend=$hsp->end('query');
my $subjectstart=$hsp->start('sbjct');
my $subjectend=$hsp->end('sbjct');
$percentidentity[$simspanid]=$hsp->percent_identity;
$hsplength[$simspanid]=$hsp->length('hit');
$simspan.="$simspanid:$querystart-$queryend:$subjectstart-$subjectend.";
$simspanid++;
}
my $sum_identical=0;
my $sum_length=0;
for (my $i=1;$i<$simspanid;$i++) {
$sum_identical+=$percentidentity[$i]*$hsplength[$i];
$sum_length+=$hsplength[$i];
}
my $percentIdent=int($sum_identical/$sum_length);
print PARSEOUT "$similarityid;$queryid;$querylen;$subjectid;$subjectlen;$pvalue;$percentIdent;$simspan\n";
} else {
print PARSEOUT "$similarityid;$queryid;$querylen;$subjectid;$subjectlen;$pvalue;0;NULL\n";
my $parse_cmd_file = "$parseoutfile.cmd";
open (CMD,">$parse_cmd_file");
foreach my $splitfile ( @a_split_files )
{
print CMD "$FindBin::RealBin/orthomcl-parseblast.pl $splitfile $parseoutfile.seqlen $parseoutfile.blast_flag $pv_cutoff \n";
}
close (CMD);
system ("(cat $parseoutfile.cmd | parallel -j $BLAST_NOCPU > $parseoutfile.cmd.out ) 2> $parseoutfile.cmd.err");
my $similarityid=1;
open (PARSEOUT,">$parseoutfile");
foreach my $splitfile ( @a_split_files )
{
open (INFILE, "$splitfile.parse");
while (my $line = <INFILE>)
{
print PARSEOUT "$similarityid;$line";
$similarityid++;
}
close (INFILE);
}
#orthomcl-parseblast.pl $splitfile $parseoutfile.seqlen $blast_flag{'format'}
write_log("Parsing blast file finished\n");
close(PARSEOUT);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment