2010年11月10日水曜日

SAMファイルをとりあえずビューワーで見る

Bowtieなどのマッピングツールでリード配列をリファレンス配列にマッピングすると、SAM (Sequence Alignment Map) というフォーマットのファイルができる。
SAMは実際、マッピングファイルの世界標準になりつつある。
サイズは、かなり大きい。 中身を見てどうこうということは無いけど、先ず、リファレンスに使った配列の情報がずらっと並ぶ。 @SQ SNから始まり、LN:の前までの、例えばgi|224589801|ref|NC_000010.10|が、リファレンス配列の名前だ。ちなみにこれはヒト10番染色体。
@HD VN:1.0 SO:unsorted
@SQ SN:gi|224589801|ref|NC_000010.10| LN:135534747
@SQ SN:gi|224589802|ref|NC_000011.9| LN:135006516
@SQ SN:gi|224589803|ref|NC_000012.11| LN:133851895
@SQ SN:gi|224589804|ref|NC_000013.10| LN:115169878
@SQ SN:gi|224589805|ref|NC_000014.8| LN:107349540
@SQ SN:gi|224589806|ref|NC_000015.9| LN:102531392
@SQ SN:gi|224589807|ref|NC_000016.9| LN:90354753
@SQ SN:gi|224589808|ref|NC_000017.10| LN:81195210
@SQ SN:gi|224589809|ref|NC_000018.9| LN:78077248
@SQ SN:gi|224589810|ref|NC_000019.9| LN:59128983
@SQ SN:gi|224589800|ref|NC_000001.10| LN:249250621
@SQ SN:gi|224589812|ref|NC_000020.10| LN:63025520
@SQ SN:gi|224589813|ref|NC_000021.8| LN:48129895
@SQ SN:gi|224589814|ref|NC_000022.10| LN:51304566
@SQ SN:gi|224589811|ref|NC_000002.11| LN:243199373
@SQ SN:gi|224589815|ref|NC_000003.11| LN:198022430
@SQ SN:gi|224589816|ref|NC_000004.11| LN:191154276
@SQ SN:gi|224589817|ref|NC_000005.9| LN:180915260
@SQ SN:gi|224589818|ref|NC_000006.11| LN:171115067
@SQ SN:gi|224589819|ref|NC_000007.13| LN:159138663
@SQ SN:gi|224589820|ref|NC_000008.10| LN:146364022
@SQ SN:gi|224589821|ref|NC_000009.11| LN:141213431
@SQ SN:gi|224589822|ref|NC_000023.10| LN:155270560
@SQ SN:gi|224589823|ref|NC_000024.9| LN:59373566
そしてBowtieのパスが書かれていて、そのあとから、各リード配列のマッピング情報がダーッと書かれるわけ。
1279_2_116_F3 4 * 0 0 * * 0 0 CCGTACGTCGTGGGTAGGGGCNGTGAGTCGCTTCGGGTCGAGGATCTGG =;?89,:4:=2,57-&*'&&/!,)&/.&),49,'1)-)'(,12+//-7) XM:i:0
これは人間が読むわけでは無いので、どうしても意味が知りたいひとはここを参照。
http://samtools.sourceforge.net/SAM1.pdf

さて、このファイルをとりあえずフリーのビューワーで見てみよう。
例として、IGV (Integrative Genomics Viewer)で見てみる。
登録(無料)してからダウンロードする。これはWindowsで動くので敷居が低い!

ところでさっきのSAMファイル、リファレンスの名前がgi|224589801|ref|NC_000010.10| のようになっているが、このままではIGVが認識しない。
ビューワーが認識できる染色体名に変換してあげなくてはいけない。
IGV なら、gi|224589801|ref|NC_000010.10| は、10 に変換する。以下24本の染色体も同じ。
やり方は色々あるだろうが、一番簡単なのは、Perlスクリプトを使って、
perl -pe "s/gi\|224589801\|ref\|NC_000010.10\|/10/g" SD_Agilent_Exome_F3.sam > temp1.sam
perl -pe "s/gi\|224589802\|ref\|NC_000011.9\|/11/g" temp1.sam > temp2.sam
perl -pe "s/gi\|224589803\|ref\|NC_000012.11\|/12/g" temp2.sam > temp.sam
perl -pe "s/gi\|224589804\|ref\|NC_000013.10\|/13/g" temp.sam > temp1.sam
のような置換コマンドを24本分繰り返すと、最終的に全ての染色体名をIGV用の染色体名に置き換えることができる。 こんな感じになった。

@SQ SN:10 LN:135534747
@SQ SN:11 LN:135006516
@SQ SN:12 LN:133851895
@SQ SN:13 LN:115169878
@SQ SN:14 LN:107349540
@SQ SN:15 LN:102531392
@SQ SN:16 LN:90354753
@SQ SN:17 LN:81195210
@SQ SN:18 LN:78077248
@SQ SN:19 LN:59128983
@SQ SN:1 LN:249250621
@SQ SN:20 LN:63025520
@SQ SN:21 LN:48129895
@SQ SN:22 LN:51304566
@SQ SN:2 LN:243199373
@SQ SN:3 LN:198022430
@SQ SN:4 LN:191154276
@SQ SN:5 LN:180915260
@SQ SN:6 LN:171115067
@SQ SN:7 LN:159138663
@SQ SN:8 LN:146364022
@SQ SN:9 LN:141213431
@SQ SN:23 LN:155270560
@SQ SN:24 LN:59373566
さて、そうしてSAMファイルをIGVでも認識できるようにしたら、インストールしたIGVを立ち上げてみよう。
そうしたら、次に、File > Run igvtools... で IGV tool を起動する。
そこで先ず、さっきのSAMファイルを染色体ごとにソートする。 CommandでSortを選び、リファレンスの名前をIGV用に置き換えたSAMファイルを指定して、Runする。

その次に、インデックスを作成する。同じくIGV tool上にて、CommandでIndexを選び、今ソートした結果ファイルを指定する。
そうすると拡張子 sai のファイルができる。
ここまできたら、さあ、IGVで取り込もう!

File > Load File で、ソート済みのSAMファイルを指定すると、自動的にインデックスファイルも認識されて取り込まれるはずだ。
しかし、最初は何も見えないだろう。 こんな風に
染色体を適当に選んでから、ズームインしてみよう!

こんな感じに見える。

まとめ: SAMファイルは、ちょっとリファレンスの名前を変更するめんどくささはあるけど、フリーのビューワーでちゃんと見ることができる。

ビューワーについてはまた今度。

0 件のコメント:

コメントを投稿