← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to examples/generate_random_seq.pl

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/generate_random_seq.pl
 
1
#!/bin/perl -w
 
2
use strict;
 
3
use vars qw($USAGE);
 
4
 
 
5
# random sequence generator #
 
6
# -c=1 option will cause prot sequences to be built 
 
7
# using vertebrate aa frequencies, 
 
8
# with option -a putting a 1st methionine residues on. Frequencies are
 
9
# calculated from the NCBI human RefSeq protein sequences 
 
10
# -c and -a only affect protein sequences
 
11
# -a only works in conjunction with -c
 
12
# -n number of random sequences, default = 1
 
13
 
 
14
use Bio::PrimarySeq;
 
15
use Bio::SeqIO;
 
16
use Getopt::Long;
 
17
my ($length,$type,$filename,$comp,$met);
 
18
 
 
19
$USAGE = 'usage: generate_random_seq.pl --length=1000 --type=dna --filename=/tmp/test.seq --number=50';
 
20
 
 
21
my %alphabets = ( 'dna' => [qw(C A G T)],
 
22
                  'rna' => [qw(C A G U)],
 
23
                  'prot'=> [qw( A C D E F G H I K L M N P Q R S T V W Y)],
 
24
              );
 
25
# make random num from 1-10000. numbers in this array reflect the frequency,
 
26
# e.g., a random number from 1.744 = A, 745-991 = C etc;
 
27
my @aa_frequencies = qw(744 991 1398 2017 2378 3104 3349 3726 4239 5273 5443 
 
28
                        5749 6410 6848 7455 8263 8760 9340 9488 9713 10000);
 
29
my $number = 1;
 
30
 
 
31
&GetOptions
 
32
  (
 
33
   'l|length:s'          => \$length,
 
34
   't|type|m|alphabet:s' => \$type,
 
35
   'f|file|filename:s'   => \$filename,
 
36
   'c|composition:s'     => \$comp,
 
37
   'a|methionine:s'      => \$met,
 
38
   'n|number:s'          => \$number
 
39
  );
 
40
 
 
41
assert ( $type && defined ($alphabets{lc $type}),
 
42
         $USAGE);
 
43
assert ( $length && $length =~ /^\d+$/, $USAGE );
 
44
 
 
45
foreach my $num (1..$number) {
 
46
   my $sequence = "";
 
47
   my $alphabet = $alphabets{lc $type};
 
48
   my $sspace = scalar @$alphabet;
 
49
   if (!$comp || $type ne 'prot') {
 
50
      foreach ( 1..$length ) {
 
51
         $sequence .= $alphabet->[ int rand($sspace) ];
 
52
      }
 
53
   }elsif ($type eq 'prot') {
 
54
      $sequence = build_seq($length, \@aa_frequencies);
 
55
   }
 
56
   my $seq =  Bio::PrimarySeq->new(-seq        => $sequence, 
 
57
                                   -display_id => 'randomseq'.$num);
 
58
   my %args = (-format => 'fasta');
 
59
   if( $filename ) { $args{-file} = ">>$filename" }
 
60
   my $seqio = Bio::SeqIO->new(%args);
 
61
   $seqio->write_seq($seq);
 
62
}
 
63
 
 
64
sub assert { die $_[1] unless( $_[0] ); }
 
65
 
 
66
sub build_seq {
 
67
   #takes seqlen and ref to frequency data as parameters
 
68
   my ($len, $pf)  = @_;
 
69
   my $str = ($met)?'M':'';
 
70
   my $i = ($met)?1:0;
 
71
   for ($i..$len-1) {
 
72
      my $aa = int(rand (10000)) ;
 
73
      my $j = 0;
 
74
      while ($pf->[$j] < $aa && $j <19) {
 
75
         $j++;
 
76
      }
 
77
      $str .= $alphabets{'prot'}[$j];
 
78
   }
 
79
   print "str is $str\n";
 
80
   return $str;
 
81
}