# 11章 アラインメントデータの操作
- `SAM`(Sequence Alignment/Mapping)と`BAM`(そのバイナリー形式)
- 大量のアラインメントを保存するための最も一般的な形式
## なぜSAM/BAMを学習するのか
1. バイオインフォマティクスの仕事の大部分はアラインメントファイルの操作
2. SAM/BAM処理に使う下記の様なスキルは、他のファイル形式でも応用できる
- バイナリファイルを扱う
- ビットフラグ(後ほど紹介)を使って情報を抽出
- APIを利用

## アラインメント形式の理解
- [githubの11章](https://github.com/vsbuffalo/bds-files/tree/master/chapter-11-alignment)
- [SAMのオリジナルの仕様書](http://samtools.github.io/hts-specs/)
- [SAM形式の原著論文](https://academic.oup.com/bioinformatics/article/25/16/2078/204688)

### SAMヘッダー

@SQ     SN:I    LN:15072434
@SQ     SN:II   LN:15279421
@SQ     SN:III  LN:13783801
@SQ     SN:IV   LN:17493829
@SQ     SN:MtDNA        LN:13794
@SQ     SN:V    LN:20924180
@SQ     SN:X    LN:17718942
@RG     ID:VB00023_L001 SM:celegans-01
@PG     ID:bwa  PN:bwa  VN:0.7.10-r789  CL:bwa mem -R @RG\tID:VB00023_L001\tSM:celegans-01 Caenorhabditis_elegans.WBcel235.dna.toplevel.fa celegans-1.fq celegans-2.fq

- @SQは、参照配列の染色体とその塩基数を示す
- @RGは、リードグループの識別IDとサンプルやシークエンシングプラットフォームのメタデータ
- @PGは、SAM/BAMファイルの作成と処理に使用されたプログラムに関するメタデータが含まれている
- 特に@RGは重要らしい。下流のプログラムでもこのメタデータを使うことが多い。
- `PL`は使ったシークエンサー(ILLUMINAやPACBIOなど)

@RG     ID:SRR035334    PL:ILLUMINA     LB:Solexa-3628  PI:0    DS:SRP000032    SM:NA12891      CN:BI

### samtools
- SAM/BAMファイルはsamtoolsというコマンドラインプログラムを介して処理を行うことが標準的
- 自分は`conda`でインストールしました
- SAMとBAMを自動検出してくれる
- samtoolsはUnixのコマンドラインと組み合わせて使用することが可能
- `samtools view -H celegans.bam | grep "^@RG"`
- 詳しくはあとで紹介される

### SAMアラインメントセクション　

I_2011868_2012306_0:0:0_0:0:0_2489      83      I       2012257 40      50M     =       2011868 -439    CAAAAAATTTTGAAAAAAAAAATTGAATAAAAATTCACGGATTTCTGGCT      22222222222222222222222222222222222222222222222222      NM:i:0  MD:Z:50 AS:i:50 XS:i:50 RG:Z:VB00023_L001       XA:Z:I,-2021713,50M,0;
I_2011868_2012306_0:0:0_0:0:0_2489      163     I       2011868 60      50M     =       2012257 439     GTGGAGACAGCGCCAAAACACCACATGGTGTATTGGCTGCGCCTACTTTT      22222222222222222222222222222222222222222222222222      NM:i:0  MD:Z:50 AS:i:50 XS:i:0  RG:Z:VB00023_L001
I_13330604_13331055_2:0:0_0:0:0_3dd5    83      I       13331006        60      50M     =       13330604        -452    CTAGCGCGCGCACCGCCCGTGTTGACCAGCGAGAGAGTGTGTGAAGAGTG      22222222222222222222222222222222222222222222222222      NM:i:0  MD:Z:50 AS:i:50 XS:i:0  RG:Z:VB00023_L001
I_13330604_13331055_2:0:0_0:0:0_3dd5    163     I       13330604        60      50M     =       13331006        452     TGGAAATAGTTCAGTTTAAAGCATAATACAAAATTAAAGCTATTTACCTT      22222222222222222222222222222222222222222222222222      NM:i:2  MD:Z:14G6A28    AS:i:40 XS:i:0  RG:Z:VB00023_L001
I_14715941_14716459_1:0:0_0:0:0_4553    99      I       14715941        53      50M     =       14716410        519     CTCAACCCGTTCCCACTTTTCTGTATTATTCCACAATAATCAATATAATT      22222222222222222222222222222222222222222222222222      NM:i:1  MD:Z:38T11      AS:i:45 XS:i:40 RG:Z:VB00023_L001XA:Z:II,+10130993,50M,2;

- 各アラインメントエントリー(それぞれの行)にどんな情報が含まれているのか
1. クエリー名(配列リードの名前)
2. ビット表現フラグ(あとで説明)
3. RNAME: リファレンス名(マッピングした染色体もしくはミトコンドリアなど)
4. クエリー配列の最初のマッピング塩基(一番左)の位置
5. マッピング品質(あとで説明)
6. CIGAR文字列(あとで説明)
7. RNEXTとPNEXTは、ペアエンドの相手リードのリファレンス名と位置を示す。`*`はRNEXTが存在せず、`=`はRNEXTがRNAMEと同じであることを示す。
9. ペアエンドリードのテンプレートの長さ
10. 元のリード配列
11. オリジナルリードの塩基品質

### ビット表現フラグ

Flags:
        0x1     PAIRED        .. paired-end (or multiple-segment) sequencing technology
        0x2     PROPER_PAIR   .. each segment properly aligned according to the aligner
        0x4     UNMAP         .. segment unmapped
        0x8     MUNMAP        .. next segment in the template unmapped
        0x10    REVERSE       .. SEQ is reverse complemented
        0x20    MREVERSE      .. SEQ of the next segment in the template is reversed
        0x40    READ1         .. the first segment in the template
        0x80    READ2         .. the last segment in the template
        0x100   SECONDARY     .. secondary alignment
        0x200   QCFAIL        .. not passing quality controls
        0x400   DUP           .. PCR or optical duplicate
        0x800   SUPPLEMENTARY .. supplementary alignment

- アラインメントに関する多くの情報をビットでコード化している
- マップされたのか
- 相補鎖にマップされたのか
- QC成功したのか
- ペアエンド配列なのか
- 例えばビットフラグ(10進数)が`147`だったとすると２進数では`1001 0011`なので16進数の`0x1`,`0x2`,`0x10`,`0x80`に対応。
- `samtools flags 147`でフラグの内容を教えてくれる

### CIGAR文字列
- ビット表現フラグに続き、アラインメントされた配列についての情報をエンコードするための特殊な方法
- アラインメントに対する塩基のマッチ/ミスマッチ、挿入/欠失、ソフトクリップ/ハードクリップなどの情報をエンコード
- ソフトクリップ：リードの一部しかアラインメントされていない場合
- ハードクリップ:SAMのSEQフィールドに格納された配列の中には存在しない?????
- ex) `43S6M1I26M`:43塩基がソフトクリップ、6塩基がマッチまたはミスマッチ、1塩基分の挿入、26塩基分のマッチが存在。合計の76塩基はSAMエントリーの配列の長さになっている。

### マッピング品質
- `Q = -10log10P`によって与えられる対数確率（塩基品質と同様)
- 下記のような時に使用する
- 不正確なアラインメントのフィルタリング
- ゲノム全体でのマッピング品質分布を評価

## samtoolsについて詳しく
### samtools viewを使ってSAMとBAMを相互変換する
- `samtools view -b celegans.sam > celegans_copy.bam`
- BAMをSAMに変換するときには、`-h`オプションをつけてヘッダーを含める
- 一般的に、ディスク利用効率を高めて処理を高速にするために、ファイルはBAMファイルで保存したほうがいい
- CRAM形式という高度に圧縮されたファイル形式も対応。
- 不可逆圧縮方法なので、元のリードの情報は失われる
- ゲノム配列が必要。EBIが作ったらしく、EBIが推しているらしい by 坊農先生

### sortサブコマンドとindexサブコマンド
- `sort`によって、BAMファイルを特定の順番に並べる
- `-m`でメモリの割り当て、`-@`でスレッド数を指定
- `index`でBAMファイルに索引をつける
- 中間ファイルなしでsamをbamに変換すると同時にsortする by 坊農先生

### アラインメントの抽出
1. index作成
2. 興味のある領域を`samtools view`で検索
3. BED形式で保存された領域が多数ある場合は`-L`オプションで、領域を指定して抽出する`samtools view -L USH2A_exons.bed NA12891_CEU_sample.bam | head -n 3`

SRR003214.11652876      163     1       215796180       60      76M     =       215796224       92      CAGAGCTGTTAAGAGTATTATTTTTTAATAATAAATAGATGATTAAGCCAGAATTAGGTCACATTCAAATATTTAT    <>=<=?A==@8<<:<::?@8?A@A@A3=@>=B5>;B;;8B>9CC=<:C>@;>2DC978>B>@=@A=;887.141)-    XT:A:U  NM:i:0  SM:i:37 AM:i:37 X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:76 RG:Z:SRR003214 OQ:Z:I6IEIIII5I48I=I-9>E28EIAID+5I=3F+3/C1H-6H.E12,C8+2G6&70*7A(.:./54/,*)=)44-,<
SRR010927.6484300       163     1       215796188       60      51M     =       215796213       76      TTAAGAGTATTATTTTTTAATAATAAATAGATGATTAAGCCAGAATTAGGT     ;?=?===>>?@>@@AAAA?A@?A@?AA@>>=?>=?@=?=;=><;=<<7554     XT:A:U  NM:i:0  SM:i:37AM:i:37  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:51 RG:Z:SRR010927  OQ:Z:IIIIIIIIIIIIIIIIIIIIIIIIIIIII2IIIIIIII:EII7IIII;)3G
SRR005667.2049283       163     1       215796190       60      51M     =       215796340       201     AAGAGTATTATTTTTTAATAATAAATAGATGATTAAGCCAGAATTAGGTCA     :?<<<;<<======>>>?>>@>>@@>?>>>?>>??@???@??A???>7<>:     XT:A:U  NM:i:0  SM:i:37AM:i:37  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:51 RG:Z:SRR005667  OQ:Z:IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIE7I9=

