R でバイナリファイルを読む練習をした。クロマトグラフィーのデータや結晶学のデータファイル(MTZファイル)など、バイナリファイルを直接読めると便利だ。
R では、file 関数でファイルを開き、seek で目的の箇所に移動して、readBin でデータを読み込む。仕事が終わったら close する。
今回は MTZ ファイルを読み込んでみる。MTZファイルを本格的に扱うときは、ccp4 のライブラリ cmtzlib を Rcpp から使うのが本道だと思われ、実際 yam_cpp 氏による実装 (http://snipt.org/zog8) もあるが、今回は全てRで実装してみる。ファイル形式は、実際に自分でパースしてみないと理解が深まらない気がする。MTZ 形式の仕様書は CCP4 Program Suite: mtz format にあるのだが、これが非常に分かりにくく、Emacs のバイナリエディタモード(M-x hexl-find-file)でファイルを眺めながら解析した。
MTZ ファイル内部のデータはすべて REAL*4 型である。Fortran には詳しくないのだが、いわゆる float 型 (32bit 浮動小数点数)らしい。このとき、readBin に与える type は double() であり、size=4 として float 型にする。double() を指定するというのが紛らわしいと感じる。なお、指数のように整数値しか取らないものでも Integer 型は使われていない。
文字列はヌル終端になっていないので、type=character() では読み込めない。raw 型として読み込んでおいて、rawToChar で変換する必要がある。
オフセットは 1-index になっているようだ。しかも、"Reflection data starting at byte 21" のように書いてあるが、実際には 21バイト目からではなく、21 word目からという意味らしく、(21 - 1) * sizeof(REAL*4) = 80 バイト目から始まる。また、ヘッダがファイルの末尾にあるという謎仕様。ヘッダを見てカラム数を特定しないとデータが読めないので、seek を使うことが必須となる。 訂正: カラム数は分からなくても、全部まとめて読んでしまえばよい。ヘッダがレコード1つにつき80バイトであることはどこにも明記していない。COLSRC レコードは最近追加されたらしく、上記仕様書では触れられていない。CRYSTAL や DATASET といった階層構造に関係するらしいのだが。その辺は今回は無視している。
と、愚痴っぽくなってしまったが、そういう謎解きも含めて楽しみつつ、以下のコードができた。最近の ccp4 や phenix が出力したファイルは読めているようだが、もちろんいいかげんな実装なので、全てのファイルとの互換性は保証しない。また、エンディアン・フィールドの値の意味の説明が見つからなかったので、エンディアンについても配慮していない。
hklin <- file("test.mtz", "rb") magic <- readBin(hklin, raw(), 4) if (rawToChar(magic) != "MTZ ") { stop("Not a MTZ file") } # TODO: Endian の違いをサポートする header_offset <- (readBin(hklin, integer(), 1) - 1) * 4 # 1-indexed seek(hklin, header_offset, origin="start") columns <- c() while (TRUE) { header_record <- rawToChar(readBin(hklin, raw(), 0x50)) if (substr(header_record, 1, 3) == "END") break; if (substr(header_record, 1, 6) == "COLUMN") { columns <- c(columns, gsub(" ", "", substr(header_record, 8, 30))) } } columns start <- (21 - 1) * 4 n_row <- (header_offset - start) / 4 / length(columns) seek(hklin, start, origin="start") data <- matrix(readBin(hklin, double(), n=n_row * length(columns), size=4), nrow=n_row, byrow=T) colnames(data) <- columns close(hklin) head(data) summary(data) table(data[, "FreeR_flag"])
