perl初试

我手头有一个脚本，用于从blastp序列比对的结果文件中，进行文本处理，

获取序列比对最优的hit记录

#!/usr/bin/perl -w
use strict;my ($blast_out) = @ARGV;
my $usage = "This script is to get the best hit from blast output file with 1 input sequence.
usage: $0 <blast_output_file>
";
die $usage if @ARGV<1;open(BLASTOUT,$blast_out)||die("open $blast_out error!\n");my $query = "";
my $score = "";
my $p_value = "";
my $hit = "";
my $flag = 0;
while(<BLASTOUT>){chomp;if($flag == 0){if(/^Query=\s*(\w+)/){$query = $1;}elsif(/^Sequences producing High-scoring Segment Pairs:/){$flag = 1;}else{next;}}else{if(/^([\w\.\-]+)\s+.+\s+([0-9]+)\s+([0-9e\-\.]+)\s+[0-9]+$/){$hit = $1;$score = $2;$p_value = $3;last;}else{next;}}
}close BLASTOUT;print "Best hit to $query is: $hit, with score $score, P-value $p_value.\n";exit;

然后我们逐行解读：
1，表头的一些注释部分

#!/usr/bin/perl -w # 表明这是一个perl脚本，-w表示打开警告
# 参考https://www.runoob.com/perl/perl-environment.htmluse strict;
# use strict表示使用严格模式，这样可以避免一些不规范的写法

2，然后是输入以及输出文件的一些准备以及注释：

我们可以使用一个参考脚本来测试一下

#!/usr/bin/perl -w
use strict;# 测试不同的参数获取方式
print "命令行参数：@ARGV\n\n";# 方式1：列表赋值 - 获取第一个元素
my ($var1) = @ARGV;
print "my (\$var1) = \@ARGV 结果：\n";
print "\$var1 = '$var1'\n\n";# 方式2：标量赋值 - 获取数组长度
my $var2 = @ARGV;
print "my \$var2 = \@ARGV 结果：\n";
print "\$var2 = '$var2'\n\n";# 方式3：列表赋值 - 获取多个元素
my ($arg1, $arg2, $arg3) = @ARGV;
print "my (\$arg1, \$arg2, \$arg3) = \@ARGV 结果：\n";
print "\$arg1 = '", (defined $arg1 ? $arg1 : "未定义"), "'\n";
print "\$arg2 = '", (defined $arg2 ? $arg2 : "未定义"), "'\n";
print "\$arg3 = '", (defined $arg3 ? $arg3 : "未定义"), "'\n";

如果不提供参数的话，效果如下：

如果只提供一个参数：

如果提供两个参数的话：

如果提供3个参数的话：

——》综上，我们可以看到效果是：

my ($var1) = @ARGV

如果是加个括号的话，获取的是数组中第1个传入的参数

my $var2 = @ARGV

如果没有加括号的话，捕获的是数组的长度

my ($arg1, $arg2, $arg3) = @ARGV

如果是这样写的话，那么就可以访问数组的下标元素

——》所以综上

my ($blast_out) = @ARGV;
# 类似于shell脚本中 outputfile=$1这种写法
# 传给脚本的命令行参数列表,参考https://www.runoob.com/perl/perl-special-variables.html
# 变量定义使用my关键字，生命期直到其所在的代码块结束或者文件的末尾
# 从 @ARGV 数组中取出第一个元素，并将其赋给 $blast_out，与my $blast_out = @ARGV不同，后者会被赋予数组长度也就是元素数量
# 与my ($arg1, $arg2, $arg3) = @ARGV;写法区分

3，还有1个是函数的帮助文档

#!/usr/bin/perl -w
use strict;# 定义使用说明
my $usage = "这个脚本演示使用说明变量的用法。
用法: $0 <参数1> [参数2] ...参数1: 必需的第一个参数参数2: 可选的第二个参数
";# 检查参数
if (@ARGV < 1) {die $usage;  # 如果参数不足，终止并显示使用说明
}# 打印获取的参数
print "成功运行！\n";
print "您提供的第一个参数是: $ARGV[0]\n";if (defined $ARGV[1]) {print "您提供的第二个参数是: $ARGV[1]\n";
}

