misa+primer3设计SSR引物（一）之修改p3_in.p

不会perl，看过两集perl语言教程视频，因为句柄这个概念我无法理解，后面的视频就看不下去了。

本文灵感来自：https://www.jianshu.com/p/ce964f0dc710 中的6.3部分

1. 修改脚本目标

以小白的眼光修改p3_in.pl ，修改目的有两个：

• 便于修改下一步primer3的引物设计常用参数：产物长度，引物长度，引物GC含量，引物TM值等。

• 便于保存设计引物参数

其实运行primer3时可以指定参数文件，但是内容过多，不宜找到目标参数，而且每次修改不利于保存。

需要注意的是：

• 循环位置；变量写法；句末的";"；换行符"\n"；空格。

修改参数依据：

• 依据来自于配置文件 primer3-2.4.0/settings_files/p3_th_settings.txt

• 配置文件参数说明参考陈连福老师博客：primer3设计引物详解

2. 原脚本及修改部分的解析

脚本运行方式

$ perl p3_in.pl filename.fa.misa
 #输出文件自动命名为filename.fa.p3in
 #表面上看是需要一个输入文件，实际上还需要fasta序列文件filename.fa。
 #就是脚本中前两个open行。

因为不懂perl，也就不给大家乱解释全部脚本了。只从第19行开始胡说八道我们需要改的部分的代码内容。

#!/usr/bin/perl -w
# Author: Thomas Thiel
# Program name: primer3_in.pl
# Description: creates a PRIMER3 input file based on SSR search results

open (IN,"<$ARGV[0]") || die ("\nError: Couldn't open misa.pl results file (*.misa) !\n\n");

my $filename = $ARGV[0];
$filename =~ s/\.misa//;
open (SRC,"<$filename") || die ("\nError: Couldn't open source file containing original FASTA sequences !\n\n");
open (OUT,">$filename.p3in");

undef $/;
$in = <IN>;
study $in;

$/= ">";

my $count;  #从这个开始是我们需要改的内容。
while (<SRC>)
  {
  next unless (my ($id,$seq) = /(.*?)\n(.*)/s);
  $seq =~ s/[\d\s>]//g;#remove digits, spaces, line breaks,...
  while ($in =~ /$id\t(\d+)\t\S+\t\S+\t(\d+)\t(\d+)/g)
    {
    my ($ssr_nr,$size,$start) = ($1,$2,$3);
    $count++;
    print OUT "SEQUENCE_ID=$id"."_$ssr_nr\nSEQUENCE_TEMPLATE=$seq\n"; 
    print OUT "PRIMER_PRODUCT_SIZE_RANGE=100-280\n";
    print OUT "TARGET=",$start-3,",",$size+6,"\n";
    print OUT "PRIMER_MAX_END_STABILITY=250\n=\n"
    };
  };
print "\n$count records created.\n"; #

网页端可见代码的行号
line28-31是循环内内容

line28：这一行的输出其实是两行。

\n 表示为linux下的换行符。

第一行为：SEQUENCE_ID=$id"."_$ssr_nr 是将filename.misa文件中的第一列ID与第二列中ssr nr进行拼接。其"." 是将两部分进行“无缝拼接”。

第二行为：SEQUENCE_TEMPLATE=$seq是上一行中$id在fasta序列文件中的整条序列。这意味当我们对较长的序列，如染色体序列设计SSR引物时，只要该染色体序列出现一次SSR（描述稍不太准，以为还可能是复合SSR），该染色体序列就会被打印一次。这也是我们设计引物时，有时filename.fa.p3in会占用超大硬盘空间的原因。

line29：设置产物长度

line30：不是很明白其内容含义，不明白没事儿，不改这个引物设计参数就可以了。

猜测应该是针对设计引物的区段。但是要注意$start与$size对应数值为filename.fa.misa中的列。也就是说，这一部分是跟每条序列有关的。

