小型门户网站模板,百度推广竞价排名,最大的房产网站,南通装饰网站建设好的#xff0c;我们可以逐步分析这个 Perl 脚本的每个部分。脚本的主要功能是基于给定的 VCF 文件和参考基因组文件#xff0c;设计引物并进行电子 PCR#xff08;e-PCR#xff09;分析。我们将从脚本的头部和初始化部分开始讲解。
第一部分#xff1a;脚本头部和初始化…好的我们可以逐步分析这个 Perl 脚本的每个部分。脚本的主要功能是基于给定的 VCF 文件和参考基因组文件设计引物并进行电子 PCRe-PCR分析。我们将从脚本的头部和初始化部分开始讲解。
第一部分脚本头部和初始化
2. 模块加载
#use strict;
use Cwd qw(abs_path getcwd);
use Getopt::Long;
use Data::Dumper;
use FindBin qw($Bin $Script);
use File::Basename qw(basename dirname);
use Config::General;use File::Path qw(make_path remove_tree);
use File::Copy;
use File::Spec;
my $BEGIN_TIMEtime();
my $version1.0;
use Bio::SeqIO;use Cwd qw(abs_path getcwd);加载 Cwd 模块用于获取当前工作目录和绝对路径。use Getopt::Long;用于处理命令行参数。use Data::Dumper;用于调试可以打印复杂数据结构。use FindBin qw($Bin $Script);获取脚本所在的目录和脚本名称。use File::Basename qw(basename dirname);用于处理文件名和路径。use Config::General;用于解析配置文件。use File::Path qw(make_path remove_tree);用于创建和删除目录。use File::Copy;用于文件复制操作。use File::Spec;用于处理文件路径。use Bio::SeqIO;用于读取和处理生物序列文件如 FASTA 格式。
3. 命令行参数解析
my ($vcf,$flank,$fa,$epcrfa,$od,$outprefix,$PRIMER_NUM_RETURN,$PRIMER_OPT_SIZE,$PRIMER_MIN_SIZE,$PRIMER_MAX_SIZE,$PRIMER_PRODUCT_SIZE_RANGE,$PRIMER_MAX_NS_ACCEPTED);GetOptions(help|? \USAGE,fas\$fa,vcfs\$vcf,flanks\$flank,PRIMER_MIN_SIZEs\$PRIMER_MIN_SIZE,PRIMER_OPT_SIZEs\$PRIMER_OPT_SIZE,PRIMER_MAX_SIZEs\$PRIMER_MAX_SIZE,PRIMER_NUM_RETURNs\$PRIMER_NUM_RETURN,PRIMER_PRODUCT_SIZE_RANGEs\$PRIMER_PRODUCT_SIZE_RANGE,epcrfas\$epcrfa,ods\$od,prefixs\$outprefix,
) or USAGE;
USAGE unless($fa and $vcf);GetOptions解析命令行参数将参数值赋给相应的变量。 -fa输入的参考基因组文件FASTA 格式。-vcf输入的 VCF 文件。-flankindel 两侧的侧翼序列长度默认为 200。-PRIMER_MIN_SIZE、-PRIMER_OPT_SIZE、-PRIMER_MAX_SIZE引物的最小、最优和最大长度。-PRIMER_NUM_RETURN返回的引物对数量。-PRIMER_PRODUCT_SIZE_RANGE引物扩增产物的长度范围。-epcrfa用于电子 PCR 的参考基因组文件。-od输出目录默认为当前工作目录。-prefix输出文件的前缀默认为 indel。 USAGE如果用户未提供必要的参数-fa 和 -vcf则调用 USAGE 子程序显示帮助信息并退出。
4. 帮助信息
sub USAGE {my $usageUSAGE;
Program: $0
Version: $version
Contact: huangls
Discription:Usage:Options:-fa file Input fa file, forced-vcf file Input indel vcf file, forced-flank 200 cut 200 ; -PRIMER_MIN_SIZE 18 Minimum acceptable length of a primer. Must be greater than 0 and less than or equal to PRIMER_MAX_SIZE-PRIMER_OPT_SIZE 20 Optimum length (in bases) of a primer. Primer3 will attempt to pick primers close to this length.-PRIMER_MAX_SIZE 23 Maximum acceptable length (in bases) of a primer. Currently this parameter cannot be larger than 35. This limit is governed by maximum oligo size for which primer3s melting-temperature is valid.-PRIMER_PRODUCT_SIZE_RANGE 100-1000 The associated values specify the lengths of the product that the userwants the primers to create, and is a space separated list of elements of the form; -PRIMER_NUM_RETURN 1 Total num primer to search -epcrfa run E-PCR to filter multi mapped primers; -od str Key of output dir,default cwd;-prefix str defoult demo;-h HelpUSAGEprint $usage;exit 1;
}USAGE 子程序定义了脚本的使用说明包括所有支持的命令行参数及其说明。
总结
这一部分主要完成了以下任务
脚本头部信息和模块加载。命令行参数的解析。提供帮助信息。 第二部分默认参数设置和全局变量初始化
1. 设置默认参数值
$PRIMER_MIN_SIZE || 18;
$PRIMER_OPT_SIZE || 20;
$PRIMER_MAX_SIZE || 23;
$PRIMER_NUM_RETURN || 1;
$PRIMER_MAX_NS_ACCEPTED || 1;
$PRIMER_PRODUCT_SIZE_RANGE || 100-1000;
$od || getcwd;
$od abs_path($od);
unless (-d $od) { mkdir $od; }
$outprefix || indel;作用为未在命令行中指定的参数设置默认值。 PRIMER_MIN_SIZE引物的最小长度默认为 18。PRIMER_OPT_SIZE引物的最优长度默认为 20。PRIMER_MAX_SIZE引物的最大长度默认为 23。PRIMER_NUM_RETURN返回的引物对数量默认为 1。PRIMER_MAX_NS_ACCEPTED引物中允许的最大 N未知碱基数量默认为 1。PRIMER_PRODUCT_SIZE_RANGE引物扩增产物的长度范围默认为 100-1000。od输出目录默认为当前工作目录。outprefix输出文件的前缀默认为 indel。 getcwd 和 abs_path获取当前工作目录的绝对路径并确保输出目录存在。
2. 初始化全局变量
my %indel_length_range (mi 2, ma 5);
my %ref ();
my %ref_len ();%indel_length_range定义了 indel 长度的范围mi 表示最小长度ma 表示最大长度。%ref 和 %ref_len用于存储参考基因组序列及其长度。
3. 读取参考基因组文件
my $in Bio::SeqIO-new(-file $fa, -format Fasta);
while (my $seq $in-next_seq()) {$ref{$seq-id} $seq;$ref_len{$seq-id} $seq-length;
}
$in-close();Bio::SeqIO用于读取 FASTA 格式的参考基因组文件。$seq-id获取序列的 ID通常是染色体编号。$seq-length获取序列的长度。作用将参考基因组序列及其长度存储到 %ref 和 %ref_len 哈希表中方便后续使用。
总结
这一部分主要完成了以下任务
为未指定的参数设置默认值。初始化全局变量用于存储参考基因组序列及其长度。读取参考基因组文件并将其存储到哈希表中。 第三部分VCF 文件的读取和处理
1. 处理 VCF 文件路径
$fa abs_path($fa);
$epcrfa || $fa;
$epcrfa abs_path($epcrfa);作用 将输入的参考基因组文件路径$fa转换为绝对路径。如果未指定用于电子 PCR 的参考基因组文件$epcrfa则默认使用 $fa。确保 $epcrfa 也是绝对路径。
2. 打开 VCF 文件
if ($vcf ~ /gz$/) {open IN, gzip -dc $vcf| or die $! $vcf;
} else {open IN, $vcf or die $! $vcf;
}作用 检查 VCF 文件是否为 gzip 压缩文件通过文件扩展名 .gz 判断。如果是压缩文件使用 gzip -dc 命令解压并读取内容。如果不是压缩文件直接打开文件。如果文件无法打开程序将终止并报错。
3. 初始化变量
$flank || 200;
my %ssr_pos ();
my %SSR ();
my %genotype ();unless (-d $od/split) { mkdir $od/split; }
open OUT, $od/split/p_in_0.txt or die $!;
open SH, $od/primer3.sh or die $!;$flank设置 indel 两侧的侧翼序列长度默认为 200。%ssr_pos 和 %SSR用于存储 SSR简单序列重复的位置信息和 SSR 序列。%genotype用于存储每个 indel 的基因型信息。mkdir $od/split创建一个目录用于存放 Primer3 的输入文件。open OUT 和 open SH分别打开两个文件句柄用于写入 Primer3 的输入文件和 Shell 脚本。
4. 读取 VCF 文件内容
while (my $line IN) {chomp $line;next if $line ~ /^##/;# next if $line ~ /Scaffold/i;my tmp split(/\t/, $line);if ($line ~ /^#C/) {{$genotype{head}} tmp[7 .. $#tmp];next;}my ($ref_len, $alt_len) (length $tmp[3], length $tmp[4]);my $indel_len abs($ref_len - $alt_len);my $indel_len_s $alt_len - $ref_len;my $ID $tmp[0]_$tmp[1]_$indel_len;$ssr_pos{$ID} $tmp[0]\t$tmp[1]\t$tmp[3]\t$tmp[4]\t$indel_len_s;$tmp[7] ~ s/^.*ANNOVAR_DATE[^;];//;{$genotype{$ID}} tmp[7 .. $#tmp];my $seq ;my $indel_len_target $indel_len;if ($alt_len 1) {$seq get_target_fa($tmp[0], $tmp[1] - $flank, $tmp[1] $flank, $flank, $tmp[4]);$indel_len_target 1;} else {$seq get_target_fa($tmp[0], $tmp[1] - $flank, $tmp[1] $flank $indel_len, $flank);}作用 逐行读取 VCF 文件内容。跳过以 ## 开头的注释行。如果是列头行以 #C 开头提取样本信息并存储到 %genotype{head} 中。对于每个 indel 变异 计算参考碱基REF和变异碱基ALT的长度。计算 indel 的长度indel_len和符号indel_len_s。构造一个唯一的 ID$ID格式为 染色体编号_位置_indel长度。将 indel 的位置信息存储到 %ssr_pos 中。提取该 indel 的基因型信息并存储到 %genotype 中。调用 get_target_fa 函数获取 indel 两侧的侧翼序列$seq。如果 ALT 是插入变异长度 1将目标长度设置为 1否则目标长度为 indel 的实际长度。
5. 检测 SSR 并生成 Primer3 输入文件
if ($seq) {my ($is_have_ssr, $ssr) has_ssr($ID, $seq);if ($is_have_ssr) {$SSR{$ID} join(,, {$ssr});} else {$SSR{$ID} --;}print OUT SEQUENCE_ID$ID\n;print OUT SEQUENCE_TEMPLATE$seq\n;printf OUT (SEQUENCE_TARGET%d,%d\n, $flank 1, $indel_len_target);print OUT PRIMER_TASKpick_detection_primers\n;print OUT PRIMER_PICK_LEFT_PRIMER1\n;print OUT PRIMER_PICK_RIGHT_PRIMER1\n;print OUT PRIMER_NUM_RETURN$PRIMER_NUM_RETURN\n;print OUT PRIMER_MAX_END_STABILITY250\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_PRODUCT_OPT_SIZE125\n;print OUT PRIMER_OPT_GC_PERCENT50\n;print OUT PRIMER_MIN_TM50\n;print OUT PRIMER_OPT_TM55\n;print OUT PRIMER_MAX_TM65\n;print OUT PRIMER_MAX_NS_ACCEPTED1\n;print OUT PRIMER_FIRST_BASE_INDEX1\n;print OUT PRIMER_MIN_GC40.0\nPRIMER_MAX_GC60.0\n;print OUT PRIMER_PRODUCT_SIZE_RANGE$PRIMER_PRODUCT_SIZE_RANGE\n;printf OUT (SEQUENCE_INTERNAL_EXCLUDED_REGION%d,%d\n, $flank 1, $indel_len_target);print OUT \n;$n;if ($n 100) {$n 0;$num;close(OUT);open OUT, $od/split/p_in_$num.txt or die $!;print SH primer3_core -p3_settings_file/share/work/biosoft/primer3/latest/default_settings.txt -default_version1 -output$od/split/${num}primers.txt $od/split/p_in_$num.txt\n;}
}作用 如果成功获取了侧翼序列$seq则调用 has_ssr 函数检查该序列是否包含 SSR。如果存在 SSR将 SSR 信息存储到 %SSR 中否则存储 --。生成 Primer3 的输入文件内容包括序列 ID、模板序列、目标区域等信息。每处理 100 个 indel将 Primer3 输入文件分批保存到不同的文件中并生成对应的 Shell 命令行。这样可以避免单个 Primer3 输入文件过大提高运行效率。
总结
这一部分主要完成了以下任务
处理 VCF 文件路径确保文件可以正确读取。初始化相关变量用于存储 SSR 和基因型信息。逐行读取 VCF 文件提取 indel 信息并获取 indel 两侧的侧翼序列。检测 SSR并生成 Primer3 的输入文件。将 Primer3 输入文件分批保存生成对应的 Shell 命令行。 第四部分运行 Primer3 和电子 PCRe-PCR
1. 关闭文件句柄并生成 Primer3 脚本
close(OUT);
close(IN);
close(SH);#system(sh $od/primer3.sh);
system(parallel -j 10 $od/primer3.sh);
system(cat $od/split/*primers.txt $od/$outprefix.primers);作用 关闭所有打开的文件句柄。使用 parallel 命令并行运行 Primer3提高效率。-j 10 表示同时运行 10 个任务。将所有分批生成的 Primer3 输出文件合并为一个文件命名为 $outprefix.primers。
2. 解析 Primer3 输出文件
$/ \n;
my %Pseq ();open IN, $od/$outprefix.primers or die $!;
open OUT, $od/$outprefix.primers.xls or die $!;
print OUT #SEQUENCE_ID\tPRIMER_LEFT_SEQUENCE\tPRIMER_RIGHT_SEQUENCE\tPRIMER_PAIR_PRODUCT_SIZE\tPRIMER_LEFT_TM\tPRIMER_RIGHT_TM\tPRIMER_LEFT_GC_PERCENT\tPRIMER_RIGHT_GC_PERCENT\tSSR\tSSR_TYPE\n;
while (my $line IN) {chomp $line;my field split(/\n/, $line);my %primers get_hash(field);next if !defined($primers{PRIMER_PAIR_NUM_RETURNED}) or $primers{PRIMER_PAIR_NUM_RETURNED} 0;if ($primers{PRIMER_PAIR_NUM_RETURNED} 1) {$Pseq{$primers{SEQUENCE_ID}} [$primers{PRIMER_LEFT_0_SEQUENCE}, $primers{PRIMER_RIGHT_0_SEQUENCE}];print OUT join(\t, ($primers{SEQUENCE_ID}, $primers{PRIMER_LEFT_0_SEQUENCE}, $primers{PRIMER_RIGHT_0_SEQUENCE}, $primers{PRIMER_PAIR_0_PRODUCT_SIZE}, $primers{PRIMER_LEFT_0_TM}, $primers{PRIMER_RIGHT_0_TM}, $primers{PRIMER_LEFT_0_GC_PERCENT}, $primers{PRIMER_RIGHT_0_GC_PERCENT}, $SSR{$primers{SEQUENCE_ID}})) . \n;} else {for my $i (0 .. ($primers{PRIMER_PAIR_NUM_RETURNED} - 1)) {$Pseq{$primers{SEQUENCE_ID}.$i} [$primers{PRIMER_LEFT_${i}_SEQUENCE}, $primers{PRIMER_RIGHT_${i}_SEQUENCE}];print OUT join(\t, ($primers{SEQUENCE_ID}.$i, $primers{PRIMER_LEFT_${i}_SEQUENCE}, $primers{PRIMER_RIGHT_${i}_SEQUENCE}, $primers{PRIMER_PAIR_${i}_PRODUCT_SIZE}, $primers{PRIMER_LEFT_${i}_TM}, $primers{PRIMER_RIGHT_${i}_TM}, $primers{PRIMER_LEFT_${i}_GC_PERCENT}, $primers{PRIMER_RIGHT_${i}_GC_PERCENT}, $SSR{$primers{SEQUENCE_ID}})) . \n;}}
}
$/ \n;
close(OUT);作用 设置输入记录分隔符 $/ 为 \n以便正确解析 Primer3 的输出文件。打开 Primer3 输出文件和结果文件。遍历 Primer3 输出文件的每一部分以 分隔解析引物设计结果。调用 get_hash 函数将 Primer3 输出转换为哈希表。如果某个序列没有设计出引物则跳过。如果设计出了一对引物将其存储到 %Pseq 中并将相关信息写入结果文件。如果设计出了多对引物将每对引物的信息分别存储和写入结果文件。最后将输入记录分隔符恢复为默认值 \n。
3. 运行电子 PCRe-PCR
print e-PCR -w9 -f 1 -m100 $od/$outprefix.primers.xls D$PRIMER_PRODUCT_SIZE_RANGE $epcrfa N0 G0 T3 $od/$outprefix.epcr.out\n;
system(e-PCR -w9 -f 1 -m100 $od/$outprefix.primers.xls D$PRIMER_PRODUCT_SIZE_RANGE $epcrfa N0 G0 T3 $od/$outprefix.epcr.out);作用 构造 e-PCR 命令行运行 e-PCR 检查引物的特异性。参数说明 -w9允许最多 9 个错配。-f 1只考虑正链。-m100最大产物长度为 100。D$PRIMER_PRODUCT_SIZE_RANGE引物产物长度范围。$epcrfa用于 e-PCR 的参考基因组文件。N0不允许未定义的核苷酸N。G0不允许间隙。T3使用 3 个核苷酸的窗口进行匹配。 将 e-PCR 的输出保存到 $od/$outprefix.epcr.out 文件中。
4. 解析 e-PCR 输出文件
open IN, $od/$outprefix.epcr.out or die $!;
open OUT, $od/$outprefix.primers.result.xls or die $!;
my %P ();
while (my $line IN) {chomp $line;my tmp split(/\t/, $line);next if $tmp[6] $tmp[7] 1;$tmp[5] ~ s#/.*$##;$P{$tmp[1]}{$tmp[1]_$tmp[3]_$tmp[4]} [$tmp[1], $ssr_pos{$tmp[1]}, $Pseq{$tmp[1]}-[0], $Pseq{$tmp[1]}-[1], $tmp[5], tmp[6 .. $#tmp], {$genotype{$tmp[1]}}];
}
close(IN);
print OUT #PRIMER_ID\tCHROM\tPOS\tREF\tALT\tINDEL_LENGTH\tPRIMER_LEFT_SEQUENCE\tPRIMER_RIGHT_SEQUENCE\tPRIMER_PAIR_PRODUCT_SIZE\tMISM\tGAPS\tPRIMER_LEFT_TM\tPRIMER_RIGHT_TM\tPRIMER_LEFT_GC_PERCENT\tPRIMER_RIGHT_GC_PERCENT\tHAS_SSR\t . join(\t, {$genotype{head}}) . \n;for my $id (sort keys %P) {my ps (keys %{$P{$id}});if (ps 1) {print OUT join(\t, {$P{$id}{$ps[0]}}) . \n;}
}
close(OUT);作用 打开 e-PCR 输出文件和最终结果文件。遍历 e-PCR 输出文件的每一行解析引物的匹配结果。跳过匹配到多个位置的引物$tmp[6] $tmp[7] 1。提取引物的相关信息并存储到 %P 哈希表中。将最终结果写入 $od/$outprefix.primers.result.xls 文件包括引物 ID、染色体位置、引物序列、产物长度、错配数、间隙数、引物退火温度、GC 含量、是否含有 SSR 以及基因型信息。
总结
这一部分主要完成了以下任务
关闭文件句柄并生成 Primer3 脚本。解析 Primer3 输出文件提取引物设计结果。运行电子 PCRe-PCR检查引物的特异性。解析 e-PCR 输出文件生成最终的引物结果文件。 好的接下来我们分析脚本的第五部分生成 README 文件和清理临时文件。
第五部分生成 README 文件和清理临时文件
1. 生成 README 文件
open OUT, $od/readme.txt or die $!;
my $readme README;
*primers.result.xls
PRIMER_ID 引物ID 命名规则indel所在来源序列_indel所在位置start_indel所在位置长度
CHROM 所在染色体
POS 所在染色体位置
REF 参考基因组碱基
ALT 变异碱基
INDEL_LENGTH indel长度负数表示缺失正数表示插入
PRIMER_LEFT_SEQUENCE-- 左引物
PRIMER_RIGHT_SEQUENCE-- 右引物
PRIMER_PAIR_PRODUCT_SIZE-- 引物产物长度
MISM-- Total number of mismatches for two primers(ePCR)
GAPS-- Total number of gaps for two primers(ePCR)
PRIMER_LEFT_TM-- 退火温度
PRIMER_RIGHT_TM-- 退火温度
PRIMER_LEFT_GC_PERCENT-- GC含量
PRIMER_RIGHT_GC_PERCENT-- GC含量
HAS_SSR-- indel周围是否含有SSR -- 表示不含有SSR
INFO vcf注释信息 参照变异检测说明FORMAT vcf注释信息 参照变异检测说明 http://www.omicsclass.com/article/6 AD,Allelic depths for the ref and alt alleles in the order listedDP,Approximate read depth (reads with MQ255 or with bad mates are filtered)GQ,Genotype QualityGT,Genotype 0表示和参考基因型相同REF1表示变异基因型ALT。PL,Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification
样品基因型 信息用:隔开与FORMAT列对应方法
1.提取indel左右各$flank bp序列用primer3设计引物[2]
2.用Epcr 在参考基因组上电子PCR去除有多处比对的引物,保证引物扩增的特异性[3]
3.检测indel附近是否有SSR[1][1] MISA:Exploiting EST databases for the development and characterization of gene-derived SSR-markers in barley (Hordeum vulgare L.).
[2] Primer3: Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3 - new capabilities and interfaces. Nucleic Acids Research 40(15):e115
[3] Schuler, G. D. (1997). Sequence Mapping by Electronic PCR. Genome Research, 7(5), 541–550.READMEprint OUT $readme;
close(OUT);作用 打开一个文件句柄准备写入 README 文件。构造一个字符串 $readme包含对输出文件格式和内容的详细说明。写入 README 文件解释每个字段的含义包括引物 ID、染色体位置、引物序列、产物长度、错配数、间隙数、退火温度、GC 含量、SSR 信息以及 VCF 注释信息。说明了脚本运行的方法和引用的工具如 Primer3 和 e-PCR。关闭 README 文件。
2. 清理临时文件
print head -1 $od/$outprefix.primers.result.xls $od/head cat $od/$outprefix.primers.result.xls|grep -v # |sort -k2,2 -k3,3n $od/$outprefix.primers.result.xls.sorted cat $od/head $od/$outprefix.primers.result.xls.sorted $od/$outprefix.primers.result.xls rm $od/head $od/$outprefix.primers.result.xls.sorted\n;
system head -1 $od/$outprefix.primers.result.xls $od/head cat $od/$outprefix.primers.result.xls|grep -v # |sort -k2,2 -k3,3n $od/$outprefix.primers.result.xls.sorted cat $od/head $od/$outprefix.primers.result.xls.sorted $od/$outprefix.primers.result.xls rm $od/head $od/$outprefix.primers.result.xls.sorted;作用 使用 head 命令提取结果文件的第一行表头。使用 grep 命令过滤掉以 # 开头的注释行。使用 sort 命令按染色体编号和位置对结果进行排序。将排序后的结果重新组合为一个完整的文件。删除临时文件如 head 文件和排序后的临时文件。
总结
这一部分主要完成了以下任务
生成 README 文件详细说明输出文件的格式和内容。清理临时文件确保结果文件是有序且整洁的。 第六部分辅助函数的定义
1. get_hash 函数
sub get_hash {my info _;my %result ();for my $i (info) {next if $i ~ /^\s*$/;my ($k, $v) split(//, $i);$result{$k} $v;}return %result;
}作用 将 Primer3 的输出以 分隔的键值对解析为一个哈希表。遍历输入的每一行跳过空行。使用 split 函数将每一行按 分隔为键和值。将键值对存储到哈希表 %result 中。返回解析后的哈希表。
2. get_target_fa 函数
sub get_target_fa {my $id shift;my $posU shift;my $posD shift;if ($posU 0) {$posU 1;}my $seqobj $ref{$id};my $len $ref_len{$id};return $seqobj-subseq($posU, $len) if $posD $len;my $tri $seqobj-subseq($posU, $posD);return $tri;
}作用 从参考基因组中提取指定区域的序列。参数 $id染色体编号。$posU起始位置。$posD结束位置。 如果起始位置 $posU 小于等于 0则将其设置为 1。使用 Bio::Seq 对象 $seqobj 提取子序列。如果结束位置 $posD 超过了染色体长度 $len则提取从 $posU 到染色体末尾的序列。否则提取从 $posU 到 $posD 的序列。返回提取的序列。
3. has_ssr 函数
sub has_ssr {my ($id, $seq) _;my SSR;my %typrep (2 5,5 4,3 4,6 4,1 8,4 4);my $amb 200; # interruptions (max_difference_between_2_SSRs)my typ sort { $a $b } keys %typrep;my $max_repeats 1; # count repeatsmy $min_repeats 1000; # count repeatsmy (%count_motif, %count_class); # countmy ($number_sequences, $size_sequences, %ssr_containing_seqs, %count_motifs); # stores number and size of all sequences examinedmy $ssr_in_compound 0;my ($nr, %start, order, %end, %motif, %repeats); # store info of all SSRs from each sequence$seq ~ s/[\d\s]//g; # remove digits, spaces, line breaks, etc.$id ~ s/^\s*//g; $id ~ s/\s*$//g;#$id ~ s/\s/_/g; # replace whitespace with _$number_sequences;$size_sequences length $seq;for (my $i 0; $i scalar(typ); $i) { # check each motif classmy $motiflen $typ[$i];my $minreps $typrep{$typ[$i]} - 1;if ($min_repeats $typrep{$typ[$i]}) { $min_repeats $typrep{$typ[$i]}; }; # count repeatsmy $search (([acgt]{$motiflen})\\2{$minreps,});while ($seq ~ /$search/ig) { # scan whole sequence for that classmy $motif uc $2;my $redundant; # reject false type motifs [e.g. (TT)6 or (ACAC)5]for (my $j $motiflen - 1; $j 0; $j--) {my $redmotif ([ACGT]{$j})\\1{ . int($motiflen / $j - 1) . };$redundant 1 if ($motif ~ m/$redmotif/);}next if $redundant;$motif{$nr} $motif;my $ssr uc $1;$repeats{$nr} length($ssr) / $motiflen;$end{$nr} pos($seq);$start{$nr} $end{$nr} - length($ssr) 1;# count repeats$count_motifs{$motif{$nr}}; # counts occurrence of individual motifs$motif{$nr}{$repeats{$nr}}; # counts occurrence of specific SSR in its appearing repeat$count_class{$typ[$i]}; # counts occurrence in each motif classif ($max_repeats $repeats{$nr}) { $max_repeats $repeats{$nr}; };}}if (!$nr) {return (0, \SSR);} # no SSRs$ssr_containing_seqs{$nr};order sort { $start{$a} $start{$b} } keys %start; # put SSRs in right ordermy $i 0;my $count_seq; # countsmy ($start, $end, $ssrseq, $ssrtype, $size);while ($i $nr) {my $space $amb 1;if (!$order[$i 1]) { # last or only SSR$count_seq;my $motiflen length($motif{$order[$i]});$ssrtype p . $motiflen;$ssrseq ($motif{$order[$i]})$repeats{$order[$i]};$start $start{$order[$i]}; $end $end{$order[$i]};next;}if (($start{$order[$i 1]} - $end{$order[$i]}) $space) {$count_seq;my $motiflen length($motif{$order[$i]});$ssrtype p . $motiflen;$ssrseq ($motif{$order[$i]})$repeats{$order[$i]};$start $start{$order[$i]}; $end $end{$order[$i]};next;}my ($interssr);if (($start{$order[$i 1]} - $end{$order[$i]}) 1) {$count_seq; $ssr_in_compound;$ssrtype c*;$ssrseq ($motif{$order[$i]})$repeats{$order[$i]}($motif{$order[$i 1]})$repeats{$order[$i 1]}*;$start $start{$order[$i]}; $end $end{$order[$i 1]};} else {$count_seq; $ssr_in_compound;$interssr lc substr($seq, $end{$order[$i]}, ($start{$order[$i 1]} - $end{$order[$i]}) - 1);$ssrtype c;$ssrseq ($motif{$order[$i]})$repeats{$order[$i]}$interssr($motif{$order[$i 1]})$repeats{$order[$i 1]};$start $start{$order[$i]}; $end $end{$order[$i 1]};#$space - length $interssr}while ($order[$i 1] and (($start{$order[$i 1]} - $end{$order[$i]}) $space)) {if (($start{$order[$i 1]} - $end{$order[$i]}) 1) {$ssr_in_compound;$ssrseq . ($motif{$order[$i 1]})$repeats{$order[$i 1]}*;$ssrtype c*;$end $end{$order[$i 1]};} else {$ssr_in_compound;$interssr lc substr($seq, $end{$order[$i]}, ($start{$order[$i 1]} - $end{$order[$i]}) - 1);$ssrseq . $interssr($motif{$order[$i 1]})$repeats{$order[$i 1]};$end $end{$order[$i 1]};#$space - length $interssr}}$i;} continue {push SSR, $ssrseq;}my $has_ssr SSR;return ($has_ssr, \SSR);
}作用 检测给定序列中是否存在简单序列重复SSR。参数 $id序列 ID。$seq待检测的 DNA 序列。 使用正则表达式匹配不同类型的 SSR如单核苷酸重复、二核苷酸重复等。检测到 SSR 后记录其位置、重复次数和类型。如果检测到多个 SSR会检查它们之间的距离是否超过阈值$amb并相应地分类为复合 SSR 或简单 SSR。返回一个布尔值是否检测到 SSR和一个包含 SSR 信息的数组引用。
总结
这一部分定义了三个辅助函数
get_hash解析 Primer3 的输出将其转换为哈希表。get_target_fa从参考基因组中提取指定区域的序列。has_ssr检测序列中是否存在 SSR并返回相关信息。
这些函数在脚本的主逻辑中被多次调用用于处理数据和生成引物设计所需的输入。
