正则实例最细精解!如何从海量序列信息中提取目标序列?

在上一期中介绍了化腐朽为神奇，编程语言不可错过之—正则表达式，结合前两期的基础内容，我们从理论上认识到了perl语言的基础数据结构，和正则表达式强大的文本处理能力。事实上，做科研的宝宝经常需要在海量的数据中提取或更改数据中的某些信息，这个时候，灵活高效的正则表达式的意义就十分重大了！

今天，我们就结合科研童鞋们经常遇到的情况，给大家精讲一个Perl正则表达式的实例。

一般来说，可以将Perl程序大概分为3个部分：Head（头）、Body（正文）和Subroutine（子程序）。Head部分的主要作用是声明（解释器路径、所调用模块、调试参数等）和定义（变量、参数选项等），Body部分就是实现程序功能的主要代码，Subroutine部分主要定义一些子程序，供前面两部分调用。

今天与大家解析的案例是，如何根据基因列表提取序列？

假设做基因组研究的小余，得到了一份目标蛋白序列文件（protein.fa），文件里包含着100条蛋白序列信息，现在小余想要根据关注的特定20条序列名称（genlist.txt），提取出对应的20条目标蛋白序列信息，并存放于名为select.fa这个文件中。小余该怎么做？

protein.fa 文件格式如下：

genlist.txt文件格式如下：

Seq1

Seq3

Seq5

Seq7

Seq9

下面，我们就来详细的看下，这个案例的的Perl程序怎么写：

#### Head ####

#!/usr/bin/perl -w

use strict;

use Getopt::Long;

#控制参数的模块

my ($infile,$outfile,$list,$help);

#定义变量

＃设置参数

GetOptions(

"i|input=s" =>\$infile,

"l|list=s" =>\$list,

"o|out=s" =>\$outfile,

"h|help:s" =>\$help,

);

($infile && -s $infile) || $help || die &Usage();

#设置帮助信息显示条件，当输入文件为空或者设置-h参数时，会在终端显示帮助信息

##### Body####

open IN, $list || die “Can’t open file $list!”;

#打开句柄，读入基因列表

my %listids;

while(<IN>){

#逐行读取基因列表，直到最后一行

chomp;

#chomp函数的作用是去掉最后的记录分隔符（默认是”\n”）

$listids{$_}=1;

#将基因列表保存到哈希，记录在案

}

close IN;

#关闭句柄

$/=">";

#特殊变量">/表示记录分隔符，这里将它设置为”>”，默认是”\n”

openIN,$infile || die $!;

#打开句柄，读入序列文件

<IN>;

#读取第一个记录分隔符”>”，由于没有序列，所以不做任何操作

openOUT,”>$outfile” || die $!;

#打开句柄，输出文件

while(<IN>){

#逐一读取基因列表，直到最后一条

chomp;

#去掉最后的记录分隔符（”>”）

s/^\s+//g;

#将序列间的所有非空字符去掉，使得序列不换行

my $id;

$id=$1 if(/^(\S+)/);

#提取序列的ID

print OUT ">$_" if($listids{$id});

#如果序列ID存在于哈希(基因列表)中，就将序列输出

}

close IN;

#读取文件完毕后，记得养成随手关闭句柄的好习惯

close OUT;

#输出文件的句柄也是一样

$/=”\n”;

#将记录分隔符改成默认的”\n”，这也是编程的好习惯之一

#### Subroutine####

＃定义获取帮助信息的子程序，一般包括主程序的作用介绍、使用方法、参数说明以及版本信息等。

sub Usage{

my $info=”

Usage: perl $0 -i <infile> -l <genlist> -o <outfile>

#使用方法说明，这里特殊变量$0储存Perl程序自身的名称

Options:

#可选参数列表及说明

-i|input <str> set input fa file

-l|list <str> inputfa id list seperated by '\\n'

-o|out <str> setoutput file

-h|help set output file

“;

print $info;

exit 0;

#正常退出

}

P.S.: 参数列表中的"|"表示或，也就是选择短参数（-i）或者长参数（-input）都表示相同意思，其后接输入文件。<str> 表示所接参数类型为字符串（string），没有则表示参数后面可以不接任何内容，比如-h参数。

通过-h参数获取程序的帮助信息：

代码写完了，来和晨光一起运行下吧：

perl test.pl -i protein.fa -l genlist.txt -o select.fa

大家快来一起试试吧...

P.S.: 我们为大家准备了一套测试数据，想获取demo data的童鞋，欢迎联系小秘书weixin：genegogo007）

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

/End.