script to calculate A,T,G,C frequency for each position in an alingment

[pony2001mx]

Dear all,
I am new in programming. I wrote a script to calculate the A,T,G,C frequency for each position in an alignment. Take the following alignment for example, in the position 1 the A frequency is 33.3% and G freqeuncy is 66.6%.

ATGCATGC
GGTCACGT
GTAGCGTA

This script is based on others' suggestions, using hash of hash. However, I am in trouble to debug it. Would you please help to indicate where are my errors? THANK YOU VERY MUCH.

#!/usr/bin/perl
use strict;
use warnings;

open SEQ, $ARGV[0];
my %hash = 0;
while (<SEQ>){
chomp;
my @seq = split (/ /, $_);
	for (my $i=0; $i<=$#seq; $i++){
	$hash {$i}{$seq[$i]} +=1;
	}
	}
close SEQ;
	foreach my $posi (keys %hash){
		foreach my $nt ('A', 'T', 'G', 'C'){
	print "$posi", "$nt", %hash{$posi}{$nt}/4, "\n";
		}
	}

[gringer]

Okay, I'll just standardise your code a little to make it a bit easier for me to understand:

#!/usr/bin/perl
use strict;
use warnings;

my %hash = (); # initialise hash
while (<>){ # read a line from standard input (or files specified on command-line)
  chomp; # delete end of line character (if any)
  my @seq = split (/ /, $_); # split line into array, delimiting by a single space
  for (my $i=0; $i<=$#seq; $i++){ # for each index of the array
    $hash {$i}{$seq[$i]} +=1; # increment the value in the hash array at position (i, <base>)
  }
}

foreach my $posi (keys %hash){ # for the keys in the first hash element (indexes)
  foreach my $nt ('A', 'T', 'G', 'C'){ # for the bases A/T/G/C
    print "$posi", "$nt", %hash{$posi}{$nt}/4, "\n"; # print the location, then the base, then the value at (<loc>,<base>) divided by 4
  }
}

Given that you've specified use warnings and use strict, Perl should tell you what you've done wrong syntax-wise. What errors are you getting when you run the code? Alternatively, why do you think it's not working?

I see three potential issues:

• You're dividing by 4 rather than the number of lines read
• The hash value does not necessarily exist when printing
• You're using %hash{a}{b} in the print, rather than $hash{a}{b} [Perl should warn you about this and suggest the correct use]

Here's my suggested modification based on theory (i.e. untested):

#!/usr/bin/perl
use strict;
use warnings;

my %hash = (); # initialise hash
my $lineCount = 0; # initialise line count
while (<>){ # read a line from standard input (or files specified on command-line)
  $lineCount++; # increment lines read
  chomp; # delete end of line character (if any)
  my @seq = split (/ /, $_); # split line into array, delimiting by a single space
  for (my $i=0; $i<=$#seq; $i++){ # for each index of the array
    $hash{$i}{$seq[$i]} +=1; # increment the value in the hash array at position (i, <base>)
  }
}

foreach my $posi (sort {$a <=> $b} (keys %hash)){ # sort the keys numerically
  foreach my $nt ('A', 'T', 'G', 'C'){ # for the bases A/T/G/C
    print("$posi", "$nt", ($hash{$posi}{$nt})/($lineCount), "\n") if(defined($hash{$posi}{$nt}));
  }
}

Yes, there are other optimisations that can be done, but your code was good enough in terms of readability and I don't think you need to worry over that for this script.

[pony2001mx]

Dear gringer,

thank you very much for your detailed corrections. I am new in perl programming, especially the perl grammer. My script is based on others' advices, but i had troubles in debugging.

Your script make me better understand perl. However, when adding a "(" before "if", i have not got any output. Would you check it please? THANKS.

my input file is as follows:
ATGCACTGACTGTATGACTG
ATGGTGACTGTGACTGACTG
ATGGACCATGACTGCATGTG
ATCCACTGTGACGTGCAACA

[pony2001mx]

The last line of your script lacks a "(". After adding it, i found no error message, but i did not get any output. Could you check it please? THANKS.

[gringer]

I changed it to 'defined'. If there's still no output, it suggests that the hash isn't being referenced appropriately, in which case you can explicitly iterate through the second-level keys:

foreach my $posi (sort {$a <=> $b} (keys %hash)){ # sort the keys numerically
  foreach my $nt (keys %{$hash{$posi}}){ # for the bases A/T/G/C
    print("$posi", "$nt", ($hash{$posi}{$nt})/($lineCount), "\n");
  }
}

If still no output, the hash probably isn't being created at all, which is a little odd....

[pony2001mx]

Hi, gringer, there still exists problems. My input file is:

ATGCACTGACTGTATGACTG
ATGGTGACTGTGACTGACTG
ATGGACCATGACTGCATGTG
ATCCACTGTGACGTGCAACA

my output file is:
0.25CACTGACTGTATGACTG
0.25GACCATGACTGCATGTG
0.25CACTGTGACGTGCAACA
0.25GTGACTGTGACTGACTG

When you have time please check again. Anyway i thank you for your help. I really learned somthing from your script and believe it is very very close. I need to stick in perl and continue learning something about hash in perl. THANKS!

[landry]

Hi, gringer, there still exists problems. My input file is:

ATGCACTGACTGTATGACTG
ATGGTGACTGTGACTGACTG
ATGGACCATGACTGCATGTG
ATCCACTGTGACGTGCAACA

my output file is:
0.25CACTGACTGTATGACTG
0.25GACCATGACTGCATGTG
0.25CACTGTGACGTGCAACA
0.25GTGACTGTGACTGACTG

When you have time please check again. Anyway i thank you for your help. I really learned somthing from your script and believe it is very very close. I need to stick in perl and continue learning something about hash in perl. THANKS!

From your output, I'm pretty sure the problem is in this statement:

my @seq = split (//, $_); #DO NOT put any space in between "//" 

[gringer]

Yes, sorry, I forgot about the space separator. That demonstrates why example input/output and a "this is wrong because..." explanation makes the whole problem solving process much easier.

[pony2001mx]

Thank you both. This script works well. Personally I will continue learning perl, which is powerful for genomic analysis. THANKS.