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 Bio::Tools::IUPAC and 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.
