[Bioperl-l] Finding possible primers regex
Hilmar Lapp
hlapp at gmx.net
Sat Aug 9 16:07:30 UTC 2008
This looks like a neat trick. Do you think it's worth including as a
SimpleAlign method (obviously w/o the printing to STDOUT)? I can
imagine that a lot of people might appreciate it.
-hilmar
On Aug 4, 2008, at 12:08 AM, Chris Fields wrote:
> On Aug 2, 2008, at 3:05 PM, Benbo wrote:
>
>>
>> Hi there,
>> I'm trying to write a perl script to scan an aligned multiple entry
>> fasta
>> file and find possible primers. So far I've produced a string which
>> contains
>> bases which match all sequences and * where they don't match e.g.
>> 1) TTAGCCTAA
>> 2) TTAGCAGAA
>> 3) TTACCCTAA
>>
>> would give TTA*C**AA.
>>
>> I want to parse this string and pull out all sequences which are
>> 18-21 bp in
>> length and have no more than 4 * in them.
>>
>> So far, I've got this:
>>
>> while($fragment_match =~ /([GTAC*]{18,21})/g){
>> print "$1\n";
>> }
>>
>> hoping to match all fragments 18-21 characters in length. However
>> even that
>> doesn't work as it has essentially chunked it into 21 char blocks,
>> rather
>> than what I hoped for of
>> 0-18
>> 0-19
>> 0-20
>> 0-21
>> 1-19
>> 1-20
>> 1-21
>> 1-22
>>
>> etc.
>>
>> Can anyone let me know if this is already possible in BioPerl, or
>> how one
>> would go about it with regex. Sadly I'm fairly new to perl and
>> getting to
>> grips with BioPerl, so please treat me gently :).
>>
>> Many thanks,
>>
>> Ben
>
> There is a trick to this which is discussed more extensively in
> 'Mastering Regular Expressions'. Essentially you have to embed code
> into the regex and trick the parser into backtracking using a
> negative lookahead. The match itself fails (i.e. no match is
> returned), but the embedded code is executed for each match attempt,
>
> The following script is a slight modification of one I used which
> checks the consensus string from the input alignment (in aligned
> FASTA format here), extracts the alignment slice using that match,
> then spit the alignment out to STDOUT in clustalw format. This
> should work for perl 5.8 and up, but it's only been tested on perl
> 5.10. You should be able to use this to fit what you want.
>
> my $in = Bio::AlignIO->new(-file => $file,
> -format => 'fasta');
> my $out = Bio::AlignIO->new(-fh => \*STDOUT,
> -format => 'clustalw');
>
> while (my $aln = $in->next_aln) {
> my $c = $aln->consensus_string(100);
> my @matches;
> $c =~ m/
> ([GTAC?]{18,21})
> (?{my $match = check_match($1);
> push @matches, [$match,
> pos(),
> length($match)]
> if defined $match;})
> (?!)
> /xig;
> for my $match (@matches) {
> my ($hit, $st, $end) = ($match->[0],
> $match->[1] - $match->[2] + 1,
> $match->[1]);
> my $newaln = $aln->slice($st, $end);
> $out->write_aln($newaln);
> }
> }
>
> sub check_match {
> my $match = shift;
> return unless $match;
> my $ct = $match =~ tr/?/?/;
> return $match if $ct <= 4;
> }
>
>
> chris
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
--
===========================================================
: Hilmar Lapp -:- Durham, NC -:- hlapp at gmx dot net :
===========================================================
More information about the Bioperl-l
mailing list