RNA-seqデータの解析パイプラインを作ろう（データの視覚化）

RNA-seqのデータ(fastqファイル)をBowtieやTopHatなどでマッピングした後、bamファイルやsamファイル（場合によってはbedファイル）として出力されます。

そのデータを用いて、Cufflinksなどのリードの集計ソフトで定量値（FRKM値やRPKM値）を算出し、各転写産物の発現量を見積もると同時に、実際にゲノム上にマッピングされたリードを確認したいケースもあると思います。

そこで、今回はbam / sam / bedファイルから視覚化データを作る方法例を示したいと思います。個人でマッピングされたリードを確認する際にはIGVが便利だと思いますが、他人とデータを共有する場合では、UCSC genome browser上で可視化させたほうが都合が良いと思います。

〈スクリプトのダウンロード〉
今回は、github上で公開されているPythonスクリプトを利用します。このスクリプトはメモリの消費が少なく、複数のファイル形式（sam / bed / bowtie）に対応しているので、使い勝手が良く個人的におすすめです。

（１）rdbio-scriptsのURLに行き、右下にある「Download ZIP」をクリックし、bed2wig.pyを含む全スクリプトをダウンロードする。もしくはbed2wig.pyのURLへ移動し、スクリプトをコピペして保存する。

bed2wig.pyのスクリプト
https://github.com/daler/rdbio-scripts/blob/master/sequenceFiles/bed2wig.py
rdbio-scripts
https://github.com/daler/rdbio-scripts

〈スクリプト使用例〉
基本的な使い方は、

bed2wig.py -h

でヘルプを見ることができます。
ここでは、bamファイルをwigファイルに変換する方法について紹介したいと思います。

（１）bamファイルをソートし、samファイルに変換する。（TopHatから出力されるbamファイルは予めソートされていたと思うのでこの作業は必要ないかも。）

samtools sort accepted_hits.bam accepted_hits_sorted.bam
samtools view -F 0x0004 accepted_hits_sorted.bam > accepted_hits_sorted.sam

⇒「-F」オプションで0x0004を指定することでマッピングされなかったリードを排除しています。

（２）samファイルからwigファイルに変換する。

bed2wig.py -i accepted_hits_sorted.sam --type sam -o  accepted_hits_sorted.wig --verbose

⇒「-i」オプションでインプットファイルを指定します。「-type」オプションでインプットファイルのファイル形式を指定します。「-o」オプションでアウトプットファイルを指定します。「--verbose」オプションで進捗状況を出力させます。

（３）wigファイルの１行目を書き換える。
lessコマンドやviでファイルの中身をみてみると、１行目は「track type=wiggle_0 alwaysZero=on」となっています。このままだと使いにくいので、下記のように、適当に書き換えておきます。（もちろん、bed2wig.pyを直接書き換えても良いと思います。たぶん、269行目を書き換えればOK.。）

track name='CTRL_TopHat' description='CTRL_TopHat' type='wiggle_0' color=255,128,128 yLineMark=0.0 yLineOnOff=on visibility=2 maxHeightPixels=40:40:20 priority=1

「track name」と「description」でデータの名称を各ファイルごとに書き換えておきましょう。UCSC genome browserにアップロードする場合は、ここで指定した名称が以前アップロードしたデータと重複していると、アップロードした時に前のデータが消えてしまうので注意が必要です。

（４）ファイルを圧縮する。
UCSC genome browserではgzip (.gz), compress (.Z), or bzip2 (.bz2)の３種類の圧縮形式に対応しており、ファイルが大きくてアップロードに時間がかかりそうな場合、ファイルを圧縮しておきましょう。

bzip2 -c accepted_hits_sorted.wig > DATA.wig.bz2

⇒「-c」オプションを指定し、標準出力でファイルを出力。

〈UCSC genome browserへのアップロード〉
（１）Homeへ移動し、左側にある「Genome Browser」をクリック。
http://genome.ucsc.edu/index.html

（２）「manage custom track」ボタンをクリック。（以前にアップロードしたファイルが残っていると、「Manage Custom Tracks」ページに移動するので、その場合は「add Custom Tracks」ボタンをクリック。）

（３）「add Custom Tracks」ページに移動するので、「Paste URLs or data:」の「Or upload:」の「ファイルを選択」をクリックし、アップロードしたいwigファイルを指定する。

（４）「Submit」ボタンをクリック。ファイルがアップロードされるので、アップロードが完了するまで放置する。（アップロードしている間、ページ移動をしないこと。）

（５）アップロード完了後、「Manage Custom Tracks」ページに移動します。アップロードしたファイルは「go to Genome browser」ボタンをクリックし、Genome browser上で確認することができます。

〈Custom tracksの維持・構成の保存〉
Custom trackをアップロードしても、一週間ぐらいすると設定がリセットされて見れなくなってしまいます。そこで、Custom tracksおよびその他のTrackの構成を保存するために、UCSC genome browserが提供している「Sesson」機能を利用しましょう。

「Sesson」では、アカウントを取得することで、Custom trackの維持やTrackの構成の保存を行うことができます。また、アップロードしたCustom trackを他人が見れるように設定することもできます。

（１）Homeへ移動し、上の「My Data」→「Sessions」をクリック。
http://genome.ucsc.edu/index.html

（２）「Creat an account」をクリックし、必要なデータを入力しアカウントを作成する。

（３）Sessionにログインした後、「Save Settings」で「Save current settings as named session:」にSessionの名称を入力。

自分のアカウントにログインしていない他人にデータが見られるようにしたい場合は「allow this session to be loaded by others」のボックスにチェックをつける。

そして、「Submit」ボタンをクリックし、Sessionの保存を行う。（もちろん、Custom trackの追加やgenome browser上の構成を自分なりにアレンジして保存したい組み合わせを事前にgenome browser上で組んでおく。）

⇒Sessionを上書きしたい場合は、上書きしたいSession名を入力した後、「Submit」ボタンをクリック。上書きしてよいか確認を求められるので、「Yes」ボタンをクリック。