pysam

Pysam - ゲノムアラインメント

ハイスループットシーケンスのパイプラインで使用します。参照ゲノムにアラインメントされた数十億の DNA 断片に効率的にアクセスできます。

使用する場面

• 次世代シーケンス（NGS）データの処理。
• ゲノム変異（SNP、インデル）の解析。
• 特定のゲノム領域からのリードの抽出。
• カスタムバイオインフォマティクスパイプラインの構築。
• シーケンスデータの品質管理。

基本原則

インデックスアクセス

ゲノム領域への効率的なランダムアクセスには、BAM ファイルのインデックス（.bai）が必要です。

座標系

ゲノム座標では位置に 0 始まり（Python 形式）を使用しますが、範囲には一部のコンテキストで 1 始まりを使用します。

リード属性

各リードには、配列、品質スコア、アラインメント位置、フラグが含まれます。

クイックリファレンス

標準インポート

import pysam

基本パターン

# 1. Open BAM file
samfile = pysam.AlignmentFile("aligned_reads.bam", "rb")

# 2. Iterate over reads in a specific genomic region
for read in samfile.fetch("chr1", 10000, 10100):
    print(f"Read: {read.query_name}, Quality: {read.mapping_quality}")
    print(f"Sequence: {read.query_sequence}")
    print(f"Position: {read.reference_start}")

# 3. Variant analysis (VCF)
vcf = pysam.VariantFile("mutations.vcf")
for rec in vcf.fetch("chr1", 10000, 10100):
    print(f"Pos: {rec.pos}, Ref: {rec.ref}, Alt: {rec.alts}")
    print(f"Genotype: {rec.samples['sample1']['GT']}")

# 4. Writing aligned reads
outfile = pysam.AlignmentFile("output.bam", "wb", template=samfile)
for read in samfile:
    if read.mapping_quality > 30:
        outfile.write(read)
outfile.close()

重要なルール

✅ 推奨事項

• 必ずインデックス付きファイルを使用する - pysam.index("file.bam") でインデックスを作成すると、高速にアクセスできます。
• リードのフラグを確認する - read.is_paired、read.is_unmapped を使ってリードをフィルタリングします。
• マッピングされていないリードを処理する - マッピングされていないリードの reference_start = -1 です。
• ファイルを明示的に閉じる - コンテキストマネージャーまたは .close() を使い、リソースリークを防ぎます。

❌ 禁止事項

• BAM 全体を反復処理しない - 効率を上げるため、領域を指定して fetch() を使用します。
• 品質スコアを無視しない - 低品質の塩基は誤った変異の原因になります。
• 座標系を混在させない - 0 始まりと 1 始まりのインデックスを一貫して使用します。

高度なパターン

遺伝子ごとのリード数のカウント

# Using a gene annotation file
genes = {}  # gene_name -> (chr, start, end)
for read in samfile.fetch():
    # Check if read overlaps any gene
    for gene, (chr, start, end) in genes.items():
        if read.reference_name == chr and start <= read.reference_start < end:
            genes[gene]['count'] += 1

変異のフィルタリング

# Filter variants by quality and depth
for rec in vcf.fetch():
    depth = rec.samples['sample1']['DP']
    qual = rec.qual
    if depth > 10 and qual > 20:
        # Process high-quality variant
        pass

Pysam は、ゲノムデータ処理に必要な低レベルのアクセスを提供し、研究者が生命そのものの生データを直接扱えるようにします。