BioPerl（三）：巧用BioPerl格式化fasta文件

Overview

在处理fasta格式序列的过程中，我们经常会发现得到的fasta格式并不是很标准，比如有一个fasta文件中有多条这样形式的序列：

>gi|28898692|ref|NP_798297.1| hypothetical protein VP1918 [Vibrio parahaemolyticus RIMD 2210633]|1
MKKTTLMSAVVATLSLVGCQSTTGSSDAQPEQTSHISQAVYEVEFHAAQSFLSQASQLEQSFADFCLAPK
NDVEPVQQQWHSTMLAWMALQGQERGPATALEQSWNVQFWPDKKNTTGRKMSALTKADKVWTVEEISTQS
VTVQGLGALEWLLYDDASTLNTNSNVCASGVAIAENLHDKAQIIANSWAENPWKSLQKTEWESEYISLLS
NQLEYSMKKLSRPLAKIGHPRPYFSESWRSETSLSNLKANLESLHQLYFANGKGLDALLRAQGKTQLADR

它的序列内容并不是一行，而是四行，尽管这是被允许的，但是有的时候出于一些目的，我们想要序列的内容是一行，而不是多行。

我自己写了一个脚本用来将每一行的换行符去掉，并保留最后一行的换行符，代码看起来很丑陋。既然BioPerl可以帮我们自动处理序列的很多信息，那是否也能帮我们格式化序列的内容，我抱着试试看的态度试了一下，真的可以...

1. `BioPerl`格式化fasta序列

其实复用了 BioPerl（二）：使用BioPerl读取fasta文件中的脚本：

#!/usr/bin/perl -w
use strict;
use warnings;
use Bio::SeqIO;
use Bio::Seq;

my $catchseq_seqio_obj = Bio::SeqIO->new(-file=>"unFormatedFastaFile.fasta", -format=>'fasta');

# 在这儿处理每个序列的信息
while(my $seq_obj = $catchseq_seqio_obj->next_seq)
{
    # ...
    my $seq = $seq_obj->seq; 
    # ... 
    # 重新拼成fasta格式
    # 注意，$display_name中是不包含fasta的标识符>的，所以拼接时要手动加上>
    my $seqContent = ">".$display_name.$desc."\n".$seq."\n";                   
}
# 省略了将 $seqContent重新写到文件..

只需要使用BioPerl把不标准的Fasta格式文件读取到，再写到新的文件中就可以了，重写之后，格式为：

>gi|28898692|ref|NP_798297.1| hypothetical protein VP1918 [Vibrio parahaemolyticus RIMD 2210633]|1
MKKTTLMSAVVATLSLVGCQSTTGSSDAQPEQTSHISQAVYEVEFHAAQSFLSQASQLEQSFADFCLAPKNDVEPVQQQWHSTMLAWMALQGQERGPATALEQSWNVQFWPDKKNTTGRKMSALTKADKVWTVEEISTQSVTVQGLGALEWLLYDDASTLNTNSNVCASGVAIAENLHDKAQIIANSWAENPWKSLQKTEWESEYISLLSNQLEYSMKKLSRPLAKIGHPRPYFSESWRSETSLSNLKANLESLHQLYFANGKGLDALLRAQGKTQLADR

可以看到序列中的换行符已经去除了。既然可以格式化序列内容，也可以顺便做点别的，比如读取一个fasta序列时，去除那些序列长度过小的序列，不然使用有些特征提取算法时会出现异常。

2. 使用`BioPerl`去除fasta文件中长度过短的序列

这里给出一个完整点的例子，写一个叫removeShortFastaEntries.pl，内容如下：

#!/usr/bin/perl -w

# usage: perl removeShortFastaEntries.pl original.fasta formatedEntries.fasta [50]

use strict;
use warnings;
use Bio::SeqIO;
use  Bio::Seq;

# 在命令行设置要处理的fasta文件名字
my $input_file = $ARGV[0] or die "Need to input fasta file on the command line\n";
# 在命令行设置处理后的fasta文件名字
my $output_file = $ARGV[1] or die "Need to input output file name on the command line\n";
# 在命令行设置序列长度阈值，只有大于这个长度的序列才保留,默认是50
my $length_threadhold = ($ARGV[2] or 50);

print "The sequence whose length is no more than $length_threadhold will be removed!\n";

my $catchseq_seqio_obj = Bio::SeqIO->new(-file=>"$input_file", -format=>'fasta');
# 创建一个文件句柄，准备写入处理后的序列
open (O,">$output_file");
while(my $seq_obj = $catchseq_seqio_obj->next_seq)
{
    # 处理每个序列的信息
    my $display_name = $seq_obj->display_name; 
    my $desc = $seq_obj->desc;                                
    my $seq = $seq_obj->seq;                                 
    my $seq_type = $seq_obj->alphabet;                   
    my $seq_length = $seq_obj->length;
    # 长度大于$length_threadhold，就写出新的文件。
    if($seq_length > $length_threadhold)
    {
        my $seqContent = ">".$display_name." ".$desc."\n".$seq."\n";
        # 打印到屏幕的提示信息
        print "Writing formated data to $output_file\n";
        print $seqContent;
        # 写入文件
        print O $seqContent;
        print "success!\n";
    }
}
#关闭文件
close O;

用法如下：

perl removeShortFastaEntries.pl original.fasta formatedEntries.fasta 60

这里的original.fasta，不是一定要后缀名是.fasta，只要文件里面的内容是fasta格式的序列就可以了，名字可以是original.out，original.txt等。最后的数字是可选参数，如果不指定，默认是50。