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ファイルの名前を指定。