pysam

pysam

概要

HTSlib がサポートするゲノム形式に低レベルでストリーミングアクセスするには、pysam を使用します。

• SAM/BAM/CRAM には AlignmentFile と AlignedSegment
• VCF/BCF には VariantFile、VariantHeader、VariantRecord
• インデックス付き FASTA には FastaFile、順次読み込む FASTA/FASTQ には FastxFile
• BGZF 圧縮され、tabix インデックスが付与された BED/GFF/GTF/カスタムテーブルには TabixFile
• ラップされたコマンドディスパッチャーには pysam.samtools と pysam.bcftools

現在の上流ベースライン：pysam 0.24.0（2026 年 4 月 27 日）。HTSlib/samtools/bcftools 1.23.1 をラップしています。バージョン固有のガイダンスを更新する前に、references/sources.md をお読みください。

インストール

再現性のある作業には、固定されたリリースを使用してください。

uv pip install "pysam==0.24.0"

ランタイムを確認します。

import pysam

print(pysam.__version__)           # 0.24.0
print(pysam.__samtools_version__)  # 1.23.1

サポート対象の macOS および Linux プラットフォーム向けに、ビルド済みホイールを利用できます。ソースからビルドするには、C コンパイラーと HTSlib のビルド依存関係が必要です。references/sources.md からリンクされている公式インストールガイドをお読みください。

まず決めること

コードを書く前に、以下を確認してください。

1. 実際の形式、圧縮方式、ソート順、および利用可能なインデックスを特定します。
2. 座標として Python の数値座標を使うか、領域文字列を使うかを決めます。両者を混在させないでください。
3. CRAM の場合は、正確な参照アセンブリと FASTA を特定します。
4. インデックス付き領域アクセスを優先し、意図がある場合に限って順次イテレーションを使用します。
5. 書き込み時はヘッダーを保持し、デフォルトでは新しいパスに書き込みます。
6. フィルタリングの意味を明示します。マッピング品質／塩基品質、フラグ、オーバーラップの処理、重複の処理、pileup の深度上限を指定してください。

不慣れなファイルの場合は、同梱の読み取り専用インスペクターから始めてください。

python scripts/inspect_hts.py sample.bam
python scripts/inspect_hts.py cohort.vcf.gz
python scripts/inspect_hts.py reference.fa

同梱スクリプト

スクリプト	用途	一般的な呼び出し方
scripts/inspect_hts.py	アライメント、バリアント、FASTA、FASTQ、tabix ファイルのメタデータのみを検査	python scripts/inspect_hts.py sample.cram --reference ref.fa
scripts/alignment_qc.py	リード／QC カウントをストリーミング集計し、JSON として出力	python scripts/alignment_qc.py sample.bam --max-records 100000
scripts/variant_summary.py	バリアント、FILTER、遺伝子型の概要をストリーミング集計し、JSON として出力	python scripts/variant_summary.py cohort.vcf.gz --region chr1:1-1000000
scripts/filter_alignments.py	レコード順を変更せずに SAM/BAM/CRAM をフィルタリング	python scripts/filter_alignments.py input.bam output.bam --exclude-secondary

すべてのスクリプトは、既存の出力を上書きしません。座標、インデックス、プライバシーに関する注意事項を確認するには、それぞれを --help 付きで実行してください。

座標の仕様

pysam API が受け取る数値座標は、0 起点で終端を含まない形式です。 数値を指定する AlignmentFile.fetch()、VariantFile.fetch()、FastaFile.fetch()、TabixFile.fetch()、および pileup() の引数がこれに該当します。

領域文字列は samtools 形式で、1 起点かつ終端を含みます。

# The same 100 bases:
bam.fetch("chr1", 99, 199)          # [99, 199)
bam.fetch(region="chr1:100-199")    # 1-based inclusive

VCF テキストでは POS は 1 起点ですが、レコードのプロパティでは両方の形式が公開されています。

record.pos    # 1-based
record.start  # 0-based inclusive
record.stop   # 0-based exclusive

形式間の変換、オーバーラップの意味、インデックスの選択、コンティグ名の確認については、references/coordinates_and_indexing.md をお読みください。

アライメントファイル

コンテキストマネージャーを使い、モードを明示してください。

import pysam

with pysam.AlignmentFile("sample.bam", "rb", threads=4) as bam:
    for read in bam.fetch("chr1", 1_000, 2_000):
        if (
            not read.is_unmapped
            and not read.is_secondary
            and not read.is_supplementary
            and read.mapping_quality >= 30
        ):
            print(read.query_name, read.reference_start, read.cigarstring)

