生信：種群歷史有效群體大小推斷 SMC++(三)
好嘛，我的文章果然有幫助到其他人！開心到爆炸！實(shí)現(xiàn)了自我價(jià)值，謝謝你陌生人?。ü奉^保命）
gen_raw_mask.pl 腳本是用于從輸入的 FASTA 文件（sm.fasta）中提取簡(jiǎn)單的重復(fù)序列（ATCGN 組成）的位置，并將結(jié)果以 BED 格式輸出。
代碼解讀如下：
use strict;?和?use Bio::SeqIO;?是引入必要的 Perl 模塊，strict?告訴 Perl 在變量使用前必須聲明，Bio::SeqIO?是處理生物序列的模塊。die "perl $0 <sm.fasta> > out.bed " unless @ARGV >= 1;?是一個(gè)錯(cuò)誤處理語(yǔ)句，如果沒(méi)有提供輸入文件名，則腳本輸出錯(cuò)誤信息，并顯示正確的用法。my $infile = shift;?從命令行參數(shù)中獲取輸入的 FASTA 文件名。my $fa = Bio::SeqIO->new(-file => $infile, -format => 'fasta');?創(chuàng)建一個(gè) Bio::SeqIO 對(duì)象，用于讀取 FASTA 文件。while (my $obj = $fa->next_seq()) { ... }?是一個(gè)循環(huán)，每次讀取 FASTA 文件中的一條序列，并執(zhí)行循環(huán)體內(nèi)的操作。my $id = $obj->id;?獲取當(dāng)前序列的標(biāo)識(shí)符（ID）。my $seq = $obj->seq;?獲取當(dāng)前序列的核苷酸序列。my $len = $obj->length;?獲取當(dāng)前序列的長(zhǎng)度。while ($seq =~ /([atcgn]+)/g) { ... }?是一個(gè)嵌套的 while 循環(huán)，用于在當(dāng)前序列中查找簡(jiǎn)單的重復(fù)序列（ATCGN 組成）。my $e = pos($seq);?獲取重復(fù)序列的匹配終止位置，位置是從 1 開始的，所以實(shí)際終止位置需要減 1。my $l = length($1);?獲取匹配到的重復(fù)序列的長(zhǎng)度。my $s = $e - $l;?計(jì)算匹配到的重復(fù)序列的起始位置，位置是從 0 開始的。print "$id\t$s\t$e\n";?將匹配到的重復(fù)序列的位置和序列標(biāo)識(shí)符以 BED 格式輸出（tab 分隔的三列）。
整個(gè)腳本的目標(biāo)是讀取輸入的 FASTA 文件，查找其中的簡(jiǎn)單重復(fù)序列（ATCGN 組成）并輸出其位置信息（BED 格式）。
perl腳本：
#!/usr/bin/env?perl
use?strict;
use?Bio::SeqIO;
#?檢查命令行參數(shù)是否提供了輸入文件名
die?"perl?$0?<sm.fasta>?>?out.bed?"?unless?@ARGV?>=?1;
#?獲取輸入的?FASTA?文件名
my?$infile?=?shift;
#?創(chuàng)建?Bio::SeqIO?對(duì)象，用于讀取?FASTA?文件
my?$fa?=?Bio::SeqIO->new(-file?=>?$infile,?-format?=>?'fasta');
#?循環(huán)讀取?FASTA?文件中的每條序列
while?(my?$obj?=?$fa->next_seq())?{
????#?獲取當(dāng)前序列的標(biāo)識(shí)符（ID）
????my?$id?=?$obj->id;
????#?獲取當(dāng)前序列的核苷酸序列
????my?$seq?=?$obj->seq;
????#?獲取當(dāng)前序列的長(zhǎng)度
????my?$len?=?$obj->length;
????#?在當(dāng)前序列中查找簡(jiǎn)單的重復(fù)序列（ATCGN?組成）
????while?($seq?=~?/([atcgn]+)/g)?{
????????#?獲取重復(fù)序列的匹配終止位置，位置是從?1?開始的，所以實(shí)際終止位置需要減?1
????????my?$e?=?pos($seq)?-?1;
????????#?獲取匹配到的重復(fù)序列的長(zhǎng)度
????????my?$l?=?length($1);
????????#?計(jì)算匹配到的重復(fù)序列的起始位置，位置是從?0?開始的
????????my?$s?=?$e?-?$l;
????????#?將匹配到的重復(fù)序列的位置和序列標(biāo)識(shí)符以?BED?格式輸出（tab?分隔的三列）
????????print?"$id\t$s\t$e\n";
????}
}
這里涉及到要安裝 Bio::SeqIO 模塊，首先需要確保已經(jīng)安裝了 Perl 和 CPAN（Comprehensive Perl Archive Network）工具。CPAN 是 Perl 社區(qū)的一個(gè)資源庫(kù)，用于下載和管理 Perl 模塊。
以下是在 Linux 和 macOS 系統(tǒng)中安裝 Bio::SeqIO 模塊的步驟：
打開終端（命令行界面）。
輸入以下命令來(lái)進(jìn)入 CPAN shell：
cpan第一次運(yùn)行 CPAN shell 時(shí)，會(huì)提示您進(jìn)行初始配置。按照提示一步一步設(shè)置即可。
在 CPAN shell 中，輸入以下命令來(lái)安裝 Bio::SeqIO 模塊：
install?Bio::SeqIOCPAN shell 將自動(dòng)下載、編譯和安裝 Bio::SeqIO 及其依賴項(xiàng)。
安裝完成后，輸入?
quit?命令退出 CPAN shell。
gen_raw_mask.pl是程序自帶的
github連接:
https://github.com/lh3/misc/tree/cc0f36a9a19f35765efb9387389d9f3a6756f08f/seq/seqbility
這段代碼是一個(gè) Perl 腳本，用于處理輸入的 BWA SAM 文件并生成一個(gè)基因組掩碼（genome mask）的輸出。
腳本首先定義了一個(gè)子例程?
main，用于處理輸入并執(zhí)行主要操作。然后通過(guò)調(diào)用?main?子例程來(lái)運(yùn)行整個(gè)腳本。在?
main?子例程中，首先檢查命令行參數(shù)是否正確，如果參數(shù)錯(cuò)誤，則輸出使用說(shuō)明并退出腳本。然后，腳本定義了一個(gè)數(shù)組?
@conv，用于保存字符轉(zhuǎn)換表。接下來(lái)，通過(guò)計(jì)算?conv?表中的元素值，來(lái)構(gòu)建一個(gè)轉(zhuǎn)換表。轉(zhuǎn)換表用于將兩個(gè)整數(shù) a0 和 a1 映射為一個(gè)字符，并將該字符附加到序列?seq?中。在核心循環(huán)中，腳本逐行讀取輸入文件（或從標(biāo)準(zhǔn)輸入讀?。?，解析 SAM 文件的內(nèi)容，并將匹配情況映射到字符。然后，將得到的字符拼接到?
seq?變量中。子例程?
print_seq?負(fù)責(zé)將生成的序列打印輸出。它根據(jù)每行最大長(zhǎng)度 60 的格式，將序列按行打印。子例程?
int_log2?用于計(jì)算給定整數(shù) v 的 log2 值，并返回結(jié)果。
綜上所述，這個(gè) Perl 腳本的目標(biāo)是將輸入的 BWA SAM 文件轉(zhuǎn)換為一個(gè)基因組掩碼，并將掩碼輸出到標(biāo)準(zhǔn)輸出或文件中。掩碼中的每個(gè)字符表示匹配情況，并通過(guò)一個(gè)轉(zhuǎn)換表將兩個(gè)整數(shù)映射為一個(gè)字符。生成的掩碼序列以每行 60 個(gè)字符的格式輸出。
