アライメントデータの操作

アライメントデータ（SAM/BAMファイル）を操作するにはsamtoolsを利用する。

samtools      # サブコマンド一覧表示
samtools view # サブコマンド使用法確認

samtools view

samtools view IN     # 入力(IN): SAM/BAM (自動検出); 出力: SAM (aln)
samtools view -S IN  # 入力: SAM (以前の仕様ではSAM入力であることを明示する必要があった)
samtools view -H IN  # 出力: SAM (header)
samtools view -h IN  # 出力: SAM (header+aln)
samtools view -b IN  # 出力: BAM/圧縮(aln)
samtools view -u IN  # 出力: BAM/非圧縮(header+aln)

# 使用例
samtools view -hb in.sam              # SAM to BAM変換
samtools view -u IN | samtools xxxxx  # samtools間パイプ接続 (時短処理)

またFLAGやMAPQの情報をもとにフィルタリングするために利用する。

samtools view -f FLAG IN  # FLAGで指定した値を有つものを抽出
samtools view -F FLAG IN  # FLAGで指定した値を有たないものを抽出
samtools view -q MAPQ IN  # MAPQで指定したMAPQ値以上のものを抽出

samtools flags

FLAG情報を確認するために利用する。

samtools flags 147  # 147(10進数)は10010011(2進数), 第1,2,5,8bitが1(True)である.
                    # 出力結果：0x93 147 PAIRED, PROPER, PAIR, REVERSE, READ2

# 逆変換
samtools flags read1,proper_pair  # 出力結果：0x42 66 PROPER_PAIR,READ1

samtools sort

samtools sort IN             # 入力(IN): SAM/BAM; 出力: CSRT BAM/標準出力
samtools sort -n IN          # 出力形式: NSRT BAM
samtools sort IN -o out.bam  # 出力先: ファイル
samtools sort IN out.prefix  # 以前の使用法

  # 計算資源量指定を追加: -m memoryNumberG -@ thread

samtools index

特定領域にアライメントされたリードに効率的にランダムアクセスできるようにするためには、 samtools indexサブコマンドによりインデックスを作成する必要がある。

samtools index CSRT.bam   # 入力: CSRT.bam; 出力: CSRT.bam.bai

  # 計算資源量指定を追加: -@ thread

インデックス作成済みBAMに対してはsamtools viewを利用して範囲データ操作が実行できる。

samtools view CSRT.bam chr:start-end  # chr:start-end領域にあるアライメントを抽出
samtools view -L target.bed CSRT.bam  # target.bedと重なる領域を抽出

samtools fixmate

mate座標やインサートサイズ情報を埋める。

samtools fixmate NSRT.bam out.bam      # 入力: NSRT.bam; 出力: out.bam
samtools fixmate -r NSRT.bam out.bam   # secondary/unmapped readsを除去

  # 計算資源量指定を追加: -@ thread

samtools merge

BAMファイルを連結する。

samtools merge out.bam in1.bam in2.bam ...
  
  # 計算資源量指定を追加: -@ thread

応用例

insert sizeの一覧取得

samtools view -f 66 $IN | cut -f9 | tr -d '-'

insert sizeによるフィルタリング(例: 150bp未満を除外)

samtools view -h $IN |
awk -F'\t' '
   /^@/
  !/^@/ {if($9 >= 150 || $9 <= -150) print}' 

ENCODE3プロジェクトにおける使用例を一部改変したものを下に示す。

# 変数定義
INF=$1
ID=$(basename $INF | sed 's/\.[^.][^.]*$//')
TMPF1=${ID}.1.bam
TMPF2=${ID}.2.bam
TMPF3=${ID}.3.bam
OUTF=${ID}_FILT.bam
# ---
JARPATH=~/.conda/envs/aln/share/picard-2.27.4-0/picard.jar
CMD="java -jar $JARPATH MarkDuplicates"

# 除去: unmapped, mate unmapped, not primary alignment, reads failing platform
# 除去: low MAPQ reads
# キープ: properly paired reads
samtools view -F 1804 -f 2 -q 30 -u $INF |

# --- この範囲の処理の意義は不明瞭 ------------------------
# 除去: orphan reads (pair was removed),
#       read pairs mapping to different chromosomes
samtools sort -n -                       |
samtools fixmate -r - $TMPF1
samtools view -F 1804 -f 2 -u $TMPF1     |
# ---------------------------------------------------------
samtools sort - > $TMPF2

# 検出: duplicates
$CMD                              \
  INPUT=$TMPF2                    \
  OUTPUT=$TMPF3                   \
  METRICS_FILE=${ID}_dupl         \
  VALIDATION_STRINGENCY=LENIENT   \
  ASSUME_SORTED=true              \
  REMOVE_DUPLICATES=false

# 除去: duplicates
samtools view -F 1804 -f 2 -hb $TMPF3 > $OUTF

rm $TMPF1 $TMPF2 $TMPF3

不明点

• samtools sortでCSRTしたらファイルサイズが小さくなった？