インデックスを必要とせず、配置されていない未マッピングリードを含むすべてのレコードをファイル順にストリーミングするには、fetch(until_eof=True) を使用します。

with pysam.AlignmentFile("sample.bam", "rb") as bam:
    for read in bam.fetch(until_eof=True):
        ...

重要な違いは以下のとおりです。

• fetch() は領域とオーバーラップするアライメントレコードを返します。
• count() はレコード数を数え、デフォルトは read_callback="nofilter" です。
• count_coverage() は A/C/G/T 塩基数を返し、デフォルトでは塩基品質 15 と read_callback="all" を使用します。
• pileup() は列ごとのリードを公開し、独自のフィルタリング、塩基品質、オーバーラップ、孤立リード、および max_depth=8000 のデフォルト値を持ちます。

正確な領域の pileup では、truncate=True と明示的なフィルターを設定してください。

with pysam.FastaFile("reference.fa") as fasta, pysam.AlignmentFile(
    "sample.bam", "rb"
) as bam:
    for column in bam.pileup(
        "chr1",
        1_000,
        2_000,
        truncate=True,
        stepper="samtools",
        fastafile=fasta,
        min_mapping_quality=20,
        min_base_quality=20,
        max_depth=100_000,
    ):
        print(column.reference_pos, column.get_num_aligned())

フラグ、CIGAR 操作、タグ、修飾塩基、レコードの書き込み、pileup の詳細、イテレーターの有効期間については、references/alignment_files.md をお読みください。

バリアントファイル

入力形式は自動検出されます。数値による取得座標は引き続き 0 始まりです。

import pysam

with pysam.VariantFile("cohort.vcf.gz", threads=4) as variants:
    for record in variants.fetch("chr1", 999_999, 2_000_000):
        print(record.contig, record.pos, record.ref, record.alts)
        for sample_name, call in record.samples.items():
            print(sample_name, call.get("GT"))

レコードを取得する前にサンプルを絞り込みます。

with pysam.VariantFile("cohort.bcf") as variants:
    variants.subset_samples(["sample_A", "sample_B"])
    for record in variants:
        ...

ヘッダーを変更する場合は、各レコードをコピーし、新たに宣言した INFO／FORMAT／FILTER フィールドを割り当てる前に、コピーしたレコードを移行先のヘッダーに変換してください。header.samples を手動で消去して再構築しないでください。

安全なヘッダー、書き込み、サンプルの絞り込み、欠損ジェノタイプ、シンボリックアレル、フィルタリング、変換、インデックス作成については、references/variant_files.md を参照してください。

FASTA、FASTQ、Tabix

インデックス付き FASTA では、数値座標は 0 始まりです。

with pysam.FastaFile("reference.fa") as fasta:
    sequence = fasta.fetch("chr1", 999, 1_099)

FastxFile は逐次処理されます。persist=False は高速ですが、反復処理が進むと取得したレコードは無効になります。

with pysam.FastxFile("reads.fastq.gz", persist=False) as reads:
    for read in reads:
        qualities = read.get_quality_array()
        ...

Tabix の入力は座標順に並べ、通常の gzip ではなく BGZF で圧縮する必要があります。元データを変更しない次の 2 段階のワークフローを使用します。

pysam.tabix_compress("regions.bed", "regions.bed.gz")
pysam.tabix_index("regions.bed.gz", preset="bed")

with pysam.TabixFile("regions.bed.gz", parser=pysam.asBed()) as tbx:
    for interval in tbx.fetch("chr1", 1_000, 2_000):
        print(interval.contig, interval.start, interval.end)

FASTA／FASTQ レコードと安全な Tabix インデックスの作成については、references/sequence_files.md を参照してください。

CRAM、リモート I/O、スレッド

pysam 0.24 では、継承される HTSlib の動作が変更されました。

• 新規に書き込まれる CRAM のデフォルトは CRAM 3.0 ではなく、CRAM 3.1 です。
• HTSlib はデフォルトで EBI リファレンスサーバーに接続しなくなりました。
• 再現性のあるローカル読み書きには、reference_filename="reference.fa" を推奨します。

with pysam.AlignmentFile(
    "sample.cram",
    "rc",
    reference_filename="reference.fa",
    threads=4,
) as cram:
    for read in cram.fetch("chr1", 1_000, 2_000):
        ...

