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.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脚本就能处理;

除了我自己的主机,

在超算上同样执行成功:

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

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

Nginx1.19.2不适配OPENSSL3.0问题

Nginx 1.19.2 是较老的版本&#xff0c;而 Nginx 1.21 版本已经适配 OpenSSL 3.0&#xff0c;所以建议 升级 Nginx 到 1.25.0 或更高版本&#xff1a; wget http://nginx.org/download/nginx-1.25.0.tar.gz tar -xzf nginx-1.25.0.tar.gz cd nginx-1.25.0 ./configure --prefix…

MySQL复合查询——通过案例讲解每个指令

0.准备工作 在开始之前可以先准备好相同的数据库 方法一&#xff1a;直接在MySQL创建相应的数据库和表 第一步&#xff1a;创建数据库并进入数据库 create database soctt_data; use soctt_data; 第二步&#xff1a;创建部门信息表 DROP TABLE IF EXISTS dept; CREATE TABL…

Kubernetes全解析:从容器编排到云原生霸主

前言 在数字化转型浪潮中&#xff0c;云原生技术已成为企业构建敏捷、弹性基础设施的核心驱动力。作为容器编排领域的“操作系统”&#xff0c;Kubernetes&#xff08;K8s&#xff09;凭借其自动化部署、弹性伸缩和跨环境一致性等能力&#xff0c;正重新定义现代应用的运维范式…

我的两个医学数据分析技术思路

我的两个医学数据分析技术思路 从临床上获得的或者公共数据库数据这种属于观察性研究&#xff0c;是对临床诊疗过程中自然产生的数据进行分析而获得疾病发生发展的规律等研究成果。再细分&#xff0c;可以分为独立危险因素鉴定和预测模型构建两种。 独立危险因素鉴定是一直以…

图像滑块对比功能的开发记录

背景介绍 最近&#xff0c;公司需要开发一款在线图像压缩工具&#xff0c;其中的一个关键功能是让用户直观地比较压缩前后的图像效果。因此&#xff0c;我们设计了一个对比组件&#xff0c;它允许用户通过拖动滑块&#xff0c;动态调整两张图像的显示区域&#xff0c;从而清晰…

迷你世界脚本UI五子棋小游戏

wzq_jm "7477124677881080183-22855"--界面id wzq_jmjxh "7477124677881080183-22855_"--界面加下划线 wzq_tc "7477124677881080183-22855_262"--退出按钮id wzq_hdlt1 "7477124677881080183-22855_267"--互动聊天按钮 快点吧&a…

大模型理论基础介绍

大模型理论基础 {docsify-ignore-all} 项目简介 本项目旨在作为一个大规模预训练语言模型的教程&#xff0c;从数据准备、模型构建、训练策略到模型评估与改进&#xff0c;以及模型在安全、隐私、环境和法律道德方面的方面来提供开源知识。 项目将以斯坦福大学大规模语言模型课…

【Spring Boot 应用开发】-04-02 自动配置-数据源-手撸一个最简持久层工具类

设计概述 有时候我们不需要太重的持久层&#xff0c;就像要一个最简的、轻量的持久层&#xff0c;便于维护和扩展&#xff0c;代码掌握在自己手里&#xff0c;那么我们可以基于springboot的自动配置&#xff0c;快速的构建一个自己的持久层轻量框架&#xff0c;不说废话&#…

MicroServer Gen8再玩之三 OCP万兆光口+12G阵列卡

前一段时间&#xff0c;做了一片双OCP的合成转接卡&#xff0c;在GEN8上用了起来&#xff0c;有些小伙伴觉得还不错&#xff0c;有些则对LSI2308这块阵列卡性能表示不甚满意。 于是乎&#xff0c;就有了后续折腾的理由。 前一段时间&#xff0c;我还不了解阵列卡有啥区别&…

PostgreSQL10 物理流复制实战:构建高可用数据库架构!

背景 PostgreSQL 10 在高可用架构中提供了物理复制&#xff0c;也称为流复制&#xff08;Streaming Replication&#xff09;&#xff0c;用于实现实例级别的数据同步。PostgreSQL 复制机制主要包括物理复制和逻辑复制&#xff1a;物理复制依赖 WAL 日志进行物理块级别的同步&…

Linux网络安全技术与实现