主要问题: die 语句会立即终止脚本执行，所以 die 之后的 print(“参数不足”) 语句永远不会被执行；

逻辑顺序: 如果您想打印错误信息，应该先打印，然后再 die

——》然后die usage的用法我们可以再仔细学习观察一下：

#!/usr/bin/perl -w
use strict;# 定义使用说明
my $usage = "当你看到这个信息的时候,表明你提供的参数有问题,你应该看看这个程序是怎么使用的！\n这个脚本该如何使用\n用法如下:\n $0 <参数1>" ;if (@ARGV < 1)
{print("参数不足!\n");die $usage; 
}
else
{print "提供参数的个数为：", scalar(@ARGV), "\n";
}

4，然后就是文件操作错误示范：
从语法规范上讲：

此处提供方法：如何在perl中查看函数帮助文档

# 查看特定函数文档
perldoc -f function_name
# 例如：
perldoc -f open
perldoc -f die# 查看 Perl 操作符
perldoc perlop# 查看特殊变量
perldoc perlvar  # 查看内置函数列表
perldoc perlfunc

此处涉及到perl中句柄的相关知识：

5，然后就是函数中的编程部分：

从循环语句开始：

然后就是正则表达式，/ /中间的部分

其实就是捕获Query=开头的行中的除开空格部分的单词字符

其实整个正则表达式的匹配和正则表达式中一个子串（也就是所谓的捕获组）概念，在我之前的博客awk的内置函数中的match出现过；

然后flag=1的部分的循环

这部分的正则表达式详解：

6，然后就是文件最后结尾的部分：

——》

综上，总体就是：
1个用于blastP输出结果的文本抓取脚本，其实完全可以使用awk或者说是shell来执行

#!/usr/bin/perl -w # 表明这是一个perl脚本，-w表示打开警告
# 参考https://www.runoob.com/perl/perl-environment.htmluse strict;
# use strict表示使用严格模式，这样可以避免一些不规范的写法my ($blast_out) = @ARGV;
# 类似于shell脚本中 outputfile=$1这种写法
# 传给脚本的命令行参数列表,参考https://www.runoob.com/perl/perl-special-variables.html
# 变量定义使用my关键字，生命期直到其所在的代码块结束或者文件的末尾
# 从 @ARGV 数组中取出第一个元素，并将其赋给 $blast_out，与my $blast_out = @ARGV不同，后者会被赋予数组长度也就是元素数量
# 与my ($arg1, $arg2, $arg3) = @ARGV;写法区分my $usage = "This script is to get the best hit from blast output file with 1 input sequence.
usage: $0 <blast_output_file>
";
# 定义一个字符串变量，用于提示用户如何使用该脚本,$0表示脚本名称die $usage if @ARGV<1;
# 如果参数个数小于1，输出提示信息，注意die语句会立即终止脚本执行open(BLASTOUT,$blast_out)||die("open $blast_out error!\n");
# 打开blast输出文件，保存到文件句柄BLASTOUT中，如果打开失败，输出错误信息# 然后是初始化变量以存储
my $query = "";  # 查询序列的标识符
my $score = ""; # 最佳匹配的得分
my $p_value = ""; # 最佳匹配的P值
my $hit = "";   # 最佳匹配的标识符
my $flag = 0;   # 标志位，用于判断是否到了匹配结果的部分
while(<BLASTOUT>){# while循环读取文件句柄BLASTOUT中的每一行chomp;# chomp函数用于去掉字符串末尾的换行符if($flag == 0){if(/^Query=\s*(\w+)/){$query = $1;}# 匹配Query=开头的行，提取查询序列的标识符，保存到$query中# 其实就是捕获Query=开头的行中的除开空格部分的单词字符，第1个捕获组elsif(/^Sequences producing High-scoring Segment Pairs:/){$flag = 1;}# 匹配到Sequences producing High-scoring Segment Pairs:这一行，将标志位设为1else{next;}# FLAG=0的状态，就是寻找查询信息和匹配结果部分的标题# 如果不是上面两种情况，继续下一次循环，而一旦找到了匹配结果部分的标题，就将标志位设为1，跳出当前循环}else{if(/^([\w\.\-]+)\s+.+\s+([0-9]+)\s+([0-9e\-\.]+)\s+[0-9]+$/){$hit = $1;$score = $2;$p_value = $3;last;}else{next;}}# FLAG=1的状态，就是寻找匹配结果部分的每一行，其实就是匹配到了Sequences producing High-scoring Segment Pairs:这一行之后的部分# 匹配到了这一行，那就可以接着循环每一行，直到提取标识符，得分和P值，保存到$hit,$score,$p_value中# 如果没有匹配到，就继续下一次循环，因为flag没有改变，所以还是在匹配结果部分，也就是当前循环的目的
}close BLASTOUT;
# 关闭blast输出文件句柄BLASTOUTprint "Best hit to $query is: $hit, with score $score, P-value $p_value.\n";
# 最后，将最佳匹配的结果输出到屏幕exit;
# 脚本执行结束

