FASTQ ファイルの処理

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

FASTQ を処理にするには FASTX Toolkit などがある。本格にデータを処理するならばそれらのツールを利用した方が正確で効率がよい。ここでは、練習の意味合いとして、いくつかの FASTQ ファイルの処理例を示す。

FASTA ファイルへの変換

FASTQ は 1 配列が 4 行によって記載されている。最初の 2 行が塩基配列に関する情報で、後の 2 行がクオリティ情報となっている。これを FASTA に変換するには、最初の 2 行を残して、後の 2 行を削除すればよい。(サンプル FASTQ

FASTQ から FASTA への変換は awk コマンド 1 行で行える。

use strict;
 
sub fq2fa {
  my $fq = shift;   # FASTQファイルの場所を受け取る
  my $i = 0;        # FASTQの何行目を読んでいるかを記録
  
  # FASTQファイルを開く
  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  # ファイルを1行ずつ読み込む
  while (my $buff = <$fh>) {
    $i++;
 
    # 現在の行数を4で割り、
    # 余りが 1 ならばヘッダー、ヘッダーならば 1 文字目を > に変更
    # 余りが 2 ならば塩基配列
    if ($i % 4 == 1) {
      $buff =~ s/^@/>/;
      print $buff;
    }
    if ($i % 4 == 2) {
      print $buff;
    }
  }
  # ファイルを閉じる
  close ($fh);
}
 
 
my $fq = "samplefq1.fq"; # FASTQファイルの場所
&fq2fa($fq);

短いリードを削除

FASTQ ファイルにある短いリードを削除する例。1 配列につき 4 行で記載されるので、ファイルを読み込むとき、まず 4 行をすべて $tmp 変数に保存する。最後に、4 行目に達した時、配列の長さを判定し出力する。(サンプル FASTQ

use strict;

sub filter_short_len {
  my $fq = shift;    # FASTQファイル
  my $len = shift;   # 許容できる最短の配列の長さ
  my $i = 0;         # FASTQの何行目を読んでいるかを記録

  # ファイルを開く
  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  my $tmp = '';      # データを一時保存するための変数
  my $accept = 0;    # 長さ条件を満たすかどうか
 
  # ファイルを1行ずつ読み込む
  while (my $buff = <$fh>) {
    $i++;

    # ファイルの内容を読み込んで$tmpに一時保存
    $tmp .= $buff;

    # 1配列(4行)が読み終えた時、長さをチェックして一時保存したデータを出力
    if ($i % 4 == 0) {
      # 改行コードを削除(1文字分とカウントされるのを回避するため)
      $buff =~ s/[\r|\n]//g;
      # 配列長が$len以上ならば一時保存したデータを出力する
      if (length($buff) >= $len) {
        print $tmp;
      }
      # $tmp変数を空にする
      $tmp = '';
    }
  }

  # ファイルと閉じる
  close ($fh);
}
 
 


my $fq = "samplefq1.fq";     #FASTQファイルの場所
&filter_short_len($fq, 60);  #60を満たないリードを削除

アダプター配列の除去

FASTQ ファイルからアダプター配列を除去するサンプルの例。すべての配列はアダプター配列を持つという前提の処理例である。また、ミスマッチを許容していない。(サンプル FASTQ

use strict;
 
sub trim {
  my $fq = shift;              # FASTQファイル
  my $adp = shift;             # アダプター配列
  my $adp_len = length($adp);  # アダプター配列の長さ
  my $i = 0;                   # FASTQファイルの何行目を読んでいるかを記録

  # FASTQファイルを開く
  open (my $fh, '<', $fq) or die "cannot open FASTQ file.";
 
  # FASTQファイルを1行ずつ読む
  while (my $buff = <$fh>) {
    $i++;

    # 各配列について2行目のデータが塩基配列であるから、2行目の塩基
    # 配列の先頭がアダプター配列と一致するならば、空文字に置換する。
    if ($i % 4 == 2) {
      $buff =~ s/^$adp//g;
    } 
    # 4行目も同様にアダプター配列の部分のクオリティ情報を削除する。
    elsif ($i % 4 == 0) {
      $buff = substr($buff, $adp_len);
    }
    print $buff;
  }
 
  # ファイルを閉じる
  close ($fh);
}
 
 
 
my $fq = "samplefq1.fq";  # FASTQファイルの場所
my $adp = "GACTACGTAC";   # アダプター配列
&trim($fq, $adp);