Data screening using Perl

Question

Background information

I've been asked to write a little Perl script that allows genomic data to be screened against reference files in order to determine locations of specific mutations.

The input file format look like this (information are tab-separated). In the script it's referred to as Funestus_NON_SYN.txt.

KB669306    9962    .   C   T   520.77  PASS    AC=13;AF=0.929;AN=14;BaseQRankSum=0.540;DP=103;Dels=0.00;EFF=SYNONYMOUS_CODING(LOW|SILENT|cgC/cgT|R15||AFUN012013|||AFUN012013-RA|1|1);FS=0.000;HaplotypeScore=0.0000;MQ0=0;MQRankSum=0.231;ReadPosRankSum=0.694;set=variant-variant2-variant3-variant5-variant6-variant7-variant10 GT:AD:DP:GQ:PL  1/1:0,18:18:42:549,42,0 1/1:0,12:12:33:399,33,0 1/1:0,13:13:39:505,39,0 ./. 0/1:9,4:13:99:104,0,297 1/1:0,17:17:45:570,45,0 1/1:0,18:18:51:624,51,0 ./. ./. 1/1:0,12:12:33:420,33,0
KB669306    10439   .   C   T   668.77  PASS    AC=15;AF=0.938;AN=16;BaseQRankSum=-0.643;DP=139;Dels=0.00;EFF=SYNONYMOUS_CODING(LOW|SILENT|gcC/gcT|A174||AFUN012013|||AFUN012013-RA|1|1);FS=0.000;MQ0=0;MQRankSum=-0.746;ReadPosRankSum=-1.106;set=variant-variant3-variant5-variant6-variant7-variant8-variant9-variant10  GT:AD:DP:GQ:PL  1/1:0,22:22:51:697,51,0 ./. 1/1:0,20:20:48:667,48,0 ./. 0/1:14,12:25:99:317,0,403   1/1:0,15:15:36:469,36,0 1/1:0,16:16:36:499,36,0 1/1:0,13:13:30:415,30,0 1/1:0,14:14:30:412,30,0 1/1:0,13:13:36:492,36,0
KB668289    903 .   T   C   577.77  PASS    AC=14;AF=0.875;AN=16;DP=173;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3412|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|ctA/ctG|L578||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2667|||AFUN005869|||AFUN005869-RA||1);MQ0=0;set=variant-variant2-variant3-variant5-variant6-variant7-variant9-variant10  GT:AD:DP:GQ:PL  1/1:0,20:20:45:606,45,0 1/1:0,19:19:42:588,42,0 1/1:0,22:22:54:720,54,0 ./. 1/1:0,29:29:63:882,63,0 0/1:16,13:28:99:401,0,391   1/1:0,13:13:30:425,30,0./.  1/1:0,26:26:57:785,57,0 0/1:8,7:15:99:223,0,243
KB668289    1224    .   A   C   572.77  PASS    AC=10;AF=0.833;AN=12;DP=122;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3091|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|cgT/cgG|R471||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2346|||AFUN005869|||AFUN005869-RA||1);HaplotypeScore=0.0000;MQ0=0;set=variant-variant3-variant6-variant7-variant8-variant9   GT:AD:DP:GQ:PL  1/1:1,20:20:42:601,42,0 ./. 1/1:0,26:26:63:865,63,0 ./. ./. 0/1:9,10:19:99:276,0,236    1/1:0,17:17:42:583,42,0 0/1:14,5:19:99:119,0,406    1/1:0,20:20:45:630,45,0 ./.

Now, I have to screen the first column against files containing genes in the following format. Referred to as gsts-funestus.gff in the script

# start gene FUNEGST002
KB668289    AUGUSTUS    gene    410926  411627  0.18    +   .   FUNEGST002; gene_name "GSTd3"; auto "g17128_modified";
KB668289    AUGUSTUS    transcript  410926  411627  0.18    +   .   FUNEGST002.t1
KB668289    AUGUSTUS    start_codon 410926  410928  .   +   0   transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289    AUGUSTUS    intron  411058  411126  0.99    +   .   transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289    AUGUSTUS    CDS 410926  411057  0.58    +   0   transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289    AUGUSTUS    CDS 411127  411627  0.74    +   0   transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289    AUGUSTUS    stop_codon  411625  411627  .   +   0   transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
# end gene FUNEGST002

