pysam

Pysam

概要

Pysam は、ゲノムデータセットの読み取り、操作、書き込みを行うための Python モジュールです。htslib に対する Pythonic なインターフェースで、SAM/BAM/CRAM アライメントファイル、VCF/BCF バリアントファイル、FASTA/FASTQ 配列の読み書きができます。tabix インデックス付きファイルへのクエリ、カバレッジ算出のためのピルアップ解析、samtools/bcftools コマンドの実行も可能です。

このスキルを使用する場合

このスキルは次の場合に使用してください:

• シーケンシングのアライメントファイル (BAM/CRAM) を扱うとき
• 遺伝的バリアント (VCF/BCF) を解析するとき
• 参照配列や遺伝子領域を抽出するとき
• 生のシーケンシングデータ (FASTQ) を処理するとき
• カバレッジやリードデプスを算出するとき
• バイオインフォマティクス解析パイプラインを実装するとき
• シーケンシングデータの品質管理を行うとき
• バリアントコーリングとアノテーションのワークフローを実行するとき

クイックスタート

インストール

uv pip install pysam

基本的な例

アライメントファイルの読み取り:

import pysam

# Open BAM file and fetch reads in region
samfile = pysam.AlignmentFile("example.bam", "rb")
for read in samfile.fetch("chr1", 1000, 2000):
    print(f"{read.query_name}: {read.reference_start}")
samfile.close()

バリアントファイルの読み取り:

# Open VCF file and iterate variants
vcf = pysam.VariantFile("variants.vcf")
for variant in vcf:
    print(f"{variant.chrom}:{variant.pos} {variant.ref}>{variant.alts}")
vcf.close()

参照配列のクエリ:

# Open FASTA and extract sequence
fasta = pysam.FastaFile("reference.fasta")
sequence = fasta.fetch("chr1", 1000, 2000)
print(sequence)
fasta.close()

主要機能

1. アライメントファイル操作 (SAM/BAM/CRAM)

AlignmentFile クラスを使用して、アライメント済みのシーケンシングリードを扱います。マッピング結果の解析、カバレッジの算出、リードの抽出、品質管理に適しています。

一般的な操作:

• BAM/SAM/CRAM ファイルの オープンと読み取り
• 特定のゲノム領域からのリードの取得
• マッピング品質、フラグ、その他の条件によるリードのフィルタリング
• フィルタリングまたは修正したアライメントの書き込み
• カバレッジ統計の算出
• ピルアップ解析 (塩基ごとのカバレッジ) の実行
• リード配列、品質スコア、アライメント情報へのアクセス

リファレンス: 以下の詳細なドキュメントは references/alignment_files.md を参照してください:

• アライメントファイルのオープンと読み取り
• AlignedSegment の属性とメソッド
• fetch() による領域指定の取得
• カバレッジのためのピルアップ解析
• BAM ファイルの書き込みと作成
• 座標系とインデックス作成
• パフォーマンス最適化のヒント

2. バリアントファイル操作 (VCF/BCF)

VariantFile クラスを使用して、バリアントコーリングのパイプラインから得られた遺伝的バリアントを扱います。バリアント解析、フィルタリング、アノテーション、集団遺伝学に適しています。

一般的な操作:

• VCF/BCF ファイルの読み書き
• 特定の領域内のバリアントのクエリ
• バリアント情報 (位置、アレル、品質) へのアクセス
• サンプルのジェノタイプデータの抽出
• 品質、アレル頻度、その他の条件によるバリアントのフィルタリング
• 追加情報によるバリアントのアノテーション
• サンプルまたは領域のサブセット化

リファレンス: 以下の詳細なドキュメントは references/variant_files.md を参照してください:

• バリアントファイルのオープンと読み取り
• VariantRecord の属性とメソッド
• INFO および FORMAT フィールドへのアクセス
• ジェノタイプとサンプルの扱い
• VCF ファイルの作成と書き込み
• バリアントのフィルタリングとサブセット化
• マルチサンプル VCF の操作

3. 配列ファイル操作 (FASTA/FASTQ)

参照配列へのランダムアクセスには FastaFile、生のシーケンシングデータの読み取りには FastxFile を使用します。遺伝子配列の抽出、参照配列に対するバリアントの検証、生のリードの処理に適しています。

一般的な操作:

• ゲノム座標による参照配列のクエリ
• 対象の遺伝子や領域の配列の抽出
• 品質スコア付き FASTQ ファイルの読み取り
• バリアントの参照アレルの検証
• 配列統計の算出
• 品質または長さによるリードのフィルタリング
• FASTA と FASTQ 形式の相互変換

リファレンス: 以下の詳細なドキュメントは references/sequence_files.md を参照してください:

• FASTA ファイルへのアクセスとインデックス作成
• 領域指定による配列の抽出
• 遺伝子に対する逆相補鎖の扱い
• FASTQ ファイルの逐次読み取り
• 品質スコアの変換とフィルタリング
• tabix インデックス付きファイル (BED, GTF, GFF) の扱い
• 一般的な配列処理のパターン

4. 統合バイオインフォマティクスワークフロー

Pysam は、複数のファイル形式を統合して包括的なゲノム解析を行うことに優れています。一般的なワークフローでは、アラインメントファイル、バリアントファイル、リファレンス配列を組み合わせます。

一般的なワークフロー:

• 特定の領域のカバレッジ統計を計算する
• アラインリードに対してバリアントを検証する
• カバレッジ情報でバリアントに注釈を付ける
• バリアント位置周辺の配列を抽出する
• 複数の条件に基づいてアラインメントまたはバリアントをフィルタリングする
• 可視化用のカバレッジトラックを生成する
• 複数のデータ型にまたがる品質管理を行う

