Commit f7a00600 authored by Sebastien's avatar Sebastien

RunOrthoMCL: Precalcul blast avec Diamond

Paralellisation par proteome vs. all
parent 7ffbba36
......@@ -5,6 +5,7 @@ use FindBin;
use lib "$FindBin::RealBin/../../lib/int/";
use lib "$FindBin::RealBin/../../lib/ext/lipmutils";
use General;
use ParamParser;
use GeneralBioinfo;
......@@ -90,6 +91,7 @@ our $WORKDIR;
our @A_OUTPUT_DATA = ();
our $OUTDIR;
our $O_CONF;
our $THREADS = 4;
MAIN:
{
......@@ -137,6 +139,11 @@ MAIN:
$rh_analysis_params->{'homologygroups'}->{'version'},
$rh_seqidspeciescode);
store $rh_homologygroups, "$WORKDIR/rh_homologygroups";
store $rh_species, "$WORKDIR/rh_species";
store $rh_seqidspeciescode, "$WORKDIR/rh_seqidspeciescode";
store $rh_analysis_params, "$WORKDIR/rh_analysis_params";
my ($rh_count_proteins_in_group_by_code, $rh_count_inparalogs_by_code, $rh_count_specific_inparalogs_by_code) =
&CheckProteomeCodesConsistency($rh_species, $rh_homologygroups);
my $rh_specific_proteins = &GetSpecificProteins($rh_homologygroups, $rh_species);
......@@ -177,11 +184,14 @@ MAIN:
sub RunOrthomcl
{
my $rh_analysis_params = shift;
RunBlast($rh_analysis_params);
my $cmd =
"( cd $WORKDIR; $ORTHOMCL"
. ' --mode 1 --fa_files '
. ' --mode 3 --fa_files '
. join('.fasta,', keys %{$rh_analysis_params->{'proteomes'}})
. '.fasta ';
. '.fasta --blast_file all.blast ';
my $parameters = $rh_analysis_params->{'homologygroups'}->{'parameters'};
$cmd .= " $parameters ";
......@@ -207,6 +217,50 @@ sub RunOrthomcl
return;
}
sub RunBlast
{
my $rh_analysis_params = shift;
#~ my $cmd =
#~ "( cd $WORKDIR; cat "
#~ . join('.fasta ', keys %{$rh_analysis_params->{'proteomes'}}) . '.fasta > all.fasta;'
#~ . "$FindBin::RealBin/../ext/usearch -makeudb_usearch $WORKDIR/all.fasta -output $WORKDIR/all.fasta.udb 2>$WORKDIR/all.fasta.makeudb.log)";
#~ my $cmd =
#~ "( cd $WORKDIR; cat "
#~ . join('.fasta ', keys %{$rh_analysis_params->{'proteomes'}}) . '.fasta > all.fasta;'
#~ . "$FindBin::RealBin/../ext/ncbi-blast/makeblastdb -in $WORKDIR/all.fasta -parse_seqids -dbtype prot 2>$WORKDIR/all.fasta.makeudb.log)";
my $cmd =
"( cd $WORKDIR; cat "
. join('.fasta ', keys %{$rh_analysis_params->{'proteomes'}}) . '.fasta > all.fasta;'
. "$FindBin::RealBin/../ext/diamond makedb --in $WORKDIR/all.fasta --db $WORKDIR/all.fasta.diamond 2>$WORKDIR/all.fasta.diamond.log)";
&System($cmd);
my $parameters = $rh_analysis_params->{'homologygroups'}->{'parameters'};
my ( $evalue ) = ($parameters =~ /pv_cutoff=(\S+)/);
my $fh_cmd_blast = &GetStreamOut("$WORKDIR/cmd.blast");
foreach my $proteome (keys %{$rh_analysis_params->{'proteomes'}})
{
#~ print $fh_cmd_blast "nice $FindBin::RealBin/../ext/usearch -ublast $WORKDIR/$proteome.fasta -db $WORKDIR/all.fasta.udb -evalue 1 -blast6out $WORKDIR/all.fasta.$proteome.blast\n";
#~ print $fh_cmd_blast "nice $FindBin::RealBin/../ext/ncbi-blast/blastp -query $WORKDIR/$proteome.fasta -db $WORKDIR/all.fasta -evalue $evalue -out $WORKDIR/all.fasta.$proteome.blast -outfmt '6 std' -num_threads 2 -max_target_seqs 1000\n";
print $fh_cmd_blast "nice $FindBin::RealBin/../ext/diamond blastp --query $WORKDIR/$proteome.fasta --db $WORKDIR/all.fasta.diamond --evalue $evalue --out $WORKDIR/all.fasta.$proteome.blast --threads 2 --max-target-seqs 1000\n";
}
$fh_cmd_blast->close();
$cmd = "(cat $WORKDIR/cmd.blast | parallel -j $THREADS > $WORKDIR/cmd.blast.out ) 2> $WORKDIR/cmd.blast.err";
&System($cmd);
$cmd = "( cd $WORKDIR; cat all.fasta.*.blast > all.blast)";
&System($cmd);
return;
}
=head2 procedure RunOrthoFinder
Title : RunOrthoFinder
......
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