Beefy Boxes and Bandwidth Generously Provided by pair Networks
No such thing as a small change

Count Occurrences of a string

by MiamiGenome (Sexton)
on Sep 09, 2003 at 22:05 UTC ( #290187=perlquestion: print w/replies, xml ) Need Help??

MiamiGenome has asked for the wisdom of the Perl Monks concerning the following question:

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 :


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) :


Thank you so much for your expertise!

Replies are listed 'Best First'.
Re: Count Occurrences of a string
by Zaxo (Archbishop) on Sep 10, 2003 at 00:27 UTC

    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, written to the specified directory, and the program will overwrite any old 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/" 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,

Re: Count Occurrences of a string
by antirice (Priest) on Sep 09, 2003 at 22:25 UTC

    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.

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

Re: Count Occurrences of a string
by asarih (Hermit) on Sep 09, 2003 at 22:40 UTC
Re: Count Occurrences of a string
by biosysadmin (Deacon) on Sep 10, 2003 at 10:57 UTC
    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.
      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 :) )

        Good luck with your research.In addition to Perl Monks, a good resource for learning BioPerl is the 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. :)
Re: Count Occurrences of a string
by Anonymous Monk on Sep 09, 2003 at 23:04 UTC

    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.

Re: Count Occurrences of a string
by Abigail-II (Bishop) on Sep 09, 2003 at 22:24 UTC
    What did you try so far?


      Homework? ;-)

Log In?

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://290187]
Approved by diotalevi
Front-paged by broquaint
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others contemplating the Monastery: (5)
As of 2022-09-25 02:11 GMT
Find Nodes?
    Voting Booth?
    I prefer my indexes to start at:

    Results (116 votes). Check out past polls.