Bio/Tools/SiRNA.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tools::SiRNA
   4 #
   5 # Cared for by Donald Jackson, donald.jackson@bms.com
   6 #
   7 # Copyright Bristol-Myers Squibb
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 SiRNA - Perl object for designing small inhibitory RNAs.
  16
  17 =head1 SYNOPSIS
  18
  19   use Bio::Tools::SiRNA;
  20
  21   my $sirna_designer = Bio::Tools::SiRNA->new( -target => $bio_seq,
  22                                                -rules  => 'saigo'
  23     );
  24   my @pairs = $sirna_designer->design;
  25
  26   foreach $pair (@pairs) {
  27       my $sense_oligo_sequence = $pair->sense->seq;
  28       my $antisense_oligo_sequence = $pair->antisense->seq;
  29
  30       # print out results
  31       print join ("\t", $pair->start, $pair->end, $pair->rank,
  32                   $sense_oligo_sequence, $antisense_oligo_sequence), "\n";
  33   }
  34
  35 =head1 DESCRIPTION
  36
  37 Package for designing siRNA reagents.
  38
  39 Input is a L<Bio::SeqI>-compliant object (the target).
  40
  41 Output is a list of Bio::SeqFeature::SiRNA::Pair objects, which are
  42 added to the feature table of the target sequence.  Each
  43 Bio::SeqFeature::SiRNA::Pair contains two subfeatures
  44 (Bio::SeqFeature::Oligo objects) which correspond to the individual
  45 oligos.  These objects provide accessors for the information on the
  46 individual reagent pairs.
  47
  48 This verion of Bio::Tools::SiRNA represents a major change in architecture.
  49 Specific 'rulesets' for siRNA selection as developed by various groups are
  50 implemented as Bio::Tools::SiRNA::Ruleset objects, which inherit from
  51 Bio::Tools::SiRNA.  This will make it easier to add new rule sets or modify
  52 existing approaches. Currently the Tuschl and Ui-Tei (2004) rules are
  53 implemented. For consistency, the Tuschl rules are implemented by default.
  54
  55 In addition, this module provides three 'extra' rules which can be added
  56 above and beyond any ruleset.
  57
  58 =over 3
  59
  60 =item 1.
  61
  62 SiRNAs that overlap known SNPs (identified as SeqFeatures with
  63 primary tag = variation) can be avoided.
  64
  65 =item 2.
  66
  67 Other regions (with primary tag = 'Excluded') can also be skipped.  I
  68 use this with Bio::Tools::Run::Mdust to avoid low-complexity regions
  69 (must be run separately), but other programs could also be used.
  70
  71 =item 3.
  72
  73 SiRNAs may also be selected in the 3 prime UTR of a gene by setting
  74 $sirna_designer-E<gt>include_3pr() to true.
  75
  76 =back
  77
  78 =head2 EXPORT
  79
  80 None.
  81
  82 =head1 SEE ALSO
  83
  84 L<Bio::Tools::Run::Mdust>, L<Bio::SeqFeature::SiRNA::Pair>,
  85 L<Bio::SeqFeature::SiRNA::Oligo>..
  86
  87 =head1 FEEDBACK
  88
  89 =head2 Mailing Lists
  90
  91 User feedback is an integral part of the evolution of this and other
  92 Bioperl modules. Send your comments and suggestions preferably to
  93 the Bioperl mailing list.  Your participation is much appreciated.
  94
  95   bioperl-l@bioperl.org                  - General discussion
  96   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  97
  98 =head2 Reporting Bugs
  99
 100 Report bugs to the Bioperl bug tracking system to help us keep track
 101 of the bugs and their resolution. Bug reports can be submitted via
 102 the web:
 103
 104   http://bugzilla.open-bio.org/
 105
 106 =head1 AUTHOR
 107
 108 Donald Jackson (donald.jackson@bms.com)
 109
 110 =head1 APPENDIX
 111
 112 The rest of the documentation details each of the object methods.
 113 Internal methods are usually preceded with a _
 114
 115 =cut
 116
 117 package Bio::Tools::SiRNA;
 118
 119 use strict;
 120 use warnings;
 121
 122 use vars qw($AUTOLOAD);
 123
 124 use Bio::Seq::RichSeq;
 125 use Bio::SeqFeature::Generic;
 126 use Bio::SeqFeature::SiRNA::Oligo;
 127 use Bio::SeqFeature::SiRNA::Pair;
 128
 129
 130 use base qw(Bio::Root::Root);
 131
 132
 133 our %COMP = ( A => 'T',
 134               T => 'A',
 135               C => 'G',
 136               G => 'C',
 137               N => 'N',
 138               );
 139
 140 our @ARGNAMES = qw(RULES START_PAD END_PAD MIN_GC CUTOFF OLIGOS AVOID_SNPS
 141                    GSTRING TMPDIR TARGET DEBUG);
 142
 143
 144 =head2 new
 145
 146  Title          : new
 147  Usage          : my $sirna_designer = Bio::Tools::SiRNA->new();
 148  Function       : Constructor for designer object
 149  Returns        : Bio::Tools::SiRNA object
 150  Args           : target - the target sequence for the SiRNAs as a Bio::Seq::RichSeq
 151                   start_pad - distance from the CDS start to skip (default 75)
 152                   end_pad - distance from the CDS end to skip (default 50)
 153                   include_3pr - set to true to include SiRNAs in the 3prime UTR (default false)
 154                   rules - rules for selecting siRNAs, currently supporting saigo and tuschl
 155                   min_gc - minimum GC fraction (NOT percent) (default 0.4)
 156                   max_gc - maximum GC fraction (NOT percent) (default 0.6)
 157                   cutoff - worst 'rank' accepted(default 3)
 158                   avoid_snps - boolean - reject oligos that overlap a variation
 159                      SeqFeature in the target (default true)
 160                   gstring - maximum allowed consecutive Gs.
 161                      Too many can cause problems in synthesis (default 4)
 162   Note          : All arguments can also be changed/accessed using autoloaded
 163                  methods such as:
 164
 165     my $start_pad = $sirna_designer->start_pad().
 166
 167 =cut
 168
 169 sub new {
 170     my ($proto, @args) = @_;
 171     my $pkg = ref($proto) || $proto;
 172
 173     my $self = {};
 174     bless ($self, $pkg);
 175
 176     my %args;
 177
 178     @args{@ARGNAMES} = $self->_rearrange(\@ARGNAMES, @args);
 179
 180     if ($args{'RULES'}) {
 181         $self->rules($args{'RULES'});
 182     }
 183
 184     $self->{'start_pad'} = $args{'START_PAD'} || 75; # nt from start to mask
 185     $self->{'end_pad'} = $args{'END_PAD'} || 50; # nt from end to mask
 186     $self->{'include_3pr'} = $args{'INCLUDE_3PR'} || 0; # look for oligos in 3prime UTR
 187     $self->{'min_gc'} = $args{'MIN_GC'} || 0.40;
 188     $self->{'max_gc'} = $args{'MAX_GC'} || 0.60;
 189     $self->{'cutoff'} = $args{'CUTOFF'} || 3; # highest (worst) rank wanted
 190     $self->{'oligos'} = [];
 191     defined($args{'AVOID_SNPS'}) ? $self->{'avoid_snps'} = $args{'AVOID_SNPS'} :
 192         $self->{'avoid_snps'} = 1; # (t/f to avoid or include reagents that cover SNPs)
 193     $self->{'gstring'} = $args{'GSTRING'} || 4; # maximum allowed consecutive Gs - too many can cause problems in oligo synthesis
 194     $self->{'tmpdir'} = $args{'TMPDIR'}  || $ENV{'TMPDIR'} || $ENV{'TMP'} || '';
 195     $self->{'debug'} = $args{'DEBUG'} || 0;
 196
 197     $self->target($args{'TARGET'}) if ($args{'TARGET'});
 198
 199     return $self;
 200 }
 201
 202
 203 =head2 target
 204
 205   Title         : target
 206   Usage         : my $target_seq = $sirna_designer->target(); # get the current target
 207                   OR
 208                   $sirna_designer->target($new_target_seq); # set a new target
 209   Function      : Set/get the target as a Bio::SeqI-compliant object
 210   Returns       : a Bio::SeqI-compliant object
 211   Args          : a Bio::SeqI-compliant object (optional)
 212
 213 =cut
 214
 215 sub target {
 216     my ($self, $target) = @_;
 217
 218     if ($target) {
 219         unless ($target->isa('Bio::SeqI')) {
 220             $self->throw(  -class => 'Bio::Root::BadParameter',
 221                            -text  => "Target must be passed as a Bio::Seq object" );
 222         }
 223         if ($target->can('molecule')) {
 224             ( grep { uc($target->molecule) eq $_ } qw(DNA MRNA CDNA)) or
 225                 $self->throw(  -class => 'Bio::Root::BadParameter',
 226                                -text  =>  "Sequences of type ". $target->molecule. " are not supported"
 227                                );
 228         }
 229         else {
 230             ($target->alphabet eq 'dna') or
 231                 $self->throw(  -class => 'Bio::Root::BadParameter',
 232                                -text  =>  "Sequences of alphabet ". $target->alphabet. " are not supported"
 233                                );
 234         }
 235
 236         $self->{'target'} = $target;
 237         return 1;
 238
 239     }
 240     elsif ($self->{'target'}) {
 241         return $self->{'target'};
 242     }
 243     else {
 244         $self->throw("Target sequence not defined");
 245     }
 246 }
 247
 248 =head2 rules
 249
 250     Title       : rules
 251     Usage       : $sirna->rules('ruleset')
 252     Purpose     : set/get ruleset to use for selecting SiRNA oligo pairs.
 253     Returns     : not sure yet
 254     Args        : a ruleset name (currently supported: Tuschl, Saigo)
 255                   or a Bio::Tools::SiRNA::RulesetI compliant object
 256
 257 =cut
 258
 259 sub rules {
 260     my ($self, $rules) = @_;
 261
 262     if ($rules) {
 263         $self->_load_ruleset($rules);
 264     }
 265     # default: use tuschl rules
 266     unless ($self->{_rules}) {
 267         $self->_load_ruleset('tuschl');
 268     }
 269     return $self->{_rules};
 270 }
 271
 272 sub _load_ruleset {
 273     my ($self, $ruleset) = @_;
 274
 275     my $rule_module = join('::', ref($self), 'Ruleset', lc($ruleset));
 276
 277     eval "require $rule_module";
 278
 279     if ($@) {
 280         #warn join("\n", '@INC contains:', @INC, undef);
 281         $self->throw("Unable to load $rule_module: $@");
 282         return;
 283     }
 284
 285     else {
 286         $self->{_rules} = $rule_module;
 287         bless($self, $rule_module); # recast as subclass
 288     }
 289
 290     return 1;
 291 }
 292
 293 =head2 design
 294
 295   Title         : design
 296   Usage         : my @pairs = $sirna_designer->design();
 297   Purpose       : Design SiRNA oligo pairs.
 298   Returns       : A list of SiRNA pairs as Bio::SeqFeature::SiRNA::Pair objects
 299   Args          : none
 300
 301 =cut
 302
 303 sub design {
 304     my ($self) = @_;
 305
 306     ($self->rules) or $self->throw('Unable to design siRNAs: no rule set specified');
 307
 308 #     unless ( grep { $_->primary_tag eq 'Target' } $self->target->top_SeqFeatures ) {
 309 #       $self->_define_target();
 310 #     }
 311
 312     my @oligos = $self->_get_oligos();
 313
 314     return ( grep { $_->isa('Bio::SeqFeature::SiRNA::Pair') } $self->target->top_SeqFeatures );
 315 }
 316
 317 sub _define_target {
 318     my ($self) = @_;
 319     my ($feat, $cds, $left, $right);
 320
 321     my $target = $self->target or
 322         $self->throw("Unable to design oligos - no target provided");
 323
 324     ($cds) = grep { $_->primary_tag eq 'CDS' } $target->top_SeqFeatures if ($target->can('top_SeqFeatures'));
 325
 326     if ($cds) {
 327         $left = $cds->start + $self->start_pad;
 328         if (!$self->include_3pr) {
 329             $right = $cds->end - $self->end_pad;
 330         }
 331         else {
 332             $right = $target->length - $self->end_pad;
 333         }
 334     }
 335     else {
 336         $left = 0 + $self->start_pad;
 337         $right = $target->length - $self->end_pad;
 338     }
 339
 340
 341     # is there anything left?
 342     if (($right - $left) < 20) {
 343         $self->throw("There isn't enough sequence to design oligos.  Please reduce start_pad and end_pad or supply more sequence");
 344     }
 345     # define target region
 346     my $targregion = Bio::SeqFeature::Generic->new( -start              => $left,
 347                                                     -end                => $right,
 348                                                     -primary            => 'Target' );
 349     $self->target->add_SeqFeature($targregion);
 350
 351     # locate excluded regions
 352     my @excluded = grep { $_->primary_tag eq 'Excluded' } $self->target->top_SeqFeatures;
 353
 354     if ($self->avoid_snps) {
 355         my @snps =  grep { $_->primary_tag eq 'variation' } $self->target->top_SeqFeatures;
 356         push(@excluded, @snps);
 357     }
 358
 359     $self->excluded(\@excluded);
 360
 361     return $targregion;
 362 }
 363
 364 sub _get_targetregion {
 365     my ($self) = @_;
 366
 367     my ($targregion) = grep { $_->primary_tag eq 'Target' } $self->target->top_SeqFeatures;
 368     $targregion ||= $self->_define_target;
 369
 370     $self->throw("Target region for SiRNA design not defined") unless ($targregion);
 371
 372     my $seq = $targregion->seq->seq;
 373     # but this way I loose start info
 374      my $targstart = $targregion->start;
 375
 376     return ($seq, $targstart);
 377 }
 378
 379 # MOVE to SiRNA::Ruleset::tuschl
 380 # sub _regex {
 381 #     my ($self, $rank) = @_;
 382 #     return $PATTERNS{$rank};
 383 # }
 384
 385 # sub _get_oligos {
 386 #     # use regular expressions to pull out oligos
 387
 388 #     my ($self, $rank) = @_;
 389 #     my $regex = $self->_regex($rank);
 390 #     my @exclude;
 391
 392
 393 #     my ($targregion) = grep { $_->primary_tag eq 'Target' } $self->target->top_SeqFeatures;
 394 #     my $seq = $targregion->seq->seq;
 395 #     # but this way I loose start info
 396 #     my $targstart = $targregion->start;
 397
 398 #     # exclude masked region
 399 #     push(@exclude, grep { $_->primary_tag eq 'Excluded' } $self->target->top_SeqFeatures);
 400
 401 #     # add SNP checking
 402 #     if ($self->avoid_snps) {
 403 #       my @snps =  grep { $_->primary_tag eq 'variation' } $self->target->top_SeqFeatures;
 404 #       push(@exclude, @snps);
 405 #     }
 406
 407 #     while ( $seq =~ /$regex/gi ) {
 408 #       my $target = $1;
 409
 410 #       # check for too many Gs (or Cs on the other strand)
 411 #       next if ( $target =~ /G{ $self->gstring,}/io );
 412 #       next if ( $target =~ /C{ $self->gstring,}/io );
 413 #       # skip Ns (for filtering)
 414 #       next if ( $target =~ /N/i);
 415
 416 #       my $start = length($`) + $targstart;
 417 #       my $stop = $start + length($target) -1;
 418
 419 #       my @gc = ( $target =~ /G|C/gi);
 420 #       my $fxGC = sprintf("%2.2f", (scalar(@gc) / length($target)));
 421 #       next if ($fxGC < $self->min_gc);
 422 #       next if ($fxGC > $self->max_gc);
 423
 424 #       my $sense = Bio::SeqFeature::SiRNA::Oligo->new( -start          => $start,
 425 #                                                       -end            => $stop,
 426 #                                                       -strand         => 1,
 427 #                                                       -seq            => _get_sense($target),
 428 #                                                       -source_tag     => ref($self),
 429 #                                                      );
 430
 431 #       my $asense = Bio::SeqFeature::SiRNA::Oligo->new( -start         => $start,
 432 #                                                        -end           => $stop,
 433 #                                                        -strand        => -1,
 434 #                                                        -seq           => _get_anti($target),
 435 #                                                        -source_tag    => ref($self),
 436 #                                                        );
 437
 438 #       my $sirna = Bio::SeqFeature::SiRNA::Pair->new( -rank            => $rank,
 439 #                                                      -fxGC            => $fxGC,
 440 #                                                      -sense           => $sense,
 441 #                                                      -antisense       => $asense,
 442 #                                                      -source_tag      => ref($self),
 443 #                                                      );
 444
 445 #       unless ($self->_has_overlap($sirna, \@exclude)) {
 446 #           $self->target->add_SeqFeature($sirna);
 447 #       }
 448 #     }
 449 # }
 450
 451 =head2 add_oligos
 452
 453   Title         : add_oligos
 454   Usage         : $sirna_designer->add_oligos($sequence, $start, $rank);
 455   Purpose       : Add SiRNA olgos to target Bio::Seq as Bio::SeqFeature::SiRNA::Pair objects
 456   Args          : Oligo sequence and start position (required), rank/score (optional)
 457
 458 =cut
 459
 460 sub add_oligos {
 461     my ($self, $seq, $start, $rank) = @_;
 462
 463     ($seq) or throw ('No sequence supplied for add_oligos');
 464     (defined $start) or throw ('No start position specified for  add_oligos');
 465
 466     my ($end) = $start + length($seq);
 467
 468     my ($sseq) = $self->_get_sense($seq);
 469     my $sense = Bio::SeqFeature::SiRNA::Oligo->new( -start              => $start,
 470                                                     -end                => ($start + length($sseq)),
 471                                                     -strand     => 1,
 472                                                     -seq                => $sseq,
 473                                                     -source_tag => ref($self),
 474                                                     );
 475
 476     my $aseq = $self->_get_anti($seq);
 477     my $asense = Bio::SeqFeature::SiRNA::Oligo->new( -start             => $end,
 478                                                      -end               => ($end - length($aseq)),
 479                                                      -strand            => -1,
 480                                                      -seq               => $aseq,
 481                                                      -source_tag        => ref($self),
 482                                                      );
 483
 484     my $sirna = Bio::SeqFeature::SiRNA::Pair->new( -rank                => $rank,
 485                                                   # -fxGC               => $fxGC,
 486                                                    -sense               => $sense,
 487                                                    -antisense   => $asense,
 488                                                    -source_tag  => ref($self),
 489                                                    );
 490
 491     unless ($self->_has_overlap($sirna, $self->excluded)) {
 492         $self->target->add_SeqFeature($sirna);
 493     }
 494 }
 495
 496 sub _has_overlap {
 497     # flag any pairs that overlap an UNDESIRED feature (eg SNP)
 498     # return true if there is overlap, false if not
 499
 500     my ($self, $test, $flist) = @_;
 501     print STDERR "Checking oligo at ", $test->start, " to ",$test->end, "\n"
 502         if ($self->debug);
 503
 504     foreach my $feat (@$flist) {
 505         if (($test->start <= $feat->end) and ($test->end >= $feat->start)) {
 506             print STDERR "Overlaps ", $feat->primary_tag, " at ",
 507             $feat->start, " to ", $feat->end, "\n" if ($self->debug);
 508             return 1;
 509         }
 510     }
 511     return 0; # default - no overlap
 512 }
 513
 514 # MOVE to SiRNA::Ruleset::tuschl
 515
 516 # sub _get_sense {
 517 #     my ($target) = @_;
 518 #     # trim off 1st 2 nt to get overhang
 519 #     $target =~ s/^..//;
 520 #     # convert T's to U's (transcribe)
 521 #     $target =~ s/T/U/gi;
 522 #     # force last 2 nt to be T's
 523 #     $target =~ s/..$/TT/;
 524
 525 #     return $target;
 526 # }
 527
 528 # sub _get_anti {
 529 #     my ($target) = @_;
 530 #     my @target = split(//, $target);
 531 #     my ($nt,@antitarget);
 532
 533 #     while ($nt = pop @target) {
 534 #       push(@antitarget, $COMP{$nt});
 535 #     }
 536 #     my $anti = join('', @antitarget);
 537 #     # trim off 1st 2 nt to get overhang
 538 #     $anti =~ s/^..//;
 539 #     # convert T's to U's
 540 #     $anti =~ s/T/U/gi;
 541 #     # convert last 2 NT's to T
 542 #     $anti =~ s/..$/TT/;
 543
 544 #     return $anti;
 545 # }
 546
 547
 548 sub AUTOLOAD {
 549     my ($self, $value) = @_;
 550     my $name = $AUTOLOAD;
 551     $name =~ s/.+:://;
 552
 553     return if ($name eq 'DESTROY');
 554
 555
 556     if (defined $value) {
 557         $self->{$name} = $value;
 558     }
 559
 560     unless (exists $self->{$name}) {
 561         $self->throw("Attribute $name not defined for ". ref($self));
 562     }
 563
 564     return $self->{$name};
 565 }
 566
 567 sub _comp {
 568     my ($self, $char) = @_;
 569
 570     return unless ($char);
 571     $char = uc($char);
 572     return $COMP{ $char };
 573 }
 574 1;