リファレンスの MD5 による検索を意図している場合に限り、REF_PATH／REF_CACHE を設定してください。CRAM が自己完結していると想定しないでください。threads= は圧縮／展開を高速化しますが、Python による解析を並列化するものではありません。

CRAM 変換、リモートアクセス、または並行反復処理を行う前に、references/cram_and_performance.md を参照してください。

samtools と bcftools のラッパー

コマンドモジュールを明示的にインポートしてください。コマンドラインの各トークンは個別の文字列として渡します。

import pysam.samtools
import pysam.bcftools

pysam.samtools.sort(
    "-@", "4", "-o", "sorted.bam", "input.bam", catch_stdout=False
)
pysam.samtools.index("-@", "4", "sorted.bam", catch_stdout=False)

pysam.bcftools.index("--csi", "variants.vcf.gz", catch_stdout=False)

ディスパッチャーはデフォルトで標準出力をキャプチャします。大容量またはバイナリの出力には、出力全体をメモリに返すのではなく、catch_stdout=False を指定してツールの -o オプションを使用するか、save_stdout=... を使用してください。

try:
    pysam.samtools.quickcheck("-v", "sample.bam")
except pysam.SamtoolsError as error:
    messages = pysam.samtools.quickcheck.get_messages()
    raise RuntimeError(messages or str(error)) from error

レコード単位の処理には Python API を使用し、sort、index、merge、view、正規化などの成熟した一括操作にはディスパッチャーを使用してください。信頼できないシェルコマンドを分割して、ディスパッチャーの引数を組み立てないでください。

書き込み規則

• 出力を開く前に、有効なヘッダーをコピーするか作成してください。
• 新しいパスに書き込み、明示的な置き換えが必要な場合を除き、force=True を使用しないでください。
• 出力にインデックスを付ける場合は、ソート順を維持してください。
• query_qualities より先に query_sequence を設定してください。
• pysam.CIGAR_OPS の列挙型メンバーを推奨します。pysam.CMATCH などのトップレベル定数は互換性のための別名であり、将来的に削除される予定です。
• アラインメントの出力は pysam.samtools.quickcheck() で検証し、バリアント／シーケンスの出力は、後続処理で使用する前に開き直してください。
• リファレンスまたは座標が従来のインデックスの上限を超える場合は、BAI／TBI ではなく CSI を使用してください。

リファレンスマップ

必要な情報	参照先
アラインメント API、フラグ、CIGAR、パイルアップ、修飾塩基	references/alignment_files.md
VCF／BCF ヘッダー、レコード、サンプル、書き込み	references/variant_files.md
FASTA／FASTQ と Tabix インデックス付きテーブル	references/sequence_files.md
座標変換とインデックスの選択	references/coordinates_and_indexing.md
CRAM リファレンス、リモート I/O、スレッド、パフォーマンス	references/cram_and_performance.md
正しい統合解析パターン	references/common_workflows.md
現行 API の簡潔なシグネチャとデフォルト値	references/api_reference.md
既存環境向けのアップグレードに関する注意事項	references/migration_to_0_24.md
公式ドキュメント、仕様、リリース情報の参照先	references/sources.md

よくある失敗例

• 数値の VariantFile.fetch() 座標を1始まりとして扱う
• BGZFとtabix/CSIが必要な場面で、通常のgzipを使う
• インデックスを使わずに領域を取得する
• fetch() に未配置の未マップアラインメントが含まれると思い込む
• 正確なパイルアップ区間を取得するための truncate=True を忘れる
• 塩基品質13、深度上限8000など、パイルアップのデフォルト設定を無視する
• アクティブなイテレーターやスレッド間でファイルハンドルを共有する
• 正確な参照配列を使わずにCRAMをデコードする
• 出力ヘッダーに宣言する前に新しいVCFフィールドを割り当てる
• 大量のsamtools/bcftools出力をメモリーに読み込む
• インデルやシンボリックアレルにSNPの塩基カウント手法を使う

関連スキル

• polars-bio：ArrowとDataFusionを基盤に構築された、Polars DataFrameでのゲノム区間操作とバイオインフォマティクスファイルの入出力に対応するPythonライブラリーpolars-bioです。
• tiledbvcf：tiledbvcf Python APIとCLIを使って、TileDB-VCFデータセットにゲノムバリアントデータを保存し、クエリします（作成、保存、エクスポート、一覧表示、統計情報の取得）。
• gtars：ゲノム区間のためのRust/Python/CLIツールキットGtarsを調査し、利用計画を立てます。BED RegionSetの集合演算（reduce、setdiff、intersec…