### アラインメントのフィルタリング
- ビット表現フラグをクエリーにしてアラインメントを抽出する
- `samtools flags`でビットを確認しながら実行
- ex1) マップされていない全てのリードを出力する方法
- unmapは10進数で`4`という数字で示される(`samtools flags unmap`)
- `samtools view -f 4 bamファイル | head -n 3`
- `-F`を使うと、unmapのリードを全て排除(マップされている全てのリードを表示)
- `samtools view -F 4 bamファイル | head -n 3`
- ex2) リードはアラインメントされていてペアエンドになっているが正しいペアではないものを抽出
- `samtools flags unmap,proper_pair`
- `0x6 6 PROPER_PAIR,UNMAP`
- `samtools view -F 6 -f 1 bamファイル | head -n 3`

## samtools tviewとintegrated genomics viewerによるアラインメントの可視化
- `samtools tview`は、アラインメントデータを可視化できるサブコマンド
- 素早くいくつかのアラインメントを確認するのに便利
- しかしBAMのアラインメント、バリアント、挿入/欠失の調査を行う場合は`Integrated Genomics Viewer(IGV)`がより適している

- `Integrated Genomics Viewer(IGV)`は、Broad Instituteのアプリケーション
- [readme](https://github.com/vsbuffalo/bds-files/blob/master/chapter-11-alignment/README.md)に従い、`brew install igv`でダウンロード
- 最近はJbrowse2がおすすめ by 坊農先生
- 画面中段ではカバレッジとアラインメント情報が記載
- 色付きの文字はリード配列とリファレンスとの間に塩基ミスマッチがあることを示す
- 例えば3種類の異なるハプロタイプが存在することを示している箇所について
- 原因1: ミスアラインメントされた
- 原因2: 誤ったバリアントコール
- 原因3: GGC配列についてIlluminaの特異的なエラー
- `samtools view -H bamファイル名 | grep PL | tail -n 1`などでIlluminaシーケンサーを使っているか確認できる

### samtools パイルアップ形式
- パイルアップ形式とは: アラインメントされたリードを積み重ねることで、染色体の各位置ごとに、リードの塩基を 要約する、プレーンテキスト形式(バリアント1を参照)
- バリアントの同定に利用
- または、サンプルの個体の遺伝型の決定にも使われる
### バリアント解析の例題として`igv`でみた領域の各バリアントについて、本当にバリアントか否かを調査する
1. `samtools mpileup`を実行(バリアント２を参照)
2. `bcftools call`では、バリアントに関する追加情報が中間VCFファイルとして渡されている。それを確認(バリアント3を参照)
- `-v`を使うことで、バリアントしととしてコールされたものだけを呼び出す。
4. `bcftools call`を`-v`フラグなしで呼び出すことで全サイトのクオリティ情報が表示される(バリアント4を参照)
- バリアントとしてコールされなかった位置のクオリティなどの情報を取得できる。
- `215906555`の位置では、`GT=1/1`,`PL=184,7,0`.　`GT=1/1`は`alt`の塩基をホモにもつという意味。 PLは左から順に`ref/ref`,`ref/alt`,`alt/alt`という順で表示されるので、この位置については最も可能性が高い遺伝型は`alt/alt`であり、次が`ref/alt`で、次が`ref/ref`とわかる。
- `samtools mpileup -u`で、`-u`オプションを入れることでBAQを調査(バリアント5を参照)
- 塩基バリアントコールが間違っている確率だけでなく、特定の塩基がミスアラインされている確率も反映する。
- これで`215906547`と`215906548`の両方がバリアントである可能性は考えられなくなるらしい。。

バリアント1
1       215906560       G       15      ,,,....,....... ><=7D>S;=8=B?7>
1       215906561       A       13      ,$,,...,......  :@?.>4>>;E@A<
1       215906562       G       15      ,,....,.......^:.       9<<A=O:?><D?@>?
1       215906563       T       13      ,,..,.......^F. ?<.71?7@<?@7>
1       215906564       C       15      ,$,...,.........        ;98317@0CF>@=6?
1       215906565       C       17      ,....,.........^].^F.   =428>6>==C>?=/?<>
1       215906566       T       18      ,...,...........^F.^],  <<8B:B=EB@??=A@??2
1       215906567       A       18      ,.T.,............,      @A4R@>::C>@=@>>;=>

