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
使用するソフトウェア
- FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - Bowtie1 (v1.1.2)
http://bowtie-bio.sourceforge.net/index.shtml - TopHat (v2.1.1)
https://ccb.jhu.edu/software/tophat/index.shtml - Samtools (v1.3.1)
http://samtools.sourceforge.net/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - bedtools (v2.26.0)
http://bedtools.readthedocs.io/en/latest/ - featureCounts (subread v1.5.0.post3)
http://subread.sourceforge.net/ - MISO (v0.5.3)
http://miso.readthedocs.io/en/fastmiso/ - DEXSeq (v1.20.1)
http://bioconductor.org/packages/release/bioc/html/DEXSeq.html - rMATS (v3.2.5)
http://rnaseq-mats.sourceforge.net/
解析ワークフロー
Alternative splicingを解析する方法として、MISO, DEXSeq, rMATSそれぞれを用いた方法を紹介したいと思います。
TopHat-MISO workflow
TopHat-featureCounts-DEXSeq workflow
TopHat-rMATS workflow
各ソフトウェアの特徴(主観的な意見です…)
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
改行コードの問題
WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。
ここではAtom
と呼ばれるテキストエディタでの設定方法を説明する。
- [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
- 左のメニューから、
Packages
をクリック。 - Installed Packagesで
line-ending-selector
と検索。 - Core Packagesの項目から、line-ending-selectorの
Settings
をクリック。 - 設定画面で、Default line endingを
LF
に変更。
これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。
各ステップの説明
0. 必要なソフトウェアのインストール
まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。
基本的なソフトのインストール
$ 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_indx
とupper_bound_indx
がint型でなく、float型であることに起因したエラーだと考えられます(Indexとして整数でない数値を指定してしまっている)。
そのため、エラーを回避するために、途中でlower_bound_indx
とupper_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
という項目にリストアップしてあります。
今回は、35サンプルのうち上6つのサンプルを使用します。
GSMから始まるIDをクリックすると、各サンプルの情報を見ることができます。 サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。
そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。
SRR4081222をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。
SRR4081222.sra
のURLをコピー。
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以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。
Per tile sequence quality
フローセルの各タイルごとのクオリティスコアを示しています。
Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。
シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。
特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。
詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf
Per sequence quality scores
各リードの全長のクオリティスコアの分布を示しています。
Per base sequence content
各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。
Per sequence GC content
リードのGC contentsの分布を示しています。
Per base N content
各塩基中に含まれるN
の含有率(塩基を読めなかった箇所)を示しています。
Sequence Length Distribution
リード長の分布を示しています。
Sequence Duplication Levels
Duplidate readsの含まれている数を示しています。
Overrepresented sequences
頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。
Adapter Content
各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。
Kmer Content
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
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
というファイルが出力されます。ファイルの中身を見るとこんな感じ。
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ボタンをクリックすると、データをみることができます。
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
でリードカウントした結果をまとめたファイルを出力します(解析に直接使うことはない)。
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つめの引数として、保存先のディレクトリを指定。
保存先のディレクトリに結果のファイルが出力されます。 ファイルの中身は以下のような感じ。
基本的に、diff
とbayes_factor
の値に着目すればOKです。
詳しくは以下のサイトを参照のこと。
http://miso.readthedocs.io/en/fastmiso/#interpreting-and-filtering-miso-output
diff
はΔ Ψを意味しています。bayes_factor
はその値が高いほど有意な差が2群間であると判断できます。
8. DEXSeqの実行(所要時間: 5hr)
アノテーション情報の準備
- Subread_to_DEXSeq
https://github.com/vivekbhr/Subread_to_DEXSeq
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")
source
でGithubからクローンを作成したSubread_to_DEXSeq
ディレクトリの直下にあるload_SubreadOutput.R
を呼び出す。data.frame
でrow.names
にBAMファイルのラベル名を、condition
にBAMファイルの属性(Control or Knockdownなど)を代入したデータフレームを作成。DEXSeqDataSetFromFeatureCounts
でfeatureCountsから出力されたファイルを指定し、flattenedfile
にDEXSeq用に用意したアノテーション情報に関するGTFファイルを代入する。write.table
でfile
に保存先のファイル名を指定する。
保存先のディレクトリに結果のファイルが出力されます。 ファイルの中身は以下のような感じ。
多重検定補正を行った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 normalizationSkipFormLen
: length of skipping form, used for normalizationIncLevel1
: inclusion level for SAMPLE_1 replicates (comma separated) calculated from normalized countsIncLevel2
: inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized countsIncLevelDifference
: average(IncLevel1) - average(IncLevel2)
以上になります。
BRIC-seqデータ解析(2群間比較, Unstranded, Single-end)
概要
BRIC-seqはゲノムワイドにmRNAとlncRNAのRNA分解速度を測定する方法です。
ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、BRIC-seqのデータ解析の詳細をStep by Stepで説明します。
ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。パソコン(もしくはスパコン)のメモリは16GB以上を確保することが望ましいです。
このTutorialを通した到達目標
- FastQCから得られたシークエンスリードのクオリティデータから、FASTQファイルの状態を確認できる。
- FASTQファイルをゲノムとトランスクリプトームにマッピングすることができる。
- BAMファイルからデータの品質チェックを行うことができる。
- UCSC genome browser上でデータを可視化することができる。
- RNAの半減期を算出することができる。
- 各遺伝子に対するRNA分解曲線(Fitting curve)を描くことができる。
- Webブラウザを介して、各遺伝子に対するRNA分解曲線のデータを共有することができる。
使用するデータ
BRIC-seq (DRA001215: Control siRNA, STAU1 siRNA)
https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA001215
使用するソフトウェア
- FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - Bowtie1 (v1.1.2)
http://bowtie-bio.sourceforge.net/index.shtml - TopHat (v2.1.1)
https://ccb.jhu.edu/software/tophat/index.shtml - Samtools (v1.3.1)
http://samtools.sourceforge.net/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - bedtools (v2.26.0)
http://bedtools.readthedocs.io/en/latest/ - Cufflinks (v2.2.1)
http://cole-trapnell-lab.github.io/cufflinks/ - BridgeR2 (v0.1.0)
https://github.com/Imamachi-n/BridgeR2
解析ワークフロー
BRIC-seqデータからRNA半減期を算出し、さらに2群間の半減期を比較する方法を紹介したいと思います。
サンプルコード
Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/BRIC-seq_Tutorial
BRIC-seq_unstr_single-end_1_mapping.sh
TopHatによるマッピング・featureCountsによるリードカウントまでを行うスクリプト。FASTQファイルを引数として指定。
$ ./BRIC-seq_unstr_single-end_1_mapping.sh hoge.fastq
11行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
12行目: indexContamFile
にrRNAなどのコンタミ配列のBowtie用のIndexファイルを指定。
13行目: indexGenomeFile
にGenome配列のBowtie用のIndexファイルを指定。
35行目: -r
オプションでHouse-keeping genesの遺伝子領域のBEDファイルを指定。
BRIC-seq_unstr_single-end_2_RPKM_calc_Control.sh/BRIC-seq_unstr_single-end_2_RPKM_calc_Knockdown.sh
Cuffnormによるリードカウントを行うスクリプト。
$ ./BRIC-seq_unstr_single-end_2_RPKM_calc_Control.sh
8行目: filename
に保存先のディレクトリ名を指定。
9行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
13-18行目: スペース()区切りでBAMファイルを時系列順に指定。
21行目: gene_list
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
BRIC-seq_unstr_single-end_3_RNA_half-life_calc.sh(BridgeR_analysis_mRNA.R/BridgeR_analysis_lncRNA.R)
BridgeR2によるRNA半減期の算出、および2群間比較を行うスクリプト。
$ ./BRIC-seq_unstr_single-end_3_RNA_half-life_calc.sh
3,4行目: 算出したmRNAとlncRNAのRNA半減期のリストを保存するディレクトリ名を指定。
9,14行目: ControlとKnockdownのリードカウント(集計)データを指定(カンマ(,
)区切り)。
BridgeR_analysis_mRNA.R/BridgeR_analysis_lncRNA.Rに関しては、
8行目: group
変数に2群のラベル名を代入。
9行目: hour
変数にタイムコースを代入。
35行目: comparisonFile
変数に2群のラベル名を代入。
37行目: “BridgeR_6_PUM1KD"を任意の名前に変更。
BridgeR_shinyapp.R
6行目: BridgeR2から得られたRNA半減期のリストを指定。
11行目: hour
変数にタイムコースを代入。
~/custom_command/<hoge.py>
となっている箇所はカスタムのPythonスクリプトが置いてあるパス・ファイル名を指定する。
シェルスクリプトを使用するときの注意点
権限の問題
実行権限をシェルスクリプトに与えた後、実行すること。
$ chmod 755 <hoge.sh>
改行コードの問題
WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。
ここではAtom
と呼ばれるテキストエディタでの設定方法を説明する。
- [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
- 左のメニューから、
Packages
をクリック。 - Installed Packagesで
line-ending-selector
と検索。 - Core Packagesの項目から、line-ending-selectorの
Settings
をクリック。 - 設定画面で、Default line endingを
LF
に変更。
これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。
各ステップの説明
0. 必要なソフトウェアのインストール
まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。
$ conda install fastqc $ conda install fastx_toolkit $ conda install bowtie=1.12 $ conda install tophat $ conda install samtools $ conda install bedtools $ conda install cufflinks
tophatによるマッピングで最新のbowtieを使うとエラーが出るので、以前のバージョンのbowtie v1.12を代わりにインストールする。
BRIC-seq解析用のRパッケージbridger2
をインストールします。
$ R > Sys.setenv(CURL_CA_BUNDLE = "/home/akimitsu/miniconda2/ssl/cacert.pem") > install.packages('bridger2', repos='http://cran.us.r-project.org')
1. SRAファイルのダウンロード
DDBJ DRAに登録されているBRIC-seqのSRAファイルをダウンロードし、FASTQファイルに変換します。
$ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012976/DRR014456/DRR014456.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012974/DRR014454/DRR014454.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012972/DRR014452/DRR014452.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012970/DRR014450/DRR014450.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012968/DRR014448/DRR014448.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012966/DRR014446/DRR014446.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012965/DRR014445/DRR014445.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012963/DRR014443/DRR014443.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012961/DRR014441/DRR014441.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012959/DRR014439/DRR014439.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012957/DRR014437/DRR014437.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012955/DRR014435/DRR014435.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012954/DRR014434/DRR014434.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012952/DRR014432/DRR014432.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012950/DRR014430/DRR014430.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012948/DRR014428/DRR014428.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012946/DRR014426/DRR014426.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012944/DRR014424/DRR014424.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012943/DRR014423/DRR014423.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012941/DRR014421/DRR014421.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012939/DRR014419/DRR014419.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012937/DRR014417/DRR014417.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012935/DRR014415/DRR014415.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX012/DRX012933/DRR014413/DRR014413.sra $ mv DRR014456.sra BRIC_siSTAU1_705min_2_DRR014456.sra $ mv DRR014434.sra BRIC_siSTAU1_705min_1_DRR014434.sra $ mv DRR014454.sra BRIC_siSTAU1_465min_2_DRR014454.sra $ mv DRR014432.sra BRIC_siSTAU1_465min_1_DRR014432.sra $ mv DRR014452.sra BRIC_siSTAU1_225min_2_DRR014452.sra $ mv DRR014430.sra BRIC_siSTAU1_225min_1_DRR014430.sra $ mv DRR014450.sra BRIC_siSTAU1_105min_2_DRR014450.sra $ mv DRR014428.sra BRIC_siSTAU1_105min_1_DRR014428.sra $ mv DRR014448.sra BRIC_siSTAU1_45min_2_DRR014448.sra $ mv DRR014426.sra BRIC_siSTAU1_45min_1_DRR014426.sra $ mv DRR014446.sra BRIC_siSTAU1_0min_2_DRR014446.sra $ mv DRR014424.sra BRIC_siSTAU1_0min_1_DRR014424.sra $ mv DRR014445.sra BRIC_siCTRL_705min_2_DRR014445.sra $ mv DRR014423.sra BRIC_siCTRL_705min_1_DRR014423.sra $ mv DRR014443.sra BRIC_siCTRL_465min_2_DRR014443.sra $ mv DRR014421.sra BRIC_siCTRL_465min_1_DRR014421.sra $ mv DRR014441.sra BRIC_siCTRL_225min_2_DRR014441.sra $ mv DRR014419.sra BRIC_siCTRL_225min_1_DRR014419.sra $ mv DRR014439.sra BRIC_siCTRL_105min_2_DRR014439.sra $ mv DRR014417.sra BRIC_siCTRL_105min_1_DRR014417.sra $ mv DRR014437.sra BRIC_siCTRL_45min_2_DRR014437.sra $ mv DRR014415.sra BRIC_siCTRL_45min_1_DRR014415.sra $ mv DRR014435.sra BRIC_siCTRL_0min_2_DRR014435.sra $ mv DRR014413.sra BRIC_siCTRL_0min_1_DRR014413.sra $ for file in *.sra do fastq-dump ${file} done $ cat BRIC_siSTAU1_0min_1_DRR014424.fastq BRIC_siSTAU1_0min_2_DRR014446.fastq > BRIC_siSTAU1_0min_DRR014424_DRR014446.fastq $ cat BRIC_siSTAU1_45min_1_DRR014426.fastq BRIC_siSTAU1_45min_2_DRR014448.fastq > BRIC_siSTAU1_45min_DRR014426_DRR014448.fastq $ cat BRIC_siSTAU1_105min_1_DRR014428.fastq BRIC_siSTAU1_105min_2_DRR014450.fastq > BRIC_siSTAU1_105min_DRR014428_DRR014450.fastq $ cat BRIC_siSTAU1_225min_1_DRR014430.fastq BRIC_siSTAU1_225min_2_DRR014452.fastq > BRIC_siSTAU1_225min_DRR014430_DRR014452.fastq $ cat BRIC_siSTAU1_465min_1_DRR014432.fastq BRIC_siSTAU1_465min_2_DRR014454.fastq > BRIC_siSTAU1_465min_DRR014432_DRR014454.fastq $ cat BRIC_siSTAU1_705min_1_DRR014434.fastq BRIC_siSTAU1_705min_2_DRR014456.fastq > BRIC_siSTAU1_705min_DRR014434_DRR014456.fastq $ cat BRIC_siCTRL_0min_1_DRR014413.fastq BRIC_siCTRL_0min_2_DRR014435.fastq > BRIC_siCTRL_0min_DRR014413_DRR014435.fastq $ cat BRIC_siCTRL_45min_1_DRR014415.fastq BRIC_siCTRL_45min_2_DRR014437.fastq > BRIC_siCTRL_45min_DRR014415_DRR014437.fastq $ cat BRIC_siCTRL_105min_1_DRR014417.fastq BRIC_siCTRL_105min_2_DRR014439.fastq > BRIC_siCTRL_105min_DRR014417_DRR014439.fastq $ cat BRIC_siCTRL_225min_1_DRR014419.fastq BRIC_siCTRL_225min_2_DRR014441.fastq > BRIC_siCTRL_225min_DRR014419_DRR014441.fastq $ cat BRIC_siCTRL_465min_1_DRR014421.fastq BRIC_siCTRL_465min_2_DRR014443.fastq > BRIC_siCTRL_465min_DRR014421_DRR014443.fastq $ cat BRIC_siCTRL_705min_1_DRR014423.fastq BRIC_siCTRL_705min_2_DRR014445.fastq > BRIC_siCTRL_705min_DRR014423_DRR014445.fastq $ rm BRIC_siSTAU1_705min_2_DRR014456.fastq $ rm BRIC_siSTAU1_705min_1_DRR014434.fastq $ rm BRIC_siSTAU1_465min_2_DRR014454.fastq $ rm BRIC_siSTAU1_465min_1_DRR014432.fastq $ rm BRIC_siSTAU1_225min_2_DRR014452.fastq $ rm BRIC_siSTAU1_225min_1_DRR014430.fastq $ rm BRIC_siSTAU1_105min_2_DRR014450.fastq $ rm BRIC_siSTAU1_105min_1_DRR014428.fastq $ rm BRIC_siSTAU1_45min_2_DRR014448.fastq $ rm BRIC_siSTAU1_45min_1_DRR014426.fastq $ rm BRIC_siSTAU1_0min_2_DRR014446.fastq $ rm BRIC_siSTAU1_0min_1_DRR014424.fastq $ rm BRIC_siCTRL_705min_2_DRR014445.fastq $ rm BRIC_siCTRL_705min_1_DRR014423.fastq $ rm BRIC_siCTRL_465min_2_DRR014443.fastq $ rm BRIC_siCTRL_465min_1_DRR014421.fastq $ rm BRIC_siCTRL_225min_2_DRR014441.fastq $ rm BRIC_siCTRL_225min_1_DRR014419.fastq $ rm BRIC_siCTRL_105min_2_DRR014439.fastq $ rm BRIC_siCTRL_105min_1_DRR014417.fastq $ rm BRIC_siCTRL_45min_2_DRR014437.fastq $ rm BRIC_siCTRL_45min_1_DRR014415.fastq $ rm BRIC_siCTRL_0min_2_DRR014435.fastq $ rm BRIC_siCTRL_0min_1_DRR014413.fastq
2. Quality check
FastQCを使って、FASTQファイルに記載されている各リードのクオリティをチェックします。
FASTQファイルはFred Quality Scoreと呼ばれる指標によって、各塩基のクオリティを記述しています。
詳しくはこちらのサイトを参照のこと。
https://en.wikipedia.org/wiki/Phred_quality_score
実行例
$ mkdir fastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446 $ fastqc -t 8 -o ./fastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446 \ ./BRIC_siSTAU1_0min_DRR014424_DRR014446.fastq -f fastq
mkdir
コマンドでfastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446
という名前のディレクトリをカレントディレクトリに作ります。
次に、fastqcを使ってBRIC_siSTAU1_0min_DRR014424_DRR014446.fastq
のデータのクオリティをチェックします。 -t
オプションで使用するCPUのコア数を指定します。コア数を増やせば単純に計算速度が向上しますが、FastQCではメモリを1コアあたり250MB消費するので、要求できるメモリ量を超えるコア数を指定するとエラーが起こります。
-o
オプションで出力先のディレクトリを指定し、-f
オプションでInputファイルのファイル形式を指定しています。
FastQCから出力されるデータが、fastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446
ディレクトリの中のBRIC_siSTAU1_0min_DRR014424_DRR014446_fastqc.html
になります。
ファイルをブラウザで開くと以下のようなグラフを確認できます。
Basic Statistics
基本的な統計量が示されています。
- 全リード数
- リード長
- GC contents
などなど
Per base sequence quality
リードの各塩基のクオリティスコアを示しています。 Phred quality scoreがだいたいグリーンの領域(Scoreが28以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。
Per tile sequence quality
フローセルの各タイルごとのクオリティスコアを示しています。
Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。
シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。
特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。
詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf
Per sequence quality scores
各リードの全長のクオリティスコアの分布を示しています。
Per base sequence content
各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。
Per sequence GC content
リードのGC contentsの分布を示しています。
Per base N content
各塩基中に含まれるN
の含有率(塩基を読めなかった箇所)を示しています。
Sequence Length Distribution
リード長の分布を示しています。
Sequence Duplication Levels
Duplidate readsの含まれている数を示しています。
Overrepresented sequences
頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。
Adapter Content
各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。
Kmer Content
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
3. Quality filtering
FASTQファイル中にクオリティスコアの低いリードが含まれていることが確認された場合、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。
実行例
$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./BRIC_siSTAU1_0min_DRR014424_DRR014446.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o ./BRIC_siSTAU1_0min_DRR014424_DRR014446_1_filtered.fastq
fastq_quality_trimmerのパラメータ
-Q
: Phred quality scoreのバージョン指定(Solexa, Illumina 1.3+, Illumina 1.5+といったかなり古いバージョンのFASTQファイルでは'64'を指定。最近のFASTQファイルは'33'と指定しておけばOK)。-t
: 指定したPhred quality scoreを下回る塩基をトリミングする。-l
: トリミング後の最低リード長。
fastq_quality_filterのパラメータ
-q 20 -q 80
: リードの80%
でPhred quality scoreが20
を下回る場合、そのリードを排除する。
4. rRNA removal
polyAセレクションやRiboZero等のrRNAのコンタミを除去しない場合、BRIC-seqではかなりの割合(50%以上)でリボソームRNAのコンタミを多く含んでいます。そのため、あらかじめrRNAのコンタミを除いたリードをGenomeやTranscriptomeにリードをマッピングしたほうが計算時間を短縮できます。ここではBowtie1を用いて、リボソームRNAにマッピングされなかったリードを抽出します。
Bowtieでマッピングするために、まず、rRNAのIndexファイルを用意する必要があります。まずは、Indexファイルを作成するための配列情報をFASTAファイルとして取ってきます。
NCBIに登録されているX12811.1| Human 5S DNA
とU13369.1|HSU13369 Human ribosomal DNA complete repeating unit
の配列情報を使います。
下記のURLからFASTA形式のデータをコピペしてFASTAファイルをつくります。
Human 5S DNA
https://www.ncbi.nlm.nih.gov/nuccore/23898
Human ribosomal DNA complete repeating unit
https://www.ncbi.nlm.nih.gov/nuccore/555853
上記のrRNAのFASTA形式のデータをcontam_Ribosomal_RNA.fa
と名前をつけてまとめます。
以下で、先ほど作成したFASTAファイルをもとに、Bowtie用のIndexファイルを作成しています。
$ bowtie-build ./contam_Ribosomal_RNA.fa contam_Ribosomal_RNA
- 1つ目の引数: 先ほど作成したrRNAのFASTAファイルを指定。
- 2つ目の引数: Indexファイル名を記載。
Indexファイルを作成した後、rRNAとRepetitive elementsへのマッピングを行います。マッピングされなかったリードについては、FASTQファイルとして出力されるようにオプションを指定します。
実行例
$ bowtie -p 8 --un ./BRIC_siSTAU1_0min_DRR014424_DRR014446_2_norrna.fastq \ ./contam_Ribosomal_RNA \ ./BRIC_siSTAU1_0min_DRR014424_DRR014446_1_filtered.fastq \ > rRNA_BRIC_siSTAU1_0min_DRR014424_DRR014446.fastq \ 2>> ./log_BRIC_siSTAU1_0min_DRR014424_DRR014446.txt
-p
: 使用するコア数。--un
: マッピングされなかったリードの出力先(FASTQで出力される)。- 1つ目の引数: マッピングさせたいFASTQファイルを指定。
- 標準出力でマッピングさせた結果は出力される。
- 標準エラー出力でログデータが出力される。ここでは、
2>>
でこの標準エラー出力の結果をtxtファイルとして追加保存している。
5. Quality check
rRNAやRepetitive elementsに由来するリードを除いた後、フィルタリング後のFASTQファイルをFastQCを用いて再確認します。
実行例
$ mkdir fastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446_filtered $ fastqc -o ./fastqc_BRIC_siSTAU1_0min_DRR014424_DRR014446_filtered ./BRIC_siSTAU1_0min_DRR014424_DRR014446_2_norrna.fastq -f fastq
6. Mapping
TopHatを用いてGenomeとTranscriptomeにリードをマッピングを行います。
注意!!
RNA-seq(BRIC-seqも同様)のデータ解析については近年、マッピングソフトとしてSTARがデファクトスタンダードになっています(あくまで個人的な意見です。ENDODEプロジェクトで使われて一気にスタンダードになった気がします)。
確かに、STARのほうがTopHatと比較してマッピングの精度が向上しているので、現時点ではSTARがファーストチョイスかと思います。
しかし、STARはメモリを大量に消費し、ヒトゲノムだと32GB以上は必要となります。メモリ不足が懸念される解析環境では、メモリ消費が少ないTopHatを選択するのが無難だと思います(昔はこれがデファクトだったので)。
Indexファイルの作成
ここでも先ほどのマッピング時と同じように、特定のGenomeのBowtie用のIndexファイルを作成する必要があります。
Bowtieのサイトから、Indexファイルはダウンロードすることもできます。右端のカラムに各ゲノムのIndexファイルが羅列されているので、任意のIndexファイルをダウンロードしてきます。
http://bowtie-bio.sourceforge.net/index.shtml
もしくは、自身でIndexファイルを作成することも可能です。以下では、その方法を説明します。
まず、Human Genome(hg19)のFASTAファイルのダウンロードしてきます。 各染色体ごとにファイルが用意されているので、catコマンドでダウンロードしたすべてのファイルをマージして、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ファイル)に対して行います。
実行例
$ tophat --bowtie1 -p 8 -o ./tophat_out_BRIC_siSTAU1_0min_DRR014424_DRR014446 \ -G ./gencode.v19.annotation.gtf \ ./hg19 ./BRIC_siSTAU1_0min_DRR014424_DRR014446_2_norrna.fastq
--bowtie1
: マッピングソフトにBowtie1
を使用。-p
: 使用するコア数を指定。-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
7. データの品質チェック
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_BRIC_siSTAU1_0min_DRR014424_DRR014446 $ samtools index ./tophat_out_BRIC_siSTAU1_0min_DRR014424_DRR014446/accepted_hits.bam
mkdir
コマンドで新しいディレクトリを作成。samtools index
でBAMファイルのIndexファイルを作成します。引数にはBAMファイルを指定するだけです。
それでは、House-keeping genesの遺伝子領域のCoverageを調べてみたいと思います。
実行例
$ geneBody_coverage.py -r ./hg19.HouseKeepingGenes.bed \ -i ./tophat_out_BRIC_siSTAU1_0min_DRR014424_DRR014446/accepted_hits.bam \ -o ./geneBody_coverage_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_RSeQC_output
-r
: 先ほどダウンロードしたHouse-keeping genesの情報が格納されているBEDファイルを指定。-i
: Coverageを調べたいBAMファイルを指定(加えて、Indexファイルを同じディレクトリにおいておく必要がある)。-o
: 出力先のディレクトリとファイル名のPrefix。
結果として、指定した名前.geneBodyCoverage.curves.pdf
というファイルが出力されます。ファイルの中身を見るとこんな感じ。
RNAの品質に問題がある(例えば、RNAが壊れている)と、3'側にリードが偏っていたり、cDNAライブラリの構築に使ったRNA量が少なく、PCR時に特定の配列が異常に増幅され(PCR bias、Duplicate readsが多い)、ガタガタしたCoverageの分布を取ることがあります。
こういったデータはRNAの定量性に問題があるため、使用しないように気をつけましょう。
このステップで問題がありそうだと判断されるデータに関しては、下記のUCSC genome browserでの可視化により、個々の遺伝子上にマッピングされたリードの分布に注意して見ていく必要があります。
8. Visualization (For UCSC genome browser)
BedGraphファイルの準備
Bedtools
を利用して、TopHatから得られたBAMファイルを、UCSC genome browser上で可視化できるファイル(Bedgraphファイル)に変換します。
まず、mkdir
コマンドでデータを保存するディレクトリを作っておきます。
$ mkdir UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446
bedtools genomecovで特定のゲノム領域でのCoverageを計算します(BAMファイルからBedgraph
ファイルを作成)。
Bedgraphファイルの詳細については、下記のURLを参照のこと。
https://genome.ucsc.edu/goldenpath/help/bedgraph.html
$ bedtools genomecov -ibam ./tophat_out_BRIC_siSTAU1_0min_DRR014424_DRR014446/accepted_hits.bam -bg -split \ > ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_4_result.bg
-ibam
: BAMファイルを入力ファイルとして指定する。-bg
: 出力ファイルの形式をBedGraphファイルに指定。-split
: Splicing juctionにマッピングされたリード(リードが分割されてマッピングされているリード)を考慮する。
標準出力でBedGraph
ファイルが出力されます。
UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。
https://genome.ucsc.edu/goldenpath/help/customTrack.html
それらの情報をヘッダー行(ファイルの一行目)に記載します。
$ echo "track type=bedGraph name=BRIC_siSTAU1_0min_DRR014424_DRR014446 description=BRIC_siSTAU1_0min_DRR014424_DRR014446 visibility=2 maxHeightPixels=40:40:20" \ > ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/tmp.txt
echo
コマンドを使って、ヘッダー行の情報をtxtファイルとして出力しています。
そのあとに、cat
コマンドを用いて、先ほど作ったBedGraphファイルと結合させます。
$ cat ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/tmp.txt ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_4_result.bg \ > ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_4_result_for_UCSC.bg
出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2ファイルに圧縮します。
$ bzip2 -c ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_4_result_for_UCSC.bg \ > ./UCSC_visual_BRIC_siSTAU1_0min_DRR014424_DRR014446/BRIC_siSTAU1_0min_DRR014424_DRR014446_4_result_for_UCSC.bg.bz2
データのアップロード
UCSC genome browserにデータをアップロードします。
ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。
Add Custom Tracksというページに移動するので、ファイルを選択をクリックし、アップロードするファイルを選択し、submitボタンをクリックします。
アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgoボタンをクリックすると、データをみることができます。
9. RPKMの算出
BRIC-seqにより得られた各時間のデータ(0h, 1h, 2h, 4h, 8h, 12h)に対して、Cuffnormを用いてRPKM値を算出すします。
Cuffnormの計算に必要なアノテーション情報(GTFファイル)は、Mapping時にダウンロードしたgencode.v19.annotation.gtf
を使用します。
# siCTRL $ cuffnorm -p 8 --compatible-hits-norm -o BRIC-seq_siCTRL1_Gencode ./gencode.v19.annotation_filtered.gtf \ tophat_out_BRIC_siCTRL_0min_DRR014413_DRR014435/accepted_hits.bam \ tophat_out_BRIC_siCTRL_45min_DRR014415_DRR014437/accepted_hits.bam \ tophat_out_BRIC_siCTRL_105min_DRR014417_DRR014439/accepted_hits.bam \ tophat_out_BRIC_siCTRL_225min_DRR014419_DRR014441/accepted_hits.bam \ tophat_out_BRIC_siCTRL_465min_DRR014421_DRR014443/accepted_hits.bam \ tophat_out_BRIC_siCTRL_705min_DRR014423_DRR014445/accepted_hits.bam # siSTAU1 $ cuffnorm -p 8 --compatible-hits-norm -o BRIC-seq_siSTAU_Gencode ./gencode.v19.annotation_filtered.gtf \ tophat_out_BRIC_siSTAU1_0min_DRR014424_DRR014446/accepted_hits.bam \ tophat_out_BRIC_siSTAU1_45min_DRR014426_DRR014448/accepted_hits.bam \ tophat_out_BRIC_siSTAU1_105min_DRR014428_DRR014450/accepted_hits.bam \ tophat_out_BRIC_siSTAU1_225min_DRR014430_DRR014452/accepted_hits.bam \ tophat_out_BRIC_siSTAU1_465min_DRR014432_DRR014454/accepted_hits.bam \ tophat_out_BRIC_siSTAU1_705min_DRR014434_DRR014456/accepted_hits.bam
-p
: 使用するコア数を指定。--compatible-hits-norm
: アノテーション情報(GTFファイル)に含まれる遺伝子群にマッピングされた総リード数からRPKM値を算出する。-o
: 出力先のディレクトリを指定。1つ目の引数: アノテーション情報のGTFファイルを指定。
- 2つ目以降の引数: BAMファイルを指定。空白区切りで、時系列順にデータを指定(例えば、0h, 1h, 2h, 4h, 8h, 12hの順)。
パラメータの詳しい解説は以下のサイトを参照のこと。
http://cole-trapnell-lab.github.io/cufflinks/cuffnorm/
10. RNA半減期の計算
カスタムRパッケージであるbridger2
を用いて、RNAの半減期を算出します。
遺伝子リストの準備
以下のURLに必要なスクリプトが置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/BRIC-seq_Tutorial
まず、extract_gene_symbol_type_from_gtf.py
スクリプトを用いて、GTFファイルから遺伝子リストを作成します(Cuffnormで集計したデータからmRNAとlncRNAに分類分けするために用います)。
$ python2 extract_gene_symbol_type_from_gtf.py gencode.v19.annotation_filtered.gtf gencode.v19.annotation_filtered_symbol_type_list.txt
- 1つ目の引数: GTFファイルを指定する。
- 2つ目の引数: 遺伝子リストが出力される。
bridger2向けのInputファイルの作成
siCTRLとsiPUM1に対するBRIC-seqデータから得られたRPKM値のリストから、mRNAとlncRNAのデータをそれぞれ抽出すします。
# siCTRL ## mRNA $ python2 ~/custom_command/BridgeR_prep.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./BRIC-seq_siCTRL_Gencode/genes.fpkm_table \ mRNA \ ./BRIC-seq_siCTRL_Gencode/BridgeR_input_file_mRNA.txt ## lncRNA $ python2 ~/custom_command/BridgeR_prep.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./BRIC-seq_siCTRL_Gencode/genes.fpkm_table \ lncRNA \ ./BRIC-seq_siCTRL_Gencode/BridgeR_input_file_lncRNA.txt # siPUM1 ## mRNA $ python2 ~/custom_command/BridgeR_prep.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./BRIC-seq_siSTAU_Gencode/genes.fpkm_table \ mRNA \ ./BRIC-seq_siSTAU_Gencode/BridgeR_input_file_mRNA.txt ## lncRNA $ python2 ~/custom_command/BridgeR_prep.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./BRIC-seq_siSTAU_Gencode/genes.fpkm_table \ lncRNA \ ./BRIC-seq_siSTAU_Gencode/BridgeR_input_file_lncRNA.txt
- 1つ目の引数: 先ほど用意した遺伝子リストのファイル (
gencode.v19.annotation_filtered_symbol_type_list.txt
)を指定する。 - 2つ目の引数: Cuffnormの出力ファイル
genes.fpkm_table
を指定する。 - 3つ目の引数:
mRNA
もしくはlncRNA
を指定。指定したRNA群を抽出する。 - 4つ目の引数: 出力ファイル名を指定する。
RNA半減期を測定する
bridger2
パッケージを用いてRNA半減期を計算します。
以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/BRIC-seq_Tutorial
mRNAとlncRNAのそれぞれに対して、まずmkdir
コマンドで保存先にディレクトリを作成します。その後、BridgeR_analysis_mRNA.R
を実行し、RNA半減期を計算します。
# mRNA half-life $ mkdir BridgeR_result_mRNA $ Rscript ./BridgeR_analysis_mRNA.R \ ./BridgeR_result_mRNA \ ./BRIC-seq_siCTRL_Gencode/BridgeR_input_file_mRNA.txt,./BRIC-seq_siSTAU_Gencode/BridgeR_input_file_mRNA.txt # lncRNA half-life $ mkdir BridgeR_result_lncRNA $ Rscript ./BridgeR_analysis_lncRNA.R \ ./BridgeR_result_lncRNA \ ./BridgeR_result_mRNA \ ./BRIC-seq_siCTRL_Gencode/BridgeR_input_file_lncRNA.txt,./BRIC-seq_siSTAU_Gencode/BridgeR_input_file_lncRNA.txt
BridgeR_analysis_mRNA.R
スクリプト
- 1つ目の引数: 作成したmRNAデータの保存先ディレクトリを指定。
- 2つ目の引数: 先ほど用意したsiCTRLとsiPUM1のデータの場所を指定。2つ以上のデータに対しては、カンマ区切りでつなげる。
BridgeR_analysis_lncRNA.R
スクリプト
- 1つ目の引数: 作成したlncRNAデータの保存先ディレクトリを指定。
- 2つ目の引数: 作成したmRNAデータの保存先ディレクトリを指定(ノーマライズに使用する重み付けデータをRスクリプト内で使用)。
- 3つ目の引数: 先ほど用意したsiCTRLとsiPUM1のデータの場所を指定。2つ以上のデータに対しては、カンマ区切りでつなげる。
BridgeR_analysis_mRNA.R
スクリプトについて
8, 9行目のgroup
とhour
の変数は、サンプルごとに変更する必要がある。
> group <- c("siCTRL", "siSTAU1") # Required > hour <- c(0, 1, 2, 4, 8, 12) # Required
BridgeR_analysis_mRNA.R
やBridgeR_analysis_lncRNA.R
の最後の引数で指定したファイルの数に応じて、groupの名前を指定する。
また、自身のサンプルのタイムコースに合わせて、hour
を変更する。
21-29行目のTimePointRemoval1
, TimePointRemoval2
, ThresholdHalfLife1
, ThresholdHalfLife2
について、タイムコースのとり方に応じて変更する必要がある。
> halflife_table <- BridgeRCore(input_matrix, group=group, hour=hour, save = TRUE, TimePointRemoval1 = c(1,2), TimePointRemoval2 = c(8,12), ThresholdHalfLife1 = 3, ThresholdHalfLife2 = 12, outputPrefix = paste(dirname, "BridgeR", sep="/"))
TimePointRemoval1
とTimePointRemoval2
では、それぞれの点を削ったときにFittingの精度が向上するかどうか判定する際に削る点を指定する。各点は、タイムコースの点と一致させる必要がある。
また、ThresholdHalfLife1
の意味は、すべての点を用いてRNA半減期を求めたときに、半減期が3時間を下回る場合、TimePointRemoval1
で指定した点を除いて計算させたときに、Fittingの精度が向上するか検討する。
ThresholdHalfLife2
の意味は、すべての点を用いてRNA半減期を求めたときに、半減期が12時間を超える場合、TimePointRemoval2
で指定した点を除いて計算させたときに、Fittingの精度が向上するか検討する。
34-37行目のcomparisonFile
とoutputPrefix
で指定する名前をサンプルごとに変更する必要がある。
> pvalue_table <- BridgeRPvalueEvaluation(halflife_table, group = group, hour = hour, comparisonFile = c("siCTRL","siSTAU1"), inforColumn = 4, CutoffTimePointNumber = 4, calibration = FALSE, save = TRUE, outputPrefix = paste(dirname, "BridgeR_6_STAU1KD", sep="/")) # Required
BridgeR_analysis_lncRNA.R
スクリプトについて
9, 10行目のgroup
とhour
の変数は、サンプルごとに変更する必要がある。
> group <- c("siCTRL", "siPUM1") # Required > hour <- c(0, 1, 2, 4, 8, 12) # Required
46-57行目のTimePointRemoval1
, TimePointRemoval2
, ThresholdHalfLife1
, ThresholdHalfLife2
について、先程と同様に、タイムコースのとり方に応じて変更する必要がある。
> halflife_table <- BridgeRHalfLifeCalcR2Select(normalized_table, group = group, hour = hour, inforColumn = 4, CutoffTimePointNumber = 4, R2_criteria = 0.9, TimePointRemoval1 = c(1,2), TimePointRemoval2 = c(8,12), ThresholdHalfLife1 = 3, ThresholdHalfLife2 = 12, save = T, outputPrefix = paste(dirname_lncRNA, "BridgeR_5", sep="/"))
59-62行目のcomparisonFile
とoutputPrefix
で指定する名前をサンプルごとに変更する必要がある。
> pvalue_table <- BridgeRPvalueEvaluation(halflife_table, group = group, hour = hour, comparisonFile = c("siCTRL", "siPUM1"), inforColumn = 4, CutoffTimePointNumber = 4, calibration = FALSE, save = TRUE, outputPrefix = paste(dirname_lncRNA, "BridgeR_6_STAU1KD", sep="/"))
11. BRIC-seqデータの可視化(RNA分解曲線を描かせる)
bridger2
パッケージのBridgeReport
関数を使うことで、RNA分解曲線を描くことができます。
実行するためには、Rの統合開発環境であるRStudio
をあらかじめ自身のパソコンにインストールしておく必要があります。
https://www.rstudio.com/products/RStudio/#Desktop
Rstudioを開き、以下のスクリプトを実行します。
> setwd("C:/Users/Naoto/Downloads") > library(bridger2) > library(data.table) > filename <- "BridgeR_6_STAU1KD_halflife_pvalue_evaluation.txt" > pvalue_table <- fread(filename, header=T) > shiny_test <- BridgeReport(pvalue_table, searchRowName = "gene_symbol", hour = c(0,1,2,4,8,12), TimePointRemoval1 = c(1, 2), TimePointRemoval2 = c(8, 12)) > shiny_test
setwd
でファイルが保存されている先を指定し、filename
にインポートしたいデータ(RNAの半減期の計算で出力されるファイル)を指定します。
表示させると、以下のようになります。
12. shinyapp.ioサーバーにBRIC-seqデータ可視化用のWebアプリを作成する
shinyapp.ioとは、Rstudioを開発する会社が提供しているShinyアプリをWeb上で公開するためのクラウドサーバーです。
無料で使えるFree版も提供されていますが、作成できるアプリ数は5つ
までで、一ヶ月に稼働できる時間は25時間
までという制限が設けられています。
shinyapp.ioのアカウントを取得する
下記のURLにアクセスし、アカウントを作成する。
https://www.shinyapps.io/admin/#/signup
Shinyアプリの作成・デプロイ
RStudioを開き、新しいプロジェクトを作成する。
保存したフォルダを開き、server.R
とui.R
を消します。代わりにapp.R
というファイルを作成してください(ここは必ずapp.R
という名前にしてください。別名にするとshinyapp.ioにデプロイできません)。
また、10. RNA半減期の計算
の最後に生成された./BridgR_result_PUM1_study_mRNA/BridgeR_6_PUM1KD_halflife_pvalue_evaluation.txt
のファイルをコピーして同じフォルダに入れます。
こんな感じになります。
上図のように、保存先のフォルダにxxxxx.Rproj
というファイルができているので、RStudioで開きます。さらに、先ほど作成したapp.R
も開きます。
app.R
スクリプトに以下の内容を書き込みます。
library(bridger2) library(shiny) library(data.table) filename <- "BridgeR_6_STAU1KD_halflife_pvalue_evaluation.txt" pvalue_table <- fread(filename, header=T) shiny_test <- BridgeReport(pvalue_table, searchRowName = "gene_symbol", hour = c(0,1,2,4,8,12), TimePointRemoval1 = c(1, 2), TimePointRemoval2 = c(8, 12)) shiny_test
bridger2
で用意されているRridgeReport
関数は、先ほど削除したserver.R
とui.R
のスクリプトをBRIC-seqデータ用にまとめたWrapperになっているので、この関数1つでShinyアプリが簡単に作れます。
元のスクリプトはこちら。
https://github.com/Imamachi-n/BridgeR2/blob/master/R/reporting.R
RStudioのConsoleから、rsconnect
パッケージをインストールします。これはshinyapp.ioのサーバーへShinyアプリをデプロイするのに必要になります。
> install.packages('rsconnect')
rsconnect
パッケージを呼び込みます。
library(rsconnect)
自身のshinyapp.ioのアカウントにログインし、サーバーとやり取りするためにトークンを取得します。
[Account] -> [Tokens]からAdd Token
をクリックし、Tokenを生成してShow
をクリックする。すると、RStudioのConsoleで実行するコマンドが表示されるので、Copy to clipboard
をクリックして内容をコピーします。
[Ctrl] + [v]でRStudioのConsoleにペーストして実行します。
> rsconnect::setAccountInfo(name='hogehoge', token='xxxxxx', secret='yyyyy')
念のため、app.R
スクリプトが動くかどうかローカル環境で確認します。
library(shiny) runApp()
実行して、11. BRIC-seqデータの可視化(RNA分解曲線を描かせる)
の時と同様にポップアップウインドウが表示され、Input gene symbol
に適当な遺伝子名を入れたときにRNAの分解曲線が表示されればOKです。
最後に、shinyapp.ioのクラウドサーバーにShinyアプリをデプロイします。
library(rsconnect) deployApp()
しばらく待ちます。初回は必要なRパッケージのインストールなどが行われるため、時間がかかります。
下記のように、Application successfully deployed to https://hogehoge.shinyapps.io/hoge
と出たら、デプロイ成功です。
表示されているURLにアクセスすると、今までRStudio上で見ていたRNAの分解曲線のデータをWeb上でも閲覧することが可能になります。
また、指定されたURLを他人に公開することで、RStudioをインストールしていないユーザーに対しても、Webブラウザを介してデータを共有することができます。
shinyapps.ioの自身のDashboardからもShinyアプリのURLを確認することができます。
Recent ApplicationsのところにアップロードしたShinyアプリが表示されており、それをクリックすると、URLなどの情報が記載されているページに飛ぶことができます。
内容は以上になります。
RNA-seqデータ解析(2群間比較, Strand-specific, Single-end)
概要
RNA-seqのデータを用いて、siRNAによるノックダウンや特定の刺激により発現量が変動したRNA群を同定する方法(2群間のデータの比較)を紹介します。
ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、RNA-seqのデータ解析の詳細をStep by Stepで説明します。
ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。 パソコン(もしくはスパコン)のメモリは16GB以上を確保することが望ましいです。
使用するデータ
RNA-seq (Single-end, 101bp; Control siRNA, UPF1 siRNA)
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148
使用するソフトウェア
- FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - Bowtie1 (v1.1.2)
http://bowtie-bio.sourceforge.net/index.shtml - TopHat (v2.1.1)
https://ccb.jhu.edu/software/tophat/index.shtml - Samtools (v1.3.1)
http://samtools.sourceforge.net/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - bedtools (v2.26.0)
http://bedtools.readthedocs.io/en/latest/ - Cufflinks (v2.2.1)
http://cole-trapnell-lab.github.io/cufflinks/ - featureCounts (subread v1.5.0.post3)
http://subread.sourceforge.net/ - edgeR (v3.14.0)
https://bioconductor.org/packages/release/bioc/html/edgeR.html
解析ワークフロー
RNA-seqデータを2群間比較する方法として、Cuffdiffを用いた方法とfeatureCount-edgeRを用いた方法を紹介したいと思います。
TopHat-Cuffdiff workflow
TopHat-featureCounts-edgeR workflow
サンプルコード
Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Strand-specific_RNA-seq_Tutorial
RNA-seq_strand-specific_single-end_1_mapping.sh
TopHatによるマッピング・featureCountsによるリードカウントまでを行うスクリプト。FASTQファイルを引数として指定。
$ ./RNA-seq_strand-specific_single-end_1_mapping.sh hoge.fastq
10行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
11行目: indexContamFile
にrRNAなどのコンタミ配列のBowtie用のIndexファイルを指定。
12行目: indexGenomeFile
にGenome配列のBowtie用のIndexファイルを指定。
31行目: -r
オプションでHouse-keeping genesの遺伝子領域のBEDファイルを指定。
RNA-seq_strand-specific_single-end_2_cuffdiff_calc.sh
Cuffdiffによる2群間比較を行うスクリプト。
$ ./RNA-seq_strand-specific_single-end_2_cuffdiff_calc.sh
8行目: filename
に保存先のディレクトリ名を指定。
9行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
13-14行目: カンマ(,
)区切りでBAMファイル(比較したい2種類)を指定。
17行目: gene_list
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
RNA-seq_strand-specific_single-end_2_edgeR_calc.sh
edgeRによる2群間比較を行うスクリプト。
$ ./RNA-seq_strand-specific_single-end_2_edgeR_calc.sh
8行目: annoList
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
11行目: saveDir
に保存先のディレクトリ名を指定。
~/custom_command/<hoge.py>
となっている箇所はカスタムのPythonスクリプトが置いてあるパス・ファイル名を指定する。
シェルスクリプトを使用するときの注意点
権限の問題
実行権限をシェルスクリプトに与えた後、実行すること。
$ chmod 755 hoge.sh
改行コードの問題
WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。
ここではAtom
と呼ばれるテキストエディタでの設定方法を説明する。
- [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
- 左のメニューから、
Packages
をクリック。 - Installed Packagesで
line-ending-selector
と検索。 - Core Packagesの項目から、line-ending-selectorの
Settings
をクリック。 - 設定画面で、Default line endingを
LF
に変更。
これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。
各ステップの説明
0. 必要なソフトウェアのインストール
まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。
$ 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 cufflinks $ conda install subread $ conda install bioconductor-edger
注意!!
tophatによるマッピングで最新のbowtieを使うとエラーが出るので、以前のバージョンのbowtie v1.12を代わりにインストールする。
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
という項目にリストアップしてあります。
今回は、35サンプルのうち上6つのサンプルを使用します。
GSMから始まるIDをクリックすると、各サンプルの情報を見ることができます。 サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。
そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。
SRR4081222をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。
SRR4081222.sra
のURLをコピー。
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以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。
Per tile sequence quality
フローセルの各タイルごとのクオリティスコアを示しています。
Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。
シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。
特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。
詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf
Per sequence quality scores
各リードの全長のクオリティスコアの分布を示しています。
Per base sequence content
各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。
Per sequence GC content
リードのGC contentsの分布を示しています。
Per base N content
各塩基中に含まれるN
の含有率(塩基を読めなかった箇所)を示しています。
Sequence Length Distribution
リード長の分布を示しています。
Sequence Duplication Levels
Duplidate readsの含まれている数を示しています。
Overrepresented sequences
頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。
Adapter Content
各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。
Kmer Content
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
3. Quality filtering
FASTQファイル中にクオリティスコアの低いリードが含まれていることが確認された場合、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。
今回のデータでは、クオリティの低いリードが含まれていないのでこの工程をパスします。
Illumina製のシークエンサーでいうと、HiSeq以降のデータはデータのクオリティが非常に高いので、いちいちクオリティフィルターにかける必要がなくなりました。
一応、実行例を示すと以下の通りになります。
実行例
$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./SRR4081222_Control_1.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o ./SRR4081222_Control_1_1_filtered.fastq
fastq_quality_trimmerのパラメータ
-Q
: Phred quality scoreのバージョン指定(Solexa, Illumina 1.3+, Illumina 1.5+といったかなり古いバージョンのFASTQファイルでは'64'を指定。最近のFASTQファイルは'33'と指定しておけばOK)。-t
: 指定したPhred quality scoreを下回る塩基をトリミングする。-l
: トリミング後の最低リード長。
fastq_quality_filterのパラメータ
-q 20 -q 80
: リードの80%
でPhred quality scoreが20
を下回る場合、そのリードを排除する。
4. Quality check
フィルタリング後のFASTQファイルをFastQCを用いて再確認します(今回のデータはフィルタリングを行わないのでこの工程もパス)。
実行例
$ mkdir fastqc_SRR4081222_Control_1_filtered $ fastqc -o ./fastqc_SRR4081222_Control_1_filtered ./SRR4081222_Control_1_1_filtered.fastq -f fastq
5. Mapping
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
6. データの品質チェック
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
というファイルが出力されます。ファイルの中身を見るとこんな感じ。
RNAの品質に問題がある(例えば、RNAが壊れている)と、3'側にリードが偏っていたり、cDNAライブラリの構築に使ったRNA量が少なく、PCR時に特定の配列が異常に増幅され(PCR bias、Duplicate readsが多い)、ガタガタしたCoverageの分布を取ることがあります。
こういったデータはRNAの定量性に問題があるため、使用しないように気をつけましょう。
このステップで問題がありそうだと判断されるデータに関しては、下記のUCSC genome browserでの可視化により、個々の遺伝子上にマッピングされたリードの分布に注意して見ていく必要があります。
7. マッピングされたリードの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で可視化させるデータを作成する際に、+と-鎖をそれぞれ分けてファイルを作成します。
8. 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ボタンをクリックすると、データをみることができます。
8. RNA-seqデータの2群間比較(Cuffdiffのケース)
Cuffdiffを用いて、siCTRLとsiUPF1の2群間のRNA-seqデータを比較します。
注意すべき点は、--library-type
オプションでサンプルがStrand-specificであることを明示する点です。
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html#library-types
cuffdiff -p 8 --multi-read-correct --library-type fr-firststrand -o ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode \ /home/akimitsu/database/gencode.v19.annotation_filtered.gtf \ ./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
-p
: 使用するコア数を指定。--multi-read-correct
: 定量の精度を上げる目的で、複数箇所にマッピングされたリードに重み付けを行う。--library-type fr-firststrand
: dUTP法を利用したStrand-specific RNA-seqサンプルに対して、この設定を使用する。-o
: 出力先のディレクトリを指定。- 1つ目の引数: アノテーション情報のGTFファイルを指定。
- 2つ目・3つ目の引数: 2種類のBAMファイルをそれぞれ指定。今回のように、siCTRLとsiUPF1の条件でn=3でデータがある場合、それぞれのBAMファイルのパスをカンマ区切りで並べる(例えば、
CTRL_1.bam,CTRL_2.bam,CTRL_3.bam KD_1.bam,KD_2.bam,KD_3.bam
)。
パラメータの詳しい解説は以下のサイトを参照のこと。
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html
最後に、Gene symbolなどの情報をCuffdiffから出力されるファイルに加えます。
遺伝子リストの準備
まず、extract_gene_symbol_type_from_gtf.py
スクリプトを用いて、GTFファイルから遺伝子リストを作成します(Cuffdiffで集計したデータからmRNAとlncRNAに分類分けするために用います)。
以下のURLに必要なスクリプトが置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
$ python2 extract_gene_symbol_type_from_gtf.py gencode.v19.annotation_filtered.gtf gencode.v19.annotation_filtered_symbol_type_list.txt
- 1つ目の引数: GTFファイルを指定する。
- 2つ目の引数: 遺伝子リストが出力される。
次に、先ほど作成した遺伝子リストをもとに、cuffdiffの出力ファイルからmRNAとlncRNAの情報を抽出し、リストアップします。
# mRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ mRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_mRNA.txt # lncRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ lncRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_lncRNA.txt
- 1つ目の引数: 先ほど用意した遺伝子リストを指定。
- 2つ目の引数: cuffdiffの出力ファイルである
gene_exp.diff
ファイルを指定。 - 3つ目の引数:
mRNA
もしくはlncRNA
を指定(指定したRNA群のデータを取得)。 - 4つ目の引数: mRNAもしくはlncRNAのリストを出力。
8. RNA-seqデータの2群間比較(featureCounts-edgeRのケース)
各遺伝子に対してのリード数のカウント
まず、featureCountsで各遺伝子領域にマッピングされたリード数をカウントします。注意すべき点は、-s
オプションでStrand-specificであることを明示する点です。
$ mkdir featureCounts_result_SRR4081222_Control_1 $ featureCounts -T 8 -t exon -g gene_id -s 2 -a ./gencode.v19.annotation_filtered.gtf \ -o featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ ./tophat_out_SRR4081222_Control_1/accepted_hits.bam
mkdir
コマンドを使って、保存先のディレクトリを作り、そのあとfeatureCountsを実行します。
-T
: 使用するコア数の指定。-t exon
: リードのカウントに使用するExonの情報を指定。ここでは、GTFファイル内で3列目のfeature
がexon
である行を使用する用に指定している。-g gene_id
: 各行のメタ情報(遺伝子ID、Transcript ID、遺伝子名など)のうち、名前として使用する情報を指定する。ここでは、gene_id
を指定することで、Gencodeで決められた遺伝子IDを使っている。-s 2
: Transcriptの配列に対してAntisense鎖にあたるリードをカウントする(dUTP法でStrand-specific RNA-seqのデータが得られているため)。-a
: アノテーション情報(GTFファイル)を指定。-o
: 出力先のパスとファイル名を指定。- 1つ目の引数: BAMファイルを指定。
featureCountsから出力されたリードの集計データのリストのうち、遺伝子IDとリード数の列だけを抽出します。
$ sed -e "1,2d" featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ | cut -f1,7 - > featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1_for_R.txt
1, 2行目はHeader行なので、sed
コマンドで除きます。次に、cut
コマンドで遺伝子IDとリード数の列だけを抽出して、データを標準出力で得ます。
edgeRを用いて2群間比較を行う
edgeR_test.R
スクリプトを使って、edgeRを用いた2群間比較を行います。
以下のURLに必要なスクリプトが置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
まず、mkdir
コマンドで保存先のディレクトリを用意します。
mkdir edgeR_result_UPF1_Knockdown
edgeR_test.R
スクリプトを用いて、2群間比較を行います。
Rscript ./edgeR_test.R \ featureCounts_result_SRR4081222_Control_1_ss/featureCounts_result_SRR4081222_Control_1_for_R.txt,\ featureCounts_result_SRR4081223_Control_2_ss/featureCounts_result_SRR4081223_Control_2_for_R.txt,\ featureCounts_result_SRR4081224_Control_3_ss/featureCounts_result_SRR4081224_Control_3_for_R.txt,\ featureCounts_result_SRR4081225_UPF1_knockdown_1_ss/featureCounts_result_SRR4081225_UPF1_knockdown_1_for_R.txt,\ featureCounts_result_SRR4081226_UPF1_knockdown_2_ss/featureCounts_result_SRR4081226_UPF1_knockdown_2_for_R.txt,\ featureCounts_result_SRR4081227_UPF1_knockdown_3_ss/featureCounts_result_SRR4081227_UPF1_knockdown_3_for_R.txt \ Control,Control,Control,Knockdown,Knockdown,Knockdown \ edgeR_result_UPF1_Knockdown \ Yes
- 1つ目の引数: 先ほど用意したリード数を集計した各データをカンマ区切りで指定する。
- 2つ目の引数: 1つ目の引数で指定したサンプルに対して、ラベル付け(例えば、
Control
とKnockdown
)を行う。データはカンマ区切りで指定する。 - 3つ目の引数: 保存先のディレクトリを指定する。
- 4つ目の引数: 複数のサンプル(n=2以上)を取っているかどうか(Duplicateのチェック)。
Yes
もしくはNo
と記入。
アノテーション情報を追加する
edgeRから出力されたデータでは、遺伝子IDしか情報がないので遺伝子名などの情報を付け加えます。
$ python2 ./annotate_gene_symbol_type.py \ ./gencode.v19.annotation_symbol_type_list.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt \ $ python2 ./split_into_each_gene_type.py \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt
annotate_gene_symbol_type.py
スクリプトでアノテーション情報を付加しています。
- 1つ目の引数: GTFファイルから作成した遺伝子リストを指定します。
- 2つ目の引数: 先ほどのスクリプトを使って得られた出力ファイルを指定します。
- 3つ目の引数: アノテーション情報を付加したデータを出力します。
split_into_each_gene_type.py
スクリプトでmRNAやlncRNAなどRNAの種類に応じてデータを分割します。
- 1つ目の引数: 上記のスクリプトから得られた出力ファイルを指定します。
以上になります。
HoloLensとパソコンの連携を簡単にする「Microsoft HoloLens」
概要
Windowsのアプリストアに「Microsoft HoloLens」という名前のHoloLensとパソコンの連携を簡単にする公式アプリがあります。
そのアプリのインストール方法とできることについてまとめたいと思います。
インストール
「Microsoft HoloLens」は日本版のアプリストアでは見つからないので、まず地域の設定を変更する必要があります。
スタートメニューから設定をクリックし、設定を開きます。 そこから、[時刻と言語] -> [地域と言語] -> [国または地域] -> [米国]を選択します。
地域を日本から米国に変更した後、Windowsのアプリで「Microsoft HoloLens」と検索します。すると、以下のようなアプリが見つかると思います。
あとは、アプリストアからこのアプリをインストールするだけです。
何ができるのか
自分の持っているHoloLensのIPアドレスを入力しておくことで、簡単にHoloLensと自分のPC間で連携ができます。Device Portalの簡易版だと考えてくれればいいと思います。
HoloLensとPCの連携
Addで自分のHoloLensを登録することができます。HoloLensのIPアドレスと入力して、Connectをクリックします。
すると、このようにOnLine状態になると思います。次回以降は、自分の設定したHoloLensのボタンをクリックするだけで、HoloLensとPCの連携が可能になります。
メニュー画面はこんな感じ。
Live Stream
ライブストリーミングを行うことができます。Device Portalのものよりも遅延が少ないそうです。誰かにHoloLens越しで見えている世界を見せたいときに使えると思います。
App manager
現在動作しているアプリの一覧が出ます。どこかの空間にアプリを開いたままになって見つからない場合に、ここから簡単に各アプリを終了することができます。
Device info
バッテリー残量やビルド番号などの情報が載っています。PC側からHoloLensをシャットダウンしたり再起動することができます。
その他
CameraでHoloLens側でのカメラ撮影や録画を制御することができます。
Photo&videosで写真や動画を管理することができます。
Virtual keyboardでPC側からキー入力することができます(一番いらない機能?)。
RNA-seqデータ解析(2群間比較, Unstranded, Single-end)
概要
RNA-seqのデータを用いて、siRNAによるノックダウンや特定の刺激により発現量が変動したRNA群を同定する方法(2群間のデータの比較)を紹介します。
ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、RNA-seqのデータ解析の詳細をStep by Stepで説明します。
ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。 パソコン(もしくはスパコン)のメモリは16GB以上を確保することが望ましいです。
使用するデータ
RNA-seq (Single-end, 101bp; Control siRNA, UPF1 siRNA)
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148
使用するソフトウェア
- FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - Bowtie1 (v1.1.2)
http://bowtie-bio.sourceforge.net/index.shtml - TopHat (v2.1.1)
https://ccb.jhu.edu/software/tophat/index.shtml - Samtools (v1.3.1)
http://samtools.sourceforge.net/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - bedtools (v2.26.0)
http://bedtools.readthedocs.io/en/latest/ - Cufflinks (v2.2.1)
http://cole-trapnell-lab.github.io/cufflinks/ - featureCounts (subread v1.5.0.post3)
http://subread.sourceforge.net/ - edgeR (v3.14.0)
https://bioconductor.org/packages/release/bioc/html/edgeR.html
解析ワークフロー
RNA-seqデータを2群間比較する方法として、Cuffdiffを用いた方法とfeatureCount-edgeRを用いた方法を紹介したいと思います。
TopHat-Cuffdiff workflow
TopHat-featureCounts-edgeR workflow
サンプルコード
Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/RNA-seq_Tutorial
RNA-seq_unstr_single-end_1_mapping.sh
TopHatによるマッピング・featureCountsによるリードカウントまでを行うスクリプト。FASTQファイルを引数として指定。
$ ./RNA-seq_unstr_single-end_1_mapping.sh hoge.fastq
10行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
11行目: indexContamFile
にrRNAなどのコンタミ配列のBowtie用のIndexファイルを指定。
12行目: indexGenomeFile
にGenome配列のBowtie用のIndexファイルを指定。
31行目: -r
オプションでHouse-keeping genesの遺伝子領域のBEDファイルを指定。
RNA-seq_unstr_single-end_2_cuffdiff_calc.sh
Cuffdiffによる2群間比較を行うスクリプト。
$ ./RNA-seq_unstr_single-end_2_cuffdiff_calc.sh
8行目: filename
に保存先のディレクトリ名を指定。
9行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
13-14行目: カンマ(,
)区切りでBAMファイル(比較したい2種類)を指定。
17行目: gene_list
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
RNA-seq_unstr_single-end_2_edgeR_calc.sh
edgeRによる2群間比較を行うスクリプト。
$ ./RNA-seq_unstr_single-end_2_edgeR_calc.sh
8行目: annoList
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
11行目: saveDir
に保存先のディレクトリ名を指定。
~/custom_command/<hoge.py>
となっている箇所はカスタムのPythonスクリプトが置いてあるパス・ファイル名を指定する。
シェルスクリプトを使用するときの注意点
権限の問題
実行権限をシェルスクリプトに与えた後、実行すること。
$ chmod 755 hoge.sh
改行コードの問題
WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。
ここではAtom
と呼ばれるテキストエディタでの設定方法を説明する。
- [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
- 左のメニューから、
Packages
をクリック。 - Installed Packagesで
line-ending-selector
と検索。 - Core Packagesの項目から、line-ending-selectorの
Settings
をクリック。 - 設定画面で、Default line endingを
LF
に変更。
これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。
各ステップの説明
0. 必要なソフトウェアのインストール
まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。
$ 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 cufflinks $ conda install subread $ conda install bioconductor-edger
注意!!
tophatによるマッピングで最新のbowtieを使うとエラーが出るので、以前のバージョンのbowtie v1.12を代わりにインストールする。
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
という項目にリストアップしてあります。
今回は、35サンプルのうち上6つのサンプルを使用します。
GSMから始まるIDをクリックすると、各サンプルの情報を見ることができます。 サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。
そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。
SRR4081222をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。
SRR4081222.sra
のURLをコピー。
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以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。
Per tile sequence quality
フローセルの各タイルごとのクオリティスコアを示しています。
Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。
シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。
特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。
詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf
Per sequence quality scores
各リードの全長のクオリティスコアの分布を示しています。
Per base sequence content
各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。
Per sequence GC content
リードのGC contentsの分布を示しています。
Per base N content
各塩基中に含まれるN
の含有率(塩基を読めなかった箇所)を示しています。
Sequence Length Distribution
リード長の分布を示しています。
Sequence Duplication Levels
Duplidate readsの含まれている数を示しています。
Overrepresented sequences
頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。
Adapter Content
各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。
Kmer Content
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
3. Quality filtering
FASTQファイル中にクオリティスコアの低いリードが含まれていることが確認された場合、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。
今回のデータでは、クオリティの低いリードが含まれていないのでこの工程をパスします。
Illumina製のシークエンサーでいうと、HiSeq以降のデータはデータのクオリティが非常に高いので、いちいちクオリティフィルターにかける必要がなくなりました。
一応、実行例を示すと以下の通りになります。
実行例
$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./SRR4081222_Control_1.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o ./SRR4081222_Control_1_1_filtered.fastq
fastq_quality_trimmerのパラメータ
-Q
: Phred quality scoreのバージョン指定(Solexa, Illumina 1.3+, Illumina 1.5+といったかなり古いバージョンのFASTQファイルでは'64'を指定。最近のFASTQファイルは'33'と指定しておけばOK)。-t
: 指定したPhred quality scoreを下回る塩基をトリミングする。-l
: トリミング後の最低リード長。
fastq_quality_filterのパラメータ
-q 20 -q 80
: リードの80%
でPhred quality scoreが20
を下回る場合、そのリードを排除する。
4. Quality check
フィルタリング後のFASTQファイルをFastQCを用いて再確認します(今回のデータはフィルタリングを行わないのでこの工程もパス)。
実行例
$ mkdir fastqc_SRR4081222_Control_1_filtered $ fastqc -o ./fastqc_SRR4081222_Control_1_filtered ./SRR4081222_Control_1_1_filtered.fastq -f fastq
5. Mapping
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ファイル)に対して行います。
実行例
$ tophat --bowtie1 -p 8 -o ./tophat_out_SRR4081222_Control_1 \ -G ./gencode.v19.annotation.gtf \ ./hg19 ./SRR4081222_Control_1_1_filtered.fastq
--bowtie1
: マッピングソフトにBowtie1
を使用。-p
: 使用するコア数を指定。-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
6. データの品質チェック
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
というファイルが出力されます。ファイルの中身を見るとこんな感じ。
RNAの品質に問題がある(例えば、RNAが壊れている)と、3'側にリードが偏っていたり、cDNAライブラリの構築に使ったRNA量が少なく、PCR時に特定の配列が異常に増幅され(PCR bias、Duplicate readsが多い)、ガタガタしたCoverageの分布を取ることがあります。
こういったデータはRNAの定量性に問題があるため、使用しないように気をつけましょう。
このステップで問題がありそうだと判断されるデータに関しては、下記のUCSC genome browserでの可視化により、個々の遺伝子上にマッピングされたリードの分布に注意して見ていく必要があります。
7. Visualization (For UCSC genome browser)
BedGraphファイルの準備
Bedtools
を利用して、TopHatから得られたBAMファイルを、UCSC genome browser上で可視化できるファイル(Bedgraphファイル)に変換します。
まず、mkdir
コマンドでデータを保存するディレクトリを作っておきます。
$ mkdir UCSC_visual_SRR4081222_Control_1
bedtools genomecovで特定のゲノム領域でのCoverageを計算します(BAMファイルからBedgraph
ファイルを作成)。
Bedgraphファイルの詳細については、下記のURLを参照のこと。
https://genome.ucsc.edu/goldenpath/help/bedgraph.html
$ bedtools genomecov -ibam ./tophat_out_SRR4081222_Control_1/accepted_hits.bam -bg -split \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result.bg
-ibam
: BAMファイルを入力ファイルとして指定する。-bg
: 出力ファイルの形式をBedGraphファイルに指定。-split
: Splicing juctionにマッピングされたリード(リードが分割されてマッピングされているリード)を考慮する。
標準出力でBedGraph
ファイルが出力されます。
UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。
https://genome.ucsc.edu/goldenpath/help/customTrack.html
それらの情報をヘッダー行(ファイルの一行目)に記載します。
$ echo "track type=bedGraph name=RNA-seq_Control_1 description=RNA-seq_Control_1 visibility=2 maxHeightPixels=40:40:20" \ > ./UCSC_visual_SRR4081222_Control_1/tmp.txt
echo
コマンドを使って、ヘッダー行の情報をtxtファイルとして出力しています。
そのあとに、cat
コマンドを用いて、先ほど作ったBedGraphファイルと結合させます。
$ cat ./UCSC_visual_SRR4081222_Control_1/tmp.txt ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result.bg \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg
出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2ファイルに圧縮します。
$ bzip2 -c ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg.bz2
データのアップロード
UCSC genome browserにデータをアップロードします。
ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。
Add Custom Tracksというページに移動するので、ファイルを選択をクリックし、アップロードするファイルを選択し、submitボタンをクリックします。
アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgoボタンをクリックすると、データをみることができます。
8. RNA-seqデータの2群間比較(Cuffdiffのケース)
Cuffdiffを用いて、siCTRLとsiUPF1の2群間のRNA-seqデータを比較します。
cuffdiff -p 8 --multi-read-correct -o ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode \ /home/akimitsu/database/gencode.v19.annotation_filtered.gtf \ ./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
-p
: 使用するコア数を指定。--multi-read-correct
: 定量の精度を上げる目的で、複数箇所にマッピングされたリードに重み付けを行う。-o
: 出力先のディレクトリを指定。- 1つ目の引数: アノテーション情報のGTFファイルを指定。
- 2つ目・3つ目の引数: 2種類のBAMファイルをそれぞれ指定。今回のように、siCTRLとsiUPF1の条件でn=3でデータがある場合、それぞれのBAMファイルのパスをカンマ区切りで並べる(例えば、
CTRL_1.bam,CTRL_2.bam,CTRL_3.bam KD_1.bam,KD_2.bam,KD_3.bam
)。
パラメータの詳しい解説は以下のサイトを参照のこと。
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html
次に、Gene symbolなどの情報をCuffdiffから出力されるファイルに加える。
遺伝子リストの準備
まず、extract_gene_symbol_type_from_gtf.py
スクリプトを用いて、GTFファイルから遺伝子リストを作成します(Cuffdiffで集計したデータからmRNAとlncRNAに分類分けするために用います)。
以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
$ python2 extract_gene_symbol_type_from_gtf.py gencode.v19.annotation_filtered.gtf gencode.v19.annotation_filtered_symbol_type_list.txt
- 1つ目の引数: GTFファイルを指定する。
- 2つ目の引数: 遺伝子リストが出力される。
次に、先ほど作成した遺伝子リストをもとに、cuffdiffの出力ファイルからmRNAとlncRNAの情報を抽出し、それぞれ別のファイルとして保存する。
# mRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ mRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_mRNA.txt # lncRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ lncRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_lncRNA.txt
- 1つ目の引数: 先ほど用意した遺伝子リストを指定。
- 2つ目の引数: cuffdiffの出力ファイルである
gene_exp.diff
ファイルを指定。 - 3つ目の引数:
mRNA
もしくはlncRNA
を指定(指定したRNA群のデータを取得)。 - 4つ目の引数: mRNAもしくはlncRNAのリストを出力。
8. RNA-seqデータの2群間比較(featureCounts-edgeRのケース)
各遺伝子に対してのリード数のカウント
まず、featureCountsで各遺伝子領域にマッピングされたリード数をカウントします。
$ mkdir featureCounts_result_SRR4081222_Control_1 $ featureCounts -T 8 -t exon -g gene_id -s 2 -a ./gencode.v19.annotation_filtered.gtf \ -o featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ ./tophat_out_SRR4081222_Control_1/accepted_hits.bam
mkdir
コマンドを使って、保存先のディレクトリを作り、そのあとfeatureCountsを実行します。
-T
: 使用するコア数の指定。-t exon
: リードのカウントに使用するExonの情報を指定。ここでは、GTFファイル内で3列目のfeature
がexon
である行を使用するように指定している。-g gene_id
: 各行のメタ情報(遺伝子ID、Transcript ID、遺伝子名など)のうち、名前として使用する情報を指定する。ここでは、gene_id
を指定することで、Gencodeで決められた遺伝子IDを使っている。-a
: アノテーション情報(GTFファイル)を指定。-o
: 出力先のパスとファイル名を指定。- 1つ目の引数: BAMファイルを指定。
featureCountsから出力されたリードの集計データのリストのうち、遺伝子IDとリード数の列だけを抽出します。
$ sed -e "1,2d" featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ | cut -f1,7 - > featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1_for_R.txt
1, 2行目はHeader行なので、sed
コマンドで除きます。次に、cut
コマンドで遺伝子IDとリード数の列だけを抽出して、データを標準出力で得ます。
edgeRを用いて2群間比較を行う
edgeR_test.R
スクリプトを使って、edgeRを用いた2群間比較を行います。
以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
まず、mkdir
コマンドで保存先のディレクトリを用意する。
mkdir edgeR_result_UPF1_Knockdown
edgeR_test.R
スクリプトを用いて、2群間比較を行う。
Rscript ./edgeR_test.R \ featureCounts_result_SRR4081222_Control_1_ss/featureCounts_result_SRR4081222_Control_1_for_R.txt,\ featureCounts_result_SRR4081223_Control_2_ss/featureCounts_result_SRR4081223_Control_2_for_R.txt,\ featureCounts_result_SRR4081224_Control_3_ss/featureCounts_result_SRR4081224_Control_3_for_R.txt,\ featureCounts_result_SRR4081225_UPF1_knockdown_1_ss/featureCounts_result_SRR4081225_UPF1_knockdown_1_for_R.txt,\ featureCounts_result_SRR4081226_UPF1_knockdown_2_ss/featureCounts_result_SRR4081226_UPF1_knockdown_2_for_R.txt,\ featureCounts_result_SRR4081227_UPF1_knockdown_3_ss/featureCounts_result_SRR4081227_UPF1_knockdown_3_for_R.txt \ Control,Control,Control,Knockdown,Knockdown,Knockdown \ edgeR_result_UPF1_Knockdown \ Yes
- 1つ目の引数: 先ほど用意したリード数を集計した各データをカンマ区切りで指定する。
- 2つ目の引数: 1つ目の引数で指定したサンプルに対して、ラベル付け(例えば、
Control
とKnockdown
)を行う。データはカンマ区切りで指定する。 - 3つ目の引数: 保存先のディレクトリを指定する。
- 4つ目の引数: 複数のサンプル(n=2以上)を取っているかどうか(Duplicateのチェック)。
Yes
もしくはNo
と記入。
アノテーション情報を追加する
edgeRから出力されたデータでは、遺伝子IDしか情報がないので遺伝子名などの情報を付け加えます。
$ python2 ./annotate_gene_symbol_type.py \ ./gencode.v19.annotation_symbol_type_list.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt \ $ python2 ./split_into_each_gene_type.py \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt
annotate_gene_symbol_type.py
スクリプトでアノテーション情報を付加しています。
- 1つ目の引数: GTFファイルから作成した遺伝子リストを指定します。
- 2つ目の引数: 先ほどのスクリプトを使って得られた出力ファイルを指定します。
- 3つ目の引数: アノテーション情報を付加したデータを出力します。
split_into_each_gene_type.py
スクリプトでmRNAやlncRNAなどRNAの種類に応じてデータを分割します。
- 1つ目の引数: 上記のスクリプトから得られた出力ファイルを指定します。
(おまけ)CuffdiffとedgeRの結果を比較する
遺伝子名(ヘッダー行を除く)でソートする。
(head -n +1 cuffdiff_result_mRNA.txt && tail -n +2 cuffdiff_result_mRNA.txt | sort -k 3 ) > cuffdiff_result_mRNA_sorted.txt
(head -n +1 edgeR_test_result_anno_plus_mRNA.txt && tail -n +2 edgeR_test_result_anno_plus_mRNA.txt | sort -k 3 ) > edgeR_test_result_anno_plus_mRNA_sorted.txt
Excel上でマージして比較する。
CuffdiffとedgeRから得られたDEG比較
CuffdiffとedgeRのFold-changeの比較
発現量が低い遺伝子については、cuffdiffとedgeRの結果にばらつきが見られる。
CuffdiffとedgeRの使い分けについて
結局どちらを使ったらよいかという話ですが、好みの問題なのでどちらでも構わないというのが私の答えです(どちらかと言えば、edgeRを勧めますが)。
edgeRを使用する場合、Rのスクリプトを少し書かなければいけないので、プログラミングがわからない人にとっては少しハードルが高いのかもしれません。
その点では、Cuffdiffはコマンドを打つだけでOKなので楽です(プログラミングゼロで結果を出せる!)。
ただ、Cuffdiffに関しては少しクセのあるソフトなので、使用するときには注意が必要だと考えています。
あと、Cuffdiffは他の手法と比較してFalse-positiveが多いという解析結果も出ています(どのデータでも同様の傾向があるかどうかはわかりませんが)。
https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-9-r95
また計算速度もCuffdiffが数時間かかるのに対して、edgeRだと数分で終了します(記事を書いてて、なんかedgeRがいい気がしてきた)。
今までの経験に基づく個人的な意見を以下にまとめました。
mRNAのみを解析対象としているならCuffdiff・edgeRを選択
mRNAのみを解析対象にする場合、Cuffdiffの結果はかなり信頼性が高い(qPCRなどで再現がほぼ取れる)と実感しています。
まぁ、edgeRを使用してもかまわないと思いますが…。
long-noncoding RNAを解析対象としているならedgeRを選択
long-noncoding RNA(lncRNA、PROMPTやeRNAを含む)を解析対象にしている場合、Cuffdiffだとおかしな結果が返ってくることが多々あります。
qPCRで再現が全く取れなかったり、そもそも集計されたリード数とUCSC genome browserで可視化したときのCovarageの具合がどうも一致していない気がする(カウントしてくれたリード数がおかしいんですが、Cuffdiffさん…)などなど…。
GC-contentsが高かったり、Repetitive elementsを含んでいたり、Transcript長が短かったり(1kb未満)するとどうもCuffdiffがうまく計算してくれない気がします(理由の1つとして、Cuffdiff内で行わるBiasの補正が問題なのではと思っていますが、詳しくはよくわかっていません…)。
なので、lncRNAを解析しているなら、CuffdiffではなくedgeRを使うようにしましょう。
Pythonパッケージ: 便利なコマンドラインパーサClickを使ってみる
概要
コマンドラインパーサとして、argparse
と呼ばれるPythonの標準パッケージが存在するが、使い勝手がいまいちだと感じています(あくまで個人的な感想ですが…)。
今回は、たまたま見つけたclick
と呼ばれるサードパーティ製のPythonパッケージを用いてコマンドラインの引数を解釈するコードを書いてみようと思います。
argparse
と比較してclick
の書き方だと、コードが書きやすく、あとから読み返しやすくなると感じます。
Clickのドキュメントはこちら。
http://click.pocoo.org/6/
内容
すでに、こちらのブログに詳しく解説されているのでここでは忘備録的な記載に留めようと思います。
インストール
AnacondaもしくはBiocondaがインストールされていることが前提。
conda install click
基本的な使い方
@click.command()
デコレータを用いて関数を修飾するだけでOK。
print
関数の代わりにclick.echo
を用いて文字列を標準出力することもできる。
@click.command() def cmd(): click.echo("Command was Done.")
続けて、@click.option()
デコレータに任意のオプションを追加していくことができる。argparseと同様に、-i
や--ifastq
などのような-
から始まるコマンドライン引数を指定することができる。
@click.command() @click.option('-i', '--ifastq') def cmd(): click.echo("Command was Done.")
また、作成したオプションに代入された値を取得するために、ifastq
のように-
を付けない形で文字列を指定することで、関数の中で値が代入された変数として扱えます。
以下のように、cmd
関数の引数にifastq
を指定しておく必要がある。
@click.command() @click.option('-i', '--ifastq', 'ifastq') def cmd(ifastq): click.echo(ifastq)
help
でHelpが呼び出されたときの説明文を指定できる。
@click.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') def cmd(ifastq): click.echo(ifastq)
実行例
#/usr/bin/env python import click @click.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def cmd(ifastq, output): click.echo(ifastq) def main(): cmd() if __name__ == '__main__': main()
実行すると、
$ python test.py -i test.fastq test.fastq $ python test.py --help Usage: test.py [OPTIONS] Options: -i, --ifastq TEXT Input file in FASTQ file. -o, --output-prefix TEXT Prefix of output file. --help Show this message and exit.
コマンドライン引数のタイプを指定したい
INT型やFLOAT型など、引数のタイプがあらかじめわかっている場合、引数のタイプを指定しておきたいケースが出てくる。
そういった場合、以下のようにargparseのときのように@click.option()
デコレータ内で指定することが可能。
type=<type>
で指定。
@click.option('-m', '--minimum-length', 'minimum_length', type=int, default=18, help='Discard trimmed reads that are shorter than LENGTH. Default: 18')
パラメータタイプは、int、float、bool、strが指定できるようです。
http://click.pocoo.org/6/parameters/
デフォルトの値を決めておきたい
特定のオプションでデフォルト値を指定しておきたいケースがある。 その場合、デフォルト値の設定は以下のように行う。
default='hogehoge'
で指定。
@click.option('-nl', '--nucleotide-trimming-list', 'nucleotide_trimming_list', default='A,N', help='Trimmed nucleotide list. Default: A,N')
指定する引数を限定したい
オプションの引数として、特定の文字列や数値しか受け付けたくない場合、以下のように指定できる引数を限定することができる。
type=click.Choice(['hoge', 'huga'])
で指定。
@click.option('-nd', '--nucleotide-trimming-direction', 'nucleotide_trimming_direction', type=click.Choice(['three', 'five']), help='Trimmed nucleotide direction. Default: three')
Helpパラメータをカスタマイズしたい
Clickのデフォルトでは、Help文は--help
を指定することで出力される。
一般的に、Help文は--help
だけでなく-h
でも出力されるケースが多い。そこで、-h
でもHelp文が出力されるように設定を変更する。
@click.command()
デコレータ内で指定。
CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @click.command(context_settings=CONTEXT_SETTINGS)
サブコマンド(後述)を利用している場合は、@click.group
デコレータ内で指定する。
CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @click.group(context_settings=CONTEXT_SETTINGS)
コールバック関数で値を検知・Helpを表示
オプションの引数を評価したい場合、コールバック関数を利用すると便利です。
@click.option()
デコレータ内でcallback=<function名>
を指定する。
また、コールバック関数内で、ctx.get_help()
でHelp文を表示させることができる。ctx.exit()
でPythonスクリプトの実行を止めることができる。
def validate_ifastq(ctx, param, value): if not value: click.echo('Error: Do not choose your FASTQ file.\n') print(ctx.get_help()) ctx.exit() return(value) @cmd.command() @click.option('-i', '--ifastq', 'ifastq', callback=validate_ifastq, help='Input file in FASTQ file. [required]')
サブコマンドを用意したい場合
@click.group()
デコレータを用いて、サブコマンドを利用することができます。
#/usr/bin/env python import click @click.group() def cmd(): pass @cmd.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def trim_polya(ifastq, output): click.echo('Start trimming polyA sequence from reads...') @cmd.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def get_pass_read(ifastq, output): click.echo('Start getting polyA site-supporting reads...') def main(): cmd() if __name__ == '__main__': main()
Help文は以下のとおり。
$ python test.py --help Usage: test.py [OPTIONS] COMMAND [ARGS]... Options: --help Show this message and exit. Commands: get_pass_read trim_polya
サブコマンドごとのオプションをHelp文で出力することもできる。
$ python test.py get_pass_read --help Usage: test.py get_pass_read [OPTIONS] Options: -i, --ifastq TEXT Input file in FASTQ file. -o, --output-prefix TEXT Prefix of output file. --help Show this message and exit.
ENCODEプロジェクトのNGSデータを活用する
概要
ENCODEプロジェクトとは、ポストゲノム研究戦略としてヒトゲノムの機能性エレメントの同定とその全体像を解明するために組織された国際プロジェクトの1つです。
RNA-seq、ChIP-seq、DNase-seq、Hi-C、ChIA-PET、eCLIP-seqなどの解析手法を用いて、毎月膨大なデータを産出しています。
ENCODEプロジェクトでは、共通のプロトコルに基づいた実験と解析によりデータが産出され、即座に一般公開される仕組みを取っています。また、プロトコルも公開されているので自身で内容を確認・再解析することもできます。
ENCODEプロジェクトのデータは、プロジェクトに参加していないラボ・研究機関であっても、ウェブサイト上で公開されているNGSデータに関しては、自由にダウンロードして解析し論文のデータとして使用することもできます。
Data Use, Software, and Analysis Release Policies – ENCODE
実際に、ENCODEのデータを利用して投稿された論文はこれまで1700本にのぼります。
Publications using ENCODE data – ENCODE
ENCODEプロジェクトのデータベースには、有用な情報やデータがいっぱい公開されているのでこれらを使わない手はありません。
そこで今回は、ENCODEプロジェクトに登録されているデータを活用して手軽にデータ解析する方法をまとめました。今回は主に、RNA結合タンパク質(RBP)が結合するRNA領域を網羅的に同定する方法であるeCLIP-seqデータについて見ていきたいと思います。
また、RBPをIPした時に使用した抗体の情報など即座に使える情報も載っているのでそれらについても紹介しようと思います。
欲しいデータを探す
キーワードで検索
まず、ENCODEプロジェクトのウェブサイトに行きます。
ENCODE: Encyclopedia of DNA Elements – ENCODE
右上の検索欄に興味のあるRBPの名前を入力します。
試しにSFPQ
というRBPで検索にかけると、eCLIP-seq、Bind-n-seq、RNA-seq等のデータや、IPに使用した抗体の情報が出てきます。
データの種類から検索
他にも、[Encyclopedia]->[about]から各種データにアクセスできます。
この他にも、いろいろな検索方法が用意されていますが割愛します。
抗体の情報
右端にAntibody
と記載されているのが、抗体の情報です。
その中でも、緑色のマークがついているものは、ENCODEでStandardな抗体として認められた抗体です(抗体の特異性など幾つかの基準をもとに選定されている)。
クリックすると、その抗体でIPしたときのWestern blottingのデータや、shRNAでノックダウンしたときに目的のバンドが減弱しているかどうかWestern blottingで確認したデータなどが公開されています。
以上の情報から、各メーカーの抗体の良し悪しを知ることができます。 ちなみに、以下がENCODEプロジェクトで確認した抗体に関する論文になります。
http://yeolab.github.io/papers/2016/sundararaman_molcell_2016.pdf
Bind-n-seqデータ
Bind-n-seqについては、下記の論文を参照のこと。方法の詳細については、ここでは割愛します。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4142047/
リコンビナントタンパク質のQCデータと、Enrichされたモチーフ配列の情報が載っています。
さらに、Processed data
の項目のところに、Enrichされた4, 5, 6, 7ntのモチーフ配列のデータがそれぞれ置いてあります(output typeがenrichment
となっているやつ)。
ダウンロードボタン(Accession
にある小さいアイコン)をクリックすると、該当のファイルをダウンロードできます。
同じページにRawデータ(FASTQファイル)とプロトコルが置いてあるので、自身の手で再解析することも可能です。
eCLIP-seqデータ
1. 必要なファイルをダウンロード
eCLIPのデータを選択し、該当するデータのページに飛ぶと下のほうにProcessed data
の項目があります。
ゲノム上にマッピングされたリードを可視化するためのデータとしてbigWig
ファイルを、各Peakの染色体座標軸に関するデータとしてbed
ファイルをそれぞれダウンロードします。
また、右上のGRCh38
(hg38)をhg19
に変更することで、異なるバージョンのヒトのリファレンスゲノムをもとにして作成したファイルにアクセスできます。
今回は、hg38 (GRCh38)
のデータをダウンロードします。
$ wget https://www.encodeproject.org/files/ENCFF823MPG/@@download/ENCFF823MPG.bigWig $ wget https://www.encodeproject.org/files/ENCFF960OTE/@@download/ENCFF960OTE.bed.gz $ gunzip ENCFF960OTE.bed.gz
ファイル名がわかりにくいのでリネームします。
$ mv ENCFF823MPG.bigWig SFPQ_eCLIP_rep1_ENCFF823MPG.bigWig $ mv ENCFF960OTE.bed SFPQ_eCLIP_rep1_ENCFF960OTE.bed
2. ファイルフォーマットの変換
bedgraphファイルの用意
bigwig
ファイルではUCSC genome browserにアップロードできないので、bedgraph
ファイルと呼ばれる形式にフォーマット変換します。
変換作業を行うために、今回はbigWigToBedGraph
というプログラムを使います。
下記のURLから、目的のプログラムのバイナリファイルを入手します。
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/
$ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigToBedGraph $ chmod 755 bigWigToBedGraph
ダウンロードした後、実行権限を与え、パスの通っているディレクトリに入れます。
続いて、bigWig
ファイルをbedgraph
ファイルにフォーマット変換します。
# bigWigToBedGraph <bigWigファイル> <bedGraphファイル> $ bigWigToBedGraph SFPQ_eCLIP_rep1_ENCFF823MPG.bigWig SFPQ_eCLIP_rep1_ENCFF823MPG.bg
- 1つ目の引数: フォーマット変換したい
bigwig
ファイルを指定。 - 2つ目の引数: 出力される
bedGraph
ファイルの名前を指定。
次に、UCSC genome brower用にデータを整えます。ENCODEのデータの中にはハプロタイプのシークエンスにもマッピングしたデータが含まれているので、UCSC genome browserでデータを可視化するために、Chr1-22, chrX, chrY, chrM
からなるCanonical genomeにマッピングされたデータのみを抽出する必要があります。
下記のような適当なスクリプトを書いて、データを整えます。
$ python ENCODE_eCLIP_mod.py SFPQ_eCLIP_rep1_ENCFF823MPG.bg SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC.bg
最後に、UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。それらの情報をヘッダー行(ファイルの一行目)に記載します。
$ echo "track type=bedGraph name=SFPQ_eCLIP_rep1_ENCFF823MPG description=SFPQ_eCLIP_rep1_ENCFF823MPG visibility=2 maxHeightPixels=40:40:20 color=255,0,0" > ./header_tmp.txt $ cat header_tmp.txt SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC.bg > SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg $ bzip2 -c SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg > SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg.bz2
echo
コマンドのところでヘッダー行の情報をtxtファイルとして出力しています。
そのあとに、cat
コマンドを用いて、先ほど作ったBedGraphファイルと結合させます。
出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2
ファイルに圧縮しています。
bedファイルの用意
ENCODEのサイトから得られるbedファイルを、bigWigファイルと同様にUCSC genome browserにアップロードできるファイルに加工します。
下記のような適当なスクリプトを書いて、データを整えます。
$ python ENCODE_eCLIP_mod_bed.py SFPQ_eCLIP_rep1_ENCFF960OTE.bed SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC.bed
最後に、UCSC genome browserへデータをアップロードするために、ファイルの種類やTrackの名前の情報をヘッダー行(ファイルの一行目)に記載します。ここも先程と同様です。
$ echo "track type=bed name=SFPQ_eCLIP_rep1_ENCFF960OTE_peak description=SFPQ_eCLIP_rep1_ENCFF960OTE_peak visibility=2 maxHeightPixels=40:40:20 color=255,0,0" > ./header_tmp.txt $ cat header_tmp.txt SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC.bed > SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed $ bzip2 -c SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed > SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed.bz2
データのアップロード
UCSC genome browserにデータをアップロードします。
ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。
Add Custom Tracksというページに移動するので、ファイルを選択
をクリックし、アップロードするファイルを選択し、submit
ボタンをクリックします。
アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgo
ボタンをクリックすると、データをみることができます。
可視化したデータをみてみると、以下のようになります。
ざっくりとですが、以上になります。