有時(shí)，我們需要判斷一對引物對序列的匹配情況，或者更進(jìn)一步，需要根據(jù)引物將匹配的序列片段切割出來?？梢杂孟旅娴拇a實(shí)現(xiàn)。

1.下載并安裝全局比對軟件mafft

軟件地址為：https://mafft.cbrc.jp/alignment/software/linux.html
下載mafft_xxx_amd64.deb文件后運(yùn)行

dpkg -i mafft_xxx_amd64.deb

2.需要的文件

a.被匹配的序列文件，fasta格式

例如，此處的測試文件Adenovirus_genome.fa

Adenovirus_genome.fa

堿基序列可以從NCBI上下載，一條堿基序列可多行書寫。

b.引物序列文件，fasta格式，同時(shí)包含兩條上下游引物

例如，此處的測試文件Adenovirus_primer.fasta

Adenovirus_primer.fasta

3.perl代碼段

a.指定預(yù)設(shè)參數(shù)。包含被匹配的序列文件$fa，引物序列文件$ref，上游引物序列號$forid，下游引物序列號$revid，上下游引物允許錯(cuò)配數(shù)$cut_limit，引物匹配切割序列向引物上下游延伸長度。

其中，允許錯(cuò)配數(shù)是指插入、缺失和錯(cuò)配的總和。當(dāng)引物實(shí)際匹配序列的錯(cuò)配數(shù)大于允許錯(cuò)配數(shù)時(shí)，該序列將被舍棄。

my $fa = "Adenovirus_genome.fa";
my $ref = "Adenovirus_primer.fasta";
my $forid = "AdV-F1";
my $revid = "AdV-R1";
my $cut_limit = 3;
my $move = 0;

b.調(diào)整fa文件格式，將每條fa的堿基序列合并為一行

$fa =~ /(.*?)\.[fasta]+$/;
my $name = $1;
open IN, "$fa";
open OUT, ">${name}_trans.fasta";
my $line1 = <IN>;
print OUT $line1;
while(my $line = <IN>){
    chomp $line;
    if($line =~ /^$/){next}
    if($line =~ /^>/){
        print OUT "\n$line\n";
    }else{
        print OUT "$line";
    }
}
close IN;
close OUT;

c.用mafft軟件分別對單一fa序列和前后引物序列進(jìn)行全局比對

open IN1, "${name}_trans.fasta";
my @ids;
unless(-d "tmp"){mkdir "tmp"}
while(my $line = <IN1>){
    chomp $line;
    if($line =~ /^>/){
        $line =~ /^>(\S+)/;
        my $id = $1;
        push @ids, $id;
        open OUT1, ">tmp/$id";
        my $fa = <IN1>;
        chomp $fa;
        print OUT1 "$line\n$fa\n";
        close OUT1;
        `cat $ref >> tmp/$id`;
        `mafft --auto tmp/$id > tmp/${id}_align`;
    }
}
close IN1;

d.整理mafft全局比對結(jié)果

if($forid !~ /^>/){$forid = ">$forid"};
if($revid !~ /^>/){$revid = ">$revid"};

foreach my $i(@ids){
    open IN2, "tmp/${i}_align";
    open OUT2, ">tmp/${i}_align_combine";
    my $line1 = <IN2>;
    print OUT2 $line1;
    while(my $line = <IN2>){
        chomp $line;
        if($line =~ /^$/){next}
        if($line =~ /^>/){
            print OUT2 "\n$line\n";
        }else{
            print OUT2 "$line";
        }
    }
    close IN2;
    close OUT2;
}

e.統(tǒng)計(jì)前后引物比對情況，然后提取并切割滿足條件的序列

