Rust-BioでGzip圧縮されたFASTAを読み込む

Rustが注目を集めているらしいので流行りにのって使ってみました。

RustにはRust-Bioというバイオインフォマティクス用パッケージが存在します。そこで、今回はこれを用いてGzip圧縮されたFASTAファイルを読み込むプログラムを書いてみます。

まずはコードの全体から。

use structopt::StructOpt;
use std::fs::File;
use bio::io::fasta::{self, FastaRead};
use flate2::bufread::MultiGzDecoder;
use std::io::BufReader;
use std::ffi::OsStr;

#[derive(StructOpt)]
struct Opt {
    fasta: std::path::PathBuf,
}

fn main() {
    // コマンドライン引数をパース
    let args = Opt::from_args();

    // FASTAファイルを開く
    let file = File::open(&args.fasta).expect("Error during opening the file");

    // FASTAのリーダーを構築する(拡張子がgzの場合にMultiGzDecoderをはさむ)
    let mut reader: Box<dyn FastaRead> = match args.fasta.extension().and_then(OsStr::to_str) {
        Some("gz") => Box::new(fasta::Reader::new(MultiGzDecoder::new(BufReader::new(file)))),
        _ => Box::new(fasta::Reader::new(file)),
    };

    // FASTAのレコード毎にループ
    let mut record = fasta::Record::new();
    loop {
        reader.read(&mut record).unwrap();
        if record.is_empty() {
            break;
        }

        // 各レコードのIDと先頭10baseを表示
        println!("{} {}", record.id(), std::str::from_utf8(&record.seq()[0..10]).unwrap());
    }

}

使用するライブラリをスコープに持ち込む

use structopt::StructOpt;
use std::fs::File;
use bio::io::fasta::{self, FastaRead};
use flate2::bufread::MultiGzDecoder;
use std::io::BufReader;
use std::ffi::OsStr;

flate2::bufread::MultiGzDecoder は圧縮ファイルを解凍しつつ読み込むために使います。 flate2にはGzDecoderも存在しますがこれはhtslibbgzipコマンドのようなBGZF ブロックに分けて圧縮するようなタイプのファイルに対応していないようです。

もう一つ重要なのはbio::io::fastaです。このモジュールが提供するReaderにより FASTAを適切にパースすることができるようになります。

FASTAリーダーを構築する

// FASTAのリーダーを構築する(拡張子がgzの場合にMultiGzDecoderをはさむ)
let mut reader: Box<dyn FastaRead> = match args.fasta.extension().and_then(OsStr::to_str) {
    Some("gz") => Box::new(fasta::Reader::new(MultiGzDecoder::new(BufReader::new(file)))),
    _ => Box::new(fasta::Reader::new(file)),
};

圧縮されていないFASTAをパースするときはfasta::Reader::new(file)gzip圧縮されているFASTAをパースするときはfasta::Reader::new(MultiGzDecoder::new(BufReader::new(file))) を用います。それぞれ場合分けしてコードを書くと冗長になるためここではBoxに入れて 1つの変数に格納できるようにしています。このとき変数の型としてFastaReadトレイトを指定することで、 FastaReadトレイトを実装するオブジェクトを格納することができるようにしています。

Rustは型に厳しいようなのでこの辺がなれるまで大変だと個人的には思いました。

レコードを1つずつ読み込む

// FASTAのレコード毎にループ
let mut record = fasta::Record::new();
loop {
    reader.read(&mut record).unwrap();
    if record.is_empty() {
        break;
    }

    // 各レコードのIDと先頭10baseを表示
    println!("{} {}", record.id(), std::str::from_utf8(&record.seq()[0..10]).unwrap());
}

fasta::Readerオブジェクトにはrecordsメソッドがあるのですが、先ほどのようなBoxでラップする様な方式をとると使えなくなるみたいなので、readメソッドをループ毎に呼び出し、レコードを1つずつ読み込みます。

最後までレコードを読み込むと、Recordオブジェクトのis_emptyメソッドがtrueを返すようになるのでループを抜け出すようにします。

RecordオブジェクトにはidメソッドによってIDに、seqメソッドによって配列にアクセスできます。他にもいくつかのメソッドがありますのでドキュメントを参照してみてください。

以上がGzip圧縮されたFASTAファイルを読み込む(非圧縮ファイルも受け付ける)プログラム解説でした。

Python等と異なりコンパイルを通すまでが一苦労ですが、安全で速い言語ということなので、選択肢の一つとしていきたいと思いました。