Perl script to count positive values doesn´t report all positive values

user_13978463 0
Viewed 90

I need some help with a Perl script. This script tells me the number of values >= 0.5 for parts. Here is the script:

use strict;
use warnings;

my $file = "229E_O.csv";
my @filearray = ();
my @array_ids = ();
my $thres = 0.5;

open (F, $file) or die;
while(my $l = <F>) {
    $l =~ s/\n//g;
    $l =~ s/\r//g;
    my @cols = split(/\s+/, $l); # divide columns for mora than one space
    next unless (scalar (@cols) == 8); ### If there aren´t 8 column, don´t add to array
    push @filearray, $l;
    my $current_id = $cols[0];
    push @array_ids, $current_id;
}

close F;
my @nr_array_ids = uniq(@array_ids);
foreach my $new_id (@nr_array_ids) { ### for each ID not redundant
    my $counter = 0;
    my $total = 0;
    foreach my $new_L (@filearray) { ### for each line in the file
        my @n_cols = split(/\s+/, $new_L);
        my $potential = $n_cols[5];
        my $idd = $n_cols[0];
        if ($new_id eq $idd) {
            ++$total;
        }

        if ( ($new_id eq $idd) and ($potential >= $thres) ) {
            ++$counter;
        }
    }

    print "$new_id\t$counter\t$total\n";
}
sub uniq {
    my %seen;
    return grep { !$seen{$_}++ } @_;
}

Here the input file:

APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        671     671     0.134197        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        672     672     0.282583        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        676     676     0.290996        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        680     680     0.376348        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        682     682     0.552045        .       .       #POSITIVE
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        688     688     0.315533        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        696     696     0.111705        .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        700     700     0.20703 .       .
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        701     701     0.284842        .       .

The input has 8 columns (0-8), column 5 has the value that is important to me. When this value is >= 0.5 means a #POSITIVE. This script works fine, however, I used for report other value (>=0.7). Now when changing the value the script doesn´t report me when there is a value >=0.5 Here is the output:

APT69890_1_NA   0   197
AFR79257_1_NA   0   198
AGT21345_1_NA   0   200
QJY77970_1_NA   0   199
QJY77962_1_NA   0   200
QEO75985_1_NA   0   199
ARK08620_1_NA   0   202

If you can see each ID, for example, APT69890_1_NA is a "part" of the output. The second column is #POSITIVE and the third is all values <0.5 Here the part APT69890_1_NA must be 1 in the second column, however, the value is 0. Here is the complete example of my real data if you want to see it: https://github.com/MauriAndresMU1313/Example_NetOGlyc/tree/main

perl
2 Answers

This line of code expects only 8 columns in each line of data or the line is ignored:

next unless (scalar (@cols) == 8); ### If there aren´t 8 column, don´t add to array

But in your data file, the one line you wanted counted (line 253 of data - 0.552045) - an additional column was added "#POSITIVE" making 9 columns.

APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        682     682     0.552045        .       .       #POSITIVE

That line is therefore rejected. Your tally of 197 is one less than the total number of entries in the data file for APT69890_1_NA. Further evidence this is a data error.

Remove that 9th column or change conditional to tolerate a 9th column.

If you wanted to document your data - then in this case you could simply remove the #POSITIVE (9th column) and move it to preceding line. It would then be ignored since it contained only one column:

#POSITIVE
APT69890_1_NA   netOGlyc-4.0.0.13       CARBOHYD        682     682     0.552045        .       .

Alternatively you could change the conditional to tolerate at least 8 columns. This has the advantage of leaving the data as is but weakens the validity check of the data.:

next unless (scalar (@cols) >= 8); ### If there are less than 8 columns, don´t add to array
Admin1

Please investigate the following script for compliance with your task.

NOTE: the script assumes consistent datafile as your provided on github

  • define $threshold
  • allocate hash to hold %data
  • read file line by line in while( <> ) loop
  • attempt to split line and take $id and $value
  • skip to next line if $id is not defined
  • increase counter for this particular $id
  • initialize positive counter to 0 if it is not defined
  • increment positive counter if $val >= $threshold
  • when all data read, output result in for loop

NOTE: run as ./script.pl datafile

use strict;
use warnings;
use feature 'say';

my $threshold = 0.5;
my %data;

while( <> ) {
    my($id,$val) = (split)[0,5];
    next unless $id;
    $data{$id}{count}++;
    $data{$id}{pos} //= 0;
    $data{$id}{pos}++ if $val >= $threshold;
}


say join("\t",$_,$data{$_}->@{qw/pos count/}) for sort keys %data;
user_123133090
Related