G'day ZWcarp,

Firstly, you're doing yourself no favours whatsoever by littering your code with meaningless variable names. While you may remember, while you're writing it and it's fresh in your mind, what $r[3] or $key2 refer to, you won't next month when you have to come back to it to make a modification.

From the code and text you've provided: the keys of %AA are gene names (if I've got that right, $gene_name would be a meaningful name); the keys of %{$site{$gene_name}} are mutation sites (again $mutation_site would be meaningful); and so on throughout your code.

I don't see any purpose to any of the sorting you're doing in your "problem section" (it's wasted processing and chews up even more memory) and I agree with AM about the [$key4].

Putting all that together, I think your "problem section" could be written as (untested):

for my $gene_name (keys %AA) {
    for (keys %{$site{$gene_name}}) {
        $site_length_catch{$gene_name}{$_} = @{$site{$gene_name}{$_}};
    }
}
[download]

And later you can access that data as:

my $mutation_count = $site_length_catch{$gene_name}{$mutation_site};
[download]

There's other parts of your code that seem dubious (e.g. the my $key4=$key3; assignment) which perhaps will become obvious to you when you apply better names. You're working with very large amounts of data and loops nested three-deep: you need to keep all the code (but especially the innermost loops) as efficient as possible: go through your code and remove unnecessary assignments, sorting and so on.

Do you really need to keep a list of $key4 ? I think you could just count them, $key4 as both a key and a value looks at bit odd, and you only seem to be counting them by keeping an array of repeated values. So, you could get rid of the useless arrays and replace them with a simple integer.

#so replace this    
    push (@{$site{$key1}{$key4}},$key4);
#with 
    $site{$key1}{$key4}++;
[download]

BTW In the code you posted as the problem area the sorts are doing nothing useful, just taking up time, but maybe you're doing other stuff in there too.

$site_length_catch{$key1}[$key4]

Did you mean for that to be an array ref instead of a hash ref for $key4? Also, what values does $key4 take?

If $key4 happens to be a large number, then you'll blow away all your memory in one step as the array attempts to grow to ludicrous size.

$key4 = the amino acid position. It is derived from the first file and edited to remove non digits so that its just a numerical position

cell $key4 is originally initialized as 0 and at this step it needs to be changed to the number of mutations at this site (the array length corresponds to the gene length of gene $key1 so each cell corresponds to an amino acid position in the gene). It is set equal to scalar ( @{$site{$key1}{$key4}} ); because this calculates the number of values at a mutation site, thus the size of this essentially opperates as a count for recurrent positions across different samples. The numbers for this are mostly 0-10 there are a few large like one 35K and a few below that.

Not sure if I fully understand but I'd like to venture a guess what you want:

use strict;
use warnings;

my %site_length_catch;
my $max = 0;
foreach (@file) {
  chomp;
  my @r = split /\t/;
  # cleaning from your second loop
  $r[13] =~ s/\D\.\D([0-9]+)\D/$1/; 
  $r[13] =~ s/(\*|\?|s\d+)//; 
  $site_length_catch{$r[0]}{$r[13]}++;
  $max = $r[13]>$max?$r[13]:$max;
}

foreach my $gene (keys %site_length_catch) {
  print $site_length_catch{$gene}{$_} // 0, "\t" for 1..$max;
  print "\n";
}
[download]

The hash %site_length_catch is a sparse matrix containing the name of the gene as the first dimension and the site of mutations as the second dimension. Each cell in the matrix contains the number of mutations at that site for that gene.

When printing the empty spaces are filled with zeros (this is what  // 0 does). I have added the regexes from your second loop as they seem to be applied to the "Mutation site". Just remove them if I have guessed wrongly.

Feedback would be appreciated, along with a few lines of your input, if possible.

Thats is extremely close to what I want, thank you so much for your time. The only difference there is another level which is patient ID ..... if two mutations are the same site but they are in the same patient its not important. What is interesting is when 2 mutations are in the same site/gene of two different patients. Pool enough of these and you can identify a drug target to treat the cancer. If the site and the Amino acid change .. ( note L64F vs L64R would be different ) then that is an even better target. Here is the input lines of the data

The file I need to read in is this structure. These 19 lines are taken from the real file which has 5,772,080 lines in the same format :

Gene Name    Patient ID    Patient Diagnosis    Ammino Acid Mutation a
+nd Sit    Protein Length 
AAK1    19679    adenocarcinoma    L661I    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19676    adenocarcinoma    L664T    21265
AAK1    19677    adenocarcinoma    L64F    21265
AAK1    19678    adenocarcinoma    L64R    21265
FKT1    101063    ER-PR-sitive_carcinoma    p.L52R    2773
FKT1    103872    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    107590    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    107600    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    1135911    NS    E17K    2773
TET3    152    chronic_lymocytic_leukaemia    p.R401H    10982
TET3    587220    adenocarcinoma    M935V    10982
TET3    587220    adenocarcinoma    R1534Q    10982
TET3    587256    adenocarcinoma    G1356R    10982
TET3    587338    adenocarcinoma    G1356W    10982
[download]

Now I need to count all positions that match in Amino Acid Site (the number but not the letters of the 4th column) but are in different samples. Note : Patient ID19679 and AA mutation L664T only corresponds to a count of 2 because all of them are in the same patient except one in patient 19676.

The out put needs to be in this format, where you have rows as genes and columns are 1-Length(the fifth column above). L is different for every gene. I've listed spans as no1.....no2 just for sake of space, but in the real file all these numbers in between have to be filled with 0's:

1-Largest Gene Length    AA site -1    AA site -2    AA site -3    4……
+……16    AA site -17    18…..51etc    AA site 52    AA site 64    65…4
+00    AA site 401    402….660    AA site 661    AA site 664    AA sit
+e 935    AA site 1356    AA site 1534
AAK1    0    0    0    0    0    0    0    2        0        1    2   
+ 0    0    0
FKT1    0    0    0    0    4    0    1    0        0        0    0   
+ 0    0    0
TET3    0    0    0    0    0    0    0    0        1        0    0   
+ 1    2    1
[download]

I'm simplifying because I also need to calculate a second table but this time with Amino acid position and mutation ( thus numbers and letters of Column 4) matching in different patients. Thats why my script is so elaborate, the $key3=$key4 is to remove the letters etc. I know i've done a poor job scripting it. Any advice would be fantastic!! Thanks so much for helping.

Here is the next iteration. This now includes some code to avoid double counting of the same site of mutation for the same patient (looks very similar to your original code...). You can chose between printing the full matrix or only the relevant ones by uncommenting one of two lines close to the end of the code. Hope this is helpful.

use strict;
use warnings;

my %site_length_catch;
my %sites;
my $maxsite = 0;

<DATA>; # skip header
foreach (<DATA>) {
  chomp;

    # split and give meaningful names
  my( $gene, $patient, $diagnosis, $mut_and_sit, $length ) = split /\s
++/;

    # clean the site
    my $sit = $mut_and_sit;
    $sit =~ s/\D//g;

    # store patient to avoid double counting
  $site_length_catch{$gene}{$sit}{$patient} = 1; 

    # store all sites with mutations
    $sites{$sit} = 1;
    $maxsite = $sit > $maxsite ? $sit : $maxsite;
}

# now remove double counted patients from the data structure
foreach my $gene ( values %site_length_catch) {
        for my $count ( values %$gene ) {
                $count = keys %$count; # in scalar context you get the
+ number of keys
        }
}

# print table in desired format
# uncomment one of the following two lines
my @sitesprinted = sort { $a <=> $b } keys %sites;   # sparse printing
#my @sitesprinted = 1..$maxsite;                      # full printing

# header first
print "Gene";
print "\tsite $_" for @sitesprinted;
print "\n";

# now the data
foreach my $gene (keys %site_length_catch) {
        print $gene;
        print "\t", $site_length_catch{$gene}{$_} // 0 for @sitesprint
+ed; 
        print "\n";
}

__DATA__
Gene Name    Patient ID    Patient Diagnosis    Ammino Acid Mutation a
+nd Sit    Protein Length 
AAK1    19679    adenocarcinoma    L661I    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19679    adenocarcinoma    L664T    21265
AAK1    19676    adenocarcinoma    L664T    21265
AAK1    19677    adenocarcinoma    L64F    21265
AAK1    19678    adenocarcinoma    L64R    21265
FKT1    101063    ER-PR-sitive_carcinoma    p.L52R    2773
FKT1    103872    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    107590    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    107600    ER-PR-sitive_carcinoma    p.E17K    2773
FKT1    1135911    NS    E17K    2773
TET3    152    chronic_lymocytic_leukaemia    p.R401H    10982
TET3    587220    adenocarcinoma    M935V    10982
TET3    587220    adenocarcinoma    R1534Q    10982
TET3    587256    adenocarcinoma    G1356R    10982
TET3    587338    adenocarcinoma    G1356W    10982
[download]

As per your second table, I am optimistic that it is a relatively simple modification only.