Perl スクリプトを利用して FASTQ ファイルを処理する

FASTQ ファイルの処理

シーケンサーを利用してシーケンシングした塩基配列は FASTQ と呼ばれる形式のファイルに保存される。このファイル形式は、1 エントリーの情報を 4 行で記録している。1 行目は @ で始まり、その配列に関する情報が記載されている。2 行目はシーケンシングした配列が記載されている。3 行目は + から始まり、その配列についてのヘッダー情報が記載されている(ほとんどの場合、最初の一文字以外、1 行目と同じ内容)。最後の 4 行目はシーケンシング時のクオリティが記載されている。

FASTQ ファイルを処理するには、sed などのコマンドや cutadaptTrimmomatic などのソフトウェアを利用すると便利である。Perl の練習であれば話は別だが、わざわざ Perl でプログラムを自作する必要がない場合が多い。

FASTA ファイルへの変換

FASTQ 形式のファイルでは、1 エントリーが 4 行で記載されている。最初の 2 行が塩基配列に関する情報で、後の 2 行がクオリティ情報となっている。最初の 2 行を残して、後の 2 行を削除すればよい。また、FASTQ 形式のファイルの各エントリーの 1 行目は @ で始まっているので、これを FASTA 形式用に > に変更しておく必要がある。

use strict;
 
sub fq2fa {
  my $fq = shift;
  my $i = 0;
  
  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  while (my $buff = <$fh>) {
    $i++;
    if ($i % 4 == 1) {
      $buff =~ s/^@/>/;
      print $buff;
    }
    if ($i % 4 == 2) {
      print $buff;
    }
  }

  close ($fh);
}
 

&fq2fa("samplefq1.fq");

短いリードを削除

FASTQ ファイルにある短いリードを削除する例。1 エントリーにつき 4 行で記載されるので、ファイルを読み込むとき、まず 4 行をすべて $tmp 変数に保存する。4 行のデータがすべて揃ってから、配列の長さを判断し、このエントリーを残すかどうかを決める。

use strict;

sub filter_short_len {
  my $fq = shift;
  my $len = shift;
  my $i = 0;

  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  my $tmp = '';
  my $accept = 0;
  while (my $buff = <$fh>) {
    $i++;
    $tmp .= $buff;
    if ($i % 4 == 0) {
      $buff =~ s/[\r|\n]//g;
      if (length($buff) >= $len) {
        print $tmp;
      }
      $tmp = '';
    }
  }

  close ($fh);
}
 

&filter_short_len("samplefq1.fq", 60);

アダプター配列の除去

FASTQ ファイルからアダプター配列を除去するサンプルの例。FASTQ ファイル中にある全配列は、同じアダプター配列を持ち、かつそのアダプター配列が途中で切れたり、ミスマッチが入ったりしていない場合に、それを削除する。

アダプターの除去は cutadaptTrimmomatic などのオープンソースソフトウェアを利用すると便利。これらのソフトウェアは、アダプター配列にミスマッチが入ってたり、途中で切れたりする場合にも対応している。

use strict;
 
sub trim {
  my $fq = shift;
  my $adp = shift;
  my $adp_len = length($adp);
  my $i = 0;

  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  while (my $buff = <$fh>) {
    $i++;
    if ($i % 4 == 2) {
      $buff =~ s/^$adp//g;
    } 
    elsif ($i % 4 == 0) {
      $buff = substr($buff, $adp_len);
    }
    print $buff;
  }
 
  close ($fh);
}
 

&trim("samplefq1.fq", "GACTACGTAC");