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
繁體地址，請注明出處：http://hk.pswp.cn/news/907299.shtml
英文地址，請注明出處：http://en.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