Bioperl: generating sequence from IUB Base codes

[Aaron J Mackey]

Try this:

(btw, if there's general interest, I can turn this into a
Bio/Tools/IUPAC.pm module ...)

-Aaron

--------------------
#!/usr/local/bin/perl -w

require 5.004;

use Bio::SeqIO;
use Bio::Seq;

my %iupac = ( A => [qw(A)],
	      C => [qw(C)],
	      G => [qw(G)],
	      T => [qw(T)],
	      U => [qw(U)],
	      M => [qw(A C)],
	      R => [qw(A G)],
	      W => [qw(A T)],
	      S => [qw(C G)],
	      Y => [qw(C T)],
	      K => [qw(G T)],
	      V => [qw(A C G)],
	      H => [qw(A C T)],
	      D => [qw(A G T)],
	      B => [qw(C G T)],
	      X => [qw(G A T C)],
	      N => [qw(G A T C)]
	    );


my $in = new Bio::SeqIO -file => 'infile',
                        -format => 'Fasta';
my $out = new Bio::SeqIO -file => 'outfile',
                         -format => 'Fasta';

my @permutations = ();

while ( my $seq = $in->next_seq() ) {
    $seq->strict(0);
    $seq->alphabet_ok() or next; # silent skipping!
    for ( get_sequences($seq) ) {
	$out->write_seq($_);
    }
}

sub get_sequences {
    my ($seq) = @_;
    @permutations = ();
    fill_permutations('', $seq->str());
    return map { Bio::Seq->new(-seq => $_,
			       -id => $seq->id(),
			       -desc => $seq->desc()
			      ) } @permutations;
}

sub fill_permutations {
    my ($parsed, $unparsed) = @_;
    return push @permutations, $parsed unless $unparsed;
    for my $letter ($iupac{uc(substr($unparsed, 0, 1))}) {
	return fill_permutations($parsed . $letter,
				 substr($unparsed, 1)
				);
    }
}

------------------------

On Mon, 12 Apr 1999, Olivier Dugas wrote:

> 	Hello,
> 
> I'm trying to write a script wich is able to generate all sequences with
> real bases (ATGC) from one sequence with IUB codes (Y,S,R etc...).
> I think it's not hard but I have spend so much time on this task without
> more result than :
> seq file - > out file
> AYRT->		A
> 		CT
> 		AG
> 		T
> There's some script wich can help me ?
> 
> Olivier DUGAS.
> =========== Bioperl Project Mailing List Message Footer =======
> Project URL: http://bio.perl.org/
> For info about how to (un)subscribe, where messages are archived, etc:
> http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
> ====================================================================
> 

-- 
 o ~   ~   ~   ~   ~   ~  o
/ Aaron J Mackey           \
\  Dr. Pearson Laboratory  / 
 \ University of Virginia  \     
 /  (804) 924-2821          \
 \  amackey@virginia.edu    /
  o ~   ~   ~   ~   ~   ~  o


=========== Bioperl Project Mailing List Message Footer =======
Project URL: http://bio.perl.org/
For info about how to (un)subscribe, where messages are archived, etc:
http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
====================================================================