open OUT3, ">${name}_cut.fasta";
foreach my $i(@ids){
    open IN3, "tmp/${i}_align_combine" or die "Cannot open tmp/${i}_align_combine: $!";
    my %id2fa;
    # 統(tǒng)計(jì)前后引物比對情況
    my($forlen, $revlen, $forend, $revbeg);
    my($forseq, $revseq, $seq, $foralign, $revalign);
    while(my $line = <IN3>){
        chomp $line;
        if($line =~ /$forid/){
            my $fa = <IN3>;
            chomp $fa;
            if($fa =~ /^(-+)/){$forlen = length($1)}else{$forlen = 0}
            if($fa =~ /^(.+?)-+$/){$forend = length($1)}else{$forend = length($fa)}
            $forseq = substr($fa, $forlen, $forend - $forlen);
        }elsif($line =~ /$revid/){
            my $fa = <IN3>;
            chomp $fa;
            if($fa =~ /(.*?)-+$/){$revlen = length($1)}else{$revlen = length($fa)}
            if($fa =~ /^(-+)/){$revbeg = length($1)}else{$revbeg = 0}
            $revseq = substr($fa, $revbeg, $revlen - $revbeg);
        }elsif($line =~ /^>/){
            my $fa = <IN3>;
            chomp $fa;
            $seq = $fa;
            $id2fa{$line} = $fa;
        }
    }
    close IN3;
    
    $foralign = substr($seq, $forlen, $forend - $forlen);
    $revalign = substr($seq, $revbeg, $revlen - $revbeg);
    my($foriden, $forgap, $formis) = &get_match_situation($foralign, $forseq);
    my($reviden, $revgap, $revmis) = &get_match_situation($revalign, $revseq);
    
    # 提取并切割滿足條件的序列
    foreach my $k(keys %id2fa){
        my $fa = $id2fa{$k};
        my $startf;
        if($forlen - $move < 0){$startf = 0}else{$startf = $forlen - $move}
        my $fan = substr($fa, $startf, $revlen-$forlen+2*$move);
        my $start = $startf + 1;
        $fan =~ s/-//g;
        my $falen = length($fan);
        my $end = $start + $falen - 1;
        if($forgap + $formis <= $cut_limit and $revgap + $revmis <= $cut_limit and $falen >= 50 + $move * 2){
            my $seqlen = $end - $start + 1;
            print OUT3 "$k ($start, $end, len: $seqlen) ($foriden, $forgap, $formis; $reviden, $revgap, $revmis)\n$fan\n";
        }
    }
}
close OUT3;

# 統(tǒng)計(jì)比對情況：插入、缺失和錯(cuò)配。注意排除簡并堿基
sub get_match_situation{
    my($seq1, $seq2) = @_[0,1];
    my $len = length($seq1);
    my $gap = 0;
    my $mismatch = 0;
    foreach my $i(0..$len-1){
        my $base1 = uc(substr($seq1, $i, 1));
        my $base2 = uc(substr($seq2, $i, 1));
        if($base1 eq "-" or $base2 eq "-"){
            $gap++;
        }elsif(($base2 =~ /[Rr]/ and $base1 =~ /[AaGg]/) or ($base2 =~ /[Yy]/ and $base1 =~ /[TtCc]/) or ($base2 =~ /[Ww]/ and $base1 =~ /[AaTt]/) or ($base2 =~ /[Ss]/ and $base1 =~ /[CcGg]/) or ($base2 =~ /[Mm]/ and $base1 =~ /[AaCc]/) or ($base2 =~ /[Kk]/ and $base1 =~ /[GgTt]/) or ($base2 =~ /[Bb]/ and $base1 =~ /[CcGgTt]/) or ($base2 =~ /[Dd]/ and $base1 =~ /[AaGgTt]/) or ($base2 =~ /[Hh]/ and $base1 =~ /[AaCcTt]/) or ($base2 =~ /[Vv]/ and $base1 =~ /[AaCcGg]/) or ($base2 =~ /[Nn]/ and $base1 =~ /[AaCcGgTt]/)){
            1;
        }elsif($base1 ne $base2){
            $mismatch++;
        }
    }
    my $match = $len - $gap - $mismatch;
    my $iden = sprintf "%.3f", $match / $len * 100;
    return ($iden, $gap, $mismatch);
}

4.結(jié)果

Adenovirus_genome_cut.fasta

“>AF083132.1 Porcine adenovirus 3 strain 6618, complete genome (20909, 21853, len: 945) (95.000, 0, 1; 85.000, 0, 3)”表示引物切割后的序列為原"AF083132.1 Porcine adenovirus 3 strain 6618, complete genome"的20909-21853號堿基片段，長度945bp；上游引物95%匹配，插入缺失堿基數(shù)為0，錯(cuò)配數(shù)為1；下游引物85%匹配，插入缺失堿基數(shù)為0，錯(cuò)配數(shù)為3。