利用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,"\n";
    print  $seq->seq,"\n";
    print  $seq->length,"\n";
    print "@{$seq->qual}","\n";
    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對象的序列質量；