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」ボタンをクリック。

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

RNA-seqデータの解析例をメモ。
あくまで一例なので、ここで取り上げた方法がすべてのデータに対して良いとは限りません。

一つ一つコマンドを入力していってもいいのですが、手間なのでシェルスクリプトを活用し複数の処理（コマンド）をまとめて行った例を示します（いわゆる、バッチ処理）。

以下では、RNA-seqデータの解析例をシェルスクリプトを示しました。

#!/bin/sh
for file in "$@"
do
#Illumina CASAVA Filter [Y]を除去
fastq_illumina_filter -N -o /path/to/"$file"_casava.fastq /path/to/"$file".fastq

#Phred quality score 20未満が80%以上のリードを除去
#Phred quality score 20未満の末端を除去
#配列長が10bp未満のリードを除去
fastq_quality_filter -Q33 -q 20 -p 80 -i /path/to/"$file"_casava.fastq | fastq_quality_trimmer -Q33 -t 20 -l 10 -o /path/to/"$file"_filtered.fastq

#Ribosomal RNA（rRNA）由来のリードの除去
bowtie -p 8 --un /path/to/"$file"_norrna.fastq /path/to/bowtie-X.XX.X/indexes/Ribosomal_RNA /path/to/"$file"_filtered.fastq > /dev/null

#Reference genome / Transcriptomeへのマッピング
tophat -p 8 -o /path/to/"$file"_out -G /path/to/Refseq_genes_May_14_2013.gtf /path/to/bowtie-X.XX.X/indexes/hg19 /path/to/"$file"_norrna.fastq
done 

  cuffdiff -p 8 -o /path/to/cuffdiff_out /path/to/Refseq_genes_May_14_2013.gtf <CTRL_rep1.bam,CTRL_rep2.bam…> <KD_rep1.bam,KD_rep2.bam…>

〈使用したソフトウェア〉
・Illumina CASAVA-1.8 FASTQ Filter
・FASTX-Toolkit 0.0.13
・Bowtie 0.12.7
・SAM tools 0.1.18
・TopHat 2.0.4
・Cufflinks 2.1.1

〈使用したファイル〉
1) Ribosomal_RNA
Ribosomal RNA（18S, 28S）のIndexファイル
⇒NCBIなどからFasta形式で配列情報を取得し、そのファイルからBowtie-buildコマンドでIndexファイルを事前に用意しておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。

2) hg19
Human genome（hg19）のIndexファイル
http://bowtie-bio.sourceforge.net/index.shtml
⇒Bowtieのサイトの右側に、「Pre-build indexes」の項目があり、様々な種のゲノム配列情報を取得することができる。Fastqファイルをゲノムにマッピング際に必要なので、事前に該当のデータをダウンロードしておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。

3) Refseq_genes_May_14_2013.gtf
Illumina igenomesの「transcript.gtf」ファイル
http://cufflinks.cbcb.umd.edu/igenomes.html

⇒上記のURLから「Homo sapiens - UCSC - hg19」を選択、該当のファイル（20GBを超える）をダウンロード。その中に「transcript.gtf」があり、このファイルの中にRefseqのアノテーション情報が記載されている。

⇒目的に応じて、gencode / Ensemblなどの他のgtfファイルを利用することも可。ただ、多くのIsoformが登録されているので、どこまで正確に各Isoformごとの定量値を見積もることができるのか疑問。（そもそも、Isoformごとの正確な定量値を必要なのかどうかにもよる。）

〈シェルスクリプトを利用する際の注意点〉

rootユーザー以外が上記のシェルスクリプトを使用する場合、
（１）一行目に「#!/bin/sh」を書く
（２）スクリプトのファイルに実行権限を与える
の２点が必要になります。

（２）に関しては、以下のcommandを実行し、事前に実行権限を付与しておきましょう。

$ chmod +x /path/to/RNA-seq_pipeline.sh

〈シェルスクリプトの内容について〉
（１）クオリティチェック～マッピング(ファイル名: RNA-seq_pipeline.sh)
順を追って説明いきます。
1行目の「#!/bin/sh」はこのファイルがシェルスクリプトであることを明示するためのものであり、指定したパス「/bin/sh」のプログラムによって該当のファイルを実行させるのに必要です。

