BamTools API はBAMファイルをRead/writeするクラスと関数を提供してくれるC++ライブラリです。日本語で言及しているエントリが見当たらなかったので、備忘録として簡単なサンプルを載せておきます。

BamToolsインストール

以下のサイトからダウンロード
https://github.com/pezmaster31/bamtools

make installすると、API用ライブラリが/usr/local/libにインストールされます。

利用例

最も簡単なケース

#include <api/BamReader.h>
#include <api/BamAlignment.h>

using namespace BamTools;

int main()
{
  std::string bamFile("sam_spec_example.bam");
  BamReader reader;
  if (!reader.Open(bamFile)) {
    std::cerr << "Failed to open the input file: " << bamFile << std::endl;
    exit(1);
  }

  uint64_t mappedReads(0);
  uint64_t unmappedReads(0);

  reader.Rewind();

  // 1行ごとにパース
  BamAlignment read;
  while (reader.GetNextAlignment(read)) {
    if (read.IsMapped()) ++mappedReads;
    else ++unmappedReads;
  }
  reader.Close();

  std::cout << "mapped reads: " << mappedReads
	    << "\nunmapped reads: " << unmappedReads
	    << std::endl;
  return 0;
}

いろいろ出力するパターン

#include <api/BamReader.h>
#include <api/BamAlignment.h>

using namespace BamTools;

int main()
{
  std::string bamFile("sam_spec_example.bam");
  BamReader reader;
  if (!reader.Open(bamFile)) {
    std::cerr << "Failed to open the input file: " << bamFile << std::endl;
    exit(1);
  }
  // 入力ファイル名
  std::cout << "filename: " << reader.GetFilename() << std::endl;
  // SAM header
  SamHeader header = reader.GetHeader();
  std::string headerStr = reader.GetHeaderText();
  std::cout << "header: " << headerStr << std::endl;
  std::cout << "header: " << header.ToString() << std::endl;
  // リファレンスfasta (主に染色体)の数
  uint64_t count(reader.GetReferenceCount());
  std::cout << "count: " << count << std::endl;

  uint64_t mappedReads(0);
  uint64_t unmappedReads(0);
  uint64_t duplicatedReads(0);
  uint64_t forwardReads(0), reverseReads(0);

  reader.Rewind();
  // リファレンス染色体情報の表示
  const RefVector &refData = reader.GetReferenceData();
  for (auto &x: refData) {
    std::cout << x.RefName << "\t" << x.RefLength << std::endl;
  }

  // 1行ごとにパース
  BamAlignment read;
  while (reader.GetNextAlignment(read)) {
    if (read.IsMapped()) {
      std::cout << "refID: " << read.RefID
                << " " << refData[read.RefID].RefName
                << "\tName: " << read.Name
                << "\tFlag: " << read.AlignmentFlag
                << "\tPosition: " << read.Position
                << "\tLength: " << read.Length
                << "\tEndPosition: " << read.GetEndPosition()
                << "\tQuery: " << read.QueryBases
                << "\tAligned: " << read.AlignedBases
                << "\tTagData: " << read.TagData
                << std::endl;

      ++mappedReads;
      if (read.IsDuplicate()) ++duplicatedReads;
      if (read.IsReverseStrand()) ++reverseReads; else ++forwardReads;
    }
    else ++unmappedReads;
  }
  reader.Close();

  std::cout << "mapped reads: " << mappedReads
	    << "\nduplicated reads: " << duplicatedReads
	    << "\n+ reads: " << forwardReads
	    << "\t- reads: " << reverseReads
	    << "\nunmapped reads: " << unmappedReads
	    << std::endl;
  return 0;
}

コンパイル

make installしてあれば、以下のようなコマンドでコンパイルできると思います。

clang++ bamtoolsAPI.cpp -I/usr/local/include/ -L/usr/local/lib -lbamtools -lz

感想

使いやすい関数が揃っているのでコードは書きやすいと思います。ただ残念ながらCRAMに対応していません。あと、実行速度はあんまり早くないようです。