htslibを使う(1) bamを読む
この記事では、SAM/BAM/CRAMを読んだり書いたりするC言語のソフトウェアhtslibの使い方を紹介したいと思います。でも実は、紹介というより筆者がhtslibを使うために勉強した事柄の忘備録という感じです。
注意：本記事の正確性は担保しておりません。あくまで筆者のパソコンとデータでは動きました。バグ・間違い・より良い方法等がありましたらお知らせください。
TLDR
今回の目的は、C++でファイルをMulti-threadingで読み込み、さらにアライメントのデータを取得することです。samtools コマンドでできること以上を実行したい中級者向けの内容です。
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <vector>
#include "htslib/sam.h"
#include "htslib/thread_pool.h"
using namespace std;
/*
usage: %prog in.bam thread_num
*/
int main(int argc, char *argv[]) {
int nthread;
if (argc <= 1) {
cout << "Please specify a bam file (usage: %prog in.bam thread_num)." << endl;
return 0;
} else if (argc <= 2) {
nthread=1;
cout << "Set nthread = 1." << endl;
} else {
nthread=atoi(argv[2]);
}
cout << "Reading " << argv[1] << " with " << nthread << " threads..." << endl;
// open bam
char *f=argv[1];
htsFile *in=hts_open(f, "r");
sam_hdr_t *h=sam_hdr_read(in);
// threading
htsThreadPool p = {NULL, 0};
p.pool = hts_tpool_init(nthread);
hts_set_opt(in, HTS_OPT_THREAD_POOL, &p);
// read bam
bam1_t *b= bam_init1();
kstring_t aux={0, 0, NULL};
while (sam_read1(in, h, b) >= 0) { // end file = -1
// details can be found in sam.c bam_read1() and sam.h
// read
char *chr = h->target_name[b->core.tid]; // chr name
hts_pos_t &start = b->core.pos; // left position, 0-based
hts_pos_t end = bam_endpos(b); // right position, 0-based
char *qname = bam_get_qname(b); // read name
int32_t &l_qseq = b->core.l_qseq; // length of read
uint16_t &l_qname = b->core.l_qname; // length of read name
uint8_t &mapq = b->core.qual; // MAPQ
uint16_t &flag = b->core.flag; // SAM flag
uint32_t *cigar = bam_get_cigar(b); // cigar, 32-bit
uint32_t &n_cigar = b->core.n_cigar; // number of CIGAR operations
hts_pos_t qlen = bam_cigar2qlen(n_cigar, cigar); // qlen
hts_pos_t rlen = bam_cigar2rlen(n_cigar, cigar); // rlen
uint16_t &bin = b->core.bin; // bin calculated by bam_reg2bin()
hts_pos_t &isize = b->core.isize; // insertion size (beween R1 and R2)
// mate
char *mchr = h->target_name[b->core.mtid]; // mate chr name
hts_pos_t &mpos = b->core.mpos; // mate left position, 0-based
// format seq
uint8_t *tmp_s = bam_get_seq(b); // seq of read, nt16
char *seq=(char *)malloc(l_qseq + 1); // seq of read, ATGCN
for (int i=0; i < l_qseq; i++) {
seq[i]=seq_nt16_str[bam_seqi(tmp_s, i)]; // get nucleotide id and convert into IUPAC id
}
seq[l_qseq]='\0';
// format cigar
string cigarstr;
for (int k = 0; k < n_cigar; ++k) {
int leng=bam_cigar_oplen(cigar[k]);
char opchr=bam_cigar_opchr(cigar[k]);
cigarstr += to_string(leng) + opchr;
}
// format query quality
uint8_t *tmp_q = bam_get_qual(b); // query quality
char *qqual;
if (tmp_q[0] == 0xff) {
qqual=(char *)malloc(2);
qqual[0]='*';
qqual[1]='\0';
} else {
qqual=(char *)malloc(l_qseq + 1);
for (int i = 0; i < l_qseq; ++i)
qqual[i]=tmp_q[i]+33;
qqual[l_qseq]='\0';
}
// format auxiliary tags
uint8_t *tmp_a = bam_get_aux(b); // auxiliary data
uint8_t *tmp_end = b->data + b->l_data;
vector<char *> auxs;
while (tmp_a && tmp_end - tmp_a >= 4) {
aux={0, 0, NULL};
tmp_a=(uint8_t *)sam_format_aux1(tmp_a, tmp_a[2], tmp_a+3, tmp_end, &aux);
auxs.push_back(aux.s);
}
// parse flag
bool is_paired = (flag & BAM_FPAIRED) > 0;
bool is_proper_pair = (flag & BAM_FPROPER_PAIR) > 0;
bool is_unmapped = (flag & BAM_FUNMAP) > 0;
bool is_mate_unmapped = (flag & BAM_FMUNMAP) > 0;
bool is_reverse = (flag & BAM_FREVERSE) > 0;
bool is_mate_reverse = (flag & BAM_FMREVERSE) > 0;
bool is_read1 = (flag & BAM_FREAD1) > 0;
bool is_read2 = (flag & BAM_FREAD2) > 0;
bool is_secondary = (flag & BAM_FSECONDARY) > 0;
bool is_qc_failed = (flag & BAM_FQCFAIL) > 0;
bool is_duplicate = (flag & BAM_FDUP) > 0;
bool is_supplementary = (flag & BAM_FSUPPLEMENTARY) > 0;
// output
cout << "chr " << chr << endl;
cout << "start " << start << endl;
cout << "end " << end << endl;
cout << "qname " << qname << endl;
cout << "l_qseq " << l_qseq << endl;
cout << "mapq " << mapq << endl;
cout << "flag " << flag << endl;
cout << "cigar " << *cigar << endl;
cout << "n_cigar " << n_cigar << endl;
cout << "cigarstr " << cigarstr.c_str() << endl;
cout << "qlen " << qlen << endl;
cout << "rlen " << rlen << endl;
cout << "bin " << bin << endl;
cout << "mchr " << mchr << endl;
cout << "mpos " << mpos << endl;
cout << "isize " << isize << endl;
cout << "seq " << seq << endl;
cout << "qqual " << qqual << endl;
for (char *cp : auxs) {
cout << "aux " << cp << endl;
}
cout << "is_paired " << is_paired << endl;
cout << "is_proper_pair " << is_proper_pair << endl;
cout << "is_unmapped " << is_unmapped << endl;
cout << "is_mate_unmapped " << is_mate_unmapped << endl;
cout << "is_reverse " << is_reverse << endl;
cout << "is_mate_reverse " << is_mate_reverse << endl;
cout << "is_read1 " << is_read1 << endl;
cout << "is_read2 " << is_read2 << endl;
cout << "is_secondary " << is_secondary << endl;
cout << "is_qc_failed " << is_qc_failed << endl;
cout << "is_duplicate " << is_duplicate << endl;
cout << "is_supplementary " << is_supplementary << endl;
cout << endl;
}
// close bam
sam_hdr_destroy(h);
sam_close(in);
bam_destroy1(b);
// close threads
hts_tpool_destroy(p.pool);
return 0;
}本スクリプトの使い方
まずはコンパイルします。
g++ -o read_bam -I /path/to/htslib-1.13 -L /path/to/htslib-1.13 this_script.cpp -lhts
使い方は、read_bam [input file] [number of threads]で、例えばこのようになります。./read_bam input_file.bam 4
htslibのバージョンは1.13、g++は7.5です。
ファイルの読み込み準備
まず、ファイルを開き htsFile オブジェクトを作成します。次にヘッダーを読み込み sam_hdr_t オブジェクトを作成します。次にアライメントを1行ずつ読み込むため、bam1_t オブジェクトを作成します。
ファイルを読み込む
アライメントを読む準備ができたら、sam_read1 関数で bam1_t オブジェクトにデータを入れていきます。基本的な情報 (e.g., ゲノムでのアライメントの開始位置)はbam1_t オブジェクトのcoreに存在します (b->core)。その他の情報もbam_endpos関数等で取得可能です。
配列の取得
少し注意が必要なのは、リードの配列が4ビットで記録されている点です。そこで、4ビットの配列情報を人が読める文字列に変換する必要があります。変換にはseq_nt16_str関数を用います (hts.h を参照)。
Multi-threadingで読み込む
ファイルをmulti-threadingで読み込むのは非常に簡単です (ありがとうhtslibの開発者！)。まずhtsThreadPool オブジェクトを作成し、hts_tpool_initでスレッド数を指定します。次にhts_set_optで読み込むファイルとスレッド数をセットします。読み込みが終わったらhts_tpool_destroyでmulti-threadingを終了します。内部ではpthreadsが動いているようです (まだこちらは勉強不足)。
こちらのmulti-threadingは読み込みのスレッド数を指定します。つまり、2スレッドで読み込んで、結果を1スレッドで処理すると、合計3スレッドの使用となります。
ファイルの読み込みの終了
ファイルを最後まで読むと、sam_read1 関数は-1を返します。そこで読み込みは終了です。ファイルを全部読み終わったら、sam_hdr_destroy、sam_close、bam_destroy1関数で作ったオブジェクトを消します。
最後に
大規模並列シーケンス時代に生きる研究者として、筆者はこれからhtslibをより深く勉強していこうと思っています。応援コメント等ありましたら是非よろしくお願いします！
筆者
小嶋将平 @ ウイルス学若手ネットワーク / RIKEN IMS
本記事はウイルス学若手ネットワークのメンバーが書いています。ウイルス学内外の若手の研究のお役に立てればと思っています。
ご意見やご感想、間違い等がありましたら以下のツイッターのDMにご連絡ください。もちろんフォローも忘れずにお願いします！！
以下、免責事項です。ご一読ください。
本ページの内容の正確性は一切保証しません。本情報は著者の所属機関を代表する見解ではありません。本情報を掲載された情報・資料を利用、使用、ダウンロードするなどの行為に関連して生じたあらゆる損害等に関しても、その理由に関わらず一切責任を負いません。本ページの内容は研究手法の解説および研究を行う上での知識共有が目的であり、それ以外の意図は一切ありません。本ページのコンテンツは、予告なく変更や公開の取り消しする場合があります。
この記事が気に入ったらサポートをしてみませんか？