PERL Beginners - Perl script

Author

Perl script

Auxence Sima

2007-07-24, 6:59 pm

Hello everyone ! I have this question: can anyone give me an idea of how i can do this. Im trying to figure out how to run Primer3 in batch mode from the Linux command line, and then process the result from a perl script ( which i have to write). any idea
how to generate primers for many sequences at once, or what the perl script would look like. i have an input file created for me to use.
Any suggestions???/

---------------------------------
Get the free Yahoo! toolbar and rest assured with the added security of spyware protection.

Chas Owens

2007-07-24, 6:59 pm

On 7/24/07, Auxence Sima <auxy25@yahoo.com> wrote:
> Hello everyone ! I have this question: can anyone give me
> an idea of how i can do this. Im trying to figure out how to
> run Primer3 in batch mode from the Linux command line,
> and then process the result from a perl script ( which i have
> to write). any idea how to generate primers for many
> sequences at once, or what the perl script would look like.
> i have an input file created for me to use.
> Any suggestions???/

The Primer3 mailing list is primer3-mail@lists.sourceforge.com. I am
sure they will be able to help you to figure out how to run it. As
for what the perl script would look like, you haven't told us what you
want to do or provided any code that has a problem, so we really can't
give you any more direction than

#!/usr/bin/perl

use strict;
use warnings;

#do stuff here

Auxence Sima

2007-07-25, 6:59 pm

Chas Owens <chas.owens@gmail.com> wrote: On 7/24/07, Auxence Sima wrote:
> Hello everyone ! I have this question: can anyone give me
> an idea of how i can do this. Im trying to figure out how to
> run Primer3 in batch mode from the Linux command line,
> and then process the result from a perl script ( which i have
> to write). any idea how to generate primers for many
> sequences at once, or what the perl script would look like.
> i have an input file created for me to use.
> Any suggestions???/

The Primer3 mailing list is primer3-mail@lists.sourceforge.com. I am
sure they will be able to help you to figure out how to run it. As
for what the perl script would look like, you haven't told us what you
want to do or provided any code that has a problem, so we really can't
give you any more direction than

#!/usr/bin/perl

use strict;
use warnings;

#do stuff here

I want to run primer3 from linux command line, and parse the outcome with a perl scrip , which im trying to write, in order to generate primers for more than one sequence at once.

here are the module i used, the sample code i wrote( but did not do what i intented to) and the input file.

#! /usr/bin/perl
#
# primer.pl
#
# Redas a list of accession numbers from the command line ...
#
#
#
#
# import required modules

use strict;
use Bio::Seq;
use Bio::DB::GenBank;
use Primer;

# set up variables;

my %sequences;
my $key;
my $feature;

#Fetch sequences from GenBank, keyed by acession
my $gbConnect = new Bio::GenBank();
my $seqIO = $gbConnect=>get_Stream_by_acc(@ARGV);
while (my $tempSeq = $seqIO=>next_seq()) {

$sequences{$tempSeq=>display_id} = $tempSeq;

}

#iterate

foreach $key (keys(%sequences)) {

#set arbitrary start point

my $start =50;
my @features = $sequences{$key}=>top_Seqfeatures();

# try to start in CDS if we can

foreach $feature (@features) {

if ($feature=>primary_tag eq `CDS`) {

if ($feature=>start() < 400) {

$start = $feature=>start();

}

last;

}
}

#search for primers

my $primer = Primer=>new($key,$sequences{$key}=>subseq($start,$start +
500));

print "$key\nLeft\t\t\tLength\tPenalty\n";

while ($primer=>getNextPair()) {

if (($primer=>penalty() > 9) && ($primer=>productGC() < 60)) {

print $primer=>leftOligo() ."\t" .
$primer=>rightOligo() . "\t".
$primer=>productLength(). "\t" .
$primer=>penalty() . "\n";
}
}
print "\n";

}

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

package primer

# import a few modules
use FileHandle;
use IPC::Open3;
# non-standard defaults: For any options not in
# this list, we use the generic primer3 default
# values.
my $primer_prog = '/usr/local/bin/primer';
my %definitions = (PRIMER_PRODUCT_SIZE_RANGE =>"200-450");

## constructor
sub new {
# setup variables
my ($class, $name, $seq) = @_;
my %results;

# create boulderIO command stream
my $command ="PRIMER_SEQUENCE_ID =$name\nSEQUENCE =$seq\n";
my $arg;
foreach $arg (keys(%definitions)){
$command .= "$arg=$definitions{$arg}\n";
}
$command .= "=\n";

# send command to primer program
my ($read,$write) = (FileHandle=>new,FileHandle=>new);
my $pid = open3($write,read, "$primer_prog");
print $write $command;
$write=>close;

# read results from primer into hash
while ($_=$ read=>getline) {
chomp;
my ($key,$value) = split("=");
$results{key}= $value;
}
return bless(\%results,$class);
}

## Iterator: set index to next primer pair
sub getnextpair {
my $self = shift;
if (exists $self=>{_CURRENT}){
# already counting
$self=>{_CURRENT}++;
if ($self=>{_CURRENT} <= $self=>count()){
# valid primer pair
return 1;
} else {
# we've run out of primer pairs
return 0;
}
} else {
# initialize counter
if ($self=>count()){
# we have primers to iterate over
$self=>{_CURRENT} = 0;
return 1;

} else {
# we have no primers to work on
return 0;
}
}
}

# Note for accessors: The first primer pair is officially
# primer pair 0. Each attribute in the hash is named with
# the primer pair number appended (or embedded) in the name
# *EXCEPT* for primer pair 0. Thus PRIMER_PAIR_PENALTY is
# the penalty for the first pair, PRIMER_PAIR_PENALTY_1 is
# the second, and so forth. This leads to the odd little
# if statement to append the underscore and number to every
# primer pair key except the first one.

