perl初试

news2025/3/6 21:24:11

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

除了我自己的主机,

在超算上同样执行成功:

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/2310706.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

VS Code C++ 开发环境配置

VS Code 是当前非常流行的开发工具. 本文讲述如何配置 VS Code 作为 C开发环境. 本文将按照如下步骤来介绍如何配置 VS Code 作为 C开发环境. 安装编译器安装插件配置工作区 第一个步骤的具体操作会因为系统不同或者方案不同而有不同的选择. 环境要求 首先需要立即 VS Code…

Web Snapshot 网页截图 模块代码详解

本文将详细解析 Web Snapshot 模块的实现原理和关键代码。这个模块主要用于捕获网页完整截图&#xff0c;特别优化了对动态加载内容的处理。 1. 模块概述 snapshot.py 是一个功能完整的网页截图工具&#xff0c;它使用 Selenium 和 Chrome WebDriver 来模拟真实浏览器行为&am…

Windows 10 下 SIBR Core (i.e. 3DGS SIBR Viewers) 的编译

本文针对在 Windows 10 上从源码编译安装3DGS &#xff08;3D Gaussian Splatting&#xff09;的Viewers 即SIBR Core及外部依赖库extlibs&#xff08;预编译的版本直接在页面https://sibr.gitlabpages.inria.fr/download.html下载&#xff09; &#xff0c;参考SIBR 的官方网站…

JavaWeb-HttpServletRequest请求域接口

文章目录 HttpServletRequest请求域接口HttpServletRequest请求域接口简介关于请求域和应用域的区别 请求域接口中的相关方法获取前端请求参数(getParameter系列方法)存储请求域名参数(Attribute系列方法)获取客户端的相关地址信息获取项目的根路径 关于转发和重定向的细致剖析…

防火墙虚拟系统实验

拓扑图 需求一 安全策略要求&#xff1a; 1、只存在一个公网IP地址&#xff0c;公司内网所有部门都需要借用同一个接口访问外网 2、财务部禁止访问Internet&#xff0c;研发部门只有部分员工可以访问Internet&#xff0c;行政部门全部可以访问互联网 3、为三个部门的虚拟系统分…

点云滤波方法:特点、作用及使用场景

点云滤波是点云数据预处理的重要步骤&#xff0c;目的是去除噪声点、离群点等异常数据&#xff0c;平滑点云或提取特定频段特征&#xff0c;为后续的特征提取、配准、曲面重建、可视化等高阶应用打下良好基础。以下是点云中几种常见滤波方法的特点、作用及使用场景&#xff1a;…

Gradle 配置 Lombok 项目并发布到私有 Maven 仓库的完整指南

Gradle 配置 Lombok 项目并发布到私有 Maven 仓库的完整指南 在 Java 项目开发中&#xff0c;使用 Lombok 可以极大地减少样板代码&#xff08;如 getter/setter 方法、构造器等&#xff09;&#xff0c;提高开发效率。然而&#xff0c;当使用 Gradle 构建工具并将项目发布到私…

ArcGIS Pro 基于基站数据生成基站扇区地图

在当今数字化的时代&#xff0c;地理信息系统&#xff08;GIS&#xff09;在各个领域都发挥着至关重要的作用。 ArcGIS Pro作为一款功能强大的GIS软件&#xff0c;为用户提供了丰富的工具和功能&#xff0c;使得数据处理、地图制作和空间分析变得更加高效和便捷。 本文将为您…

【Python · Pytorch】Conda介绍 DGL-cuda安装

本文仅涉及DGL库介绍与cuda配置&#xff0c;不包含神经网络及其训练测试。 起因&#xff1a;博主电脑安装了 CUDA 12.4 版本&#xff0c;但DGL疑似没有版本支持该CUDA版本。随即想到可利用Conda创建CUDA12.1版本的虚拟环境。 1. Conda环境 1.1 Conda环境简介 Conda&#xff1…

leetcode:2965. 找出缺失和重复的数字(python3解法)

难度&#xff1a;简单 给你一个下标从 0 开始的二维整数矩阵 grid&#xff0c;大小为 n * n &#xff0c;其中的值在 [1, n2] 范围内。除了 a 出现 两次&#xff0c;b 缺失 之外&#xff0c;每个整数都 恰好出现一次 。 任务是找出重复的数字a 和缺失的数字 b 。 返回一个下标从…

