利用Bioperl的SeqIO模块解析fastq文件

　　测序数据中经常会接触到fastq格式的文件，比如说拿到fastq格式的原始数据后希望查看测序碱基的质量并去除低质量碱基。一般而言大家都是用现有的工具，比如说fastqc这个Java写的小程序，确实很好用，运行速度快，检查的项目也多。有时候我们也需要对这些数据进行个性化的分析，那么这个时候这些小工具就不能胜任了，需要我们自己写程序（脚本）来处理。本人目前才疏学浅，因此只有一下三种方案：

1.完全自己写脚本，读取每一行，手动解析，然后实现个性化分析。（显然这个比较慢，相当于重造了一个转速很慢的轮子）

2.利用前人写好的工具，找出源码里面的核心解析程序，然后加以改进，实现个性化。（当然这个只能自己私下用用。这里推荐一个C语言的库 kseq.h，可以用来解析fasta/fastq格式的文件，底层语言运行速度非常快！）

3.利用Bioperl或者Biopython里面的工具解析文件，然后再写脚本个性化分析。（鉴于python的速度，这里推荐Bioperl）

　　下面具体介绍如何使用Bioperl的SeqIO模块解析fastq格式文件。

　　首先是安装Bioperl。

sudo apt-get install perlbrew

sudo perlbrew install-cpanm

sudo /path/cpanm   Bio::Perl

　　解析 head10000.fastq 文件的前四行（第一条序列）。

#!/bin/perl -w
use Bio::SeqIO;
use Bio::Seq::Quality;
my $in = Bio::SeqIO->new(-format => "fastq",
                         -variant => "sanger",
                         -file => 'head10000.fastq' );
while(my $seq = $in->next_seq){

    print  $seq->id,"
";
    print  $seq->seq,"
";
    print  $seq->length,"
";
    print "@{$seq->qual}","
";
    last；
}

　　运行结果如下：

E00552:20:HHCM5ALXX:2:1101:9404:1573
NTCGAAACGGCGGATCATGCCAGGCTGCAACTGCAGCTGGCCTACAACTGGCACTTTGAGGTGAATGACCGGAAGGACCCCCAAGAGACGGCCAAGCTCGTTTCAGTGCCAGACTTTGTAGGTGATGCCTGCAAAGCCATCGCATCCCGG
150
2 32 32 27 37 41 41 41 37 41 41 41 41 37 41 41 41 41 41 41 41 41 41 37 41 41 41 41 27 12 37 27 41 37 22 41 41 41 41 41 41 41 41 41 41 37 27 41 37 22 37 41 37 32 41 37 41 41 41 41 41 41 27 37 12 37 41 37 32 41 41 37 41 37 41 37 32 41 41 41 41 41 41 41 32 32 37 37 27 32 41 37 22 32 37 41 41 32 37 12 32 41 41 41 37 41 41 41 41 41 41 41 41 37 41 41 37 27 32 32 37 41 41 37 41 41 41 37 32 37 37 32 32 37 37 41 12 32 37 41 41 32 27 12 22 37 8 12 27 22

　　与此同时，我们查看 head10000.fastq 文件的前四行：

@E00552:20:HHCM5ALXX:2:1101:9404:1573 1:N:0:TACAGCAT
NTCGAAACGGCGGATCATGCCAGGCTGCAACTGCAGCTGGCCTACAACTGGCACTTTGAGGTGAATGACCGGAAGGACCCCCAAGAGACGGCCAAGCTCGTTTCAGTGCCAGACTTTGTAGGTGATGCCTGCAAAGCCATCGCATCCCGG
+
#AA<FJJJFJJJJFJJJJJJJJJFJJJJ<-F<JF7JJJJJJJJJJF<JF7FJFAJFJJJJJJ<F-FJFAJJFJFJFAJJJJJJJAAFF<AJF7AFJJAF-AJJJFJJJJJJJJFJJF<AAFJJFJJJFAFFAAFFJ-AFJJA<-7F)-<7

　　对照ascii码表可以发现，运行结果的最后一行，即为原始文件的第四行的ascii码对应的十进制数值减去33。例如 "#"(35) - 33 = 2；"A"(65) - 33 = 32；“<”(60) - 33 = 27。也就是说这里的碱基质量用的是phred33。

　　最后解释一下这几行命令的意思：

    print  $seq->id;              #打印$seq对象的序列ID；
    print  $seq->seq;　　　　　　　 #打印$seq对象的序列碱基；
    print  $seq->length;          #打印$seq对象的序列长度；
    print "@{$seq->qual}";　　　　 #打印$seq对象的序列质量；