バリアント2
#CHROM  POS     REF     ALT     QUAL    INFO    FORMAT  NA12891
1       215906546       G       C,<*>   0       DP=23   PL      0,49,239,66,242,246
1       215906547       C       G,<*>   0       DP=22   PL      123,0,103,144,127,233
1       215906548       G       C,<*>   0       DP=22   PL      23,0,163,68,175,207
1       215906549       G       <*>     0       DP=22   PL      0,66,248
1       215906550       G       <*>     0       DP=22   PL      0,63,247
1       215906551       G       <*>     0       DP=22   PL      0,54,242
1       215906552       G       <*>     0       DP=21   PL      0,57,239
1       215906553       G       <*>     0       DP=20   PL      0,48,222
1       215906554       C       G,A,<*> 0       DP=18   PL      0,25,192,27,192,195,39,195,198,198
1       215906555       G       A,<*>   0       DP=19   PL      184,7,0,190,42,204

バリアント3
#CHROM  POS     REF     ALT     QUAL    INFO    FORMAT  NA12891
1       215906547       C       G       90      DP=22   GT:PL   0/1:123,0,103
1       215906555       G       A       157     DP=19   GT:PL   1/1:184,7,0

バリアント4
1       215906546       G       .       29.9006 DP=23   GT      0/0
1       215906547       C       G       90      DP=22   GT:PL   0/1:123,0,103
1       215906548       G       .       11.796  DP=22   GT      0/0
1       215906549       G       .       278     DP=22   GT      0/0
1       215906550       G       .       277     DP=22   GT      0/0
1       215906551       G       .       272     DP=22   GT      0/0
1       215906552       G       .       269     DP=21   GT      0/0
1       215906553       G       .       252     DP=20   GT      0/0
1       215906554       C       .       26.9824 DP=18   GT      0/0
1       215906555       G       A       157     DP=19   GT:PL   1/1:184,7,0
1       215906556       G       .       282     DP=18   GT      0/0

