生信：种群历史有效群体大小推断 SMC++(三)

好嘛，我的文章果然有帮助到其他人！开心到爆炸！实现了自我价值，谢谢你陌生人！（狗头保命）

gen_raw_mask.pl 脚本是用于从输入的 FASTA 文件（sm.fasta）中提取简单的重复序列（ATCGN 组成）的位置，并将结果以 BED 格式输出。

代码解读如下：

1. use strict; 和 use Bio::SeqIO; 是引入必要的 Perl 模块，strict 告诉 Perl 在变量使用前必须声明，Bio::SeqIO 是处理生物序列的模块。

2. die "perl $0 <sm.fasta> > out.bed " unless @ARGV >= 1; 是一个错误处理语句，如果没有提供输入文件名，则脚本输出错误信息，并显示正确的用法。

3. my $infile = shift; 从命令行参数中获取输入的 FASTA 文件名。

4. my $fa = Bio::SeqIO->new(-file => $infile, -format => 'fasta'); 创建一个 Bio::SeqIO 对象，用于读取 FASTA 文件。

5. while (my $obj = $fa->next_seq()) { ... } 是一个循环，每次读取 FASTA 文件中的一条序列，并执行循环体内的操作。

6. my $id = $obj->id; 获取当前序列的标识符（ID）。

7. my $seq = $obj->seq; 获取当前序列的核苷酸序列。

8. my $len = $obj->length; 获取当前序列的长度。

9. while ($seq =~ /([atcgn]+)/g) { ... } 是一个嵌套的 while 循环，用于在当前序列中查找简单的重复序列（ATCGN 组成）。

10. my $e = pos($seq); 获取重复序列的匹配终止位置，位置是从 1 开始的，所以实际终止位置需要减 1。

11. my $l = length($1); 获取匹配到的重复序列的长度。

12. my $s = $e - $l; 计算匹配到的重复序列的起始位置，位置是从 0 开始的。

13. print "$id\t$s\t$e\n"; 将匹配到的重复序列的位置和序列标识符以 BED 格式输出（tab 分隔的三列）。

整个脚本的目标是读取输入的 FASTA 文件，查找其中的简单重复序列（ATCGN 组成）并输出其位置信息（BED 格式）。

perl脚本：

#!/usr/bin/env perl

use strict;
use Bio::SeqIO;

# 检查命令行参数是否提供了输入文件名
die "perl $0 <sm.fasta> > out.bed " unless @ARGV >= 1;

# 获取输入的 FASTA 文件名
my $infile = shift;

# 创建 Bio::SeqIO 对象，用于读取 FASTA 文件
my $fa = Bio::SeqIO->new(-file => $infile, -format => 'fasta');

# 循环读取 FASTA 文件中的每条序列
while (my $obj = $fa->next_seq()) {

    # 获取当前序列的标识符（ID）
    my $id = $obj->id;

    # 获取当前序列的核苷酸序列
    my $seq = $obj->seq;

    # 获取当前序列的长度
    my $len = $obj->length;

    # 在当前序列中查找简单的重复序列（ATCGN 组成）
    while ($seq =~ /([atcgn]+)/g) {

        # 获取重复序列的匹配终止位置，位置是从 1 开始的，所以实际终止位置需要减 1
        my $e = pos($seq) - 1;

        # 获取匹配到的重复序列的长度
        my $l = length($1);

        # 计算匹配到的重复序列的起始位置，位置是从 0 开始的
        my $s = $e - $l;

        # 将匹配到的重复序列的位置和序列标识符以 BED 格式输出（tab 分隔的三列）
        print "$id\t$s\t$e\n";
    }
}


这里涉及到要安装 Bio::SeqIO 模块，首先需要确保已经安装了 Perl 和 CPAN（Comprehensive Perl Archive Network）工具。CPAN 是 Perl 社区的一个资源库，用于下载和管理 Perl 模块。

以下是在 Linux 和 macOS 系统中安装 Bio::SeqIO 模块的步骤：

1. 打开终端（命令行界面）。

2. 输入以下命令来进入 CPAN shell：

   cpan


3. 第一次运行 CPAN shell 时，会提示您进行初始配置。按照提示一步一步设置即可。

4. 在 CPAN shell 中，输入以下命令来安装 Bio::SeqIO 模块：

   install Bio::SeqIO


5. CPAN shell 将自动下载、编译和安装 Bio::SeqIO 及其依赖项。

6. 安装完成后，输入 quit 命令退出 CPAN shell。

gen_raw_mask.pl是程序自带的

github连接:
https://github.com/lh3/misc/tree/cc0f36a9a19f35765efb9387389d9f3a6756f08f/seq/seqbility


这段代码是一个 Perl 脚本，用于处理输入的 BWA SAM 文件并生成一个基因组掩码（genome mask）的输出。

• 脚本首先定义了一个子例程 main，用于处理输入并执行主要操作。然后通过调用 main 子例程来运行整个脚本。

• 在 main 子例程中，首先检查命令行参数是否正确，如果参数错误，则输出使用说明并退出脚本。

• 然后，脚本定义了一个数组 @conv，用于保存字符转换表。接下来，通过计算 conv 表中的元素值，来构建一个转换表。转换表用于将两个整数 a0 和 a1 映射为一个字符，并将该字符附加到序列 seq 中。

• 在核心循环中，脚本逐行读取输入文件（或从标准输入读取），解析 SAM 文件的内容，并将匹配情况映射到字符。然后，将得到的字符拼接到 seq 变量中。

• 子例程 print_seq 负责将生成的序列打印输出。它根据每行最大长度 60 的格式，将序列按行打印。

• 子例程 int_log2 用于计算给定整数 v 的 log2 值，并返回结果。

综上所述，这个 Perl 脚本的目标是将输入的 BWA SAM 文件转换为一个基因组掩码，并将掩码输出到标准输出或文件中。掩码中的每个字符表示匹配情况，并通过一个转换表将两个整数映射为一个字符。生成的掩码序列以每行 60 个字符的格式输出。