有时，我们需要判断一对引物对序列的匹配情况，或者更进一步，需要根据引物将匹配的序列片段切割出来。可以用下面的代码实现。
1.下载并安装全局比对软件mafft
软件地址为：https://mafft.cbrc.jp/alignment/software/linux.html
下载mafft_xxx_amd64.deb文件后运行
dpkg -i mafft_xxx_amd64.deb
2.需要的文件
a.被匹配的序列文件，fasta格式
例如，此处的测试文件Adenovirus_genome.fa
Adenovirus_genome.fa
碱基序列可以从NCBI上下载，一条碱基序列可多行书写。
b.引物序列文件，fasta格式，同时包含两条上下游引物
例如，此处的测试文件Adenovirus_primer.fasta
Adenovirus_primer.fasta
3.perl代码段
a.指定预设参数。包含被匹配的序列文件$fa，引物序列文件$ref，上游引物序列号$forid，下游引物序列号$revid，上下游引物允许错配数$cut_limit，引物匹配切割序列向引物上下游延伸长度。
其中，允许错配数是指插入、缺失和错配的总和。当引物实际匹配序列的错配数大于允许错配数时，该序列将被舍弃。
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.调整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序列和前后引物序列进行全局比对
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全局比对结果
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.统计前后引物比对情况，然后提取并切割满足条件的序列
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;
# 统计前后引物比对情况
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;
# 统计比对情况：插入、缺失和错配。注意排除简并碱基
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.结果
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%匹配，插入缺失碱基数为0，错配数为1；下游引物85%匹配，插入缺失碱基数为0，错配数为3。