Basically I have to extract the first line containing gene in the 3rd column and print a hybrid output in the following format.

GSTd3   KB668289    903 T=>C    AC=14;AF=0.875;AN=16;DP=173;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3412|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|ctA/ctG|L578||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2667|||AFUN005869|||AFUN005869-RA||1);MQ0=0;set=variant-variant2-variant3-variant5-variant6-variant7-variant9-variant10  GT:AD:DP:GQ:PL  1/1:0,20:20:45:606,45,0 1/1:0,19:19:42:588,42,0 1/1:0,22:22:54:720,54,0 ./. 1/1:0,29:29:63:882,63,0 0/1:16,13:28:99:401,0,391   1/1:0,13:13:30:425,30,0./.  1/1:0,26:26:57:785,57,0 0/1:8,7:15:99:223,0,243
GSTd3   KB668289    1224    A=>C    AC=10;AF=0.833;AN=12;DP=122;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3091|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|cgT/cgG|R471||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2346|||AFUN005869|||AFUN005869-RA||1);HaplotypeScore=0.0000;MQ0=0;set=variant-variant3-variant6-variant7-variant8-variant9   GT:AD:DP:GQ:PL  1/1:1,20:20:42:601,42,0 ./. 1/1:0,26:26:63:865,63,0 ./. ./. 0/1:9,10:19:99:276,0,236    1/1:0,17:17:42:583,42,0 0/1:14,5:19:99:119,0,406    1/1:0,20:20:45:630,45,0 ./.

The Script

My script contains this part to read and process all lines in the input file.

#!/usr/bin/perl -w

use strict;
use warnings;

print "Enter your input file containing all SNPs [e.g. Funestus_NON_SYN.txt]:\t";
chomp (my $input_file = <STDIN>);

# Open and import the information found in the files #
my @input = &InputFileHandler($input_file);

# Store the gene information in a hash with the key as scaffold and location #
my %genes = ();
foreach (@input)
{
    next if $_ !~ /^KB.*/;
    my($scaffold,$location,$dot,$start,$end,$unknown,$pass,$rest) = split(/\t/, $_, 8);
    my $mutation = "$start=>$end";
    my $key = $scaffold."_".$location;
    my $value = "$scaffold\t$location\t$mutation\t$rest";
    $genes{$key} = $value;
}

The next part allows the user to define the number of reference files and screens out input file against it. Furthermore, it prints the results in a tab-separated output file.

# Open the output file
my $output_file = &GenerateFilename();
open my $fh_out,">",$output_file.".tsv" or die "Cannot open $output_file: $!\n";

# Print the header in the CSV file
print $fh_out "Gene\tSNP\tLocation\tMutation\tRemaining information\n###\n";

print "Enter the number of reference files [e.g. gff files]:\t";
chomp (my $i = <STDIN>);

# Generic loop that let's the user enter as many reference files as required #
for (my $j = 1; $j <= $i;$j++)
{
    print "Enter filename $j [e.g. gsts-funestus.gff]:\t";
    chomp (my $file = <STDIN>);
    my @file_raw = &InputFileHandler($file);
    my @file_genes = &ScanArray(@file_raw);
    my @file_genes_print = &ScanPrint(@file_genes);
    print $fh_out "$_\n" foreach (@file_genes_print);
    print $fh_out "\n###\n";
}

Finally, my subroutines are the following.

sub GenerateFilename    # Generates a filename based on the current date and time.
{
    my @timeData = localtime(time);
    my $filename = join('_', @timeData);
    return $filename;
}

sub InputFileHandler    # Opens a file and saves its content in an array
{
    my $file = shift;
    my @lines;
    open my $fh,"<",$file or die "Cannot open $file: $!\n";
    while (my $line = <$fh>)
    {
        chomp($line);
        push (@lines,$line);
    }
    close ($fh);
    return @lines;
}

sub Information         # Separates the mutation from the remaining information
{
    my $rest = shift;
    my @rest = split (/\t/,$rest);
    my $aa_change = "$rest[1]=>$rest[2]";
    my $return = "$aa_change,$rest[5]";
    return $return;
}