2行目はループ関数である「for」を用いて、複数のFastqファイルに対して、doからdoneの間の処理を繰り返し行わせています。ここで、「$@」は1つ以上のコマンドライン引数（SRRxxxxxxのようなファイル名）を受け取り、「file」に代入しています。

5行目では、「Illumina CASAVA-1.8 FASTQ Filter」のソフトを利用し、CASAVA filterをパスしなかったクオリティの低いリードを除去します（Illumina以外のプラットフォームから得られたFastqファイルやNCBI GEOやDDBJ SRAからダウンロードしてきたFastqファイルに対しては不要）。

HiseqなどのIlluminaのシークエンサーから得られたFastqファイルには、CASAVA filterをパスしたかどうかの目印が各リードに与えられていますので（フィルターをパスしていれば「N」、フィルターをパスしていなければ「Y」）、ソフトを利用しフィルターをパスしたリードのみ抽出します。

「-N」オプションによりフィルターをパスしたリードのみを出力させ、「-o」オプションにより、指定したディレクトリに出力したFastqファイルを保存させます。

10行目では、さらにクオリティの低いリードのフィルタリング、及びクオリティの低い塩基のトリミングを行なっています。具体的には、Phred quality score 20未満が80%以上のリードを除去、Phred quality score 20未満の末端を除去、配列長が10bp未満のリードを除去を行なっています。

ここで、Phred quality scoreとは、Fastqファイル中に記載されているリードのクオリティを示す指標のことです。Phred quality scoreは以下の式によって表されます。

Phred quality score (Q) = -10 log10P

Pは「Possibility」、つまりはシークエンサーから読まれた各塩基の誤差率を示しています。例えば、P = 20であれば、誤差率は1%であり、100回に1回読み間違いが起こることを意味しています。そのため、リードのクオリティチェックでは、リードの各塩基のクオリティが「P >= 20」もしくは「P >= 30」ぐらいになるようにフィルターをかけるとよいのではないかと思います（シークエンサーで読んだリード長に応じてクライテリアを設定）。

13行目では、Ribosomal RNA由来のリードを排除します。「--un」オプションは指定したディレクトリにマッピングされなかったリードをFastqファイルで出力させます。また、Bowtieでは、コマンドライン引数は「2つ」指定する必要があり、１つ目は、Fastqファイル中のリードをマッピングするリファレンス配列（ゲノム配列・その他の配列）を指定し、2つ目は、マッピングしたいFastqファイルを指定します。上記では、コマンドライン引数の1つ目が「Ribosomal_RNA」、2つ目が「"$file"_filtered.fastq」に該当します。また、標準出力としてマッピングされたリード（Ribosomal RNA由来のリード）を「/dev/null」（ゴミ箱のようなもの）へリダイレクトするように指定しています。

「-p」オプションでは、パソコンのコア数を指定しています。デフォルトでは「-p 2」となっており、2コア分で計算を行なっています。CPUのコア数が2よりも多い場合は、自分のパソコンのコア数に応じて指定したほうがよいです。その分、計算速度も向上します。

16行目で、いよいよリファレンス配列へのマッピングを行なっていきます。対象がヒトである場合、パラメータはいじくらなくてもヒト用に最適化されている（はず）なので、デフォルトでとりあえずマッピングしてます。

「-o」オプションで出力先のディレクトリを指定します（ディレクトリがない場合は生成してくれます）。

「-G」オプションで使いたいアノテーション情報（gtfファイル）を指定します。また、コマンドライン引数についてはBowtieの時と同様です。

（２）FRPM値(RPKM値)の算出
Cuffilnksのコマンドである「Cuffdiff」を用いてFPKM値（RPKM値）の算出と発現変動遺伝子の同定を行います。

コマンドライン引数は3つ指定する必要があります。1つ目は、アノテーション情報（gtfファイル）を指定し、2つ目と3つ目には比較したいデータ（bamファイル）をそれぞれ指定します。Duplicateのデータ（n >= 2）がある場合は、各データ（bamファイル）を「,」カンマで区切り指定することができます（比較するデータをn = 1でしか取っていない場合、有意差を計算することができません。そのため、この場合Fold-change（例えば2以上）をもとに発現変動遺伝子を絞り込む必要があります）。

〈シェルスクリプト実行例〉
/path/to/RNA-seq_pipeline.sh SRRxxxxxx SRRyyyyyy SRRzzzzzz
⇒シェルスクリプトを呼び出し、コマンドライン引数として処理したいFastqファイルの名前を指定。