&#x1f345; 点击文末小卡片 &#xff0c;免费获取网络安全全套资料&#xff0c;资料在手&#xff0c;涨薪更快 Linux 网络安全和优化 Jephe Wu 翻译整理 简介 网络安全是一个非常重要的课题,基本上你运行的服务后台越多,你就可能打开更多的安全漏洞.如果配置的恰当的话,Li…

[黑马点评]关于原子性,锁的笔记

不得不说&#xff0c;黑马点评是一个非常不错的课程&#xff0c;对于线程安全方面的讲解十分详细且明朗&#xff0c;故写下这篇笔记方便复习及帮助后人&#xff08;&#xff09; 目标 我们的目标是对于大量对于优惠劵的访问时&#xff0c;要防止超卖问题以及一人多单问题。 单J…

mapbox高阶,结合threejs(threebox)添加三维球体

👨‍⚕️ 主页: gis分享者 👨‍⚕️ 感谢各位大佬 点赞👍 收藏⭐ 留言📝 加关注✅! 👨‍⚕️ 收录于专栏:mapbox 从入门到精通 文章目录 一、🍀前言1.1 ☘️mapboxgl.Map 地图对象1.2 ☘️mapboxgl.Map style属性1.3 ☘️threebox Sphere静态对象二、🍀使用t…

MAC 本地搭建部署 dify(含 github访问超时+Docker镜像源拉取超时解决方案)

目录 一、什么是 dify&#xff1f; 二、安装 docker 1. 什么是 docker&#xff1f; 2. docker下载地址 三、安装 dify 1. dify下载地址 2.可能遇到问题一&#xff1a; github访问超时 3.下载后完成解压 4.进入到 cmd 终端环境&#xff0c;执行下面三个命令 5.可能遇到…

Pytorch xpu环境配置 Pytorch使用Intel集成显卡

1、硬件集显要为Intel ARC并安装正确驱动 2、安装Intel oneAPI Base Toolkit &#xff08;https://www.intel.cn/content/www/cn/zh/developer/tools/oneapi/base-toolkit-download.html&#xff09;安装后大约20G左右&#xff0c;注意安装路径 3、安装Visual Studio Build To…

若依前后端分离版使用Electron打包前端Vue为Exe文件

1.前言 本文详细介绍如何使用electron将若依框架前后端分离版的前端Vue页面打包为Exe文件&#xff0c;并且包括如何实现应用更新。使用若依基础代码体现不出打包功能&#xff0c;因此我使用开发的文件管理系统&#xff0c;介绍上述过程&#xff0c;具体可以查看我的文章《若依…

docker:Dockerfile案例之自定义centos7镜像

1 案例需求 自定义centos7镜像。要求&#xff1a; 默认登录路径为 /usr可以使用vim 2 实施步骤 编写dockerfile脚本 vim centos_dockerfile 内容如下&#xff1a; #定义父镜像 FROM centos:7#定义作者信息 MAINTAINER handsome <handsomehandsome.com># 设置阿里云…

SpringBoot校园管理系统设计与实现

在现代校园管理中&#xff0c;一个高效、灵活的管理系统是不可或缺的。本文将详细介绍基于SpringBoot的校园管理系统的设计与实现&#xff0c;涵盖管理员、用户和院校管理员三大功能模块&#xff0c;以及系统的部署步骤和数据库配置。 管理员功能模块 管理员是系统的核心管理…

[项目]基于FreeRTOS的STM32四轴飞行器: 四.LED控制

基于FreeRTOS的STM32四轴飞行器: 四.LED控制 一.配置Com层二.编写驱动 一.配置Com层 先在Com_Config.h中定义灯位置的枚举类型&#xff1a; 之后定义Led的结构体&#xff1a; 定义飞行器状态&#xff1a; 在Com_Config.c中初始化四个灯&#xff1a; 在Com_Config.h外部声明…

Linux部署java项目

前言 Xshell下载地址 点击连接 常见命令 ls ls:显示当前目录下的文件 ll:可以显示隐藏文件和非隐藏文件与ls -l一样 ls -a -l这两个掌握就可以了 ls --help就可以知道这个后面可以跟什么 ls -al还可以这样 cd cd&#xff1a;进入文件夹 cd后面可以跟相对路径&#xff0…