sub ScanArray           # Scans a file array for lines containing SNPs, discards rest
{
    my @array = @_;
    my @return;
    foreach (@array)
    {
        next if $_ !~ /^KB.*\t.*\t[Gg][Ee][Nn][Ee]/;
        my ($scaffold,$second,$type,$start,$end,$sixth,$seventh,$eigth,$info) = split (/\t/,$_);

        my ($unknown,$gene_name,$auto) = split (/; /,$info,3);
        $gene_name =~ s/.*"(.*?)".*/$1/s;
        my $gene = $scaffold."_".$gene_name;
        push (@return,$gene);
    }
    return @return;
}

sub ScanPrint           # Filters for genes that contain SNPs
{
    my @array = @_;
    my @return;
    foreach (@array)
    {
        my ($scaffold,$gene_name) = split (/_/,$_,2);
        foreach my $key (keys %genes)
        {
            next if $key !~ /$scaffold/;
            my $line = "$gene_name\t$genes{$key}";
            push (@return, $line);
        }
    }
    return @return;
}

Question

Do you people have any suggestions or ways of improving this script. It all works and no warnings have been returned with all search queries I've done (I used different files of the same format).

EDIT 1

As pointed out to me by @Edward: Yes, my output should be the two lines. I have adjusted it in the main part of the question.

Answer 1

Style Guides

The Perl community does not have a single universal style – there is more than one way to do it. But sometimes, consistency isn't a bad thing either. Here are the three important cornerstones of Perl style:

• A perlstyle manpage exists which explains a sensible core of a style guide.
• Damian Conway published the Perl Best Practices book with many valuable tips, although I'd view his preferences as over-zealous.
• The perlcritic tool exists to automatically check your code for a set of style violations. It's inspired by Perl Best Practices. You can install it via the Perl::Critic module. You can layout your code automatically with perltidy, installable via Perl::Tidy.

Based on those references, I make the recommendations below.

Objective Style Issues

Now lets look at what you've written.

• Don't use the -w command line flag, only use warnings, and use it always. They are fundamentally the same (they activate warnings), but -w has an unfortunate global effect. In your case it doesn't matter, but -w is a bad habit to get into.

• Don't call functions with an ampersand, like &foo(). This has a specific meaning (it disables so-called prototypes), but it hasn't been necessary to use it for over 20 years. Simply invoke your functions like foo().

