广告时间：如果觉得这篇文章对你有用，可以看下我的资料介绍，关注下。
背景
有时会需要针对整个基因组seqs.fa设计SSR引物，但是由于p3_in.pl 脚本提取序列的特点，每个SSR位点都会将整个序列存储一次，就连基因组稍小的拟南芥都需要超大的硬盘空间（测试的时候一条染色体序列没测试完，文件实在太大了，时间也有点久）。
有推文是利用misa.pl 文件提取SSR位点左右延长150bp 的bed文件，然后将SSR位点及侧翼序列用bedtools提取出来 bed.seqs.fa。再用MISA+Primer3的标准流程对bed.seqs.fa 文件进行SSR位点搜索与引物设计。这样，大大缩减了p3_in.pl 文件大小，节省了硬盘空间。
这里我写了个perl 脚本bug，以提取SSR位点及侧翼。
碎碎念
但是这样的方式（就常用的参数100bp-300bp产物长度，100bp内SSR视为复合SSR而言），提取150bp侧翼序列，可能会在额外的50bp处提取到重复的SSR位点，并可能涉及处100bp+的产物。所以，个人不建议用这个方法来统计SSR位点与引物设计成功率的问题。只是设计引物还是没问题的，注意先e-PCR验证即可。
脚本与使用
脚本需要指定seqs.fa.msia 文件和侧翼长度，侧翼长度可以不指定，默认150 bp。结果文件是一个四列的bed文件seqs.fa.bed 。
perl exact_bed.pl <seqs.fa.misa> <flank length(int,default:150)>
seqs.fa.bed 第一列为原序列ID；第二列、第三列为SSR及侧翼两端start和end，且考虑到了不超过原序列的边界；第四列是seqs.fa.msia 中前两列组合作为唯一name。
实际上使用了seqs.fa.misa 和seqs.fa 文件，不过其中seqs.fa 在参数中不直接体现，且需要两个文件位于同文件中（还是按原脚本的习惯写的）。
打印到屏幕的内容，第二列是考虑了序列起始长度后对第一列的矫正，第四列是考虑了序列总长度后对第三列的矫正（结果文件seqs.fa.bed 是考虑了序列长度的）。
有人可能注意到了，既然使用了seqs.fa 辅助产生bed文件，那为什么不直接提取出序列呢？因为……我现在还不会写。所以，还必需配合bedtools来跑流程。
下面是perl bug：（对于一个新手来说，听说写脚本和写bug是一样的）
#!/usr/bin/perl -w
# Program name: exact_bed.pl
## Description: creates a bed format file for fasta exactraction based on seqs.fa and seqs.fa.misa
## Usage: perl exact_bed.pl <seqs.fa.misa> <flank length(int,default:150)>
if (@ARGV == 0)
{
open (IN,"<$0");
while (<IN>) {if (/^\#\# (.*)/) {$message .= "$1\n"}};
close (IN);
die $message;
};
# Check if help is required #
if ($ARGV[0] =~ /-help|-h/i)
{
open (IN,"<$0");
while (<IN>) {if (/^\#\#\#(.*)/) {$message .= "$1\n"}};
close (IN);
die $message;
};
#open seqs.fa.misa
open (IN,"<$ARGV[0]") ||die ("\nError: Couldn't open misa.pl results file (*.misa) !\n\n");
#open seqs.fa
my $filename = $ARGV[0];
$filename =~ s/\.misa//;
open (FA, "<$filename") || die ("\nError: Couldn't open source file containing original FASTA sequences !\n\n");
#$filename = basename "$ARGV[1]";
#$filename =~ s/\.misa//;
#create output file: bed file
open (OUT, ">$filename.bed")|| die("\n$filename.bed Couldn't create! \n");
undef $/;
$in = <IN>;
study $in;
$/=">";
while(<FA>){
next unless (my ($id,$seq) = /(.*?)\n(.*)/s);
# print $id, "\n";
$seq =~ s/[\d\s]//g;
chomp $seq;
my $fa_len = length($seq)-1; #整理获取序列长度
# print "$id\t$fa_len\n";
while ($in =~ /$id\t(\d+)\t\S+\t\S+\t\d+\t(\d+)\t(\d+)/g){
my ($ssr_nr,$start,$end) = ($1,$2,$3);
my $flank = $ARGV[1]||150;
my $margin_left=$2-$flank-1;
print $margin_left,"\t";
if ($margin_left < 0){
$margin_left = 0;
}else{
$margin_left=$2-$flank-1;
};
print $margin_left,"\t";
my $margin_right = $3+$flank-1;
print $margin_right,"\t";
if ($margin_right > $fa_len){
$margin_right = $fa_len;
}else{
$margin_right = $3+$flank-1;
};
print $margin_right,"\n";
print OUT "$id\t$margin_left\t$margin_right\t$id"."_$ssr_nr\n";
}
}
close(IN);
close(FA);
close(OUT);
print "DONE: SSRs bed of $filename !\n";
参考：
小麦misa+primer3设计SSR引物,并用e-PCR进行筛选 :https://www.jianshu.com/p/d9a9b0f3f09f
网友评论