alterlab-pysam

Pysam

概要

Pysam は、ゲノムデータセットの読み取り、操作、書き込みを行うための Python モジュールです。htslib の Pythonic なインターフェースを使用して、SAM/BAM/CRAM アライメントファイル、VCF/BCF バリアントファイル、FASTA/FASTQ シーケンスの読み書きができます。tabix でインデックス化されたファイルのクエリ、カバレッジのパイルアップ解析、samtools/bcftools コマンドの実行も可能です。

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

次のような場合に、このスキルを使用してください。

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

対象外

シナリオ	代わりに使用
バリアントコーリングパイプライン全体を実行する（アライメント -> GATK/DeepVariant -> アノテーション）	alterlab-nf-core-sarek
多数のサンプルのバリアントを、クエリ可能な配列として保存・検索する	alterlab-tiledbvcf
カバレッジトラック、bigWig の生成、deepTools 形式の BAM サマリー	alterlab-deeptools
シーケンスレコードの解析、翻訳、アライメントオブジェクト（BAM/VCF は対象外）	alterlab-biopython
RNA-seq FASTQ からのトランスクリプト単位の定量	alterlab-rnaseq-quant

クイックスタート

インストール

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 に書き込む

CRAM：pysam 0.24 / htslib 1.22 以降の 2 つの動作変更

これらの変更は気づきにくいため、データに問題があると判断する前に確認してください。

1. デフォルトの出力 CRAM バージョンは 3.0 ではなく 3.1 になりました。 "wc" で書き込んだファイルは、古い samtools や、古い htslib に固定された後続ツールでは読み取れない場合があります。利用側の環境を管理できない場合は、明示的に 3.0 を指定してください。

   out = pysam.AlignmentFile("out.cram", "wc", header=src.header,
                             reference_filename="ref.fa",
                             format_options=["version=3.0"])
   

   format_options は str のリストを受け取ります（以前の pysam では bytes も受け付けていました）。

2. CRAM の参照配列は、EBI から自動取得されなくなりました。 CRAM では参照配列に対する相対形式でリードを保存するため、一致する FASTA を指定せずに読み込むと、暗黙にダウンロードされるのではなく失敗します。開くときに reference_filename= を渡すか、REF_PATH / REF_CACHE 環境変数にローカルキャッシュを設定してください。

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

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

よくある落とし穴

1. 座標の混同： コンテキストによって 0 始まりと 1 始まりの方式が異なることに注意してください
2. インデックスの欠落： 多くの操作にはインデックスファイルが必要です。先に作成してください
3. 部分的なオーバーラップ： fetch() は領域内に完全に収まるリードだけではなく、領域の境界と重なるリードも返します
4. 範囲指定のない pileup： pileup(chrom, start, stop) は、重なるリードがまたがるすべての位置の列を生成し、start..stop の範囲に限られません。要求した範囲に限定するには truncate=True を渡してください
5. 塩基品質による意図しない除外： count_coverage() のデフォルトは quality_threshold=15 のため、品質の低い塩基は除外されます。すべてをカウントするには quality_threshold=0 を設定してください
6. イテレーターのスコープ： 「PileupProxy accessed after iterator finished」エラーを避けるため、pileup イテレーターへの参照を保持してください
7. 品質スコアの編集： query_sequence を変更した後に query_qualities をその場で変更することはできません。先にコピーを作成してください
8. ストリームの制限： ストリーミングでサポートされるのは stdin/stdout のみで、任意の Python ファイルオブジェクトには対応していません
9. スレッドセーフ性： I/O 中は GIL が解放されますが、包括的なスレッドセーフ性は完全には検証されていません
10. 参照配列なしの CRAM： 上記の CRAM に関する注意を参照してください。現在は reference_filename= または REF_PATH が必要です

コマンドラインツール

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 操作、pileup 解析、アラインメントの書き込みを含みます

• 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/