Question: Perl script review
0
gravatar for Mostafa Mahmoud
3 months ago by
Mostafa Mahmoud0 wrote:

I wrote Perl script to run spades commands but the script doesn't process the cmd

 #!/usr/bin/env perl
use strict;
use warnings;
use Getopt::Long qw(:config no_ignore_case bundling);
use FindBin;
use Cwd;

######################################################
### Set to base directory :
#my $BASEDIR = ".";
#######################################################

my $usage = <<__EOUSAGE__;
##########################################################################################################
##
##  Required:
##
##  --kmer  <string>                K-mer size
##
##  --samples_file <string>         tab-delimited text file indicating biological replicate relationships.
##                                   ex.
##                                        cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq
##                                        cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq
##                                        cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq
##                                        cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq
##
##                                   
## --CPU=6
## --max_memory=10G
##                               
##
##
##
##  Optional:
##
##  -I                              Interactive mode, waits between commands.
##
############################################################################################################


__EOUSAGE__

    ;

my $help_flag;
my $read_samples_descr_file;
my $PAUSE_EACH_STEP = 0;
my $mem;
my $kmer;
my $cpu;
&GetOptions ( 'CPU=s' => \$cpu,
          'kmer=s' => \$kmer,
              'max_memory=s' => \$mem,
              'h' => \$help_flag,
              'samples_file=s' => \$read_samples_descr_file,
              'I' => \$PAUSE_EACH_STEP,

);

if ($help_flag) {
    die $usage;
}

unless ( $read_samples_descr_file) {
    die $usage;
}

{

## Check for required software
my @needed_tools = qw spades.py); 
    my $missing_flag = 0;
    foreach my $prog (@needed_tools) {
        my $path = `which $prog`;
        unless ($path =~ /\w/) {
            print STDERR "\n** ERROR, cannot find path to required software: $prog **\n";
            $missing_flag = 1;
}
    }
    if ($missing_flag) {
        die "\nError, at least one required software tool could not be found. Please install tools and/or adjust your PATH settings before retryi
ng.\n";
    }
}


my $workdir = cwd();
my %conditions_to_read_info = &parse_sample_descriptions($read_samples_descr_file);
my @conditions = sort keys %conditions_to_read_info;
foreach my $condition (@conditions) {
my $replicates_href = $conditions_to_read_info{$condition};
my @replicates = keys %$replicates_href;
foreach my $replicate (@replicates) {
my ($left_fq_file, $right_fq_file) = @{$replicates_href->{$replicate}};




        my $cmd = "spades.py --rna   "
                  ."-k $kmer" 
                  ."-t $cpu"
                  ."-1 $$left_fq_file"
                  ."-2 $right_fq_file"
                  ."-m $mem"
                  ."-o spades/$replicate"
            ;
if ($left_fq_file && $right_fq_file) {
            $cmd .= " --left $left_fq_file --right $right_fq_file ";
        }
        elsif ($left_fq_file) {
            $cmd .= " --single $left_fq_file ";
        }
        else {
            die "Error, no left or right read files";
        }



        ;&process_cmd($cmd,"Running assemply  on sample $replicate"),}}

sub process_cmd {
    my ($cmd, $msg) = @_;


    if ($msg) {
        print "\n\n";
        print "#################################################################\n";
        print "$msg\n";
        print "#################################################################\n";
    }

   print "CMD: $cmd\n";
    if ($PAUSE_EACH_STEP) {
        print STDERR "\n\n-WAITING, PRESS RETURN TO CONTINUE ...";
        my $wait = <STDIN>;
        print STDERR "executing cmd.\n\n";

    }


    my $time_start = time();

    my $ret = system($cmd);
    my $time_end = time();

    if ($ret) {
        die "Error, CMD: $cmd died with ret $ret";
    }

    my $number_minutes = sprintf("%.1f", ($time_end - $time_start) / 60);

    print "TIME: $number_minutes min. for $cmd\n";


    return;
}

