龙空技术网

实用案例精讲!如何用perl写一个截序列的脚本?

快乐的山丘hc 165

前言:

而今姐妹们对“linux运行perl脚本”大致比较看重,各位老铁们都想要分析一些“linux运行perl脚本”的相关知识。那么小编同时在网上网罗了一些有关“linux运行perl脚本””的相关内容,希望朋友们能喜欢,咱们一起来了解一下吧!

上一节正则表达式的实操课中,我们介绍了Perl程序主要包含3个部分:Head(头)、Body(正文)和Subroutine(子程序)。今天这一期内容,我们再来结合大家日常分析中经常遇到的问题来学习下字符串的操作。

我们时常发现,二代测序得到的原始序列的末端碱基质量一般比较差,这时往往需要将这些质量不好的碱基截掉。那么如何截取指定长度的序列呢?今天的实操课就告诉你答案。

假设,做基因组研究的小王,有一批reads的序列长度为60bp,现在想要截掉末端的10bp,用前50bp用于后续数据分析,小王该怎么做?

输入文件(in.fq)格式如下:

下面,我们就来详细的看下,这个案例下的Perl脚本怎么写:

#### Head ####

#!/usr/bin/perl- w

use strict;

use PerlIO::gzip

#调用压缩文件读取相关的模块

@ARGV==3 ||die"usage: perl $0 <in.fq> <out.fq> <cutlen>\n";

#设置帮助信息显示条件,当输入文件的数目不等于3个时,会终止程序并在终端显示帮助信息

PS: @ARGV是Perl默认用来接收参数的数组,Perl程序可以有多个参数,$ARGV[0]是表示接收到的第一个参数,$ARGV[1]表示第二个,以此类推。当程序出错或者满足指定条件时,可以用die函数终止程序运行,并将错误信息输出到屏幕,后面一般可以接$!(系统报错信息)或指定的报错信息。

#### Body ####

fq_lencut(@ARGV);

#正文只有一个子程序

#### Subroutine####

sub fq_lencut{

my ($infq,$outfq,$lcut) = @_;

$infq =~ /\.gz$/ ? open IFQ,"<:gzip",$infq || die $! :open IFQ,$infq || die $!;

#读取原始序列文件(压缩或非压缩)

$outfq =~ /\.gz$/ ? open OFQ,">:gzip",$outfq || die $!: open OFQ,">$outfq" || die$!; #输出截取后的序列文件(压缩或非压缩)

while(my $head = <IFQ>){

#读取序列的ID

my $seq = <IFQ>;

#读取序列的内容

my $stand = <IFQ>;

#读取序列的名称,一般用"+"代替

my $qual = <IFQ>;

#读取序列的碱基质量值

chomp($seq,$qual);

my $len = length $seq;

if($len > $lcut){

substr($seq,$lcut) = "";

#用substr函数从原始序列截取指定长度

substr($qual,$lcut) = "";

#从原始序列的碱基质量值截取指定长度

}

print OFQ $head,$seq,"\n",$stand,$qual,"\n";

#输出截取后的序列信息,包括序列的ID、截取后的序列、序列的名称("+")以及截取后的碱基质量值

}

close IFQ;

#关闭句柄,不要忘了~

close OFQ;

#这个也是

}

PS: 这里采用三目操作符(条件?语句1:语句2),判断输入的fastq文件是否.gz文件 (是否为压缩文件),从而选择打开句柄的方式。三目操作符,满足?前条件执行:前的语句1,否则执行:后的语句2。

如何获取帮助信息?

PS: 因为程序的参数不等于3(实际为0),所以会终止程序,并输出帮助信息。

如何运行?

perl fqcut.pl in.fq out.fq 50

PS: 其中in.fq是截取前的原始序列(输入文件),50是截取的长度,截取后的序列指定输出到out.fq(输出文件)。

out.fq文件格式:

以上的代码主要就是实现序列截取的功能,我们结合部分代码(带注释)介绍了Perl语言的一些常用函数的基本用法。因篇幅有限,故不能面面俱到,大家可以结合之前所讲的课程内容自己进行扩展练习,从而到达学以致用的目的,也祝大家能够有所收获!

标签: #linux运行perl脚本 #linux运行perl脚本正则