Count Occurrences of a string

[MiamiGenome]

Thanks for the scripting tips!

I have a 'counting' issue which I need to quickly resolve. A typical sequence input file (5 - 700 bases) looks like :

AGTAGTCGATCATNATANCTANTACNACTACTAACTATGCTAGNNAATATAAAAAAAAANAAA

I have over 500 files, named *.seq. I would like to create a script which :

a. runs through all the files,

b. counts the length of the 'poly A' tail (defined as the longest stretch containing either A or N)

c. sends the output to a file, eg.

25 1.seq

87 2.seq

13 3.seq

Example valid poly A tails (I desire the full length reported) :

AAAANANANANAAANNAAAAAA
AAAAAAAAAAAAAA
NNNNNNNNNNNNN
AAANNNNNNNNNNNAAAAAAAAA

Thank you so much for your expertise!

[Zaxo]

Let's write this step by step. We will write a program called ctpA which takes a directory path as argument. The outpur file will be called pA.cat, written to the specified directory, and the program will overwrite any old pA.cat file.

We want to run this as an interpreter script, so we start with a shebang line and shift in the directory argument. If there is none, we use the current working directory.

#!/usr/bin/perl -w
use strict;

die 'Usage: ctpA [DIR]' if @ARGV > 1;
my $dir = shift || '.';

The || operator brings in the default if no argument is given. To be careful, you should check $dir here for mischievous values.
We need to look at all the files matching a naming pattern, so glob is the easy way to get a list of those names. Since we'll want the file names stripped out of the path for the summary we will write, we make the File::Basename module available. We also open the output file for writing.

use File::Basename;
open local(*OUT), "> $dir/pA.cat" or die $!;

We begin to look at the files one at a time by starting a for loop over their names. for (glob "${dir}/*.seq" ) { This makes $_ contain the file path for each in turn. We set the input record seperator, $/, to undef to slurp the entire file, open, read, and close,

    local $/;
    open local(*IN), "< $_" or die $!;
    my $seq = <IN>;
    close IN or die $!;

We get rid of all whitespace, including newlines, just in case there are some.     $seq =~ s/\s//g; and do an anchored match at the tail of the string. This regex will match any $seq, since the putative tail may be of zero length.

    $seq =~ /([AN]*)$/;
    printf OUT "%d %s.seq\n", length $1, basename $_;
}

That regular expression makes A and N a character class and looks for the longest run of them at the end of the string. The parentheses capture the result to $1, used on the output line next. If there is no tail, the match still succeeds and length $1 is zero. That concludes the work to be done for a file, so we close the loop.

All that remains is to close the output file. close OUT or die $!; This is untested - I don't have that sort of data around - but it is all straightforward stuff.

After Compline,
Zaxo

[antirice]

Not too difficult:

my $string = "AAACANANANTANAACATNNCAAAAAA";
my @matches = $string =~ /([AN]+)/g;
$\=$,=$/;
print @matches;
__DATA__
output:
AAA
ANANAN
ANAA
A
NN
AAAAAA

As for everything else: opendir, readdir, closedir, glob, open, perlopentut, and of course the Schwartzian Transform.

Hope this helps.

antirice    
The first rule of Perl club is - use Perl
The ith rule of Perl club is - follow rule i - 1 for i > 1

[asarih]

For your long-term benefit, one of a growing number of Perl books with bioinformatic slant, such as ISBN 0596000804, ISBN 0596003072, ISBN 052180177X, ISBN 1565926641 may be beneficial to you. (I am not qualified to comment on merit or demerit on any of these.)

A CPAN search on bioinformatics turns up a few modules, and bioperl.org might be of intereste to you.

[biosysadmin]

Here's a solution using Bioperl:

#!/usr/bin/perl -w

use strict;
use Bio::SeqIO;

my @files   = <*.seq>;
my $pattern = '([AN]*)';

foreach my $file (@files) {
   my $seqio = Bio::SeqIO->new( -file => $file );
   while ( my $seqobj = $seqio->next_seq() ) {
       my $raw_seq = $seqobj->seq();
       my @length_of_matches;

       while ( $raw_seq =~ /$pattern/og ) {
           push @length_of_matches, length $1;
       }

       @length_of_matches = sort { $b <=> $a } @length_of_matches;
       print $length_of_matches[ 0 ], ' ', $file, "\n";
   }
}

Using Bioperl for easy scripts likes this is a good idea, because it allows you to very easily port your program to almost any sequence file format known to man.

[MiamiGenome]

Dear Perl Monks,

Thank you all for your constructive and instructive help.

I am dealing with 2 conflicting, yet synergistic issues in my research, namely 1. getting my analyses completed, and 2. learning BioPerl/Perl.

For nearly all my scripting, I am whip out quick jiffies to answer the questions of our genome center, but as my BioPerl comes up to speed, I sporadically need more guidance.

I again thank you all.

(No, my question wasn't homework :) )

[biosysadmin]

Good luck with your research.In addition to Perl Monks, a good resource for learning BioPerl is the doc.bioperl.org website, and the BioPerl mailing lists.

Your name suggests that we may have a mutual friend in common, do you know Ken Rudd at the University of Miami? If so, tell him that Tex Thompson from the Rochester Institute of Technology says hello. :)

[Anonymous Monk]

Your questions are building on top of one another. Did you try using BioPerl to loop through the files in the directory, and then adding the reg-ex to it? That's 90% of the battle, and if you did you'd probably get more helpful and enthusiastic responses.

People don't mind helping, but try to show that you're learning and making an attempt. A nice way to prove that, as Abigail said, is to show either a code fragment or a list of what you've tried.

[Abigail-II]

What did you try so far?

Abigail

[Marcello]

Homework? ;-)