samtools-bam-processing

samtools — SAM/BAM/CRAM アラインメントツールキット

概要

samtools は、SAM、BAM、CRAM 形式のシーケンスアラインメントファイルを処理するための標準的なコマンドラインツールキットです。形式変換、座標順ソート、インデックス作成、品質管理統計、リードのフィルタリング、重複のマーキング、複数ファイルのマージなど、アラインメントファイルのライフサイクル全体に対応します。samtools は、アラインメント（STAR、BWA）から下流解析（バリアントコール、ピークコール）までの NGS パイプラインで、ほぼ不可欠な構成要素です。

使用する場面

• アラインメント後に BAM ファイルを座標順にソートする場合（インデックス作成前に必須）
• ソート済み BAM ファイルにインデックスを付け、ランダムアクセスや領域クエリを可能にする場合
• ストレージを節約するために SAM、BAM、CRAM 形式を相互変換する場合
• アラインメントの QC 指標（マッピング率、インサートサイズ、染色体ごとの統計など）を生成する場合
• マッピング品質、FLAG ビット、またはゲノム領域に基づいてリードをフィルタリングする場合
• バリアントコール前に PCR 重複をマーキングまたは除去する場合
• 異なるレーンやサンプルの複数の BAM ファイルをマージする場合
• ターゲット領域の塩基ごとの深度またはカバレッジ率を計算する場合
• カスタムスクリプトで Python ネイティブの BAM 操作が必要な場合は、代わりに pysam を使用してください
• 正規化済み bigWig カバレッジトラックが必要な場合は、代わりに deeptools bamCoverage を使用してください
• 全ゲノムの塩基ごとの深度を計算する場合は、代わりに mosdepth を使用してください（より高速で並列化されています）

前提条件

• インストール：samtools 1.17 以降を推奨
• 入力要件：SAM/BAM/CRAM ファイル。CRAM には FASTA リファレンスが必要
• 関連ツール：FASTA のインデックス作成には samtools faidx、samtools index の前には samtools sort

インストール前に確認してください：現在の環境（たとえば pixi / conda 環境内）で、すでにツールを利用できる場合があります。まず command -v samtools を実行し、パスが返された場合は、以下のインストールコマンドをスキップしてください。pixi プロジェクト内で実行する場合は、samtools を直接実行するのではなく、pixi run samtools を使ってください。

# Bioconda (recommended)
conda install -c bioconda samtools

# Homebrew (macOS)
brew install samtools

# Verify
samtools --version | head -1

クイックスタート

# Typical post-alignment workflow: sort → index → QC
samtools sort -@ 8 -o sorted.bam input.bam
samtools index sorted.bam
samtools flagstat sorted.bam

主要 API

モジュール 1：BAM/SAM の入出力と形式変換

SAM/BAM/CRAM 形式を相互変換し、サブセットを抽出します。

# SAM → BAM (saves ~75% disk space)
samtools view -b -h input.sam -o output.bam

# BAM → CRAM (saves additional 40-50%)
samtools view -C -T reference.fa input.bam -o output.cram

# Filter: mapping quality ≥20, exclude unmapped (-F 4)
samtools view -q 20 -F 4 input.bam -o filtered.bam

# Extract specific region (requires index)
samtools view -h sorted.bam "chr1:1000000-2000000" -o region.bam

# Count reads matching filter
samtools view -c -F 4 input.bam
# Output: 45231923 (number of mapped reads)

# Extract reads as FASTQ (for realignment or de novo assembly)
samtools fastq -@ 4 -1 R1.fastq.gz -2 R2.fastq.gz -0 unpaired.fastq.gz input.bam

# Extract reads as FASTA
samtools fasta input.bam > reads.fasta

# Filter by read group
samtools view -r SAMPLE_001 multi_rg.bam -o sample001.bam

モジュール 2：ソートとインデックス作成

効率的なランダムアクセスができるように BAM ファイルを整理します。

# Sort by coordinate (required before indexing)
samtools sort -@ 8 -m 2G input.bam -o sorted.bam

# Sort by read name (required for fixmate/markdup)
samtools sort -n -@ 8 input.bam -o namesorted.bam

