【pysamの使い方】SAM/BAMファイルを変換する

次世代シーケンサーの解析において、SAM/BAM形式のファイル操作は非常に重要となります。

そこで本記事では、SAM/BAMファイルの読み込みや変換ができるPythonのライブラリであるpysamを使い方を解説します。

SAMファイルとは

シーケンスのアライメント結果を扱う際に広く利用される形式がSAM（Sequence Alignment/Map）ファイルです。

SAMファイルは大きく「ヘッダー部」と「アライメント部」の2つのセクションに分かれています。

ヘッダー部

ヘッダー部はファイルの先頭にあり、各行が「@」で始まります。主なヘッダー行の種類は以下の通りです。

• @HD (Header): ファイル全体の情報（フォーマットバージョン、ソート順など）
• @SQ (Sequence): 参照配列の情報。各参照配列（例：染色体）の名前や長さを記述
• @PG (Program): 使用プログラム

アライメント部

ヘッダー部の後に続くのがアライメント部で、各行が1つのリードのアライメント情報を示します。

以下の表は、SAMファイルのアライメント部における必須11項目のフィールド名、説明、および例をまとめたものです。

1. FLAG値の解釈

具体例を用いて、説明します。

例1: FLAG = 99

99 (10進数) は、次のビットの和として表されます。

• 0x1 (1)

• 0x2 (2)

• 0x20 (32)

• 0x40 (64)

つまり、FLAG=99 の場合、以下の情報が含まれています：

• ペアエンドリードである（0x1）

• 正しくペアになっている（0x2）

• このリードがテンプレート内の最初のリードである（0x40）

• ペアのリードはリバースコンプリメントでマッピングされている（0x20）
  また、FLAG 99 には 0x10 が立っていないため、このリード自体は正方向にマッピングされていると解釈されます。

例2: FLAG = 147

147 は次のビットの和です。

• 0x1 (1)

• 0x2 (2)

• 0x10 (16)

• 0x80 (128)

これにより、FLAG=147 の場合は：

• ペアエンドリードである（0x1）

• 正しくペアになっている（0x2）

• このリードがテンプレート内の最後のリードである（0x80）

• このリードがリバースコンプリメントでマッピングされている（0x10）

このように、FLAG値は各ビットの組み合わせによって、リードの状態（マッピングの有無、向き、ペアの情報など）をコンパクトに表現しています。

2. CIGAR文字列の解釈

CIGER文字列についても、具体例をみながら解説していきます。

• 例: “10M1I39M”

  • 意味:

    • 最初の10塩基が参照と整列（マッチまたはミスマッチ）

    • 続いて、リード側に1塩基の挿入がある（参照側には対応なし）

    • その後、39塩基が整列

  • 解説: リードの長さは10 + 1 + 39 = 50塩基ですが、参照配列上では10 + 39 = 49塩基分の領域にマッピングされます。

• 例: “5S45M”

  • 意味:

    • 最初の5塩基はソフトクリッピングされ、整列には含まれないが、リードのSEQには含まれる

    • その後、45塩基が整列

  • 解説: リード全体は50塩基ですが、アライメントとしては45塩基が参照に対してマッピングされます。ソフトクリップされた部分は整列から除外されるため、解析対象外となります。リードの両端など、整列に使われなかった部分を「ソフトクリップ」として扱います。この部分はリードのシーケンス（SEQフィールド）には残っており、後で品質チェックや再解析のために利用できる場合がありますが、アライメント解析時には無視されます。

これらの表を参照することで、SAMファイルのFLAG値やCIGAR文字列が示す情報を正確に解釈し、アライメントの状態やリードの特徴を理解することができます。

SAM/BAM/CRAMのファイルの違い

SAM, BAM, CRAMのファイルの違いについて、以下に示します。

• SAM: テキスト形式で記述されているため、内容が直接確認できるが、可読性が高い反面、ファイルサイズが大きい。