バリアント5
1       215906545       G       <*>     0       DP=22   PL      0,39,215
1       215906546       G       <*>     0       DP=23   PL      0,42,216
1       215906547       C       <*>     0       DP=22   PL      0,21,141
1       215906548       G       <*>     0       DP=22   PL      0,42,200
1       215906549       G       <*>     0       DP=22   PL      0,45,221
1       215906550       G       <*>     0       DP=22   PL      0,48,228

- ファイル形式の詳細
- パイルアップファイルの形式
- BAQ: Base Alignment Qualityの略
- `.`は参照配列の正鎖にマッチ、`,`は参照配列の相補鎖にマッチしたことを示す。
- `大文字の塩基`は正鎖でのミスマッチ、`小文字の塩基`は相補鎖でのミスマッチを示す
- `^`と`$`はリードの開始と終了
- 挿入は`+`、欠失は`-`
- 右の方にある項目は塩基のクオリティを示す
- mpileupの情報は`Variant Call Format (VCF)`として保存できる
- VCF形式
- DPはリードの深さ
- PLはPhred-scaled genotype likelyhoods(遺伝型尤度)
- 上記のような定義はファイルの最初に表示されるヘッダーに書いてある

## Pysamで独自のSAM/BAM処理ツールを作成する
- SAM/BAMを使うためのpythonライブラリー
- [わかりやすい説明ページ](https://bi.biopapyrus.jp/python/module/pysam.html)

### 例1:BAMファイルを開き、領域からアラインメントを取り出し、リードに対する操作を繰り返す
- `pysam.AlignmentFile`でファイルを開く。SAMなら`r`で開き、BAMなら`rb`で開く
- 開いたファイルを`for read in bamfile`とかで1行ずつ読み込んで処理するのが基本

### 例2: AlignmentFileオブジェクトからのSAM/BAMヘッダー情報の抽出
- `bafile.header.keys()`でDictとして保管されているヘッダー情報を確認
- 例えば`RG`セクションの3番目のレコードが見たければ`bamfile.header['RG'][3]`で確認できる
- 参照配列名とその長さは`bamfile.reference`、`bamfile.lengths`

### 例3: AlignedSegmentオブジェクトを操作する
- 例えばソフトクリップされたリードを見つけてプログラムを終了させる方法
- forループでリードを順にみていって、ソフトクリップを表すCIGARの`S`が現れたらループをBreakする方法を使う

for read in bamfile:
    if `S` in read.cigarstring:
        break
read.cigarstring
`35M16S`

### 例4: アラインメント統計を記録するプログラムの作成
- 内容は本を参照
- 注意点としては、結果が正しいものか検証することが大事
- 今回の場合は`samtools flagstat`というコマンドで統計値を検証している

### その他のPysam機能、およびその他のSAM/BAM API
- `AlignmentFile.pileup()`パイルアップデータを生成
- `pysam.sort`,`pysam.view`など
- `Picard`はJava APIを提供
- SAMファイルのヘッダーがおかしい場合などはPicardが有効らしい
- 発音方法はピカード
#### その他のアラインメントパッケージ
- `Rsamtools`
- `GenomicAlignments`