title: Counting nucleotide frequency using perl script
body: <p>I have this perl script below to calculate sequence length and their frequency along with nucleotide frequency(<code>A</code>,<code>T</code>,<code>G</code> and <code>C</code>). This script works fine for a file with large number of sequences, but it does not give the right result for a file of small size like this: </p>
  
  <p><code>infile.fasta</code></p>
  
  <pre><code>&gt;NC_013116_1051_1114
  TTGTCCCTTTGAGTCTCTGG
  &gt;NC_013116_1051_1114
  GCGCAGCCGATATGGATAA
  &gt;NC_013116_1051_1114
  TCGAGACTTTGTAATGTTTGGG
  &gt;NC_013116_1051_1114
  TATTCCACGTCAGGTGCTTTT
  &gt;NC_013116_1051_1114
  TAGAGCCGATTCCAGACTGTTCC
  &gt;NC_013116_1051_1114
  TACAGGACCAAGCTCTTCACTC
  &gt;NC_013116_1051_1114
  CCGTCAAGTTCAGCTCCAATAA
  &gt;NC_013116_6_301
  CCACGCAACGGACAATCAAACA
  &gt;NC_013116_6_301
  GGACACTTCCAACTATAAATA
  &gt;NC_013116_6_301
  CCACGCAACGGACAATCAAACA
  &gt;NC_013116_1051_1114
  GCTCTTCACTCTTCCTCGTCT
  &gt;NC_013116_1051_1114
  TTGGGAAAAAGAAGTTGCTGCAGC
  &gt;NC_013116_1051_1114
  TCGCAGTATCTCTGAAGTTG
  </code></pre>
  
  <p>count.pl</p>
  
  <pre><code>#!/usr/bin/perl -w
  
  #usage ./count.pl infile min_length max_length
  #usage ./count.pl infile 18 34
  
  my $min_len = $ARGV[1];
  my $max_len = $ARGV[2];
  my $read_len = 0;
  my @lines = ("header1","sequence","header2","quality");
  my @lray = ();
  my $count = 0;
  my $total = 0;
  my $i = 0;
  
  my @Aray = ();
  my @Cray = ();
  my @Gray = ();
  my @Tray = ();
  
  my$FN = "";
  
  for ($i=$min_len; $i&lt;=$max_len; $i++){
     $lray[$i] = 0;
  }
  
  open (INFILE, "&lt;$ARGV[0]") || die "couldn't open input file!";
     while (&lt;INFILE&gt;) {
        $lines[$count] = $_;
        chomp($lines[$count]);
        $count++;
        if($count eq 4){
           $read_len = length($lines[1]); 
  #         print "$read_len $lines[1]\n";
           $FN = substr $lines[1], 0, 1;  
           $lray[$read_len]++;
           if ($FN eq "T") { $Tray[$read_len]++;} 
           else {        
              if ($FN eq "A"){ $Aray[$read_len]++;}
              else {
                 if ($FN eq "C"){ $Cray[$read_len]++;}
                 else {
                   if ($FN eq "G"){ $Gray[$read_len]++;}
                 }   
              }
           }           
           $count = 0;
        }
     }
  print "length\tnumber\tA\tC\tG\tT\n";
  for ($i=$min_len; $i&lt;=$max_len; $i++){
     print "$i\t$lray[$i]\t$Aray[$i]\t$Cray[$i]\t$Gray[$i]\t$Tray[$i]\n";
  }
  exit;
  </code></pre>
  
  <p>This is the type of result I get from a big file with many sequences. </p>
  
  <pre><code>length  number  A   C   G   T
  18  4473    542 710 471 2750
  19  12647   990 1680    1103    8874
  20  31194   3010    3354    2743    22087
  21  61214   6288    7196    5784    41946
  22  128642  14596   11902   12518   89626
  23  65190   6859    6525    7773    44033
  24  10012   611 1401    1112    6888
  25  1406    231 192 435 548
  26  661 169 91  105 296
  27  407 126 81  65  135
  28  602 391 49  68  94
  29  520 54  30  370 66
  30  175 26  93  18  38
  31  156 35  28  29  64
  32  106 22  16  24  44
  33  97  45  17  16  19
  34  0           
  </code></pre>
  
  <p>I would really appreciate if you could help me correct this code. Thanks    </p>