# Index sorted BAM (creates sorted.bam.bai)
samtools index sorted.bam

# For chromosomes > 512 Mbp: use CSI index instead
samtools index -c sorted.bam

# Group reads by name (fast, for fixmate — no full sort needed)
samtools collate -o collated.bam input.bam

モジュール 3：品質管理と統計

アラインメントの QC 指標とカバレッジレポートを生成します。

# Quick summary: total, mapped, paired, properly paired
samtools flagstat sorted.bam
# Example output:
# 50000000 + 0 in total (QC-passed reads + QC-failed reads)
# 48523111 + 0 mapped (97.05% : N/A)
# 50000000 + 0 paired in sequencing
# 48490234 + 0 properly paired (96.98% : N/A)

# Per-chromosome mapped/unmapped read counts
samtools idxstats sorted.bam
# chr1  248956422  12345678  0
# chr2  242193529  11234567  0

# Comprehensive stats (insert sizes, GC content, base quality)
samtools stats -r reference.fa sorted.bam > full_stats.txt
grep "^SN" full_stats.txt | cut -f2,3  # Summary Numbers only

# Coverage report (min/max/mean per region/chromosome)
samtools coverage sorted.bam

# Per-base read depth for specific regions
samtools depth -b target_regions.bed sorted.bam > depth.txt
# Output: chr  pos  depth (e.g., chr1  1000  45)

# Statistics split by read group
samtools stats -S RG sorted.bam > per_rg_stats.txt

モジュール 4：リードのフィルタリングと FLAG 操作

SAM FLAG ビットを使って、特定のサブセットにリードを絞り込みます。

# FLAG reference — common masks:
# 1    = paired         4  = unmapped
# 2    = proper pair    8  = mate unmapped
# 16   = reverse strand 64 = R1 (first in pair)
# 128  = R2             256= secondary alignment
# 1024 = PCR duplicate  2048= supplementary

# Extract properly paired, mapped reads (FLAG 2 set, 4 unset)
samtools view -f 2 -F 4 sorted.bam -o proper_pairs.bam

# Extract R1 reads only
samtools view -f 64 sorted.bam -o R1.bam

# Remove secondary and supplementary alignments
samtools view -F 2304 sorted.bam -o primary.bam

# Extract reads from BED file regions
samtools view -L regions.bed -b sorted.bam -o regions.bam

モジュール 5：重複の処理

バリアントコール前に PCR 重複をマーキングまたは除去します。

# Full duplicate marking workflow (collate → fixmate → sort → markdup)
samtools collate -@ 8 -o collated.bam input.bam
samtools fixmate -m -@ 8 collated.bam fixmated.bam
samtools sort -@ 8 -o sorted.bam fixmated.bam
samtools markdup -@ 8 sorted.bam marked.bam
samtools index marked.bam

# Check duplication rate
samtools flagstat marked.bam | grep "duplicates"
# Output: 2345678 + 0 duplicates (4.83%)

# NovaSeq optical duplicate detection (2500 pixel distance)
samtools markdup -d 2500 sorted.bam marked_novaseq.bam

# Remove duplicates instead of marking
samtools markdup -r sorted.bam deduped.bam

# Get duplication stats without writing output
samtools markdup -s sorted.bam /dev/null

モジュール 6：複数ファイル操作と領域解析

BAM ファイルをマージし、領域レベルの解析を実行します。

# Merge multiple BAM files (all must be sorted)
samtools merge -@ 8 merged.bam lane1.bam lane2.bam lane3.bam

# Merge files listed in a text file (one per line)
samtools merge -b bam_list.txt -@ 8 merged.bam

# Merge with read group tags from filenames
samtools merge -r merged.bam sample1.bam sample2.bam

# Extract specific chromosome region from merged output
samtools view -h merged.bam chr1 -b -o chr1.bam

主な概念

SAM FLAG ビット

FLAG は、ビット値の合計としてリードの属性を表します。よく使われるフィルターパターンは次のとおりです。

