我手头有一个脚本,用于从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.html
use 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.html
use 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:这一行,将标志位设为1
else{
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输出文件句柄BLASTOUT
print "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:这一行,将标志位设为1
else{
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输出文件句柄BLASTOUT
print "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脚本就能处理;
除了我自己的主机,
在超算上同样执行成功:
——》其实最主要修改部分就是正则表达式那部分