indel_snp_ssr_primer

好的，我们可以逐步分析这个 Perl 脚本的每个部分。脚本的主要功能是基于给定的 VCF 文件和参考基因组文件，设计引物并进行电子 PCR（e-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_TIME=time();
my $version="1.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,"fa=s"=>\$fa,"vcf=s"=>\$vcf,"flank=s"=>\$flank,"PRIMER_MIN_SIZE=s"=>\$PRIMER_MIN_SIZE,"PRIMER_OPT_SIZE=s"=>\$PRIMER_OPT_SIZE,"PRIMER_MAX_SIZE=s"=>\$PRIMER_MAX_SIZE,"PRIMER_NUM_RETURN=s"=>\$PRIMER_NUM_RETURN,"PRIMER_PRODUCT_SIZE_RANGE=s"=>\$PRIMER_PRODUCT_SIZE_RANGE,"epcrfa=s"=>\$epcrfa,"od=s"=>\$od,"prefix=s"=>\$outprefix,
) or &USAGE;
&USAGE unless($fa and $vcf);
  • GetOptions：解析命令行参数，将参数值赋给相应的变量。
    • -fa：输入的参考基因组文件（FASTA 格式）。
    • -vcf：输入的 VCF 文件。
    • -flank：indel 两侧的侧翼序列长度，默认为 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 $usage=<<"USAGE";
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 primer3's 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. 脚本头部信息和模块加载。
  2. 命令行参数的解析。
  3. 提供帮助信息。

第二部分：默认参数设置和全局变量初始化

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 哈希表中，方便后续使用。

总结

这一部分主要完成了以下任务：

  1. 为未指定的参数设置默认值。
  2. 初始化全局变量，用于存储参考基因组序列及其长度。
  3. 读取参考基因组文件，并将其存储到哈希表中。

第三部分：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_TASK=pick_detection_primers\n";print OUT "PRIMER_PICK_LEFT_PRIMER=1\n";print OUT "PRIMER_PICK_RIGHT_PRIMER=1\n";print OUT "PRIMER_NUM_RETURN=$PRIMER_NUM_RETURN\n";print OUT "PRIMER_MAX_END_STABILITY=250\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_SIZE=125\n";print OUT "PRIMER_OPT_GC_PERCENT=50\n";print OUT "PRIMER_MIN_TM=50\n";print OUT "PRIMER_OPT_TM=55\n";print OUT "PRIMER_MAX_TM=65\n";print OUT "PRIMER_MAX_NS_ACCEPTED=1\n";print OUT "PRIMER_FIRST_BASE_INDEX=1\n";print OUT "PRIMER_MIN_GC=40.0\nPRIMER_MAX_GC=60.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_version=1 -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 输入文件过大，提高运行效率。

总结

这一部分主要完成了以下任务：

  1. 处理 VCF 文件路径，确保文件可以正确读取。
  2. 初始化相关变量，用于存储 SSR 和基因型信息。
  3. 逐行读取 VCF 文件，提取 indel 信息，并获取 indel 两侧的侧翼序列。
  4. 检测 SSR，并生成 Primer3 的输入文件。
  5. 将 Primer3 输入文件分批保存，生成对应的 Shell 命令行。

第四部分：运行 Primer3 和电子 PCR（e-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. 运行电子 PCR（e-PCR）
print "e-PCR -w9 -f 1 -m100 $od/$outprefix.primers.xls D=$PRIMER_PRODUCT_SIZE_RANGE $epcrfa N=0 G=0 T=3 > $od/$outprefix.epcr.out\n";
system("e-PCR -w9 -f 1 -m100 $od/$outprefix.primers.xls D=$PRIMER_PRODUCT_SIZE_RANGE $epcrfa N=0 G=0 T=3 > $od/$outprefix.epcr.out");
  • 作用：
    • 构造 e-PCR 命令行，运行 e-PCR 检查引物的特异性。
    • 参数说明：
      • -w9：允许最多 9 个错配。
      • -f 1：只考虑正链。
      • -m100：最大产物长度为 100。
      • D=$PRIMER_PRODUCT_SIZE_RANGE：引物产物长度范围。
      • $epcrfa：用于 e-PCR 的参考基因组文件。
      • N=0：不允许未定义的核苷酸（N）。
      • G=0：不允许间隙。
      • T=3：使用 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 以及基因型信息。

总结

这一部分主要完成了以下任务：

  1. 关闭文件句柄并生成 Primer3 脚本。
  2. 解析 Primer3 输出文件，提取引物设计结果。
  3. 运行电子 PCR（e-PCR），检查引物的特异性。
  4. 解析 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 listed"DP,"Approximate read depth (reads with MQ=255 or with bad mates are filtered)"GQ,"Genotype Quality"GT,"Genotype"   0表示和参考基因型相同（REF），1表示变异基因型（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 文件和排序后的临时文件）。

总结

