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
も存在しますがこれはhtslib
のbgzip
コマンドのような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等と異なりコンパイルを通すまでが一苦労ですが、安全で速い言語ということなので、選択肢の一つとしていきたいと思いました。