一般的なフィルター	-f（必須）	-F（除外）	選択対象
マッピング済みリード	—	4	アラインメント済みのすべてのリード
適切なペア	2	—	両方がマッピングされた適切なペア
ユニークなプライマリー	—	2308	セカンダリー／補助／重複なし
R1 のみ	64	—	ペアの 1 番目のリード
未マッピング	4	—	アラインメントに失敗したリード

CRAM と BAM、SAM の比較

形式	サイズ	速度	必要なもの
SAM	BAM の約 10 倍	I/O が遅い	なし
BAM	1 倍	高速	ランダムアクセス用の .bai インデックス
CRAM	BAM の約 0.6 倍	やや遅い	リファレンス FASTA とインデックス

長期保存には CRAM、進行中の解析には BAM を使用してください。

よくあるワークフロー

ワークフロー 1：アラインメント後の QC と準備

目的：アライナーの出力を QC 指標付きの解析対応 BAM に変換します。

#!/bin/bash
SAMPLE="sample_001"
REF="reference.fa"
THREADS=8

# 1. Sort and index (aligner often outputs unsorted SAM/BAM)
samtools sort -@ $THREADS -o ${SAMPLE}.sorted.bam ${SAMPLE}.bam
samtools index ${SAMPLE}.sorted.bam

# 2. QC metrics
samtools flagstat ${SAMPLE}.sorted.bam > ${SAMPLE}.flagstat.txt
samtools stats -r $REF ${SAMPLE}.sorted.bam > ${SAMPLE}.stats.txt
samtools coverage ${SAMPLE}.sorted.bam > ${SAMPLE}.coverage.txt

# 3. Per-chromosome stats
samtools idxstats ${SAMPLE}.sorted.bam > ${SAMPLE}.idxstats.txt

echo "QC complete: $(grep 'mapped (' ${SAMPLE}.flagstat.txt | head -1)"

ワークフロー 2：重複マーキングの全工程

目的：重複除去済みの入力を必要とする GATK などのバリアントコーラー向けに BAM を準備します。

#!/bin/bash
INPUT="aligned.bam"
FINAL="deduped.bam"
THREADS=8

# Collate → fixmate → sort → markdup
samtools collate -@ $THREADS -o collated.bam $INPUT
samtools fixmate -m -@ $THREADS collated.bam fixmated.bam
samtools sort -@ $THREADS -o sorted.bam fixmated.bam
samtools markdup -@ $THREADS -s sorted.bam $FINAL

# Clean up intermediates
rm collated.bam fixmated.bam sorted.bam

# Index and verify
samtools index $FINAL
samtools flagstat $FINAL | grep "duplic"
# Expected: 3-15% duplicates (WGS); 10-30% for amplicon

主要なパラメーター

パラメーター	コマンド	デフォルト	範囲／オプション	効果
-@	Most	0	1～N コア	圧縮／I/O 用スレッドを追加します
-m	sort	768M	例：2G、4G	ソート用のスレッドあたりのメモリ
-q	view	0	0～60	マッピング品質の最小フィルター
-f	view	0	FLAG ビット	すべてのビットが設定されたリードを含めます
-F	view	0	FLAG ビット	いずれかのビットが設定されたリードを除外します
-b	view	—	フラグ	BAM 形式で出力します
-C	view	—	フラグ	CRAM を出力します（-T が必要です）
-T	view	—	FASTA パス	CRAM 出力用のリファレンス
-d	markdup	0	0～2500	光学的重複のピクセル距離
-r	markdup	—	フラグ	重複を削除します（マークのみではありません）
-n	sort	—	フラグ	位置ではなくリード名でソートします
-c	index	—	フラグ	CSI インデックスを作成します（chr が 512 Mb を超える場合に必要です）

ベストプラクティス

1. インデックス作成の前に必ずソートします：samtools index には座標順にソートされた入力が必要です。ソートされていない BAM のインデックス作成を試みると、失敗するか、誤った結果が生成されます。

2. すべての本番実行で -@ を使用します：ほとんどの samtools コマンドは I/O がボトルネックになります。-@ 8 を追加すると、オーバーヘッドを最小限に抑えながら、圧縮／展開をほぼ線形に高速化できます。