这一部分主要完成了以下任务：

  1. 生成 README 文件，详细说明输出文件的格式和内容。
  2. 清理临时文件，确保结果文件是有序且整洁的。

第六部分：辅助函数的定义

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 信息的数组引用。

总结

这一部分定义了三个辅助函数：

  1. get_hash：解析 Primer3 的输出，将其转换为哈希表。
  2. get_target_fa：从参考基因组中提取指定区域的序列。
  3. has_ssr：检测序列中是否存在 SSR，并返回相关信息。

这些函数在脚本的主逻辑中被多次调用，用于处理数据和生成引物设计所需的输入。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.pswp.cn/news/907299.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

2.4GHz 射频前端芯片AT2401C

2.4GHz 射频前端芯片AT2401C

射频前端芯片作为无线通信系统的核心组件&#xff0c;涵盖功率放大器&#xff08;PA&#xff09;、滤波器、开关、低噪声放大器&#xff08;LNA&#xff09;等关键器件&#xff0c;其性能直接影响通信质量、功耗及信号稳定性。 AT2401C是一款面向 Zigbee&#xff0c;无线传感网…
阅读更多...
Batch Normalization[[

Batch Normalization[[

error surface如果很崎岖,那么就代表比较难train,我们有没有办法去改变这个landscape呢 可以用batch normalization. 如果 ( x_1 ) 的取值范围很小&#xff08;如 1, 2&#xff09;&#xff0c;而 ( x_2 ) 的取值范围很大&#xff08;如 100, 200&#xff09;&#xff0c;那么…
阅读更多...
c++结构化绑定

c++结构化绑定

author: hjjdebug date: 2025年 05月 28日 星期三 15:57:58 CST descrip: c结构化绑定: 结构化绑定: 名称辨析: 名称叫绑定好还是叫解绑好&#xff1f; 解绑意思是原来是一个整体,现在被分成了若干个部分,所以叫解. 绑定强调的意思是. 被分解的某个变量,绑定到了整体的某个变量…
阅读更多...
大数据治理：理论、实践与未来展望(一）

大数据治理：理论、实践与未来展望(一）

文章目录 一、大数据治理的定义与重要性&#xff08;一&#xff09;定义&#xff08;二&#xff09;重要性 二、大数据治理的应用场景&#xff08;一&#xff09;金融行业&#xff08;二&#xff09;医疗行业&#xff08;三&#xff09;制造业&#xff08;四&#xff09;零售行…
阅读更多...
AI系统化学习月计划6月计划

AI系统化学习月计划6月计划

以下是为技术总监设计的 AI系统化学习月计划&#xff08;每天投入2小时&#xff0c;共30天&#xff09;&#xff0c;结合战略思维、技术基础、实战应用和行业趋势&#xff0c;帮助您快速掌握AI的核心知识&#xff0c;并转化为业务决策能力。 第一周&#xff1a;AI基础与战略思维…
阅读更多...
详解MySQL调优

详解MySQL调优

目录 1. SQL 语句优 1.1 避免低效查询 1.2 索引优化 1.3 分析执行计划 2. 数据库配置优化 2.1 核心参数调整 2.2 表结构与存储引擎 2.3 存储引擎选择 3. 事务与锁优化 3.1 事务控制 3.2 锁机制优化 3.3 批量操作优化 4. 其他优化手段 4.1 监控与分析工具 4.2 读写…
阅读更多...
VScode单双引号、分号格式

VScode单双引号、分号格式

1、settings.json中添加&#xff1a; 1 2 3 "prettier.semi": false, // 取消自动加分号 "prettier.singleQuote": true, // 保持单引号&#xff0c;不自动变双引号 "prettier.trailingComma": "none" // 去掉结尾的逗号 2、如上一步…
阅读更多...
自动驾驶规划控制教程——不确定环境下的决策规划

自动驾驶规划控制教程——不确定环境下的决策规划

引言：驾驭未知——不确定性下的自动驾驶决策挑战 自动驾驶汽车 (Autonomous Vehicles, AVs) 的愿景是彻底改变交通运输的面貌，提高道路安全、提升交通效率、改善驾乘体验。然而，要将这一愿景安全可靠地付诸实践，自动驾驶系统必须能够在复杂、动态且充满不确定性的真实世界…
阅读更多...
电缆中性点概念

电缆中性点概念

电缆中性点概念 电缆中性点（也称“中性点”或“中性线”）是电力系统和电气设备中一个非常重要的概念，尤其在三相电系统中。下面是对中性点概念的系统性解释。 1. 基本定义 中性点：三相电缆（A/B/C相）的电压矢量交汇点，理想情况下三相平衡时该点电压为零。对于星形（Y形…
阅读更多...
MyBatis 动态 SQL 详解：灵活构建强大查询

MyBatis 动态 SQL 详解：灵活构建强大查询

MyBatis 的动态 SQL 功能是其最强大的特性之一&#xff0c;它允许开发者根据不同条件动态生成 SQL 语句&#xff0c;极大地提高了 SQL 的灵活性和复用性。本文将深入探讨 MyBatis 的动态 SQL 功能&#xff0c;包括 OGNL 表达式的使用以及各种动态 SQL 元素&#xff08;如 if、c…
阅读更多...
嵌入式自学第三十天（5.28）

嵌入式自学第三十天（5.28）

&#xff08;1&#xff09;多线程资源竞争问题&#xff1a; 互斥&#xff1a;在多线程中对临界资源的排他性访问。 解决方案&#xff1a;互斥锁 mutex互斥锁在进程pcb块&#xff0c;ret 为0说明别人在用&#xff0c;1说明空闲。 阻塞锁 man pthread_mutex_init man pthread_…
阅读更多...
【HW系列】—web常规漏洞（SQL注入与XSS）

【HW系列】—web常规漏洞（SQL注入与XSS）

SQL注入与XSS攻防解析&#xff08;安全防御指南&#xff09; 一、SQL注入基础&#xff08;防御视角&#xff09; ​​1. 简介​​ SQL注入是一种通过构造非预期SQL语句操纵数据库的攻击技术。作为开发者&#xff0c;需重点关注输入验证与查询安全&#xff0c;建立全流量监测…
阅读更多...
Accelerate 2025北亚巡展正式启航！AI智御全球·引领安全新时代

Accelerate 2025北亚巡展正式启航！AI智御全球·引领安全新时代

近日&#xff0c;网络安全行业年度盛会Accelerate 2025北亚巡展正式在深圳启航&#xff01;智库专家、产业领袖及Fortinet高管、产品技术团队和300余位行业客户齐聚一堂&#xff0c;围绕“AI智御全球引领安全新时代”主题&#xff0c;共同探讨AI时代网络安全新范式。大会聚焦三…
阅读更多...
RAG系统构建之嵌入模型性能优化完整指南

RAG系统构建之嵌入模型性能优化完整指南

导读&#xff1a;在企业级RAG系统的实际部署中&#xff0c;您是否遇到过这样的困扰&#xff1a;嵌入计算成本不断攀升&#xff0c;API调用频繁触及限制&#xff0c;而系统响应速度却始终达不到用户期望&#xff1f;这些看似分散的问题&#xff0c;实际上都指向同一个技术核心&a…
阅读更多...
python 自动生成不同行高的word

python 自动生成不同行高的word

python 自动生成不同行高的word # -*- coding: utf-8 -*- from docx import Document from docx.shared import Cm, Pt, Inches from docx.oxml import OxmlElement from docx.oxml.ns import qn from docx.enum.text import WD_ALIGN_PARAGRAPHclass DynamicTableGenerator:d…
阅读更多...
如何训练意志力

如何训练意志力

设定清晰的目标 目标需要是具体的&#xff0c;可实现的&#xff0c;有时间限制的。比如不要说“我要锻炼”&#xff0c;而是改成“每周跑步3次&#xff0c;每次30分钟”。 从小事开始 起步通常都是困难的&#xff0c;一开始定一个很大很复杂的任务也超出了自己的能力&#x…
阅读更多...
FastAPI 依赖注入

FastAPI 依赖注入

依赖注入常用于以下场景&#xff1a; 共享业务逻辑&#xff08;复用相同的代码逻辑&#xff09; 共享数据库连接 实现安全、验证、角色权限 等…… 上述场景均可以使用依赖注入&#xff0c;将代码重复最小化。 创建依赖项 依赖项就是一个函数&#xff0c;且可以使用与路…
阅读更多...
接口幂等性原理与方案总结

接口幂等性原理与方案总结

文章目录 接口幂等概念典型场景核心解决方案一锁二判三更新 方案选型对比 接口幂等概念 定义&#xff1a;无论调用接口多少次&#xff0c;对系统的影响与单次调用一样 范畴&#xff1a;在后端开发中&#xff0c;通常更关注写接口的幂等&#xff0c;因为写接口才会对系统数据造…
阅读更多...
【已解决】windows gitbash 出现CondaError: Run ‘conda init‘ before ‘conda activate‘

【已解决】windows gitbash 出现CondaError: Run ‘conda init‘ before ‘conda activate‘

在 Git Bash 中执行&#xff1a; source /c/Users/你的用户名/miniconda3/etc/profile.d/conda.sh # 注意填入你自己的路径 conda init bash关闭并重新打开 Git Bash 终端。测试激活环境&#xff1a; conda activate your_env_name注意事项 要把上述命令中的 你的用户名 替…
阅读更多...
软件包管理系统的架构与生态机制

软件包管理系统的架构与生态机制

文章目录 前言一、总结二、如何上传自己的软件包 前言 在日常软件开发中&#xff0c;我们经常使用诸如apt install, pip install, npm install之类的命令&#xff0c;但有一个问题是&#xff0c;这些下载命令是从哪里下载的这些软件包&#xff0c;以及我们是否能上传自己的代码…
阅读更多...
最新文章

Copyright @ 2022~2023