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

Posted 2021-04-17 生信者言

上一节正则表达式的实操课中，我们介绍了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> "; 

#设置帮助信息显示条件，当输入文件的数目不等于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," ",$stand,$qual," ";

        #输出截取后的序列信息，包括序列的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语言的一些常用函数的基本用法。因篇幅有限，故不能面面俱到，大家可以结合之前所讲的课程内容自己进行扩展练习，从而到达学以致用的目的，也祝大家能够有所收获！

/End.

推荐阅读

猛戳下方链接即可阅读

扫码关注，获取更多精彩内容

我

是

彩

蛋

喜马拉雅FM搜索并订阅：生信者言；收听内容：

《一分钟听懂NGS基础概念》，让生信分析不再遥不可及

《亲爱的姑娘，你值得被温柔以待》，11个真实的人物故事

《众病之王：癌症传》，一起聆听人类对抗癌症的斗争史

回复文字：果然科学，看一篇好玩的科普文。