sub parse_sample_descriptions {
    my ($read_samples_descr_file) = @_;

    my %samples_descr;


    open (my $fh, $read_samples_descr_file) or die $!;
    while (<$fh>) {
        if (/^\#/) { next; }
        unless (/\w/) { next; }
        s/^\s+|\s+$//g;
        chomp;
        my @x = split(/\t/);
        if (/=/ && scalar(@x) == 1) {
            my ($key, $val) = split(/=/, $x[0]);

        }
        else {


            my ($condition, $replicate, $reads_left, $reads_right) = @x;

            #remove gzip extension, will convert to gunzipped version later 
                        $reads_left =~ s/\.gz$//;
                                    $reads_right =~ s/\.gz$// if $reads_right;

      }
    }

    close $fh;

    return(%samples_descr);

}
spades perl • 299 views
ADD COMMENTlink modified 3 months ago by Michael Dondrup45k • written 3 months ago by Mostafa Mahmoud0

I wrote Perl script to run spades commands but the script doesn't process the cmd

what is the question ?

ADD REPLYlink written 3 months ago by Pierre Lindenbaum116k

actually, I am new to Perl I wrote this script to run spades by reading the samples one by one from a text file but when I run it shows nothing and I can't figure out where is the problem ???

ADD REPLYlink written 3 months ago by Mostafa Mahmoud0

Please use ADD REPLY for comments or edit your question by incorporating your comment, but not the answer field.

ADD REPLYlink written 3 months ago by ATpoint13k

And, what happens? This is a rather complex perl wrapper to run a python tool and it just maps the parameters over and reads some of them from a file, I would normally wrap the command call into a shell script for one-off calls. First, make sure the posted code doesn't contain a syntax error, if you want to continue using the script. Then, please post example input files and the command call as well as the output of the call.

perl -c test-script.pl
Bareword "py" not allowed while "strict subs" in use at test-script.pl line 70.
syntax error at test-script.pl line 70, near "py)"
test-script.pl had compilation errors.
ADD REPLYlink modified 3 months ago • written 3 months ago by Michael Dondrup45k

it just maps the parameters and read the sample and run every in sample file line one by one

cat sample file



cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq
cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq
cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq
cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq
ADD REPLYlink modified 3 months ago • written 3 months ago by Mostafa Mahmoud0

Yes, I can see that it is supposed to do that. Please fix the syntax error first if it is not just a copy paste problem.

ADD REPLYlink written 3 months ago by Michael Dondrup45k
0
gravatar for Michael Dondrup
3 months ago by
Bergen, Norway
Michael Dondrup45k wrote:

The problem is in the function below. You want to return a hash %samples_descr but you never assign anything to it. Therefore, there are no keys in the hash and your outer loop is never entered.

 sub parse_sample_descriptions {
    my ($read_samples_descr_file) = @_;

    my %samples_descr = (); #initialize empty hash


    open (my $fh, $read_samples_descr_file) or die $!;
    while (<$fh>) {
        if (/^\#/) { next; }
        unless (/\w/) { next; }
        s/^\s+|\s+$//g;
        chomp;
        my @x = split(/\t/);
        if (/=/ && scalar(@x) == 1) {
            my ($key, $val) = split(/=/, $x[0]);
          # ignoring this for now, don't understand which format this is for
        }
        else {


            my ($condition, $replicate, $reads_left, $reads_right) = @x;

            #remove gzip extension, will convert to gunzipped version later 
                        $reads_left =~ s/\.gz$//;
                        $reads_right =~ s/\.gz$// if $reads_right;
                     ## EDIT: try something like
                     $samples_descr{$condition}->{$replicate} = [$reads_left, $reads_right]; 


      }
    }

    close $fh;

    return(%samples_descr);

}
ADD COMMENTlink modified 3 months ago • written 3 months ago by Michael Dondrup45k

i run it but it gives wrong mapping

#################################################################
Running assemply  on sample repA2
#################################################################
CMD: spades.py --rna   -k 25 -t 6 -1 reads2.left.fq.gz     reads2.right.fq -2  -m 10G -o spades/repA2 -1 reads2.left.fq.gz     reads2.right.fq 


== Error ==  Please specify option (e.g. -1, -2, -s, etc) for the following paths: reads2.right.fq, 10G, reads2.right.fq


In case you have troubles running SPAdes, you can write to spades.support@cab.spbu.ru
or report an issue on our GitHub repository github.com/ablab/spades
Please provide us with params.txt and spades.log files from the output directory.
Error, CMD: spades.py --rna   -k 25 -t 6 -1 reads2.left.fq.gz     reads2.right.fq -2  -m 10G -o spades/repA2 -1 reads2.left.fq.gz     reads2.right.fq  died with ret 256 at ./spades.pl line 147.
ADD REPLYlink modified 3 months ago • written 3 months ago by Mostafa Mahmoud0
1

You are obviously building your command line wrong. Are you sure that $left_fq_file and $right_fq_file are being assigned correctly? And that you want to be outputting $$left_fq_file,?

ADD REPLYlink written 3 months ago by swbarnes24.8k
1

I think there are more syntax problems in the script, but OP has not responded to the request to fix them first. I suggest that we do not invest more time into fixing this, until the syntax is corrected.

ADD REPLYlink written 3 months ago by Michael Dondrup45k
1

So, the first fix changed the output? Please invest a little more effort in proper feedback or we will refuse to assist you further.

ADD REPLYlink written 3 months ago by Michael Dondrup45k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 699 users visited in the last hour