リファレンス: 以下の詳細な例については references/common_workflows.md を参照してください:

• 品質管理ワークフロー(BAM の統計、リファレンスの整合性)
• カバレッジ解析(塩基ごとのカバレッジ、低カバレッジの検出)
• バリアント解析(注釈、リードのサポート数によるフィルタリング)
• 配列抽出(バリアントのコンテキスト、遺伝子配列)
• リードのフィルタリングとサブセット化
• 統合パターン(BAM+VCF、VCF+BED など)
• 複雑なワークフローのパフォーマンス最適化

主要な概念

座標系

重要: Pysam は 0 始まり・半開区間 の座標を使用します(Python の慣例):

• 開始位置は 0 始まりです(最初の塩基は位置 0)
• 終了位置は排他的です(範囲に含まれません)
• 領域 1000-2000 は塩基 1000～1999 を含みます(合計 1000 塩基)

例外: fetch() の領域文字列は samtools の慣例(1 始まり)に従います:

samfile.fetch("chr1", 999, 2000)      # 0-based: positions 999-1999
samfile.fetch("chr1:1000-2000")       # 1-based string: positions 1000-2000

VCF ファイル: ファイル形式では 1 始まりの座標を使用しますが、VariantRecord.start は 0 始まりです。

インデックスの要件

特定のゲノム領域へのランダムアクセスにはインデックスファイルが必要です:

• BAM ファイル: .bai インデックスが必要です(pysam.index() で作成)
• CRAM ファイル: .crai インデックスが必要です
• FASTA ファイル: .fai インデックスが必要です(pysam.faidx() で作成)
• VCF.gz ファイル: .tbi tabix インデックスが必要です(pysam.tabix_index() で作成)
• BCF ファイル: .csi インデックスが必要です

インデックスがない場合は、順次読み取りのために fetch(until_eof=True) を使用してください。

ファイルモード

ファイルを開く際は形式を指定します:

• "rb" - BAM の読み取り(バイナリ)
• "r" - SAM の読み取り(テキスト)
• "rc" - CRAM の読み取り
• "wb" - BAM の書き込み
• "w" - SAM の書き込み
• "wc" - CRAM の書き込み

パフォーマンスに関する考慮事項

1. ランダムアクセス操作には、必ずインデックス付きファイルを使用してください
2. 繰り返しの fetch 操作の代わりに、列単位の解析には pileup() を使用してください
3. 手動で反復してカウントするのではなく、カウントには count() を使用してください
4. 独立したゲノム領域を解析する場合は、領域を並列に処理してください
5. リソースを解放するために、ファイルは明示的にクローズしてください
6. インデックスなしで順次処理する場合は、until_eof=True を使用してください
7. 必要な場合を除き、複数のイテレーターの使用は避けてください(必要な場合は multiple_iterators=True を使用)

よくある落とし穴

1. 座標の混同: 文脈によって 0 始まりと 1 始まりの座標系が異なることに注意してください
2. インデックスの欠如: 多くの操作にはインデックスファイルが必要です。先に作成してください
3. 部分的な重なり: fetch() は、領域内に完全に含まれるリードだけでなく、領域の境界と重なるリードも返します
4. イテレーターのスコープ: "PileupProxy accessed after iterator finished" エラーを避けるため、pileup イテレーターの参照は保持したままにしてください
5. 品質スコアの編集: query_sequence を変更した後に query_qualities をその場で変更することはできません。先にコピーを作成してください
6. ストリームの制限: ストリーミングでサポートされるのは stdin/stdout のみで、任意の Python ファイルオブジェクトはサポートされません
7. スレッドセーフティ: I/O 中は GIL が解放されますが、包括的なスレッドセーフティは十分に検証されていません

コマンドラインツール

Pysam は samtools と bcftools のコマンドへのアクセスを提供します:

# Sort BAM file
pysam.samtools.sort("-o", "sorted.bam", "input.bam")

# Index BAM
pysam.samtools.index("sorted.bam")

# View specific region
pysam.samtools.view("-b", "-o", "region.bam", "input.bam", "chr1:1000-2000")

# BCF tools
pysam.bcftools.view("-O", "z", "-o", "output.vcf.gz", "input.vcf")

エラー処理:

try:
    pysam.samtools.sort("-o", "output.bam", "input.bam")
except pysam.SamtoolsError as e:
    print(f"Error: {e}")

リソース

references/

各主要機能の詳細なドキュメント:

• alignment_files.md - SAM/BAM/CRAM の操作に関する完全なガイドです。AlignmentFile クラス、AlignedSegment の属性、fetch 操作、ピルアップ解析、アライメントの書き出しを扱います

• variant_files.md - VCF/BCF の操作に関する完全なガイドです。VariantFile クラス、VariantRecord の属性、ジェノタイプの扱い、INFO/FORMAT フィールド、マルチサンプルの操作を扱います

• sequence_files.md - FASTA/FASTQ の操作に関する完全なガイドです。FastaFile と FastxFile クラス、配列の抽出、品質スコアの扱い、tabix インデックス付きファイルへのアクセスを扱います

• common_workflows.md - 複数のファイル形式を組み合わせた、統合バイオインフォマティクスワークフローの実践的な例です。品質管理、カバレッジ解析、バリアントの検証、配列の抽出を扱います

ヘルプの参照先

特定の操作の詳細は、該当するリファレンスドキュメントを参照してください:

• BAM ファイルの操作やカバレッジの計算 → alignment_files.md
• バリアントやジェノタイプの解析 → variant_files.md
• 配列の抽出や FASTQ の処理 → sequence_files.md
• 複数のファイル形式を統合する複雑なワークフロー → common_workflows.md

公式ドキュメント: https://pysam.readthedocs.io/