二，使用ssearch

 /lustre/share/class/BIO8402/tools/fasta-36.2.6/bin/ssearch36 -q u1_human.fa /lustre/share/class/BIO8402/C.elegans/Proteome/ws_215.protein.fa > new.out

现在我们改一下工具，我们使用fasta也就是S-W经典算法中的工具实现ssearch，来进行序列比对，

然后同样要求使用perl脚本进行提取记录实现。

#!/usr/bin/perl -w use strict;my ($ssearch_out) = @ARGV;my $usage = "This script is to get the best hit from ssearch output file with 1 input sequence.
usage: $0 <ssearch_output_file>
";die $usage if @ARGV<1;open(SSEARCHOUT,$ssearch_out)||die("open $ssearch_out error!\n");my $query = "";  
my $s_w = "";
my $bits = "";	
my $e_value = ""; 
my $hit = "";	
my $flag = 0;	
while(<SSEARCHOUT>){chomp;if($flag == 0){if(/^Query:\s*(\w+)/){$query = $1;}# 匹配Query:开头的行，提取查询序列的标识符，保存到$query中# 其实就是捕获Query=开头的行中的除开空格部分的单词字符，第1个捕获组elsif(/^The best scores are:/){$flag = 1;}# 匹配到The best scores are:这一行，将标志位设为1else{next;}# FLAG=0的状态，就是寻找查询信息和匹配结果部分的标题# 如果不是上面两种情况，继续下一次循环，而一旦找到了匹配结果部分的标题，就将标志位设为1，跳出当前循环}else{if(/^([\w\.\-]+)\s+.+\s+([0-9]+)\s+([0-9\.]+)\s+([0-9e\-\.]+)+$/){$hit = $1;$s_w = $2;$bits = $3;$e_value = $4;last;}else{next;}}# FLAG=1的状态，就是寻找匹配结果部分的每一行，其实就是匹配到了The best scores are:这一行之后的部分# 匹配到了这一行，那就可以接着循环每一行，直到提取标识符，得分和e值，保存到对应变量中# 如果没有匹配到，就继续下一次循环，因为flag没有改变，所以还是在匹配结果部分，也就是当前循环的目的
}close SSEARCHOUT;
# 关闭blast输出文件句柄BLASTOUTprint "Best hit to $query is: $hit, with s-w $s_w, bites $bits, e-value $e_value.\n";
# 最后，将最佳匹配的结果输出到屏幕exit;

实际上就是这么一小条序列

使用 fasta 软件包中的 ssearch36 工具，将 u1_human.fa 文件中的蛋白质序列与 ws_215.protein.fa 文件中的蛋白质序列进行比对，并将比对结果保存到 ssearch.out 文件中；

我们得到的结果也只是针对这个蛋白的分析

按理来说得到的也应该是这个序列：也就是K08D10.3

修改部分其实很简单，就是文本处理，其使用三剑客awk、sed、grep写shell脚本就能处理；

除了我自己的主机，

在超算上同样执行成功：

——》其实最主要修改部分就是正则表达式那部分