RNA-seqデータ解析(Alternative splicingの解析, Strand-specific, Single-end)

概要

RNA-seqのデータを用いて、Alternative splicingのイベント(後述)の変化を検出する方法(2群間のデータの比較)を紹介します。

ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、Alternative splicingのデータ解析の詳細をStep by Stepで説明します。

ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。 パソコン(もしくはスパコン)のメモリは16-32GB以上を確保することが望ましいです。

使用するデータ

RNA-seq (Single-end, 101bp; Control siRNA, UPF1 siRNA)

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148

使用するソフトウェア

解析ワークフロー

Alternative splicingを解析する方法として、MISO, DEXSeq, rMATSそれぞれを用いた方法を紹介したいと思います。

TopHat-MISO workflow

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-03-28 14.34.30.png (48.7 kB)

TopHat-featureCounts-DEXSeq workflow

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-03-28 14.34.56.png (53.5 kB)

TopHat-rMATS workflow

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-03-28 14.34.42.png (48.6 kB)

各ソフトウェアの特徴(主観的な意見です…)

MISO
  • 最も引用数が多く、精度も高い(個人的な意見です)。Alternative splicingの解析では、デファクトスタンダードのソフト。
  • skipped exons (SE), mutually exclusive exons (MXE), alternative 3' splice sites (A3SS) , alternative 5' splice sites (A5SS), retained introns (RI)のパターンごとに解析してくれる。
  • n=1のデータしか扱えない(n=2同士のデータを比較することができない)。
  • 扱えるアノテーション情報が限定されている(UCSC gene, RefSeq, Gencodeの組み合わせ。UCSC genome browserにアップロードされている最新のデータにしか対応していない)。
  • 結果ファイルの中身がSplicingのイベントごとに評価されているため、どの遺伝子でSplicingのパターンの変化が起こっているかわかりにくい。
  • p-valueによる有意差検定でなく、Bayes factorを計算しているので、どこでThresholdを引いたらいいかわかりにくい。
DEXSeq
  • MISOについで有名なソフト。精度自体はイマイチなイメージ(個人的な意見です)。
  • Duplicatesに対応しているので、n>=2のデータ同士の比較が可能(内部的にはDESeqのアルゴリズムで2群間比較している)。
  • アノテーション情報は自由に選択できる。
  • 遺伝子のIDとTranscriptのIDがタグ付けされているのでリストが見やすい。
  • Exonブロックごとに評価する(この方法がイマイチな印象)ので、SE, MXE, A3SS, A5SS, RIのどのパターンの変化か判別しづらい。
rMATS
  • あまり引用されてないソフト。
  • Duplicatesに対応しているので、n>=2のデータ同士の比較が可能。
  • Novel splicing sitesの検出も可能。
  • アノテーション情報を自由に選択できる。
  • 遺伝子のIDとGene symbolがタグ付けされているのでリストが見やすい。
  • SE, MXE, A3SS, A5SS, RIごとに解析してくれる。
  • 3つのソフトの中で、最も計算時間がかかる(長いと1日半くらいかかる)。

サンプルコード

Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/RNA-seq_for_alternative_splicing_Tutorial

RNA-seq_strand-specific_single-end_1_mapping_for_AS.sh

TopHatによるを行うスクリプト。FASTQファイルを引数として指定。

$ ./RNA-seq_strand-specific_single-end_1_mapping_for_AS.sh  hoge.fastq
  • 10行目: gtfFileにGTFファイル(アノテーション情報)を指定。
  • 11行目: indexContamFileにrRNAなどのコンタミ配列のBowtie用のIndexファイルを指定。
  • 12行目: indexGenomeFileにGenome配列のBowtie用のIndexファイルを指定。
  • 24行目: -rオプションでHouse-keeping genesの遺伝子領域のBEDファイルを指定。
  • 29行目: -rオプションでGencodeアノテーション情報に関するBEDファイルを指定。
MISO_test.sh

MISOでリードカウントを行うスクリプト。FASTQファイルを引数として指定。

$ ./MISO_test.sh hoge.fastq  
  • 14行目: MISO_setting_fileにMISOの設定ファイルを指定。
MISO_comparison.sh

MISOで2群間比較を行うスクリプト。引数として、2群のFASTQファイルを指定する。

$ ./MISO_comparison.sh control.fastq knockdown.fastq
DEXSeq_test.sh

featureCountsでリードカウントを行い、DEXSeqで2群間比較を行うスクリプト

$ ./DEXSeq_test.sh  
  • 8行目: gtfFileにDEXSeq用に用意したアノテーション情報に関するGTFファイルを指定。
  • 9行目: resultDirに保存先のディレクトリ名を指定。
  • 10行目: featureCountFileにfeatureCountsの出力ファイルの名前を指定。
  • 11行目: resultFileにDEXSeqの出力ファイルの名前を指定。
  • 18-23行: リードカウントするBAMファイルをスペース区切りで指定。
  • 30行目: 18-23行で指定したBAMファイルのラベル名をカンマ(,)区切りで指定。
  • 31行目: 18-23行で指定したBAMファイルの属性(Control or Knockdownなど)をカンマ(,)区切りで指定。
rMATS_test.sh

rMATSで2群間比較を行うスクリプト

$ ./rMATS_test.sh  
  • 8行目: gtfFileにGencodeなどのアノテーション情報に関するGTFファイルを指定。
  • 10行目: -b1にControlのBAMファイルをカンマ(,)区切りで指定。
  • 11行目: -b2にKnockdownのBAMファイルをカンマ(,)区切りで指定。

~/custom_command/<hoge.py>となっている箇所はカスタムのPythonスクリプトが置いてあるパス・ファイル名を指定する。

シェルスクリプトを使用するときの注意点

権限の問題

実行権限をシェルスクリプトに与えた後、実行すること。

$ chmod 755 hoge.sh

改行コードの問題

WindowsLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。

