Как я могу получить правильный счет для точных совпадений FASTA?

Question

Я должен искать шаблон wTTTAYRTTTW, где W знак равно A или же T, Y знак равно C или же T, R знак равно A или же R, в файле FASTA последовательности генома. Должно быть разрешено некоторое несоответствие, то есть строка с точным соответствием и их положение. Мой подход:

#!user/bin/perl -w
use strict;
my @name = qw( NC_004314.2_10 );
for (my $i = 0 ; $i< scalar(@name);$i++)
{
my $fname = $name[$i];
print "$fname\n";
my $read_pat= "WTTTAYRTTTW";
#print"\nPlease enter how many mismatch is allowd : ";
my $m =<STDIN>;
chomp $m;
unless(open(fh1, "$fname")){
    print "Cannot open file \"$fname\"\n\n";
    exit;
                }
my @fh=<fh1>;
close fh1;
if ($fh[0] !~ /^>/)
    {
        print "not fasta file\n";
        exit;
    }
my $seq='';
foreach my $line(@fh)
    {
        if($line =~ /^>/)
        {
        next;
        }
        else
        {
        $seq=$seq.$line;
        }
    }
sub trans_pat
    {
     my $pat=shift;
     $pat=~s/R/\[CG\]/g;
        $pat=~s/W/\[AT\]/g;
        $pat=~s/Y/\[AG\]/g;
    return $pat;
    }
open(FH1,">$fname.csv");
sub find_pat
{
my ($pat,$seq) = (@_);
#print FH1 "Looking for pattern $pat\n";
}

find_pat (trans_pat($read_pat),$seq);

# Allowing for a single mismatch

my $pat=trans_pat($read_pat);
print FH1 "Looking for pattern $pat\n";
while ($seq=~m/(?=$pat)/g)
{
print FH1"match at\t$-[0]\t$&\n"
}
foreach my $i (1..(length $read_pat)-($m-1))
{
my $mis_pat = $read_pat;
substr($mis_pat,$i-1,$m)=".{$m}";
my $pat1=trans_pat($mis_pat);
print FH1 "Looking for pattern $pat1\n";
while ($seq=~m/(=?$pat1)/g)

{
print FH1 "match at\t$-[0]\t$&\n";
}
#print FH1"$& \n";
}
close FH1;

Результат, найденный этим кодом, неверен для данной последовательности NC_004314.2 в файле FASTA для точного соответствия, общее количество совпадений должно быть 829. Может ли кто-нибудь исправить этот код?

Answer 1

Я ничего не знаю о поиске генома и не могу понять вашу программу быстро, но как совет:

Вы действительно должны сканировать:

• http://www.bioperl.org/ и особенно
• http://www.bioperl.org/wiki/HOWTO:Getting_Genomic_Sequences

Может быть, помогает, и вы можете получить некоторые идеи, как делать вещи.

Perl - это язык биоинформатики.