line31：上下游引物5bp 3' 碱基最大的稳定性。

不是很理解数值是如何确定碱基的，还是不修改就可以了。

必须注意的是末尾的\n=\n ，意思是打印完250 就换行，然后再打印= 后换行。也就是说这个小循环是以"="结尾的，否则下一步会报错。

line34在循环外：

line34：统计处理的记录，输出到标准输出，可以理解为打印到屏幕上。

都说到这儿了，确定不试着自己改改？

3. 第一次修改：设定引物设计参数

完成设定引物设计常用的引物长度范围，GC含量范围，TM值范围。

修改范围：line28-31。

    print OUT "SEQUENCE_ID=$id"."_$ssr_nr\nSEQUENCE_TEMPLATE=$seq\n";
    print OUT "PRIMER_PRODUCT_SIZE_RANGE=100-280\n"; #产物长度
    print OUT "PRIMER_MIN_SIZE=18\n"; #引物长度
    print OUT "PRIMER_OPT_SIZE=20\n";
    print OUT "PRIMER_MAX_SIZE=23\n";
    print OUT "PRIMER_MIN_TM=57.0\n"; #引物TM范围
    print OUT "PRIMER_OPT_TM=59.0\n";
    print OUT "PRIMER_MAX_TM=62.0\n";
    print OUT "PRIMER_MIN_GC=30.0\n"; #引物GC含量
    print OUT "PRIMER_MAX_GC=70.0\n";
    print OUT "TARGET=",$start-3,",",$size+6,"\n";
    print OUT "PRIMER_MAX_END_STABILITY=250\n=\n"

注意：引物参数来自文件p3_th_settings.txt 。

需要注意的是：

• 该参数文件中所有 = 两侧均无空格。所以，我们在改的时候需要注意这点。(不过这步我没测试= 两侧加上空格会不会有报错)。

• 写法模仿原脚本line28中产物长度的写法，注意转行与";"

• 另，再次强调每个小循环需要以= 结尾，否则会报错。

这样改完后，输出的filename.fa.p3in 文件中，每个需要设计SSR位点的位置都会紧跟一遍引物设定参数。

但是，此时还未能将设置的参数打印到标准输出。

4. 将引物设定的参数打印到标准输出

产物长度范围，引物长度范围，引物TM范围，引物GC含量范围，整体来说是固定的。没必要每个引物设计次打印一遍，所以要在循环外打印设置的参数。同时为了以后修改参数方便，我们在循环内的参数（即第一次修改内容）之前定义好参数值，循环外打印参数到标准输出。

从原脚本的line 19行开始修改。