ここではAtomと呼ばれるテキストエディタでの設定方法を説明する。

  • [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
  • 左のメニューから、Packagesをクリック。
  • Installed Packagesでline-ending-selectorと検索。
  • Core Packagesの項目から、line-ending-selectorのSettingsをクリック。
  • 設定画面で、Default line endingをLFに変更。

Settings — C__Users_Naoto_OneDrive_NGS解析_NGS-Tutorial_RNA-seq_Tutorial — Atom 2017-01-31 11.20.21.png (250.6 kB)

これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。

各ステップの説明

0. 必要なソフトウェアのインストール

まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。

imamachi-n.hatenablog.com

基本的なソフトのインストール

$ conda install fastqc
$ conda install fastx_toolkit
$ conda install bowtie=1.12
$ conda install tophat
$ conda install samtools
$ conda install rseqc
$ conda install bedtools
$ conda install subread
$ conda install bioconductor-dexseq
注意!!

tophatによるマッピングで最新のbowtieを使うとエラーが出るので、以前のバージョンのbowtie v1.12を代わりにインストールする。

MISOとrMATSのインストールは一筋縄ではいかないので、個別にインストール方法について説明します(あくまで、自分の解析環境にインストールしたときのやり方で、一般的な方法ではありません)。

MISOのインストール

インストールするもの
  • fastmiso (MISO)
  • rnaseqlib
インストールの方法

詳細は下記のURLを参照のこと。
http://miso.readthedocs.io/en/fastmiso/
http://rnaseqlib.readthedocs.io/en/clip/

まず、biocondaで新しい仮想環境を作ります。ここでは、python2.7をベースに仮想環境を構築します。

$ conda create -n miso python=2.7
  • -nで仮想環境の名前を指定します。
  • pythonに2.7を指定することで、python2.7を仮想環境にインストール

source activateに続けて、先ほど作成した仮想環境の名称を指定することで、仮想環境に入ることができます。

$ source activate miso

conda info --envsで現在どの環境にいる確認できるので、前後でどの環境にいるか確認してみましょう。

$ conda info --envs
  # conda environments:
  #
  miso                     /home/akimitsu/miniconda2/envs/miso
  root                  *  /home/akimitsu/miniconda2

$ source activate miso
$ conda info --envs
  # conda environments:
  #
  miso                 *  /home/akimitsu/miniconda2/envs/miso
  root                     /home/akimitsu/miniconda2

最初はrootという環境にいたのが、source activate misoの後、misoという仮想環境に移動していることがわかると思います。

MISOの実行は、このmisoと名付けた仮想環境上で行うこととします。また、その他の必要なソフトウェアに関しても同様の環境上にインストールしていきます。

それでは、fastmiso, rnaseqlibに必要なパッケージやソフトウェアをBiocondaを利用してあらかじめインストールしておきましょう。

$ conda install matplotlib
$ conda install scipy 
$ conda install pysam
$ conda install pandas
$ conda install pybedtools
$ conda install cutadapt
$ conda install samtools=0.1.19
$ conda install bedtools
$ conda install gffutils
注意!!

MISOは内部でSamtoolsを利用していますが、バージョンはv1.0未満のものに対応しています。なので、v1.0以上だとエラーが発生してMISOが動きません。

conda install samtools=0.1.19
fastmiso (MISO)のインストール

下記のURLにアクセスして、
http://genes.mit.edu/burgelab/miso/software.html
MISO from the latest GitHub repository can be downloaded as a zipfile without going through GitHub here:
* C version of MISO (fastmiso)

をクリックしてサイトから直接ソースコードをダウンロードします。 PyPyからインストールしたものだとうまく機能しないことがありました。
http://genes.mit.edu/burgelab/miso/software.html

ダウンロード後、ファイルを解凍します。

$ unzip rden-MISO-09b2a92.zip
yarden-MISO-09b2a92/misopy/credible_intervals.pyスクリプトの修正

現行のソースコードでは、MISOの実行中にcredible_intervals.pyでエラーが発生してしまいます。エラーを回避するために、面倒ですがソースコードを訂正する必要があります。

46行目:
lower_bound_indx = round((alpha/2)*num_samples) - 1から
lower_bound_indx = int(round((alpha/2)*num_samples) - 1)に変更

49行目:
upper_bound_indx = round((1-alpha/2)*num_samples) - 1から
upper_bound_indx = int(round((1-alpha/2)*num_samples) - 1)に変更

ちなみに元のソースコードだと、54行目で以下のようなエラーが発生します。

cred_interval = [samples[lower_bound_indx], samples[upper_bound_indx]]
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

これは、lower_bound_indxupper_bound_indxがint型でなく、float型であることに起因したエラーだと考えられます(Indexとして整数でない数値を指定してしまっている)。

そのため、エラーを回避するために、途中でlower_bound_indxupper_bound_indxをint型に型変換を行うように書き換える必要があったということです。

ソースコードを修正し終わったら、MISOを現在の仮想環境にインストールします。

$ cd rden-MISO-09b2a92
$ python setup.py install
インストールのチェック

MISO側で、必要なパッケージがちゃんとインストールされているかどうかのチェックを行うコマンドを用意してます。実行してエラーが出ないことを確認します。

$ module_availability
Checking availability of Python modules for MISO
Looking for required Python modules..
Checking for availability of: numpy
Checking for availability of: scipy
Checking for availability of: json
Checking for availability of: matplotlib
Checking for availability of: pysam
All modules are available!
Looking for required executables..
Checking if samtools is available
  - samtools is available
Checking if bedtools is available
  - bedtools is available

また、ユニットテストも用意されているのであわせて実行します。

$ python -m unittest discover misopy
Testing fr-unstranded...
Checking read  f_read  against  +
Checking read  f_read  against  -
Checking read  r_read  against  +
Checking read  r_read  against  -
Testing fr-firststrand...
.Testing conversion of SAM to BAM...
Executing: sam_to_bam --convert /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-data/sam-data/c2c12.Atp2b1.sam /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output
Converting SAM to BAM...
  - Executing: samtools view -Sbh /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-data/sam-data/c2c12.Atp2b1.sam  > /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.bam
[samopen] SAM header is present: 35 sequences.
Sorting BAM file...
  - Executing: samtools sort /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.bam /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted
Indexing BAM...
  - Executing: samtools index /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted.bam
Conversion took 0.00 minutes.
.Testing gene-level Psi...
Testing GFF indexing of: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1.mm9.gff
Executing: index_gff --index /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1.mm9.gff /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed
Indexing GFF...
/sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1.mm9.gff appears to already be indexed. Aborting.
Executing: miso --run /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted.bam --output-dir /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/gene-psi-output --read-len 36
MISO (Mixture of Isoforms model)
Probabilistic analysis of RNA-Seq data for detecting differential isoforms
Use --help argument to view options.

Using MISO settings file: /home/akimitsu/miniconda2/envs/miso/lib/python2.7/site-packages/misopy-0.5.3-py2.7-linux-x86_64.egg/misopy/settings/miso_settings.txt
Computing Psi values...
  - GFF index: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed
  - BAM: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted.bam
  - Read length: 36
  - Output directory: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/gene-psi-output
Checking your GFF annotation and BAM for mismatches...
Checking if BAM has mixed read lengths...
Found reads of length 36 in BAM.
Mapping genes to their indexed GFF representation, using /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed
Searching for /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/genes_to_filenames.shelve..
  - File not found.
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/compressed_ids_to_genes.shelve.bak
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/genes.gff
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/genes_to_filenames.shelve.bak
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/compressed_ids_to_genes.shelve.dat
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/genes_to_filenames.shelve.dir
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/genes_to_filenames.shelve.dat
Skipping: /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1/indexed/compressed_ids_to_genes.shelve.dir
Preparing to run 1 batches of jobs...
Running batch of 1 genes..
  - Executing: python /home/akimitsu/miniconda2/envs/miso/lib/python2.7/site-packages/misopy-0.5.3-py2.7-linux-x86_64.egg/misopy/run_miso.py --compute-genes-from-file "/sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/gene-psi-output/batch-genes/batch-0_genes.txt" /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted.bam /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/gene-psi-output --read-len 36  --overhang-len 1 --settings-filename /home/akimitsu/miniconda2/envs/miso/lib/python2.7/site-packages/misopy-0.5.3-py2.7-linux-x86_64.egg/misopy/settings/miso_settings.txt
  - Submitted thread batch-0
Waiting on 1 threads...
.Testing gene-level Psi...
Executing: python /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/index_gff.py --index /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1.mm9.gff /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/indexed
Indexing GFF...
/sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/genes/Atp2b1.mm9.gff appears to already be indexed. Aborting.
Executing: python /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/run_events_analysis.py  --compute-genes-psi /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/gff-events/mm9/indexed /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/sam-output/c2c12.Atp2b1.sorted.bam --output-dir /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/gene-psi-output --read-len 36  --paired-end 250 30 --use-cluster
MISO (Mixture of Isoforms model)
To run MISO, please use "miso" instead.
.Testing single-end SE event interface...
Executing: python /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/run_events_analysis.py  --compute-events-psi se-sample /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-data/se-counts/se_test.counts --output-dir /sshare4/home_mig/akimitsu/software/yarden-MISO-09b2a92/misopy/test-output/SE-output --read-len 35 --overhang-len 4  --event-type SE --use-cluster
MISO (Mixture of Isoforms model)
To run MISO, please use "miso" instead.
.
----------------------------------------------------------------------
Ran 5 tests in 9.273s

OK

エラーが出ずに、最後にOKが出力されればインストールが適切に行われたと判断します。

rnaseqlibのインストール

次に、MISO用のアノテーション情報を用意するためのスクリプト群をまとめたrnaseqlibをインストールします。Github上にソースコードがあるので、クローンを自身の環境につくります。

git clone https://github.com/yarden/rnaseqlib.git

しかし、現在のバージョンのソースコードではインストールできないので注意してください。こちらも先ほどと同様に、ソースコードを訂正する必要があります。

現行のソースコードでは、setuptoolsをsetup.pyでインポートしていないことによるエラーが発生してしまいます。

問題点の詳細については、下記を参照。
http://stackoverflow.com/questions/12767023/python-packaging

そのため、setup.pyに下記の一文を書き加える必要があります。パッケージをインポートしている序盤の行に書き加えます。

from setuptools import setup

スクリプトを訂正した後、rnaseqlibをインストールします。

# rnaseqlib/
python setup.py install

rMATSのインストール

BiocondaにすでにrMATSは登録されているのですが、現行のレシピでは正常に動作しません。そのため、ソースコードからインストールを行います。
http://rnaseq-mats.sourceforge.net/

まず、ソースコードをダウンロードします。

wget https://sourceforge.net/projects/rnaseq-mats/files/MATS/rMATS.3.2.5.tgz
tar zxvf rMATS.3.2.5.tgz

rMATS.3.2.5ディレクトリにパスを通します。下記のように、ホームディレクトリにあるbashrcにパスを書き加えます。

export PATH=/home/akimitsu/software/rMATS.3.2.5:${PATH}

RNASeq-MATS.pyの1行目を書き換えます。rMATS.3.2.5/RNASeq-MATS.pyの一行目に、下記のようにPythonのパスを書いておくことで、どこからでも実行できるようにします(ここでは、Bioconda2をインストールしているケースを示している)。

#!/home/<ユーザーの名前>/miniconda2/bin/python

1. SRAファイルのダウンロード

NGSデータはシークエンスされたリード(塩基配列)の情報とクオリティスコアのデータをまとめたFASTQ形式のファイルとして保存されています。

FASTQファイルの詳細についてはこちらを参照。
https://en.wikipedia.org/wiki/FASTQ_format

FASTQデータは非常に大きいファイルなので、NCBI GEOでは、FASTQファイルを圧縮したSRAファイルと呼ばれるファイル形式でデータを公開しています。

まず、下記のURLにアクセスすると、
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148

GSEと呼ばれるIDで管理された研究プロジェクトごとのページで、研究の概要や引用されている論文情報を見ることができます。

また、そのプロジェクトで解析したNGSデータがSamplesという項目にリストアップしてあります。

GEO Accession viewer - Google Chrome 2017-02-07 16.27.03.png (236.8 kB)

今回は、35サンプルのうち上6つのサンプルを使用します。

GSMから始まるIDをクリックすると、各サンプルの情報を見ることができます。 サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。

そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。

GEO Accession viewer - Google Chrome 2017-02-07 16.34.35.png (280.0 kB)

SRR4081222をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。

SRR4081222.sraのURLをコピー。

_sra_sra-instant_reads_ByExp_sra_SRX_SRX205_SRX2059672 のインデックス - Google Chrome 2017-02-07 17.13.13.png (61.2 kB)

wgetで該当するSRAファイルをダウンロードします。

$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081227
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081226
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081225
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081224
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081223
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081222

SRRのIDだけだとわかりにくいので、ファイル名をリネームします。

$ mv SRR4081222 SRR4081222_Control_1.sra
$ mv SRR4081223 SRR4081223_Control_2.sra
$ mv SRR4081224 SRR4081224_Control_3.sra
$ mv SRR4081225 SRR4081225_UPF1_knockdown_1.sra
$ mv SRR4081226 SRR4081226_UPF1_knockdown_2.sra
$ mv SRR4081227 SRR4081227_UPF1_knockdown_3.sra

SRAファイルからFASTQファイルに変換

先ほどインストールしたSRA Toolkitを使って、NCBI GEOからダウンロードしたSRAファイルをFASTQファイルに展開します。

$ fastq-dump SRR4081222_Control_1.sra
$ fastq-dump SRR4081223_Control_2.sra
$ fastq-dump SRR4081224_Control_3.sra
$ fastq-dump SRR4081225_UPF1_knockdown_1.sra
$ fastq-dump SRR4081226_UPF1_knockdown_2.sra
$ fastq-dump SRR4081227_UPF1_knockdown_3.sra

2. Quality check

FastQCを使って、FASTQファイルに記載されている各リードのクオリティをチェックします。

FASTQファイルはPhred Quality Scoreと呼ばれる指標によって、各塩基のクオリティを記述しています。

詳しくはこちらのサイトを参照のこと。
https://en.wikipedia.org/wiki/Phred_quality_score

実行例

$ mkdir fastqc_SRR4081222_Control_1
$ fastqc -t 8 -o ./fastqc_SRR4081222_Control_1 \
  ./SRR4081222_Control_1.fastq -f fastq

mkdirコマンドでfastqc_SRR4081222_Control_1という名前のディレクトリをカレントディレクトリに作ります。

次に、fastqcを使ってSRR4081222_Control_1.fastqのデータのクオリティをチェックします。 -tオプションで使用するCPUのコア数を指定します。コア数を増やせば単純に計算速度が向上しますが、FastQCではメモリを1コアあたり250MB消費するので、要求できるメモリ量を超えるコア数を指定するとエラーが起こります。

-oオプションで出力先のディレクトリを指定し、-fオプションでInputファイルのファイル形式を指定しています。

FastQCから出力されるデータが、fastqc_SRR4081222_Control_1ディレクトリの中のSRR4081222_Control_1_fastqc.htmlになります。 ファイルをブラウザで開くと以下のようなグラフを確認できます。

Basic Statistics

基本的な統計量が示されています。

  • 全リード数
  • リード長
  • GC contents

などなど

Per base sequence quality

リードの各塩基のクオリティスコアを示しています。 Phred quality scoreがだいたいグリーンの領域(Scoreが28以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。

ダウンロード.png (8.5 kB)

Per tile sequence quality

フローセルの各タイルごとのクオリティスコアを示しています。

Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。

シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。

特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。

詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf

ダウンロード (1).png (8.6 kB)

Per sequence quality scores

各リードの全長のクオリティスコアの分布を示しています。 ダウンロード (2).png (16.2 kB)

Per base sequence content

各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。

ダウンロード (3).png (33.0 kB)

Per sequence GC content

リードのGC contentsの分布を示しています。 ダウンロード (4).png (32.4 kB)

Per base N content

各塩基中に含まれるNの含有率(塩基を読めなかった箇所)を示しています。 ダウンロード (5).png (7.2 kB)

Sequence Length Distribution

リード長の分布を示しています。 ダウンロード (10).png (20.3 kB)

Sequence Duplication Levels

Duplidate readsの含まれている数を示しています。 ダウンロード (6).png (22.9 kB)

Overrepresented sequences

頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。

Adapter Content

各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。 ダウンロード (7).png (8.8 kB)

Kmer Content

特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。 ダウンロード (8).png (60.2 kB)

3. Mapping

注意!!

マッピングするリード長は必ず同じ長さになるようにしてください。今回使用するAlternative splicingの解析を行うソフトは、異なるリード長が混在するデータを受け付けてくれません。リードのトリミングは避けてください。

TopHatを用いてGenomeとTranscriptomeにリードをマッピングを行います。

注意!!

RNA-seqのデータ解析については近年、マッピングソフトとしてSTARがデファクトスタンダードになっています(あくまで個人的な意見です。ENDODEプロジェクトで使われて一気にスタンダードになった気がします)。

確かに、STARのほうがTopHatと比較してマッピングの精度が向上しているので、現時点ではSTARがファーストチョイスかと思います。

しかし、STARはメモリを大量に消費し、ヒトゲノムだと32GB以上は必要となります。メモリ不足が懸念される解析環境では、メモリ消費が少ないTopHatを選択するのが無難だと思います(昔はこれがデファクトだったので)。

Indexファイルの作成

TopHatの内部で使用されているBowtieでゲノムへマッピングするために、ヒトゲノムに対するIndexファイルを用意する必要があります。

Bowtieのサイトから、Indexファイルはダウンロードすることもできます。右端のカラムに各ゲノムのIndexファイルが羅列されているので、任意のIndexファイルをダウンロードしてきます。 http://bowtie-bio.sourceforge.net/index.shtml

もしくは、自身でIndexファイルを作成することも可能です。以下では、その方法を説明します。

まず、Human Genome(hg19)のFASTAファイルのダウンロードしてきます。 次に、catコマンドで各染色体ごとのFASTAファイルをマージして、hg19.faという名前で出力します。

$ for file in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X Y M
  do
    wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/chr${file}.fa.gz
    gunzip chr${file}.fa.gz
  done

$ cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chr20.fa chr21.fa chr22.fa chrX.fa chrY.fa chrM.fa > hg19.fa

最後に、hg19.faファイルからBowtie用のIndexファイルを作成します。

$ bowtie-build ./hg19.fa hg19
  • 1つ目の引数: リファレンスゲノムのFASTAファイルを指定。
  • 2つ目の引数: Indexファイルの名前を指定します。ここでは、Indexの名前をhg19とします。

アノテーション情報(GTFファイル)の準備

TopHatでは、GTF形式のファイルのTranscriptomeのアノテーション情報を使って、Splice junctionサイトに対してもマッピングを行うことができます。

今回はTranscriptomeの情報として、Gencodeのアノテーション情報を利用します。
https://www.gencodegenes.org/releases/19.html

上記のサイトにアクセスすると、さまざまなデータが置いてあるが基本的にgencode.v19.annotation.gtfというメインのアノテーション情報を利用すればOKです。
ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

wgetコマンドで該当のファイルをダウンロードし解凍します。

$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
$ gunzip ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

カレントディレクトリにgencode.v19.annotation.gtfが保存されます。

GenomeとTranscriptomeへのマッピング

ここまでで、ゲノムのIndexファイルと、アノテーション情報の入ったGTFファイルが用意できたので、いよいよTopHatを用いて各シークエンスリードをゲノムとトランスクリプトームにマッピングします。

以下のコマンドを各サンプル(FASTQファイル)に対して行います。注意すべき点は、--library-typeオプションでサンプルがStrand-specificであることを明示する点です。
https://ccb.jhu.edu/software/tophat/manual.shtml

実行例

$ tophat --bowtie1 -p 8 --library-type fr-firststrand -o ./tophat_out_SRR4081222_Control_1_ss \
  -G ./gencode.v19.annotation.gtf \
  ./hg19 ./SRR4081222_Control_1_1_filtered.fastq
  • --bowtie1: マッピングソフトにBowtie1を使用。
  • -p: 使用するコア数を指定。
  • --library-type fr-firststrand: dUTP法を利用したStrand-specific RNA-seqサンプルに対して、この設定を使用する。
  • -o: 出力先のディレクトリ名を指定(ディレクトリは作成していなくても自動的に作成してくれる)。
  • -G: アノテーション情報として、GTFファイルを指定。
  • 1つ目の引数: 先ほど作成したヒトゲノム(hg19)のBowtie1のIndexファイルを指定。
  • 2つ目の引数: マッピングしたいFASTQファイルを指定。

TopHatによるマッピングが終わると、

  • align_summary.txt
  • accepted_hits.bam

というファイルが出力されます。

align_summary.txtマッピング結果をまとめたテキストファイルで、中身を見てみると、

Reads:
          Input     :  13282516
           Mapped   :  12196063 (91.8% of input)
            of these:   1887950 (15.5%) have multiple alignments (985 have >20)
91.8% overall read mapping rate.

上記のように、マッピングされたリードの数・割合や、Multi-hitしたリードの数・割合などが記載されています。

accepted_hits.bamマッピングされたリードの情報、具体的には、各リードがマッピングされた染色体座標軸、アライメントスコア、リファレンスゲノムと比較したときのIndelやMutationの有無などの情報が記載されています。

ただし、BAMファイルはバイナリファイルなので、ファイルを開いて直接中身を見ることができません。

ファイルの中身を確認したい場合はsamtoolsを用いることで、SAM形式に変換してファイルの中身を確認することができます。

$ samtools view ./tophat_out_SRR4081222_Control_1/accepted_hits.bam | less

SAMファイルの中身。

HWI-ST1075L:319:C42GNACXX:1:1315:16479:46448    256     chr1    12137   0       36M     *       0       0       CCTGCATGTAACTTAATACCACAACCAGGCATAGGG    @@@FFFFFHHHHHJJJJJJIJJJJJGIGEHIGGIEG    XA:i:1  MD:Z:0T35       NM:i:1  XS:A:+  NH:i:6  CC:Z:chr15      CP:i:102518998  HI:i:0
HWI-ST1075L:319:C42GNACXX:1:1214:15035:37072    272     chr1    12333   0       36M     *       0       0       GGCTGTGACTGCTCAGACCAGCCGGCTGGAGGGAGG    JIJJIIIIIJJIGIIIIJJJJJIHHHGHFFFFFCCC    XA:i:0  MD:Z:36 NM:i:0  XS:A:+  NH:i:6  CC:Z:chr15      CP:i:102518802  HI:i:0
HWI-ST1075L:319:C42GNACXX:1:1301:7343:9547      272     chr1    12693   1       36M     *       0       0       CTCTTGTGAGTGTCCCCAGTGTTGCAGAGGTGAGAG    JJIHIJJJJJJIGJJJJJJJJJJHHHHHFFFFFCCC    XA:i:0  MD:Z:36 NM:i:0  XS:A:+  NH:i:4  CC:Z:chr15      CP:i:102518442  HI:i:0
HWI-ST1075L:319:C42GNACXX:1:2204:8186:34501     16      chr1    12693   1       36M     *       0       0       CTCTTGTGAGTGTCCCCAGTGTTGCAGAGGTGAGAG    JJJJJJJJIJJJIIJJJJJJIJJHHHHHFFFFFCCC    XA:i:0  MD:Z:36 NM:i:0  XS:A:+  NH:i:4  CC:Z:chr15      CP:i:102518442  HI:i:0
HWI-ST1075L:319:C42GNACXX:1:1305:5797:18762     272     chr1    12957   0       36M     *       0       0       CATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCAGG    DEHIGGIFF?2F;GGFAIIHIGGHDFHDDFDFD@@@    XA:i:0  MD:Z:36 NM:i:0  XS:A:+  NH:i:6  CC:Z:chr12      CP:i:92621      HI:i:0
HWI-ST1075L:319:C42GNACXX:1:2213:9952:39080     272     chr1    13024   0       36M     *       0       0       GTGCATGAAGGCTGTCAACCAGTCCATAGGCAAGCC    JJJJJJJJJJIIJIIHIGJJJJJHHGHHFFFFFCCC    XA:i:0  MD:Z:36 NM:i:0  XS:A:+  NH:i:7  CC:Z:chr12      CP:i:92554      HI:i:0

BAM/SAMファイルについての詳細を知りたい場合は、下記のURLを参照のこと。
https://samtools.github.io/hts-specs/
https://samtools.github.io/hts-specs/SAMv1.pdf
https://samtools.github.io/hts-specs/SAMtags.pdf

4. データの品質チェック

RNA-seq用のRNAサンプルはバイオアナライザーなどの機械であらかじめRNAの品質をチェックしてからcDNAライブラリを作成しています。

ただ、古いデータの中にはRNAの品質に問題があるものもあり、遺伝子領域を見たときに、例えば、5'側から3'側でのCoverageがガタガタするケースがあります。

その他にも、PolyAセレクションを行ったサンプルの場合、RNAの品質が悪いと3'側に極端にリードが偏る傾向にあることが知られています。

そこで、マッピングされたリードから遺伝子領域上のシークエンスリードのCoverageを調べることによって、大まかにRNAの品質も問題がないか確認を行います。

今回は、RSeQCスクリプトの1つであるgeneBody_coverage.pyを用いて、遺伝子領域上でのCoverageを調べてみようと思います。

このスクリプトの実行には、BAMファイルとそのIndexファイル、BED形式で記述された遺伝子領域の情報(アノテーション情報)が必要になります。

まず、アノテーション情報に関しては、基本的にどの細胞でも発現しているHouse-keeping genesのみに絞って解析したほうがよいです。

そこで、RSeQCの開発元が配布しているHouse-keeping genesのBEDファイルのリストを入手してきます。

$ wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19.HouseKeepingGenes.bed

次に、結果の保存先のディレクトリの用意と、BAMファイルのIndexファイルを作成しておきます。

$ mkdir geneBody_coverage_SRR4081222_Control_1
$ samtools index ./tophat_out_SRR4081222_Control_1/accepted_hits.bam

mkdirコマンドで新しいディレクトリを作成。samtools indexでBAMファイルのIndexファイルを作成します。引数にはBAMファイルを指定するだけです。

それでは、House-keeping genesの遺伝子領域のCoverageを調べてみたいと思います。

実行例

$ geneBody_coverage.py -r ./hg19.HouseKeepingGenes.bed \
  -i ./tophat_out_SRR4081222_Control_1/accepted_hits.bam  \
  -o ./geneBody_coverage_SRR4081222_Control_1/SRR4081222_Control_1_RSeQC_output
  • -r: 先ほどダウンロードしたHouse-keeping genesの情報が格納されているBEDファイルを指定。
  • -i: Coverageを調べたいBAMファイルを指定(加えて、Indexファイルを同じディレクトリにおいておく必要がある)。
  • -o: 出力先のディレクトリとファイル名のPrefix。

結果として、指定した名前.geneBodyCoverage.curves.pdfというファイルが出力されます。ファイルの中身を見るとこんな感じ。

image.png (73.1 kB)

RNAの品質に問題がある(例えば、RNAが壊れている)と、3'側にリードが偏っていたり、cDNAライブラリの構築に使ったRNA量が少なく、PCR時に特定の配列が異常に増幅され(PCR bias、Duplicate readsが多い)、ガタガタしたCoverageの分布を取ることがあります。

こういったデータはRNA定量性に問題があるため、使用しないように気をつけましょう。

このステップで問題がありそうだと判断されるデータに関しては、下記のUCSC genome browserでの可視化により、個々の遺伝子上にマッピングされたリードの分布に注意して見ていく必要があります。

5. マッピングされたリードのStrandnessの確認

RNA-seqの各シークエンスリードがStrand-specificに読まれているかどうか(リードのStrandness)を確認します。ここでは、RSeQCのinfer_experiment.pyというスクリプトを用いました。

アノテーション情報としてBEDファイルを用意

下記のURLに置いてあるgtf2bed.plを使用します。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

使用するスクリプトアノテーション情報としてBEDファイルしか受け付けないので、マッピングの際に使用したGTFファイルgencode.v19.annotation.gtfをここで一旦BEDファイルに変換します。

$ perl gtf2bed.pl gencode.v19.annotation.gtf gencode.v19.annotation.bed
  • 1つ目の引数: 変換したいGTFファイルを指定
  • 2つ目の引数: 出力するBEDファイルを指定。

Strandnessの確認

マッピングしたデータであるBAMファイルと、アノテーション情報であるBEDファイルを指定して実行します。

$ infer_experiment.py -i ./tophat_out_SRR4081222_Control_1_ss/accepted_hits.bam \
  -r ./gencode.v19.annotation.bed -s 400000 \
  >  ./tophat_out_SRR4081222_Control_1_ss/SRR4081222_Control_1_strand_stat.txt
  • -i: Strandnessを調べるBAMファイルを指定。
  • -r: 先ほど用意したBEDファイルのアノテーション情報を指定。
  • s: 解析に使用するリードを指定。ここでは、40万リードを抽出してStrandnessを調べています。指定する数値を大きくすれば計算精度も上がると思いますが、その分時間もかかります。

下記がその結果になります。"++“や”+-“のように2文字で表現されていますが、1文字目がアノテーション情報でのStrandの向き、2文字目がマッピングされたリードのStrandの向きを表しています。

unstranded RNA-seqの場合、Fraction of reads explained by “++,–"とFraction of reads explained by ”+-,-+“の割合が、50%に近くなりますが、今回のようにStrand-specificに読まれているサンプルの場合、下記のように偏りが見られます。

dUTP法でシークエンスされたRNA-seqデータでは、RNAに対してAntisense鎖が読まれます。そのため、+鎖のRNA配列は-鎖として、-鎖のRNA配列は+鎖としてそれぞれ読まれているはずです。

This is SingleEnd Data
Fraction of reads failed to determine: 0.1910
Fraction of reads explained by "++,--": 0.1045
Fraction of reads explained by "+-,-+": 0.7045

例えば上記のように、予想したとおり+-,-+の組み合わせのリードが70%程度大半を占めていることがわかります。

次に、この情報をもとに、UCSC genome browserで可視化させるデータを作成する際に、+と-鎖をそれぞれ分けてファイルを作成します。

6. Visualization (For UCSC genome browser)

BedGraphファイルの準備

Bedtoolsを利用して、TopHatから得られたBAMファイルを、UCSC genome browser上で可視化できるファイル(Bedgraphファイル)に変換します。

まず、mkdirコマンドでデータを保存するディレクトリを作っておきます。

$ mkdir UCSC_visual_SRR4081222_Control_1_ss

bedtools genomecovで特定のゲノム領域でのCoverageを計算する(BAMファイルからBedgraphファイルを作成)。

Bedgraphファイルの詳細については、下記のURLを参照のこと。
https://genome.ucsc.edu/goldenpath/help/bedgraph.html

$ bedtools genomecov -ibam ./SRR4081222_Control_1_ss/accepted_hits.bam \
  -bg -split -strand - \
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_Forward.bg
$ bedtools genomecov -ibam ./SRR4081222_Control_1_ss/accepted_hits.bam \
  -bg -split -strand + \
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_Reverse.bg
  • -ibam: BAMファイルを入力ファイルとして指定する。
  • -bg: 出力ファイルの形式をBedGraphファイルに指定。
  • -split: Splicing juctionにマッピングされたリード(リードが分割されてマッピングされているリード)を考慮する。
  • -strand: +もしくは-を指定することで、ゲノムに対して+鎖もしくは-鎖にあたるリードを選択的に抽出することができる。

標準出力でBedGraphファイルが出力されます。

UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。それらの情報をヘッダー行(ファイルの一行目)に記載します。

$ echo "track type=bedGraph name=SRR4081222_Control_1_Fw description=SRR4081222_Control_1_Fw visibility=2 maxHeightPixels=40:40:20 color=255,0,0" \
  > ./UCSC_visual_SRR4081222_Control_1_ss/tmp_SRR4081222_Control_1_forward.txt
$ echo "track type=bedGraph name=SRR4081222_Control_1_Re description=SRR4081222_Control_1_Re visibility=2 maxHeightPixels=40:40:20 color=0,0,255" \
  > ./UCSC_visual_SRR4081222_Control_1_ss/tmp_SRR4081222_Control_1_reverse.txt

echoコマンドを使って、ヘッダー行の情報をtxtファイルとして出力しています。

そのあとに、catコマンドを用いて、先ほど作ったBedGraphファイルと結合させます。

$ cat ./UCSC_visual_SRR4081222_Control_1_ss/tmp_SRR4081222_Control_1_forward.txt ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_Forward.bg \
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_forward_for_UCSC.bg
$ cat ./UCSC_visual_SRR4081222_Control_1_ss/tmp_SRR4081222_Control_1_reverse.txt ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_Reverse.bg 
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_reverse_for_UCSC.bg

出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2ファイルに圧縮しています。echoコマンドのところでヘッダー行の情報をtxtファイルとして出力しています。

$ bzip2 -c ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_forward_for_UCSC.bg \
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_forward_for_UCSC.bg.bz2
$ bzip2 -c ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_reverse_for_UCSC.bg \
  > ./UCSC_visual_SRR4081222_Control_1_ss/SRR4081222_Control_1_reverse_for_UCSC.bg.bz2

データのアップロード

UCSC genome browserにデータをアップロードします。

ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。

Add Custom Tracksというページに移動するので、ファイルを選択をクリックし、アップロードするファイルを選択し、submitボタンをクリックします。

アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgoボタンをクリックすると、データをみることができます。

20170205212758.png (1.7 MB)

7. MISOの実行(所要時間: 1hr)

MISO用のアノテーション情報のIndexファイル作成

MISOの実行には、MISO用のアノテーション情報(Alternative splicingのパターン)を用意する必要があります。使用するアノテーション情報が決まっており、UCSCのknownGene、RefSeqのrefGene、GencodeのensGeneの3つを使用します。

まず、必要なファイルのダウンロードし、解凍します。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/knownGene.txt.gz
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/ensGene.txt.gz
$ gunzip *.gz

GFFファイルの作成

skipped exons (SE), mutually exclusive exons (MXE), alternative 3' splice sites (A3SS) , alternative 5' splice sites (A5SS), retained introns (RI)のそれぞれのパターンをGFFファイルにまとめます。

gff_make_annotation ./ ./gff --flanking-rule commonshortest --genome-label hg19_Gencode_v19
  • --genome-labelでGFFファイルの名称を指定。

Indexファイルの作成

GFFファイルを作成したら、それぞれのパターンに対して、Indexファイルを作成します。

index_gff --index A3SS.hg19_Gencode_v19.gff3 indexed_A3SS_events/
index_gff --index A5SS.hg19_Gencode_v19.gff3 indexed_A5SS_events/
index_gff --index MXE.hg19_Gencode_v19.gff3 indexed_MXE_events/
index_gff --index SE.hg19_Gencode_v19.gff3 indexed_SE_events/
index_gff --index RI.hg19_Gencode_v19.gff3 indexed_RI_events/

MISOの設定ファイルの準備

miso_settings.txtというファイル名のMISOの設定ファイルを作成します。 以下のような内容を記述しておきます。

[data]
filter_results = True
min_event_reads = 20
strand = fr-firststrand

[sampler]
burn_in = 500
lag = 10
num_iters = 5000
num_processors = 8
  • filter_results: 低Coverageだった場合、解析対象から外すかどうか
  • min_event_reads: 最低リード数を指定。
  • strand: Strandの情報について指定。
  • num_processors: 使用するコア数を指定。

strandに関しては以下のサイトを参照のこと。
http://onetipperday.sterding.com/2012/07/how-to-tell-which-library-type-to-use.html

MISOを実行する

5つそれぞれのパターン(A3SS, A5SS, MXE, RI, SE)についてMISOで解析を行います。

$ for type in A3SS A5SS MXE RI SE
  do
      gff_index="./MISO_index/gff/commonshortest/indexed_${type}_events"
      outputDir="MISO_output_SRR4081222_Control_1_${type}_events/"

      miso --run ${gff_index} ./tophat_out_SRR4081222_Control_1/accepted_hits.bam --output-dir ${outputDir} \
      --read-len 101 --settings-filename ./miso_settings.txt

      summarize_miso --summarize-samples ${outputDir} ${outputDir}
  done
miso

2群間比較する前に、各サンプルについてSplicingのパターンごとにリードカウントを行います。

  • --runで先ほど用意したGFFのIndexファイルを指定します。
  • --output-dirで保存先のディレクトリを指定。
  • --read-lenでリード長を指定。
  • --settings-filenameで先ほど用意したMISOの設定ファイルを指定。
  • 1つめの引数として、解析するBAMファイルを指定。
miso summary

misoでリードカウントした結果をまとめたファイルを出力します(解析に直接使うことはない)。

  • --summarize-samplesmisoで保存先として指定したディレクトリを指定。
  • 1つめの引数として、保存先のディレクトリを指定。
miso_compare

Alternative splicingのパターンの変化を2群間で比較します。

$ for type in A3SS A5SS MXE RI SE
  do
      outputDir="./MISO_comparison_${type}_events"
      MISOFile1="./MISO_output_SRR4081222_Control_1_${type}_events"
      MISOFile2="./MISO_output_SRR4081225_UPF1_knockdown_1_${type}_events"

      compare_miso --compare-samples ${MISOFile1} ${MISOFile2} ${outputDir}
  done
  • --compare-samplesで比較したい2群のデータを指定します(スペース区切り)。このとき指定するのは、misoで保存先に指定したディレクトリです。
  • 1つめの引数として、保存先のディレクトリを指定。

保存先のディレクトリに結果のファイルが出力されます。 ファイルの中身は以下のような感じ。

Book1 - Excel 2017-03-28 22.35.39.png (38.0 kB)

基本的に、diffbayes_factorの値に着目すればOKです。 詳しくは以下のサイトを参照のこと。
http://miso.readthedocs.io/en/fastmiso/#interpreting-and-filtering-miso-output

diffはΔ Ψを意味しています。bayes_factorはその値が高いほど有意な差が2群間であると判断できます。

8. DEXSeqの実行(所要時間: 5hr)

アノテーション情報の準備

必要なスクリプトGithubから入手します。

git clone https://github.com/vivekbhr/Subread_to_DEXSeq.git

Subread_to_DEXSeqディレクトリにパスを通します。下記のように、ホームディレクトリにあるbashrcにパスを書き加えます。

export PATH=/home/akimitsu/software/Subread_to_DEXSeq:${PATH}

次に、アノテーション情報を用意します。

dexseq_prepare_annotation2.py -f gencode.v19.annotation_filtered_DEXSeq.gtf gencode.v19.annotation_filtered.gtf gencode.v19.annotation_filtered_DEXSeq.gff
  • -fオプションでfeatureCountsで集計する用のGTFファイルの出力先を指定。
  • 1つ目の引数で、アノテーションに使用するGTFファイルを指定。
  • 2つ目の引数で、DEXSeq用のGTFファイルの出力先を指定。

featureCountsで集計する

$ mkdir DEXSeq_output_test
$ featureCounts -T 8 -f -O -s 2 \
  -F GTF -t exon -a ./gencode.v19.annotation_filtered_DEXSeq.gtf \
  -o DEXSeq_output_test \
  ./tophat_out_SRR4081222_Control_1/accepted_hits.bam \
  ./tophat_out_SRR4081223_Control_2/accepted_hits.bam \
  ./tophat_out_SRR4081224_Control_3/accepted_hits.bam \
  ./tophat_out_SRR4081225_UPF1_knockdown_1/accepted_hits.bam \
  ./tophat_out_SRR4081226_UPF1_knockdown_2/accepted_hits.bam \
  ./tophat_out_SRR4081227_UPF1_knockdown_3/accepted_hits.bam

mkdirコマンドを使って、保存先のディレクトリを作り、そのあとfeatureCountsを実行します。

  • -T: 使用するコア数の指定。
  • -f: 特徴(Exonなど)ごとにリードをカウントする。
  • -O: 複数の領域にOverlapしているリードもカウントする。
  • -s 2: Transcriptの配列に対してAntisense鎖にあたるリードをカウントする(dUTP法でStrand-specific RNA-seqのデータが得られているため)。
  • -F: GTF/GFFファイルをInputファイルとしていることを明示。
  • -t exon: リードのカウントに使用するExonの情報を指定。ここでは、GTFファイル内で3列目のfeatureがexonである行を使用する用に指定している。
  • -a: アノテーション情報(GTFファイル)を指定。
  • -o: 出力先のパスとファイル名を指定。
  • 1つ目の引数: BAMファイルを指定。

DEXSeqでExonごとのUsageを調べる

内部的には、DESeq2の方法でDispersionを求めてGeneごとでなく、Exonごとに有意差検定を行っています(MISOと計算方法が異なる)。

Rで以下のコードを実行します(もしくはスクリプトにまとめて実行)。

source("~/software/Subread_to_DEXSeq/load_SubreadOutput.R")
suppressPackageStartupMessages({
  require(dplyr)
})

# Load featureCounts result file
samp <- data.frame(row.names = c("cont_1","cont_2","test_1","test_2"), 
                   condition = c("Control","Control","Treatment","Treatment"))
dxd.fc <- DEXSeqDataSetFromFeatureCounts("featureCounts_result_test.txt",
                                         flattenedfile = "/home/akimitsu/database/gencode.v19.annotation_filtered_DEXSeq.gtf",
                                         sampleData = samp)

# Adjust for coverage biases among samples
dxd <- estimateSizeFactors(dxd.fc)

# Estimate variance or dispersion parameters individually exon by exon
dxd <- estimateDispersions(dxd)

png(filename = "test.png", width = 600, height = 600)
plotDispEsts(dxd)
dev.off()

# Test for differential exon usage in each gene
dxd <- testForDEU(dxd)

# Estimate relative exon usage fold changes
dxd <- estimateExonFoldChanges(dxd, fitExpToVar="condition")

# Output DEXSeq result
dxr1 <- DEXSeqResults(dxd)
write.table(dxr1, quote = F, sep = "\t", file = "DEXSeq_test_result.txt")
  • sourceGithubからクローンを作成したSubread_to_DEXSeqディレクトリの直下にあるload_SubreadOutput.Rを呼び出す。
  • data.framerow.namesにBAMファイルのラベル名を、conditionにBAMファイルの属性(Control or Knockdownなど)を代入したデータフレームを作成。
  • DEXSeqDataSetFromFeatureCountsでfeatureCountsから出力されたファイルを指定し、flattenedfileにDEXSeq用に用意したアノテーション情報に関するGTFファイルを代入する。
  • write.tablefileに保存先のファイル名を指定する。

保存先のディレクトリに結果のファイルが出力されます。 ファイルの中身は以下のような感じ。

Book1 - Excel 2017-03-28 23.12.34.png (51.0 kB)

多重検定補正を行ったp-value値であるpadjと、発現量の変化を示したlog2fold_Treatment_Controlに着目すればよい。

9. rMATSの実行(所要時間: 1.5 days)

RNASeq-MATS.py \
-b1 ./tophat_out_SRR4081222_Control_1/accepted_hits.bam,./tophat_out_SRR4081223_Control_2/accepted_hits.bam,./tophat_out_SRR4081224_Control_3/accepted_hits.bam \
-b2 ./tophat_out_SRR4081225_UPF1_knockdown_1/accepted_hits.bam,./tophat_out_SRR4081226_UPF1_knockdown_2/accepted_hits.bam,./tophat_out_SRR4081227_UPF1_knockdown_3/accepted_hits.bam \
-gtf /home/akimitsu/database/gencode.v19.annotation_filtered.gtf \
-o rMATS_test -t single -len 101 -c 0.0001 -analysis U -libType fr-firststrand -novelSS 1
  • -b1: ControlのBAMファイルをカンマ(,)区切りで指定。
  • -b2: KnockdownのBAMファイルをカンマ(,)区切りで指定。
  • -gtf: Gencodeなどのアノテーション情報に関するGTFファイルを指定。
  • -o: 保存先のディレクトリを指定。
  • -t: Single-endかPaired-endか明示する。
  • -len: リード長を指定。
  • -c: Splicingパターンの変化のcutoff値を指定。0.0001の場合、Transcript全体から見て0.01%以上の変化を検出する。
  • -analysis: paired analysisかunpaired analysisか指定する。
  • -libType: Strandの向きの指定(MISOのときと同じ)。
  • -novelSS: Noval splicing sitesを検出するかどうか。1でYes、0でNo。

オプションの詳細はこちらのサイトを参照のこと。
http://rnaseq-mats.sourceforge.net/user_guide.htm

保存先のディレクトリに結果のファイルが出力されます。 ファイルの中身は以下のような感じ。

AS_Event.MATS.JunctionCountOnly.txt

Junctionをまたぐリードだけを集計。

AS_Event.MATS.ReadsOnTargetAndJunctionCounts.txt

Junctionリードを含むExon上にマッピングされたすべてのリードを集計。

  • IncFormLen: length of inclusion form, used for normalization
  • SkipFormLen: length of skipping form, used for normalization
  • IncLevel1: inclusion level for SAMPLE_1 replicates (comma separated) calculated from normalized counts
  • IncLevel2: inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized counts
  • IncLevelDifference: average(IncLevel1) - average(IncLevel2)

以上になります。