Question: Using Bioperl To Match Primer Sites. Difficulty With Iupac.Pm
0
gravatar for Daniel
7.1 years ago by
Daniel3.8k
Cardiff University
Daniel3.8k wrote:

Hi All

I have been trying to write a script to take IUPAC compliment primers and pull out the corresponding subsequence from a full length (Aiming to select sequence regions from an amplicon).

The code works for direct matches but I am having trouble implementing the IUPAC module for generating a regexp on the fly. The code below only includes generating the forward primer regex at the moment. I also want to change the trimming portion to a substring argument but that's for when I've worked out how to find the character position of the match.

I have been mainly using http://doc.bioperl.org/bioperl-live/Bio/Tools/IUPAC.html and http://www.bioperl.org/wiki/Degenerate_primers but I have hit a brick wall, cant find any similar example scripts and would really appreciate any advice.

Thanks.

This is where I am so far:

#!/usr/bin/perl
use strict;
use warnings;

use Bio::Seq;
use Bio::SeqIO;
use Bio::Tools::IUPAC;

my $infile = $ARGV[0];
my $primer_f = $ARGV[1];
my $primer_r = $ARGV[2];

my $in = Bio::SeqIO->new(-file => "$infile", -format => 'fasta');
my $out = Bio::SeqIO->new(-file => ">out.fasta", -format => 'fasta');

#####
##Generate forward matching regex 
#####
my $ambiseq_f = Bio::Seq->new(-seq => "$primer_f", -alphabet => 'dna');

# Create all possible non-degenerate sequences
my $iupac = Bio::Tools::IUPAC->new(-seq => $ambiseq_f);
while (my $uniqueseq = $iupac->next_seq()) {
     # process the unique Bio::Seq object.    <---This is from the module documentation, Dont know what referring to.
        print $uniqueseq->seq . "\n";
}

my $regexp = $iupac->regexp();
#print "\n$regexp\n";

#####
##Trim to matching region
#####
while ( my $seq = $in->next_seq() ){
        my $trim5 = 0;
        my $trim53 = 0;

        if ( $seq->seq =~ m/^.*$primer_f(.*)/i ){
         $trim5 = $1;
        }

        if ($trim5 =~ m/^(.*)$primer_r.*/i ){
                $trim53 = $1;
        }
        #print ">" . $seq->id . "\n" . $trim53 . "\n";
        #$out->write_seq($trim53);
}

EDIT: removed an array - pretty sure that's not right but forgot to remove it before.

bioperl • 3.1k views
ADD COMMENTlink modified 5.4 years ago by Biostar ♦♦ 20 • written 7.1 years ago by Daniel3.8k
4
gravatar for Daniel
7.1 years ago by
Daniel3.8k
Cardiff University
Daniel3.8k wrote:

The issue was that I was reading the manual for the bioperl-live development package, then trying to use the stable version without realising.

For reference the development IUPAC.pm has a great method for generating degenerate primer regexs. Came in very useful.

ADD COMMENTlink written 7.1 years ago by Daniel3.8k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2177 users visited in the last hour