Android U 分屏——SystemUI侧处理

WMShell相关的dump命令 手机分屏启动应用后运行命令&#xff1a;adb shell dumpsys activity service SystemUIService WMShell 我们可以找到其中分屏的部分&#xff0c;如下图所示&#xff1a; 分屏的组成 简图 分屏是由上分屏(SideStage)、下分屏(MainStage)以及分割线组…

flink集成tidb cdc

Flink TiDB CDC 详解 1. TiDB CDC 简介 1.1 TiDB CDC 的核心概念 TiDB CDC 是 TiDB 提供的变更数据捕获工具&#xff0c;能够实时捕获 TiDB 集群中的数据变更&#xff08;如 INSERT、UPDATE、DELETE 操作&#xff09;&#xff0c;并将这些变更以事件流的形式输出。TiDB CDC 的…

推荐1款OCR的扫描仪软件,无需安装,打开即用!

聊一聊 现在日常办公&#xff0c;很多时候还是需要扫描仪配合。 很多时候需要将文件搜索成PDF再传输。 今天给大家分享一款OCR扫描仪软件。 软件介绍 OCR的扫描仪软件 支持扫描仪共享。 支持WIA、TWAIN、SANE和ESCL驱动程序。 还可以批量多扫描仪配置扫描&#xff0c;支持…

SpringBoot为什么默认使用CGLIB?

大家好&#xff0c;我是锋哥。今天分享关于【SpringBoot为什么默认使用CGLIB?】面试题。希望对大家有帮助&#xff1b; SpringBoot为什么默认使用CGLIB? 1000道 互联网大厂Java工程师 精选面试题-Java资源分享网 Spring Boot 默认使用 CGLIB&#xff08;Code Generation Li…

神经网络|(十三)|SOM神经网络

【1】引言 前序已经对神经网络有了基础认识&#xff0c;今天先学习SOM神经网络。 前序学习文章链接包括且不限于&#xff1a; 神经网络|(十一)|神经元和神经网络-CSDN博客 神经网络|(十二)|常见激活函数-CSDN博客 【2】SOM神经网络 SOM神经网络是一种结构比较简单、但是理…

IP协议、DNS协议、DHCP协议、Telent协议的记忆总结

首先记忆一下几个协议的端口号 HTTP&#xff1a;超文本传输协议 80 HTTPS&#xff1a;安全传输协议 443 DHCP&#xff1a;动态主机配置协议 67/68 DNS&#xff1a;域名解析协议 53 FTP&#xff1a;文件传输协议 20/21 TFTP&#xff1a;简单文件传输协议 69 TELENT&#xff1a;远…

Pico 4 Enterprise(企业版)与Unity的交互-有线串流调试篇

入手了Pico 4 E做VR开发&#xff0c;谁知入了天坑...根据官方文档&#xff0c;尝试了串流助手、企业串流、PICO Developer Center&#xff0c;陷入了各种版本问题、环境问题的陷阱。而且Pico4E的OS自24年12开始就不再更新&#xff0c;头盔中预装的企业串流版本也较低&#xff0…

DeepSeek-R1:使用KTransformers实现高效部署指南

KTransformers作为一个开源框架&#xff0c;专门为优化大规模语言模型的推理过程而设计。它支持GPU/CPU异构计算&#xff0c;并针对MoE架构的稀疏性进行了特别优化&#xff0c;可以有效降低硬件要求&#xff0c;允许用户在有限的资源下运行像DeepSeek-R1这样庞大的模型。 硬件…

任务9:交换机基础及配置

CSDN 原创主页&#xff1a;不羁https://blog.csdn.net/2303_76492156?typeblog 一、交换机基础 交换机的概念&#xff1a;交换机是一种网络设备&#xff0c;用于连接多台计算机或网络设备&#xff0c;实现数据包在局域网内的快速交换。交换机基于MAC地址来转发数据包&#x…

Notepad++ 8.6.7 安装与配置全攻略(Windows平台)

一、软件定位与核心优势 Notepad 是开源免费的代码/文本编辑器&#xff0c;支持超过80种编程语言的高亮显示&#xff0c;相比系统自带记事本具有以下优势&#xff1a; 轻量高效&#xff1a;启动速度比同类软件快30%插件扩展&#xff1a;支持NppExec、JSON Viewer等200插件跨文…