3. 解析の前に flagstat を実行します：samtools flagstat は数秒で完了し、後続の処理に時間を費やす前に、アライメントの失敗（マッピング率の低さ、想定外のペアエンド率）を検出できます。

4. collate → fixmate → sort → markdup のパイプラインを使用します：fixmate を実行せずに座標順の BAM に対して samtools markdup を直接実行すると、重複検出が不正確になります。fixmate -m で追加される mate 情報は不可欠です。

5. アーカイブには CRAM を推奨します：CRAM は品質を損なうことなく、BAM と比べてストレージを 40～50% 削減します。CRAM ファイルと一緒にリファレンス FASTA を必ず保存してください。

6. 対象を絞った解析には -L bed_file を使用します：BED で定義した対象領域（WES キャプチャ、アンプリコン）に samtools view を限定すると、後続処理の I/O を大幅に削減できます。

よく使うレシピ

レシピ：複数サンプルの Flagstat を一括実行する

# Process all BAM files in directory
for bam in *.sorted.bam; do
    echo "=== $bam ==="
    samtools flagstat $bam | grep -E "mapped|properly paired|duplicates"
done

レシピ：De Novo アセンブリ用に未マッピングリードを抽出する

# Pull both unmapped reads (useful for pathogen detection)
samtools view -f 4 -b input.bam -o unmapped.bam
samtools fastq -@ 4 -1 unmapped_R1.fastq -2 unmapped_R2.fastq unmapped.bam
echo "Unmapped pairs ready for de novo assembly"

レシピ：BAM を目標カバレッジにダウンサンプリングする

# Estimate current depth, then subsample to ~30×
TOTAL=$(samtools flagstat input.bam | grep "mapped (" | head -1 | awk '{print $1}')
GENOME_SIZE=3100000000  # hg38
READ_LEN=150
CURRENT_COV=$(echo "scale=1; $TOTAL * $READ_LEN / $GENOME_SIZE" | bc)
TARGET_FRAC=$(echo "scale=3; 30 / $CURRENT_COV" | bc)
echo "Current: ${CURRENT_COV}×; subsample fraction: $TARGET_FRAC"
samtools view -b -s $TARGET_FRAC input.bam -o downsampled.bam
samtools index downsampled.bam

トラブルシューティング

問題	原因	解決策
[bam_index_build2] fail to index	BAM が座標順にソートされていません	先にソートします：samtools sort -o sorted.bam input.bam
BAI index too large for chromosome	染色体が 512 Mbp を超えています	CSI インデックスを使用します：samtools index -c input.bam
CRAM: reference not found	リファレンス FASTA がないか、正しくありません	REF_PATH 環境変数を設定するか、-T ref.fa を使用します
重複のマーキングが不正確です	fixmate の手順が省略されています	パイプライン全体を実行します：collate → fixmate → sort → markdup
flagstat で properly paired が 0% と表示されます	ペアエンド BAM に mate 情報がありません	samtools fixmate を実行して mate 座標を設定します
ソートが非常に遅いです	スレッドあたりのメモリが少なすぎます	-m 4G に増やします。メモリが制限されている場合は -@ を減らします
領域クエリで何も返されません	BAM にインデックスがないか、座標が誤っています	samtools index を実行します。1 始まりの座標を使用します：chr1:1000-2000
[E::hts_open_format] fail to open	ファイルパスが誤っているか、BAM が破損しています	パスを確認し、samtools quickcheck file.bam でテストします

関連スキル

• deeptools-ngs-analysis — samtools の後続処理として、正規化された bigWig カバレッジトラックと ChIP-seq の可視化を行います
• pysam-genomic-files — カスタムスクリプトで BAM を操作するための Python API
• bedtools-genomic-intervals — samtools が生成した BAM／BED ファイルに対するゲノム区間操作

参考資料

• samtools のドキュメント — 公式マニュアルページとコマンドリファレンス
• GitHub: samtools/samtools — ソース、リリース、課題トラッカー
• Danecek et al. (2021)「Twelve years of SAMtools and BCFtools」— GigaScience 10(2)
• SAM 形式仕様 — FLAG ビット、CIGAR 文字列、オプションタグ