faidxと言えばsamtoolsに実装されている機能の1つで、FASTAやFASTQファイルのインデックスを作ったり部分配列を取り出したりするときに便利なプログラムです。今回はこれについて理解を深めるため簡単な実験をしてみました。
以下のようなテストデータを用意しました。
>a AGGAC CAGGG GGGGG GG >abcd AGGGCC AGGGCA NNAGGA
上記テストデータについて、samtools faidxを実行すると以下の内容を記録したfaidxファイルが生成されます。
a 17 3 5 6 abcd 18 30 6 7
これはタブ区切りの表となっていて各列については以下の通りです。(https://www.htslib.org/doc/faidx.htmlより引用) +
NAME 配列の名前 LENGTH 塩基単位で数えた配列の長さ
OFFSET FASTA/FASTQファイルにおける配列の最初の塩基のオフセット
LINEBASES １行に記載されている塩基数
LINEWIDTH 改行コードを含む１行に記載されるバイト数
QUALOFFSET FASTQファイルにおいて最初のクオリティー値のオフセット
NAMEとLENGTHはわかりやすいかと思います。OFFSETというのはその配列がファイルの先頭から何バイト目から始まるかを示しており、先の例では配列aは0-indexで3バイト目から始まっていることがわかります。
実際0-indexで3バイト目から5バイト分のデータ、30バイト目から5バイト分をodコマンドで出力してみると以下のようになります。
$ od -c -j 3 -N 5 x.fa 0000003 A G G A C 0000010 $ od -c -j 30 -N 5 x.fa 0000036 A G G G C 0000043
これらFAIDXファイルに記載される値が事前に与えられると配列の任意の部分配列がバイトで数えてどの位置にあるかがわかります。
それでは上記式に基づいてa:3-8とabcd:14-18を得てみましょう。
samtoolsの場合は以下のようになります。
$ samtools faidx x.fa a:3-8 >a:3-8 GACCAG $ samtools faidx x.fa abcd:14-18 >abcd:14-18 NAGGA
上記式で計算されるインデックスはそれぞれ(start, end) = (6, 12), (46, 51)でありこれは1-indexedなのでodコマンドのオプション引数はそれぞれ-j 5 -N 7、-j 45 -N 6となり、以下のような結果となります。
$ od -c -j 5 -N 7 x.fa 0000005 G A C \n C A G 0000014 $ od -c -j 45 -N 6 x.fa 0000055 N A G G A \n 0000063
ちゃんと欲しい部分配列が取れていると思います。LINEWIDTHが必要な理由としては上記結果でも\nとして見えていますが、改行コード等の見えない文字を加味してバイトインデックスを計算するためです。
そして部分配列のバイトインデックスがわかると何がいいかというと部分配列の抽出が高速にできるということです。
このことを実感するためにまず、インデックスを使用しない場合の部分配列の抽出法を考えてみます。おそらく最もナイーブな方法はファイルに記録された全塩基配列を文字列として読み込み各種プログラミング言語で提供されている部分文字列へのアクセス機能を利用することです。例えばC言語では以下のようになると思います。
void subseq_without_faidx(char filename[], size_t start[], size_t end[], size_t n, int debug) { char* seq2; size_t seq_size = 10000000000; seq2 = malloc(seq_size); if (seq2 == NULL) { fprintf(stderr, "Out of memmory\n"); exit(8); } FILE *fp = fopen(filename, "r"); if (fp == NULL) { fprintf(stderr, "Unable to open the file.\n"); exit(8); } int ch; char line[1001]; char seqid[1001]; int j = 0; while (fgets(line, sizeof(line), fp) != NULL) { if (line[0] == '>') { for (int i = 1; i < strlen(line)-1; ++i) { seqid[i-1] = line[i]; } seqid[strlen(line)-2] = '\0'; continue; } for (size_t i = 0; i < sizeof(line); ++i) { if (line[i] == '\n') { continue; } else if (line[i] == '\0') { break; } else { seq2[j] = line[i]; j++; if (j+1 > seq_size) { fprintf(stderr, "Out of memmory\n"); fclose(fp); free(seq2); exit(8); } } } } seq2[j] = '\0'; fclose(fp); char region[30]; char seq[100000]; for (int i = 0; i < n; ++i) { sprintf(region, "%s:%zu-%zu", seqid, start[i], end[i]); sprintf(seq, "%.*s", (int)(end[i]-start[i]+1), seq2+start[i]-1); if (debug) { printf("%s\n", region); printf("%s\n", seq); } } free(seq2); }
上記はあくまで実装の1例で、配列数や配列長に暗黙の制限があるなど、いろいろとツッコミどころはあると思いますが、これでFASTAファイル名とスタートおよびエンドのインデックスを配列で与えればn個の部分配列を得ることができます。
ただ、この実装では一度ファイル全体を読み込んでから処理を行なっています。コンピュータによる計算において外部記憶装置（HDDやSSD）からメインメモリにデータを読み込むという作業はかなりコストがかかるので時間がかかりそうです。
一方でFAIDXを用いて部分配列を得る場合は以下のような実装になると思います。ここではhtslibを用いています。
void subseq_with_faidx(char filename[], size_t start[], size_t end[], size_t n, int debug) { faidx_t *fai = fai_load(filename); hts_pos_t seq_len; char region[30]; for (int i = 0; i < n; ++i) { sprintf(region, "random:%zu-%zu", start[i], end[i]); char *seq = fai_fetch64(fai, region, &seq_len); if (debug) { printf("%s\n", region); printf("%s\n", seq); } } fai_destroy(fai); }
上記関数を実行するために以下のようなmain関数を実装しました。
int main(int argc, char* argv[]) { if (argc < 2) { fprintf(stderr, "Please specify fasta file.\n"); exit(8); } if (argc == 2 || argv[2][0] == '0') { // debug mode size_t start[] = {1, 10, 100}, end[] = {10, 40, 200}; clock_t start_clock, end_clock; printf("with faidx\n"); subseq_with_faidx(argv[1], start, end, 3, 1); printf("without faidx\n"); subseq_without_faidx(argv[1], start, end, 3, 1); } else { // ランダムにn回部分配列を取り出し時間を計測 int n = atoi(argv[2]); srand((unsigned int)time(NULL)); size_t start[n], end[n]; for (int i = 0; i < n; ++i) { // 文字列長,部分文字列配列要素数は1000000000,100000なのでそれを超えないようにする start[i] = rand() % (1000000000 - 200000) + 1; end[i] = rand() % (99000) + start[i]; } clock_t start_clock, end_clock; start_clock = clock(); subseq_with_faidx(argv[1], start, end, n, 0); end_clock = clock(); printf("clock (with faidx): %f\n", (double)(end_clock - start_clock) / CLOCKS_PER_SEC); start_clock = clock(); subseq_without_faidx(argv[1], start, end, n, 0); end_clock = clock(); printf("clock (without faidx): %f\n", (double)(end_clock - start_clock) / CLOCKS_PER_SEC); } return 0; }
この関数を第3引数なしで実行してみると以下の出力が得られます。
$ ./a.out test.fa with faidx random:1-10 CTTGAAGGGC random:10-40 CTCCTTTCGAACGTTAGGATATACGTCAATA random:100-200 AATGATCCTATGCCGGACCCCCCGGGTTGCATCTATGTGACTCCCACAATCGTAGGTAGATGATCGGTGAAGCCCGCATGCCCATTACGATGCGAGTTCAG without faidx random:1-10 CTTGAAGGGC random:10-40 CTCCTTTCGAACGTTAGGATATACGTCAATA random:100-200 AATGATCCTATGCCGGACCCCCCGGGTTGCATCTATGTGACTCCCACAATCGTAGGTAGATGATCGGTGAAGCCCGCATGCCCATTACGATGCGAGTTCAG
無事FAIDX使用バージョンと不使用バージョンで結果が一致しているので、実際にどのくらいの時間がかかるのかを測ってみたいと思います。テスト用のファイルとして1Gbのランダム配列1つが格納されたFASTAファイルtest.faを用意しました。1000回、ランダムな位置についての部分配列を取得する操作を行なったときにかかる時間を測ってみます。
$ ./a.out test.fa 1000 clock (with faidx): 2.560794 clock (without faidx): 1.334881 $ ./a.out test.fa 1000 clock (with faidx): 0.155571 clock (without faidx): 1.341478
プログラムを2回連続して実行しています。面白いことに1回目はFAIDXを使用した方が遅いのに、2回目では格段にはやくなっています。これは1回目では.faiファイルが作成されていないので、 FAIDX使用バージョンもFASTAファイル全体を一旦読み込んでインデックスを作成する必要があるためです。2回目以降は作成済みのインデックスをもとに部分配列のみを対象にメモリにロードするので、ファイル全体をロードする場合より高速に処理ができています。
それでは先ほど1000回の部分配列抽出の速度を比較しましたが、100000回の場合も見てみましょう。今回はtest.fa.faiファイルが存在する状態で実行します。
$ ./a.out test.fa 100000 clock (with faidx): 13.303650 clock (without faidx): 1.738566
すると今度はインデックス存在下でもFAIDX使用バージョンの方が遅くなりました。これは100000と部分配列アクセス回数が増えたことで、ファイル全体を読み込む以上の量のデータ読み込みが発生したためと考えられます。このように大量の部分配列にアクセスする場合はいっそファイル全体をメインメモリに読み込んでメインメモリ上で何回もアクセスする方がコストが低くなります。
最後に使用したCプログラムの全体を載せて締めたいと思います。
#include <stdio.h> #include <stdlib.h> #include <time.h> #include <string.h> #include <htslib/faidx.h> void subseq_with_faidx(char filename[], size_t start[], size_t end[], size_t n, int debug) { faidx_t *fai = fai_load(filename); hts_pos_t seq_len; char region[30]; for (int i = 0; i < n; ++i) { sprintf(region, "random:%zu-%zu", start[i], end[i]); char *seq = fai_fetch64(fai, region, &seq_len); if (debug) { printf("%s\n", region); printf("%s\n", seq); } } fai_destroy(fai); } void subseq_without_faidx(char filename[], size_t start[], size_t end[], size_t n, int debug) { char* seq2; size_t seq_size = 10000000000; seq2 = malloc(seq_size); if (seq2 == NULL) { fprintf(stderr, "Out of memmory\n"); exit(8); } FILE *fp = fopen(filename, "r"); if (fp == NULL) { fprintf(stderr, "Unable to open the file.\n"); exit(8); } int ch; char line[1001]; char seqid[1001]; int j = 0; while (fgets(line, sizeof(line), fp) != NULL) { if (line[0] == '>') { for (int i = 1; i < strlen(line)-1; ++i) { seqid[i-1] = line[i]; } seqid[strlen(line)-2] = '\0'; continue; } for (size_t i = 0; i < sizeof(line); ++i) { if (line[i] == '\n') { continue; } else if (line[i] == '\0') { break; } else { seq2[j] = line[i]; j++; if (j+1 > seq_size) { fprintf(stderr, "Out of memmory\n"); fclose(fp); free(seq2); exit(8); } } } } seq2[j] = '\0'; fclose(fp); char region[30]; char seq[100000]; for (int i = 0; i < n; ++i) { sprintf(region, "%s:%zu-%zu", seqid, start[i], end[i]); sprintf(seq, "%.*s", (int)(end[i]-start[i]+1), seq2+start[i]-1); if (debug) { printf("%s\n", region); printf("%s\n", seq); } } free(seq2); } int main(int argc, char* argv[]) { if (argc < 2) { fprintf(stderr, "Please specify fasta file.\n"); exit(8); } if (argc == 2 || argv[2][0] == '0') { // debug mode size_t start[] = {1, 10, 100}, end[] = {10, 40, 200}; clock_t start_clock, end_clock; printf("with faidx\n"); subseq_with_faidx(argv[1], start, end, 3, 1); printf("without faidx\n"); subseq_without_faidx(argv[1], start, end, 3, 1); } else { // ランダムにn回部分配列を取り出し時間を計測 int n = atoi(argv[2]); srand((unsigned int)time(NULL)); size_t start[n], end[n]; for (int i = 0; i < n; ++i) { // 文字列長,部分文字列配列要素数は1000000000,100000なのでそれを超えないようにする start[i] = rand() % (1000000000 - 200000) + 1; end[i] = rand() % (99000) + start[i]; } clock_t start_clock, end_clock; start_clock = clock(); subseq_with_faidx(argv[1], start, end, n, 0); end_clock = clock(); printf("clock (with faidx): %f\n", (double)(end_clock - start_clock) / CLOCKS_PER_SEC); start_clock = clock(); subseq_without_faidx(argv[1], start, end, n, 0); end_clock = clock(); printf("clock (without faidx): %f\n", (double)(end_clock - start_clock) / CLOCKS_PER_SEC); } return 0; }
