当前位置： 首页 > news >正文

perl初试

news 来源：原创 2025/6/9 22:08:40

我手头有一个脚本，用于从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脚本就能处理；

除了我自己的主机，

在超算上同样执行成功：

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

相关文章：

  • 网络服务之SSH协议
  • 【计算机视觉】手势识别
  • DeepSeek R1大语言模型实战工作坊02：deepseek发展演进
  • linux nginx 安装后，发现SSL模块未安装，如何处理？
  • AGI 之 【Dify】 之 使用 Docker 在 Windows 端本地部署 Dify 大语言模型（LLM）应用开发平台
  • 基于物联网技术的电动车防盗系统设计（论文+源码）
  • 【星云 Orbit • STM32F4】07. 用判断数据尾来接收据的串口通用程序框架
  • linux服务器根据内核架构下载各种软件依赖插件（例子：Anolis服务器ARM64架构内核Nginx依赖插件下载）
  • golang反射
  • cenos7网络安全检查
  • 机器学习——回归树
  • 前端基础之动画效果
  • 信贷风控系统架构设计
  • opencompass框架测试Deepseek使用教程
  • 【ORACLE】ORACLE19C在19.13版本前的一个严重BUG-24761824
  • js操作字符串的常用方法
  • 【万字长文】基于大模型的数据合成(增强）及标注
  • Pytorch的一小步，昇腾芯片的一大步
  • 【Elasticsearch】reindex
  • Pythonweb开发框架—Flask工程创建和@app.route使用详解
  • 网站建设报价分析/网络推广一个月的收入
  • qq教程网站源码/网页制作用什么软件做
  • wordpress和站点/企业推广软文
  • 网站建设 企业/百度运营平台
  • 阿里云上怎么做网页网站/互联网销售包括哪些
  • 电商货源网站/今日国内新闻最新消息大事