Perl script review

Question

[Mostafa Mahmoud]

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);

}

[Pierre Lindenbaum]

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

what is the question ?

[Mostafa Mahmoud]

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 ???

[ATpoint]

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

[Michael]

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.

[Mostafa Mahmoud]

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

[Michael]

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.

Answer 1

[Michael]

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);

}

[Mostafa Mahmoud]

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 [e-mail]
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.

[swbarnes2]

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,?

[Michael]

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.

[Michael]

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