## Accessor: Number of primer pairs
sub count {
my $self = shift;
return int(keys(%{$self})/20 - 1);
}

## Method: calculate GC content of product
sub productGC {
my $self = shift;

# find the start position of the current left primer
my ($leftKey) = "PRIMER_LEFT";
if ($self=>{_CURRENT}) {
$leftKey .= "_self=>{_CURRENT}";
}

# create substring for product
my ($start) = split(',', $self=>{$leftKey});
my $length = $self=>productlength();
my $prodSeq = substr($self=>{SEQUENCE}, $start, $length);

# calculate GC percentage rounded to two decimal places
my $GC = $prodseq =~ tr/GC/GC/;
return int(($GC/$length)*10000/100;
}

## Accessory: penalty
sub penalty {
my $self = shift;
# create hash key for current pair
my $key = "PRIMER_PAIR_PENALTY";
if ($self=>{_CURRENT}) {
$key .= " _$self=>{_CURRENT}";
}
return $self=>{$key};
}

## Accessor: product length
sub productlength {
my $self = shift;

# create hash key for current primer pair
my $key = "PRIMER_PRODUCT_SIZE";
if ($self=>{_CURRENT}) {
$key .= "_$self=>{_CURRENT}";
}
return $self=>{$key};
}

## Accessor: left oligo
sub leftOligo {
my $self = shift;

# create hash key for current primer pair
my $key = "PRIMER_LEFT";
if ($self=>{_CURRENT})";
$key .= "_self=>{_CURRENT}";
}
$key .= "_SEQUENCE";
return $self=>{$key};
}

# Accessor: right oligo
sub rightOligo {
my $self = shift;

# create hash key for current primer pair
my $key = "PRIMER_RIGHT";
if ($self=>{_CURRENT}){
$key .= "_$self=>{_CURRENT}";
}
$key .= "_SEQUENCE";
return $self=>{$key};
}
1;

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

PRIMER_SEQUENCE_ID=9622(KLK4) SEQUENCE=cctcctccctcagacccaggagtccagcccc
tcctccctcagacccaggagtccagaccccccagcccctc
ctccctcagacccaggggtccaggcccccaacccctcctc
cctcagactcagaggtccaagcccccaacccctccttccc
cagacccagaggtccaggtcccagcccctcctccctcaga
cccagcggtccaatgccacctaga
ctctccctgtacacagtgcccccttgtggcacgttgaccc
aaccttaccagttggtttttcattttttgtccctttcccc
tagatccagaaataaagtctaagagaagcgccaggttctg
tttgtagtgtcctccttgcttccccttccattcaagagat
tccagcaaatggctttccctccatgggtcgcatttcccta
tctgggaaatgcgagctttgaccaggagggccacagagtt
tcagagggaggggcc
attccctggatagttaagagaaaccagtgccatggatgtt
gtccagtttctgaaaa TARGET=100,377 PRIMER_COMMENT=Complete, or mostly-complete (> 300 bp), PSR Hit. PSR hit 5'-flush (within 100 bp of target) with target--100 5' bp of target allocated for forward primer. PSR in
formation used to specify selection of primers. PRIMER_PRODUCT_SIZE_RANGE=424-526 PRIMER_OPT_SIZE=21 PRIMER_MIN_SIZE=17 PRIMER_MAX_SIZE=24 PRIMER_OPT_TM=50.55 PRIMER_MIN_TM=34.4 PRIMER_MAX_TM=66.7
PRIMER_MISPRIMING_LIBRARY=/home/mpauley/bin/lib/human_mispriming_library.fas PRIMER_FILE_FLAG=0 PRIMER_EXPLAIN_FLAG=1 PRIMER_FIRST_BASE_INDEX=1 = PRIMER_SEQUENCE_ID=4633(MYL2) SEQUENCE=aggttgaccagatgttcgccgccttcccccc
tgacgtgactggaaacttggactacaagaacc
tggtgcacatcatcacccacggagaagagaaggactagga
gggggctcgctcgtggaccctgggctcgtctttgcagagt
ggtccctgccctcatctctctcccccgagtaccgcctctg
tccctaccttgtctgttagccatgtggctgccccatttat
ccacctccatcttctttgcagcctgggtggctatgggtac
ttcgtggccgcacatcctacagttggaaatccatccagag
gccatgttccaataa
acaggaggtcgtgtatttggtcacgacatttctctgacaa
acagtgccttgtgtaaggaaggcaaaggcagatgagctag
tgacagccttgcagtcattattattgaatcttggtgactc
tttatttaaggtagacctttgcagccactcagaactgttt
ctatcaaggatcacgttttcctggaatgttacctggatga
ccaaagccctctcttccaaagagctcaagcatttaga PRIMER_COMMENT=I
ncomplete (< 300 bp) or absent PSR hit. PSR information not used to specify selection of primers. PRIMER_PRODUCT_SIZE_RANGE=524-555 PRIMER_OPT_SIZE=21 PRIMER_MIN_SIZE=17 PRIMER_MAX_SIZE=24 PRIMER_OPT_TM=50.55
PRIMER_MIN_TM=34.4 PRIMER_MAX_TM=66.7 PRIMER_MISPRIMING_LIBRARY=/home/mpauley/bin/lib/human_mispriming_library.fas PRIMER_FILE_FLAG=0 PRIMER_EXPLAIN_FLAG=1 PRIMER_FIRST_BASE_INDEX=1

---------------------------------
Ready for the edge of your seat? Check out tonight's top picks on Yahoo! TV.

Sponsored Links