• When looping over a range, do not use the C-style for (my $j = 1; $j <= $i; $j++). Instead, use the range operator, and a foreach loop:

  for my $j (1 .. $i) {
  

• Since perl 5.10.0 (released 2007), you can use feature 'say' to enable the say function. It behaves exactly like print but always appends a newline, making it more convenient to use in most cases.

• The keywords for and foreach are absolutely equivalent, so we'd use the shorter form. When used as a statement modifier, the parens around the iterated values are optional:

  ... for @file_genes_print;
  

• In a foreach loop, always name the iteration variable rather than using $_. This is not possible in the statement modifier use case.

  for my $iteration_variable (@list) {
      ...
  }
  

• Hash variables are initialized by declaration. Assigning them the empty list is unnecessary – my %hash; rather than my %hash = ();.

• For readability, any comma should be followed by a space.

  open my $fh_out, ">", "$output_file.tsv" or die ...
  

• Prefer interpolation over concatenation. $scaffold."_".$location could be "${scaffold}_${location}" and $output_file.".tsv" could be "$output_file.tsv".

• Sometimes, join is more preferable than both interpolation or concatenation. For example, "$scaffold\t$location\t$mutation\t$rest" becomes more readable as

  join "\t", $scaffold, $location, $mutation, $rest
  

• The print (and say) syntax is horrible because of the difference between print $foo $bar and print $foo, $bar – the former prints $bar to the filehandle $foo, whereas the latter prints "$foo$bar" to the currently selected filehandle (usually STDOUT).

  • We could use the object-oriented interface.

    use IO::File;
    $foo->print($bar);
    

  • We can use curly braces for the indirect object notation

    print { $foo } $bar
    

  I'd recommend to always use the curly braces. You only really need them when the file handle is something else than a simple variable, but it's a good visual distinction.

• The ScanPrint uses the %genes hash which is therefore essentially a global variable. Instead, pass it in as an argument:

  scan_print(\%genes, @file_genes);
  
  ...
  
  sub scan_print {
      my ($genes, @array) = @_;
      ...
      for my $key (keys %$genes) {
  

• Don't load a whole file into memory unless you have to. If you can, process files line-by-line. Especially your use of InputFileHandler violates this. On the other hand this may be absoluetly OK for smaller files, or when you are going to use that amount of memory anyway.

Subjective Style Issues

• Variable and function names should generally be snake_case: lowercase, and words separated by underscores. You apply this for variables, but not for subroutine names.

• Consider using the K&R style/egyptian brackets. This means that the opening curly brace is on the same line as it's introducing keyword and not on a line of its own:

  for (...) {
      ...
  }
  

  Why? For no other reason than that most Perl programmers do this.

• It is customary to call builtins without parentheses. This is partly because builtins are actually operators, not functions, and partly because it looks nicer. So let's do close $fh rather than close ($fh).

• In regexes, put symbols into character classes to emphasize them: /; / might become /[;][ ]/.

Performance Issues, Bugs, and Other Notes

• You are using a regex /^KB.*/. The trailing .* is entirely unnecessary as this always matches any string. Remove it: /^KB/.

• Do not prompt for file names. Instead, take these from the command line arguments. The arguments are in the @ARGV array. So you could do:

  my ($input_file, @reference_files) = @ARGV;
  

  and invoke your script like perl my_script.pl input_file.tsv ref1 ref2.

• $key !~ /$scaffold/ interprets $scaffold as a regex. To test that $key does not contain the $scaffold, use the more efficient -1 == index $key, $scaffold. To test that they are not equal, use $key ne $scaffold, in which case you could reduce the surrounding loop to

  if (exists $genes{$scaffold}) {
      push @return, "$gene_name\t$genes{$scaffold}";
  }
  

• The substitution s/.*"(.*?)".*/$1/s is iffy. If you want to match a quoted string, "[^"]*" is simpler for the regex engine to handle than ".*?". Of course, it also doesn't support escapes.

  If you want to replace the whole string by the thing in the quotes, just use a normal match rather than a substitution, e.g:

  if ($gene_name =~ /["]([^"]*)["]/) {
      $gene_name = $1;
  }
  

  which is more straightforward.

• Creating a time-based filename is fragile and unintuitive for a user. Instead, let the user specify the output file name manually. Alternatively, generate a more unique filename (e.g. by using the process ID $$ in the filename).

  Another approach, which I have taken below, is to print to standard output, and let the user write the results to a file using shell redirection.

• The Information subroutine isn't used.

Suggested Rewrite

This rewrite doesn't use any prompting (instead it uses command line arguments and shell redirection), and it doesn't use any subroutines. The result should be shorter, simpler, and more efficient.

This code has not been tested.

#!/usr/bin/perl

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

# we take the input file as STDIN
# and the output file as STDOUT
# the user can redirect this to desired files:
# $ perl my_script.pl ref1 ref2 <input >output

my @reference_files = @ARGV;

# early input validation
if (not @reference_files) {
    die "You have to specify at least one reference file as argument\n";
}
for my $file (@reference_files) {
    die "The reference file '$file' does not exist\n" if not -e $file;
}

# parse the input file
my %genes;
while (my $line = <STDIN>) {
    chomp $line;
    next if not $line =~ /^KB/;
    my ($scaffold, $location, $dot, $start, $end, $unknown, $pass, $rest) = split /\t/, $line, 8;
    $genes{$scaffold}{$location} = join "\t", $scaffold, $location, "$start=>$end", $rest;
}

# print the header
say join "\t", 'Gene', 'SNP', 'Location', 'Mutation', 'Remaining Information';
say "###";

# go through each reference file
for my $file (@reference_files) {
    open my $ref, "<", $file or die "Can't open reference file '$file': $!";
    while (my $line = <$ref>) {
        chomp $line;
        my ($scaffold, undef, $type, $start, $end, undef, undef, undef, $info) = split /\t/, $line;
        next if not $scaffold =~ /^KB/;
        next if not $type =~ /^GENE/i;

        my ($transcript_id, $gene_name, $auto) = split /[;][ ]/, $info;
        $gene_name = $1 if $gene_name =~ /["]([^"]*)["]/;

        if (my matching_genes = $genes{$scaffold}) {
            say join "\t", $gene_name, $_ for values %$matching_genes;
        }
    }

    say "";
    say "###";
}

Answer 2

Scan file directly into data structures

You have a routine InputFileHandler that puts every line of an input file into an array. That array is then subsequently dissected into either the %genes hash or @file_genes array. What would be more Perl-ish would be to avoid the intermediate step and scan directly into the desired data structures. That would eliminate InputFileHandler and modify and rename ScanArray.

Remove unused Information subroutine

It's possible that it's needed for some other processing that you intend to do, however. If that's the case, I'd suggest renaming it to be a little more suggestive of its purpose. The comment suggests that a better name might be ExtractMutation.

Use undef rather than dummy variable names for split

When you're splitting variables, use undef instead of dummy variable names. It will be clearer to those who read the code which fields are not of interest and your program will run faster and with less memory as a result.

Prefer command line arguments to interactive I/O

It's easier and more clear if you accept command-line arguments as filenames rather than prompting for them. Since you have one input file and then one or more gff files, you can treat the first argument as the txt file and all subsequent args as gff files.

GenerateFilename might be improved

The current GenerateFilename function generates filenames from timestamps, which might not be too bad, but it does so in a way that will make it difficult to parse the resulting files because the seconds are first, followed by minutes, hours, etc. Better would be to generate file names in a way that they can be rationally sorted (e.g. with year, month, day first) or to omit generation of output filenames entirely and assume that the user will redirect output to an appropriately named file. This is what I would choose, but your application may dictate some other choice, which is why I left it in place in the example code below.

Consolidate reading and output

Once the %genes memory structure is populated, there seems to be little use in storing the subsequent gff file contents in memory. Better would be to process each gff line as received and correlate it to the matching %genes for immediate output. This will save a lot of runtime memory and be faster.

Results:

With all of that done, my rewrite looks like this:

#!/usr/bin/perl 

use strict;
use warnings;

my $input_file = shift;

# Store the gene information in a hash with the key as scaffold and location #
my %genes;
open my $fh,"<",$input_file or die "Cannot open $input_file: $!\n";
while (<$fh>)
{
    next if $_ !~ /^KB.*/;
    chomp;
    my($scaffold,$location,undef,$start,$end,undef,undef,$rest) = split(/\t/, $_, 8);
    my $mutation = "$start=>$end";
    my $key = "$scaffold\_$location";
    my $value = "$scaffold\t$location\t$mutation\t$rest";
    $genes{$key} = $value;
}
close $fh;

# Open the output file
my $output_file = &GenerateFilename();
open my $fh_out,">",$output_file.".tsv" or die "Cannot open $output_file: $!\n";

# Print the header in the CSV file
print $fh_out "Gene\tSNP\tLocation\tMutation\tRemaining information\n###\n";

# Generic loop that lets the user enter as many reference files as required
while (my $file = shift) {
    &ReadSNPFile($fh_out, $file);
}
print $fh_out "\n###\n";

sub GenerateFilename    # Generates a filename based on the current date and time.
{
    my @timeData = localtime(time);
    my $filename = join('_', @timeData);
    return $filename;
}

# passed an output file name, and an SNP filename, 
# this scans the file for lines containing SNPs
# and emits matches.
#
sub ReadSNPFile           
{
    my $out_file = shift;
    my $file = shift;
    open my $fh,"<",$file or die "Cannot open $file: $!\n";
    while (<$fh>)
    {
        next if $_ !~ /^KB.*\t.*\t[Gg][Ee][Nn][Ee]/;
        my ($scaffold,undef,undef,undef,undef,undef,undef,undef,$info) = split (/\t/,$_);
        my ($unknown,$gene_name,$auto) = split (/; /,$info,3);
        $gene_name =~ s/.*"(.*?)".*/$1/s;
        foreach my $key (keys %genes)
        {
            next if $key !~ /$scaffold/;
            print $out_file "$gene_name\t$genes{$key}\n";
        }
    }
    close ($fh);
}

I've tested it and it produces results identical to your original but in about half the number of lines of code.