• BAM: バイナリ形式で、直接内容を確認するには専用ツール（例：samtools view）が必要。

• CRAM: BAM同様のバイナリ形式ですが、参照ベースの圧縮を利用しており、さらに高い圧縮率を実現。

pysamとは

pysamは、HTSlibの機能をPythonから利用できるようにしたライブラリです。

• 主な機能
  • SAM/BAM/CRAMファイルの読み込み、書き出し
  • インデックスファイルの作成と利用
  • アライメント情報やシーケンスデータの操作・抽出

これにより、独自のパイプライン構築や解析ツールの作成が容易になります。

pysamのインストール方法

pysamはpipを使って簡単にインストールできます。ターミナルで以下のコマンドを実行してください。

pip install pysam

これで、pysamをインストールできました。

ファイルの読み込み

UCSCからサンプルのBAMファイルをダウンロードします。

下記のコードをpythonで実行して、BAMファイルをダウンロードします。

import urllib
src_file = "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam"
dest_file = "input.bam"
urllib.request.urlretrieve(src_file, dest_file)

import pysam
samfile = pysam.AlignmentFile("input.bam", "rb")
num = 0
for read in samfile.fetch(until_eof=True):
    if num > 5:
        break
    print(read.query_name, read.reference_name, read.mapping_quality, read.flag)

samfile.close()

SOLEXA-1GA-2_2_FC20EMB:5:251:979:328 chr1 25 0 
SOLEXA-1GA-2_2_FC20EMB:5:102:214:278 chr1 25 0
SOLEXA-1GA-2_2_FC20EMB:5:195:284:685 chr1 25 16
SOLEXA-1GA-2_2_FC20EMB:5:35:583:827 chr1 25 0
SOLEXA-1GA-2_2_FC20EMB:5:248:130:724 chr1 25 0
SOLEXA-1GA-2_2_FC20EMB:5:236:644:107 chr1 25 16

SAM/BAMファイル変換

BAMファイルからSAMファイルに変換してファイルを出力します。

# BAM ファイルを読み込み (バイナリモード)
with pysam.AlignmentFile("input.bam", "rb") as bamfile:
    # 同じヘッダー情報を使って SAM ファイルをテキストモードで書き出す
    with pysam.AlignmentFile("output.sam", "w", header=bamfile.header) as samfile:
        # ファイル内の全リードを反復して書き出す
        for read in bamfile.fetch(until_eof=True):
            samfile.write(read)

これで、output.samというSAMファイルに変換されました。SAMファイルはテキスト形式なので読むことができます。

先ほど、変換したSAMファイルをBAMファイルに変換します。

# SAM ファイルを読み込み
with pysam.AlignmentFile("output.sam", "r") as samfile:
    # 同じヘッダー情報を使って SAM ファイルをテキストモードで書き出す
    with pysam.AlignmentFile("output.bam", "wb", header=samfile.header) as bamfile:
        # ファイル内の全リードを反復して書き出す
        for read in samfile.fetch(until_eof=True):
            bamfile.write(read)

うまくいけば、output.bamファイルが作成されます。

インデックスの作成

BAMファイルのインデックスを作成すると、より効率的にファイルを扱うことができます。

bam_path = "input.bam"
pysam.index(bam_path)

コマンドを実行すると、input.bam.baiというインデックスファイルが作成されます。

インデックスがあれば、下記のように特定の領域を効率よく取り出すことができます。

samfile = pysam.AlignmentFile("input.bam", "rb")
for read in samfile.fetch('chr1', 10000, 10200):
    print(read)

samfile.close()

まとめ

本記事では、pysamの基本的な使い方から、特定領域のデータ抽出、インデックス作成、そしてファイルの書き出しまで、バイオインフォマティクスにおけるSAM/BAMファイルの操作に役立つ機能を解説しました。pysamはシンプルながらも強力なツールであり、独自の解析パイプライン構築に大いに役立ちます。ぜひ、実際のデータを用いて試してみてください。