my $count;
 my $PRIMER_PRODUCT_SIZE_RANGE_left = 100; #产物长度
 my $PRIMER_PRODUCT_SIZE_RANGE_right = 280;
 my $PRIMER_MIN_SIZE = 18;#引物大小
 my $PRIMER_OPT_SIZE = 20;
 my $PRIMER_MAX_SIZE = 23;
 my $PRIMER_MIN_TM = 57.0;#引物TM范围
 my $PRIMER_OPT_TM = 59.0;
 my $PRIMER_MAX_TM = 62.0;
 my $PRIMER_MIN_GC = 30.0;#引物GC含量
 my $PRIMER_MAX_GC = 70.0;
 while (<SRC>)
  {
  next unless (my ($id,$seq) = /(.*?)\n(.*)/s);
  $seq =~ s/[\d\s>]//g;#remove digits, spaces, line breaks,...
  while ($in =~ /$id\t(\d+)\t\S+\t\S+\t(\d+)\t(\d+)/g)
  {
  my ($ssr_nr,$size,$start) = ($1,$2,$3);
  $count++;
  print OUT "SEQUENCE_ID=$id"."_$ssr_nr\nSEQUENCE_TEMPLATE=$seq\n";
  print OUT "PRIMER_PRODUCT_SIZE_RANGE=$PRIMER_PRODUCT_SIZE_RANGE_left"."-$PRIMER_PRODUCT_SIZE_RANGE_right\n"; #产物长度
  print OUT "PRIMER_MIN_SIZE=$PRIMER_MIN_SIZE\n"; #引物大小
  print OUT "PRIMER_OPT_SIZE=$PRIMER_OPT_SIZE\n";
  print OUT "PRIMER_MAX_SIZE=$PRIMER_MAX_SIZE\n";
  print OUT "PRIMER_MIN_TM=$PRIMER_MIN_TM\n"; #引物TM范围
  print OUT "PRIMER_OPT_TM=$PRIMER_OPT_TM\n";
  print OUT "PRIMER_MAX_TM=$PRIMER_MAX_TM\n";
  print OUT "PRIMER_MIN_GC=$PRIMER_MIN_GC\n"; #引物GC含量
  print OUT "PRIMER_MAX_GC=$PRIMER_MAX_GC\n";
  print OUT "TARGET=",$start-3,",",$size+6,"\n";
  print OUT "PRIMER_MAX_END_STABILITY=250\n=\n"
  };
  };
 print "\n$count records created.\n"; #以下内容是为了可以将参数输出到标准输出，利于重定向保存操作的参数
 print "PRIMER_PRODUCT_SIZE_RANGE=$PRIMER_PRODUCT_SIZE_RANGE_left"."-$PRIMER_PRODUCT_SIZE_RANGE_right\n"; #产物长度
 print "PRIMER_MIN_SIZE=$PRIMER_MIN_SIZE\n"; #引物大小
 print "PRIMER_OPT_SIZE=$PRIMER_OPT_SIZE\n";
 print "PRIMER_MAX_SIZE=$PRIMER_MAX_SIZE\n";
 print "PRIMER_MIN_TM=$PRIMER_MIN_TM\n"; #引物TM范围
 print "PRIMER_OPT_TM=$PRIMER_OPT_TM\n";
 print "PRIMER_MAX_TM=$PRIMER_MAX_TM\n";
 print "PRIMER_MIN_GC=$PRIMER_MIN_GC\n"; #引物GC含量
 print "PRIMER_MAX_GC=$PRIMER_MAX_GC\n";
 #print "TARGET=",$start-3,",",$size+6,"\n"; #跟每一条序列有关，是循环内变量，不打印。
 print "PRIMER_MAX_END_STABILITY=250\n"

修改需要注意的是：

• 产物长度定义，如果直接写为my $PRIMER_PRODUCT_SIZE_RANGE = 100-280 ，在脚本中，会让其误以为这是个减号，输出文件filename.fa.p3in 会是-180 。

  所以这里的书写方式模仿原脚本中line28行，用"."拼接两端范围值。

• print OUT "TARGET=",$start-3,",",$size+6,"\n"; 行跟每个序列有关，位于循环内，意义不大，不做记录。

• 此外，打印标准输出时，可以尝试用"."分别 拼接引物大小，TM范围，GC含量的范围值，使其更加易读。

• 此时运行方式：

$ perl p3_in.pl filename.fa.misa > p3_in.log
 #将引物设计条件进行重定向保存。

5. 小结与讨论

在运行该脚本时修改、记录参数更加方便了。但是缺点也显而易见，输出文件filename.fa.p3in 会更加占用硬盘空间了。

这一缺点，对于短序列还不明显，对长序列较为明显。这一缺点也可以通过上游用bedtools提取序列减轻硬盘负担。这一步，为什么网上教程取左右各150长度，是为了更好的设计引物？如果p3_in.pl中TARGET 是设定引物的区间，那么将150改为3，6不行？另外，这一步中如果序列范围溢出该如何操作？还需要注意，fasta序列中1-100的序列位置，在bed格式文件中要写为0-99。这一步暂时还用不到，等忙完小论文再回来思考、摸索这个问题。