0

Please help improve the following code . I am not able to print the sequence in one single line. Would like to have output printed in four lines each with nucleotide frequency of one of the four characters. Thanks in advance.enter code here

#!/usr/bin/perl
use strict;
use warnings;
my $A;    
my $T;
my $G;
my $C;
my $fileIN;
my $fileOUT;

my $seq ;
open ($fileIN, "basecount.nfasta") or die "can't open file ";
open ($fileOUT, ">basecount.out") or die "can't open file ";

while (<$fileIN>)
{

             if ($_ =~/^>/)  #ignore header line
             {next;}

             else
                   {
                    $seq  = $_; #copy the all line with only nucleotide characters ATGC
                   }
            $seq  =~ s/\n//g; #create one single line containing all ATGC characters

             print "$seq\n"; # verify previous step

             my @dna = split ("",$seq); #create an array to include each nucleotide as array element

             foreach my $element (@dna)

            {
            if ($element =~/A/) # match nucleotide pattern and countstrong text
                            {
                             $A++;
                            }
             if ($element =~/T/)
                            {
                             $T++;
                            }
             if ($element =~/G/)
                            {
                             $G++;
                            }
             if ($element =~/C/)
                            {
                             $C++;
                            }

            }

            print $fileOUT "A=$A\n";
            print $fileOUT "T=$T\n";
            print $fileOUT "G=$G\n";
            print $fileOUT "C=$C\n";
}

close ($fileIN);
close ($fileOUT);
4
  • What is exactly the problem? If you want to get the frequency for all lines together, just move the print statements outside the loop Jan 13, 2014 at 12:29
  • thanks, moving print statements out of the loop worked.My concern is ,although I have replaced \n character, why is that I still see nucleotides in multiple line instead of being in one single line ?
    – harsh
    Jan 13, 2014 at 15:25
  • Because you are reading the file line by line with the while statement. If you want to read all the lines in a scalar you need to make it different. For example use file::slurp module (there are more ways just look for them in google) Jan 13, 2014 at 19:18
  • It appears that you just want a total nucleotide frequency for all sequences in you fasta file. Here's a one-liner that will do this: perl -076 -ne 'chomp; s/.+?\n//; $h{$1}++ while /(A|T|C|G)/g; END{print "$_=$h{$_}\n" for keys %h}' inFile.fa [>outFile.out] The last, optional parameter directs output to a file.
    – Kenosis
    Jan 13, 2014 at 20:50

1 Answer 1

1

At first, i would use some shortcuts. Its easier to read:

use strict;
use warnings;
use feature 'say';
my $A;
my $T;
my $G;
my $C;
my $fileIN;
my $fileOUT;

open $fileIN,  '<',"basecount.nfasta" or die "can't open file basecount.nfasta for reading";
open $fileOUT, '>','basecount.out' or die "can't open file basecount.out for writing";

while ( my $seq = <$fileIN> ) {

  next if $seq =~ /^>/;
  $seq =~ s/\n//g;
  say $seq;

  my @dna = split //, $seq;

  foreach my $element ( @dna ) {
    $A++ if $element =~ m/A/;
    $T++ if $element =~ m/T/;
    $G++ if $element =~ m/G/;
    $C++ if $element =~ m/C/;
  }

  say $fileOUT "A=$A";
  say $fileOUT "T=$T";
  say $fileOUT "G=$G";
  say $fileOUT "C=$C";
}

close $fileIN;
close $fileOUT;

Using the 3 statement open is also recommended ( and a good die warning as well ).

EDIT: I used use feature 'say' here because all of your prints end with a newline. say does exactly the same like print, just with adding newlines at the end.

1
  • thanks for improvement; particularly for 'say'; its helpful.
    – harsh
    Jan 13, 2014 at 15:28

Not the answer you're looking for? Browse other questions tagged or ask your own question.