Ribo-seqデータ解析 (RiboTaper, short ORF prediction)
概要
Ribo-seq (Ribosome Profiling)とは、リボソームが結合しているRNA領域をシークエンスする手法です。この方法を利用することで、mRNAのどの領域が翻訳されているか、つまりはORFを予測することができます。また、特定の刺激などによる翻訳活性の変化を捉えることができます。
ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、Ribo-seqのデータ解析の詳細をStep by Stepで説明します。ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。
使用するデータ
SRA160745 (Gao X et al. Nat. Methods 2015)
https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=SRA160745
使用するソフトウェア
- SRA Toolkit (v.2.8.0)
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software - FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - Cutadapt (v1.12)
http://cutadapt.readthedocs.io/en/stable/guide.html - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - STAR (v2.5.2b)
https://github.com/alexdobin/STAR - samtools (v0.1.19)
samtools.sourceforge.net/ - bedtools (v2.17.0)
http://bedtools.readthedocs.io/en/latest/ - RiboTaper (v1.3)
https://ohlerlab.mdc-berlin.de/software/RiboTaper_126/
解析ワークフロー
今回の解析の流れを以下にまとめました。
![PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-03-30 10.21.48.png (140.1 kB)](https://img.esa.io/uploads/production/attachments/3909/2017/03/30/11755/24ea6ea0-4d5d-4ae1-b1b5-89aa5e4c0089.png)
サンプルコード
Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Ribo-seq_RiboTaper
シェルスクリプトを使用するときの注意点
権限の問題
実行権限をシェルスクリプトに与えた後、実行すること。
$ 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のインストールに関しては、下記の記事を参照のこと。
まず、Biocondaで仮想環境をつくります。ここでは、python2.7をベースに仮想環境を構築します。
conda create -n ribotaper python=2.7
-nで仮想環境の名前を指定します。pythonに2.7を指定することで、python2.7を仮想環境にインストール。
source activateに続けて、先ほど作成した仮想環境の名称を指定することで、仮想環境に入ることができます。
source activate ribotaper
conda info --envsで現在どの環境にいる確認できるので、前後でどの環境にいるか確認してみましょう。
$ conda info --envs # conda environments: # ribotaper /home/akimitsu/miniconda2/envs/ribotaper root * /home/akimitsu/miniconda2 $ source activate miso $ conda info --envs # conda environments: # ribotaper * /home/akimitsu/miniconda2/envs/ribotaper root /home/akimitsu/miniconda2
最初はrootという環境にいたのが、source activate ribotaperの後、ribotaperという仮想環境に移動していることがわかると思います。
RiboTaperの実行は、このribotaperと名付けた仮想環境上で行うこととします。また、その他の必要なソフトウェアに関しても同様の環境上にインストールしていきます。
それでは、RiboTaperに必要なパッケージやソフトウェアをBiocondaを利用してあらかじめインストールしておきましょう。
注意!!
samtoolsとbedtoolsはバージョンに依存しているので、最新のものをインストールしてはならない。
$ conda install fastqc $ conda install cutadapt $ conda install fastx_toolkit $ conda install star $ conda install samtools $ conda install rseqc $ conda install bedtools $ conda install subread # For RiboTaper $ conda install samtools=0.1.19 $ conda install bedtools=2.17.0 $ conda install r-seqinr # conda install r-ade4 $ conda install r-multitaper $ conda install r-xnomial # DOWNGRADED due to dependency conflicts (R(v3.2.2-0)) $ conda install r-domc # conda install r-iterators # conda install r-foreach
RiboTaperのソースコードをダウンロードして、ファイルを解凍します。
$ wget https://ohlerlab.mdc-berlin.de/files/RiboTaper/RiboTaper_v1.3.tar.gz $ tar zxvf RiboTaper_v1.3.tar.gz
RiboTaper_v1.3のディレクトリにパスを通します。下記のように、ホームディレクトリにある.bashrcにパスを書き加えます。
export PATH=/home/akimitsu/software/RiboTaper_v1.3/scripts:${PATH}
1. SRAダウンロード
NGSデータはシークエンスされたリード（塩基配列）の情報とクオリティスコアのデータをまとめたFASTQ形式のファイルとして保存されています。
FASTQファイルの詳細についてはこちらを参照。
https://en.wikipedia.org/wiki/FASTQ_format
FASTQデータは非常に大きいファイルなので、NCBI GEOでは、FASTQファイルを圧縮したSRAファイルと呼ばれるファイル形式でデータを公開しています。
まず、下記のURLにアクセスすると、
https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=SRA160745
DDBJ DRAのサイトはこんな感じ。NCBI GEOと比較してIDが羅列してあるだけなので、どのIDがどのサンプルに対応しているかぱっと見わかりません…。
下記のサイトでDRAに登録されているデータの詳細をリスト化してみることができます。
DBCLS SRA: Survey of Read Archives
http://sra.dbcls.jp/v1/
早速、SRA#の欄にSRA160745と入力して、検索してみましょう。
Expsの項目をクリックします。
すると、各IDに対してタイトルなどの情報を見ることができます。
今回使用するのは、No16とNo19のデータになります。Exp#の項目のSRX740748とSRX740751をそれぞれクリックします。
DDBJ DRAのページに移動します。NavigationのRunのSRAボタンをクリックします。
SRR1630831のURLをコピーします、
wgetで該当するSRAファイルをダウンロードします。
$ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/SRX/SRX740/SRX740748/SRR1630831/SRR1630831.sra $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/SRX/SRX740/SRX740751/SRR1630838/SRR1630838.sra
2. SRAファイルからFASTQファイルに変換
まず、SRA IDだけだったファイル名にサンプルの情報を付け加えます（SRRのIDだけだと何のファイルかわからなくなるので）。mvコマンドを使ってリネームを行います。
$ mv SRR1630831.sra Ribo-Seq_HEK293_Cell_Control_SRR1630831.sra $ mv SRR1630838.sra RNA-Seq_HEK293_Cell_SRR1630838.sra
リネーム後、先ほどインストールしたSRA Toolkitを使って、NCBI GEOからダウンロードしたSRAファイルをFASTQファイルに展開します。
$ fastq-dump Ribo-Seq_HEK293_Cell_Control_SRR1630831.sra $ fastq-dump RNA-Seq_HEK293_Cell_SRR1630838.sra
カレントディレクトリにFASTQファイルが展開されているはずです。
3. Quality check
FastQCを使って、FASTQファイルに記載されている各リードのクオリティをチェックします。
FASTQファイルはPhred Quality Scoreと呼ばれる指標によって、各塩基のクオリティを記述しています。
詳しくはこちらのサイトを参照のこと。
https://en.wikipedia.org/wiki/Phred_quality_score
実行例
$ mkdir fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831 $ fastqc -t 8 -o ./fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831 ./Ribo-Seq_HEK293_Cell_Control_SRR1630831.fastq -f fastq
mkdirコマンドでfastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831という名前のディレクトリをカレントディレクトリに作ります。
次に、fastqcを使ってRibo-Seq_HEK293_Cell_Control_SRR1630831.fastqのデータのクオリティをチェックします。
-tオプションで使用するCPUのコア数を指定します。コア数を増やせば単純に計算速度が向上しますが、FastQCではメモリを1コアあたり250MB消費するので、要求できるメモリ量を超えるコア数を指定するとエラーが起こります。
-oオプションで出力先のディレクトリを指定し、-fオプションでInputファイルのファイル形式を指定しています。
FastQCから出力されるデータが、fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831ディレクトリの中のRibo-Seq_HEK293_Cell_Control_SRR1630831_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
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
4. Adapter trimming
PAR-CLIPのリードはアダプター配列の一部が含まれているので、Cutadaptを利用してそれを除きます。 アダプター配列の設定はサンプルに応じて変更が必要です。
どのようなアダプター配列がリードに混入しているかどうか調べるには、ダウンロード時に参照したNCBI GEOの各サンプルのページから得るか、もしくは論文のMaterial & Methodsを参照して得ます。
今回は、論文から情報を得ました。Ribo-seqで用いられるアダプター配列はほとんど同じで、下記の論文に記載されている、TCGTATGCCGTCTTCTGCTTGという配列です。
Hafner, et al. Methods, 44:3-12 (2008);
https://www.ncbi.nlm.nih.gov/pubmed/18158127
まず、TCGTATGCCGTCTTCTGCTTGというアダプター配列がそもそもFASTQファイル内に保存されている各リードに含まれているかどうか、目で確認したいと思います。
$ less Ribo-Seq_HEK293_Cell_Control_SRR1630831.fastq
ここでは、lessコマンドでFASTQファイルの中身を見ています。lessで表示したとき、十字キーでさらに下の情報も見ることができます。
また、/に続けて文字列を入力してEnterを押すと、その文字列を含む列を検索することができます。アダプター配列全長を検索にかけるとヒットしにくくなるので、5'end側から6-7bpぐらいの配列で検索にかけるとよいかもしれません（TCGTATGとか）。
検索した文字列と一致する箇所に、このようにハイライトがつきます。
上図では、今回のサンプルとは異なるサンプルを見ています。
実行例
$ cutadapt -m 10 -a TCGTATGCCGTCTTCTGCTTG \ Ribo-Seq_HEK293_Cell_Control_SRR1630831.fastq > Ribo-Seq_HEK293_Cell_Control_SRR1630831_1_trimmed_adapter.fastq 2>> ./log_Ribo-Seq_HEK293_Cell_Control_SRR1630831.txt
-f: フォーマットの指定。-m: トリミング後のリードの最低配列長。例えば、10と指定すると、トリミング後にリード長が10未満のリードは排除する。-a: 3'側にライゲーションされている、指定したアダプター配列を除去する。
5. Quality filtering
FASTQファイル中にクオリティスコアの低いリードが含まれている場合、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。
コード例
$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_1_trimmed_adapter.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o Ribo-Seq_HEK293_Cell_Control_SRR1630831_2_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を下回る場合、そのリードを排除する。
6. リボソームRNA由来のリードを除く
リボソームをIPしている関係で、Ribo-seqのFASTQファイルにはある程度リボソームRNA由来のコンタミしています。そのため、ゲノムにマッピングする前に、リボソームRNA由来のリードを除去します。
まず、リボソームRNAの配列を用意します。
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としてまとめます。
rRNAを除去するために、FASTQファイルに格納されている各リードをそれらの配列にアライメント（マッピング）して、それらの配列にマッピングされないリードを抽出します。こうすることで、rRNAに由来しないリードのみをもとのFASTQファイルから抽出することができます。
STARと呼ばれるソフトでマッピングさせるために、まず、rRNAのIndexファイルを用意する必要があります。以下で、先ほど作成したFASTAファイルをもとに、STAR用のIndexファイルを作成しています。
実行例
$ mkdir ./STAR_Index_rRNA_contam
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./STAR_Index_rRNA_contam \
--genomeFastaFiles ./contam_Ribosomal_RNA.fa
まず、Indexファイルを保存するディレクトリをmkdirコマンドで作成します。ここでは、STAR_Index_rRNA_contamという名前のディレクトリを用意しました。
STARのパラメータ。
--runThreadN: 使用するコア数を指定。--runMode:genomeGenerateを指定して、Indexファイルを作成するモードに設定。--genomeDir: 出力先のディレクトリを指定。--genomeFastaFiles: Indexを作成したいFASTAファイルを指定。ここでは、先ほど用意したリピート配列とrRNAの配列データ（FASTAファイル）を指定する。
rRNAを除く
Indexファイルを作成した後、rRNAへのマッピングを行います。マッピングされなかったリードについては、FASTQファイルとして出力されるようにオプションを指定します。
実行例
$ STAR --runMode alignReads --runThreadN 8 --genomeDir ./STAR_Index_rRNA_contam \ --readFilesIn ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_2_filtered.fastq \ --outSAMunmapped Within --outFilterMultimapNmax 30 --outFilterMultimapScoreRange 1 \ --outFileNamePrefix ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq \ --outSAMattributes All --outStd BAM_Unsorted --outSAMtype BAM Unsorted \ --outFilterType BySJout --outReadsUnmapped Fastx --outFilterScoreMin 10 \ --alignEndsType EndToEnd > ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_repbase_rrna_comtam.bam $ mv ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastqUnmapped.out.mate1 ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq
STARのパラメータ。
--runMode:alignReadsを指定して、アライメントを行うモードに設定。--runThreadN: 使用するコア数を指定。--genomeDir: Indexファイルの保存されているディレクトリを指定。--readFilesIn: マッピングさせたいFASTQファイルを指定。--outSAMunmapped Within: 出力されるBAMファイルにマッピングされなかったリード情報も含める。--outFilterMultimapNmax: 複数箇所にマッピング（Multi-mapping）されることをどこまで許容するか設定する。指定した数を超える箇所にマッピングされたリードは排除する。--outFilterMultimapScoreRange: Multi-mappingアライメントのスコアの範囲。複数箇所マッピングされてもアライメントスコアが高ければ残す。--outFileNamePrefix: 出力ファイルのPrefixを指定。-outSAMattributes: 基本情報以外に、BAMファイルに追加する情報を指定。--outStd BAM_Unsorted: 出力されるBAMファイルをソートしない。--outSAMtype BAM Unsorted: SAMファイルでなく、BAMファイルとしてマッピングされた情報を直接出力する。ただし、ソートはしない。--outFilterType BySJout:SJ.out.tabファイルとして出力される信頼性の高いExon junctionをまたぐリードのみを残す。Filteringをクリアできないリードは排除する。--outReadsUnmapped Fastx: マッピングされなかったリードをFASTQファイルとして出力する。--outFilterScoreMin 10: アライメントスコアが10を下回るリードは排除する。--alignEndsType EndToEnd: End-to-Endアライメントを行う（Localアライメントでない）。
7. Quality check
rRNAに由来するリードを除いた後、フィルタリング後のFASTQファイルをFastQCを用いて再確認します。
実行例
$ mkdir fastqc_${file}_filtered
$ fastqc -o ./fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831_filtered ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq -f fastq
8. Mapping
Indexファイルの作成
ここでも先ほどのマッピング時と同じように、GenomeとTranscriptomeのSTAR用のIndexファイルを作成する必要があります。このIndexファイルの作成には、50G程度のメモリが必要になるので注意してください。
まず、Human Genome（hg38）のFASTAファイルのダウンロードしてきます。
各染色体ごとにファイルが用意されているので、catコマンドでダウンロードしたすべてのファイルをマージして、hg38.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/hg38/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
Transcriptomeの情報として、Gencodeのアノテーション情報を利用します。
https://www.gencodegenes.org/releases/24.html
上記のサイトにアクセスすると、さまざまなデータが置いてあるが基本的にgencode.v19.annotation.gtfというメインのアノテーション情報を利用すればOKです。
ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_24/gencode.v24.annotation.gtf.gz
wgetコマンドで該当のファイルをダウンロードし解凍します。
$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_24/gencode.v24.annotation.gtf.gz $ tar zxvf ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_24/gencode.v24.annotation.gtf.gz
カレントディレクトリにgencode.v24.annotation.gtfが保存されます。
最後に、STARを使ってhg24+Gencode_v24のIndexファイルを作成します。
ここまでで、ゲノムとトランスクリプトームのIndexファイル用意できたので、いよいよSTARを用いて各シークエンスリードをマッピングします。
実行例
$ mkdir hg38_Genome_Gencode_v24
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir hg38_Genome_Gencode_v24 \
--genomeFastaFiles ./hg38.fa --sjdbGTFfile ./gencode.v24.annotation.gtf --sjdbOverhang 100
STARのパラメータ。
--runThreadN: 使用するコア数を指定。--runMode:genomeGenerateを指定して、Indexファイルを作成するモードに設定。--genomeDir: 出力先のディレクトリを指定。--genomeFastaFiles: Indexを作成したいFASTAファイルを指定。ここでは、先ほど用意したヒトゲノムのデータ（hg19、FASTAファイル）を指定する。--sjdbGTFfile: STARではアノテーション情報のGTFファイルからExon junctionの情報を抽出する。ここでは、先ほどダウンロードしたgencode.v19.annotation.gtfを使用する。--sjdbOverhang 100: Exon junctionの周囲、100塩基のゲノム配列をIndex化する。マッピングに使用するFASTQファイルのリード長と設定すると良い。公式のおすすめは100。
ゲノムとトランスクリプトームへのマッピング
Indexファイルが用意できたので、ゲノム（hg38）とトランスクリプトーム（Gencode v24）にリードをマッピングします。
実行例
$ mkdir STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831
$ STAR --runMode alignReads --runThreadN 8 --genomeDir ${indexFile} \
--readFilesIn ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq \
--outFilterMultimapNmax 8 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --sjdbScore 1 \
--outFilterMismatchNmax 4 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 \
--outFileNamePrefix ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_ \
--outSAMattributes All --outSAMtype BAM SortedByCoordinate --limitBAMsortRAM ${maxRAM} \
--outFilterType BySJout --outReadsUnmapped Fastx
STARのパラメータ。
--outSAMtype BAM SortedByCoordinate: SAMファイルでなく、BAMファイルとしてマッピングされた情報を直接出力する。出力ファイルは染色体の座標でソートする。--limitBAMsortRAM 32000000000: BAMファイルをソートするときに割り当てられるメモリの上限を設定。
9. RiboTaperによるORF予測
80Sリボソームと結合するRNA領域の長さは28-30nt程度であることが知られており、実際にRibo-seqのリード長も同じ長さを取っています。
また、コドンは3塩基ごとに読まれていくので、マッピングされたシークエンスリードの分布には、ORFの読み枠に合わせてバイアスがかかっていることが知られています。
例えば、RibosomeのP-siteにmRNAのAUG配列がロードされている場合、シークエンスリードの5'末から12nt下流の位置を1とすると、AUGのAがそれに該当します（13-15nt目がリボソームのP-siteの位置になります）。
------------AUG--------------- <-- 12nt -->1<----15/16nt---->
Ribo-seqのデータでは、先ほど定義したP-siteの相対的な位置（1: 5'末から13nt目）に、コドンの最初の塩基が来るケースが多いです（これがマッピングされたリードの偏り）。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4193932/
RiboTaperでは、上記のようなRibo-seqにおけるシークエンスリードの分布の偏りから、リードがマッピングされた領域がORFとしてリボソームに読まれて翻訳されている領域かどうか判断しています。
また、RiboTaperでは、このシークエンスリードのマッピングバイアスを事前にチェックするプログラムが用意されているので、それを使ってデータのクオリティに問題がないか最初に確認しましょう。
# RiboTaper用のアノテーション情報の準備
まず、ORFなどのアノテーション情報をRiboTaper用に用意する必要があります。RiboTaper側で用意されているcreate_annotations_files.bashシェルスクリプトを用いて、GencodeなどのGTFファイルからRiboTaper用のアノテーション情報を作成します。
$ create_annotations_files.bash \ ./gencode.v24.basic.annotation.gtf \ ./hg38.fa \ true \ true \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \ /home/akimitsu/software/RiboTaper_v1.3/scripts
- 1つ目の引数: アノテーション情報のGTFファイルを指定。
- 2つ目の引数: ヒトゲノム（hg38）のFASTAファイルを指定。
- 3つ目の引数: CCDS annotationを利用するかどうか（Gencodeのデータを推奨）。
- 4つ目の引数: Appris annotationを利用するかどうか（Gencodeのデータを推奨）。
- 5つ目の引数: 保存先のディレクトリを指定。
- 6つ目の引数: 使用するBedtoolsのパスを指定。
- 7つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。
CCDSとApprisについては下記のサイトを参照のこと。 https://www.gencodegenes.org/gencode_tags.html
# Ribo-seqデータのチェック
冒頭で説明したとおり、まず、マッピングしたRibo-seqのBAMファイルをチェックします。
実行例
$ create_metaplots.bash \ ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24/start_stops_FAR.bed \ metaplots_result \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir/ \ /home/akimitsu/software/RiboTaper_v1.3/scripts/
- 1つ目の引数: マッピングしたRibo-seqのBAMファイルを指定。
- 2つ目の引数: 先ほど作成したRiboTaper用のアノテーション情報の
start_stops_FAR.bedというファイルを指定。 - 3つ目の引数: 保存するファイルの名称を指定。
- 4つ目の引数: 使用するBedtoolsのパスを指定。
- 5つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。
実行すると、リード長ごとに、開始コドンと終止コドンの周辺のリードの分布を図にしてくれます。
緑色のバーがORFの読み枠（例えば、AUGのA）と一致しているリードの数になります。結果の図を見たときに、良い例のように緑色のバーにリードの分布が偏っているのが使えるシークエンスデータになります。
良い例
一方で、悪い例のようにリードの分布にバイアスがイマイチかかっていないサンプルは、ノイズの大きいデータになるので、RiboTaperによる解析には使えません。
悪い例
データのクオリティに問題がないと判断できたら、RiboTaperのシェルスクリプトを実行します。
RiboTaperの実行
RiboTaperには計算上のボトルネックが存在したり、プログラムが不安定で途中で処理が止まってしまったりします。
今回は、スクリプトを修正したり、処理を細分化して実行時にエラーが起こらないようにしたファイルを用意したので、そちらを使用した場合の実行例を示したいと思います。
とりあえず、下記のサイトにおいてあるscriptsディレクトリを、自分の環境のRiboTaper_v1.3の直下に上書きすればOKです。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Ribo-seq_RiboTaper
（注: RiboTaperのソースコードはGNU General Public Licenseのもとに配布されています。）
参考（スクリプトの修正箇所）
念のため、スクリプトの修正箇所について言及しておきます。 修正を加えている点はすべて、Rスクリプトを実行している箇所です。
もとのRスクリプトでは、１行目が、
#!/usr/bin/Rscript
となっており、スクリプトの実行はRscriptを利用したコマンドライン上での実行を行っていません。
例えばこんな感じ。
./hoge.R
そのため、Rscriptが/usr/bin/Rscriptに存在する環境下では機能しますが、Biocondaなどでローカル環境にRをインストールした場合には機能しません。
そこで、以下のように書き換えます。
Rscript ./hoge.R
こうすることで、/usr/bin/RscriptにRscriptがなくても、Rスクリプトが機能します。
修正箇所は、以下のとおりです。
CCDS_orf_finder.R (821行目)
# syst_scr<-paste(scr,"bed_tocheck_ccds",sep = " ")
scr2<-paste("Rscript",scr,sep=" ")
syst_scr<-paste(scr2,"bed_tocheck_ccds",sep = " ")
NONCCDS_orf_finder.R (738行目)
# syst_scr<-paste(scr,"bed_tocheck_nonccds",sep = " ")
scr2<-paste("Rscript",scr,sep=" ")
syst_scr<-paste(scr2,"bed_tocheck_nonccds",sep = " ")
Ribotaper_ORF_find.R (105-121行目)
# $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores Rscript $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores echo "NONCCDS ORF finding..." # $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores Rscript $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores # Groups ORFs and creates BED files + protein fasta database echo "Grouping ORFs and creating protein fasta database..." # $scripts_dir"/create_protein_db.R" Rscript $scripts_dir"/create_protein_db.R" # makes summary plot for the found ORFs echo "Summarizing ORF finding results" # $scripts_dir"/ORF_final_results.R" Rscript $scripts_dir"/ORF_final_results.R"
create_annotations_files.bash (220行目, 225行目, 235行目)
# $scripts_dir_full"/write_startstops.R" Rscript $scripts_dir_full"/write_startstops.R"
# $scripts_dir_full"/gtf_to_start_stop_tr.R" Rscript $scripts_dir_full"/gtf_to_start_stop_tr.R"
# $scripts_dir_full"/genes_coor.R" Rscript $scripts_dir_full"/genes_coor.R"
create_metaplots.bash (68行目)
# $scripts_dir"metag.R" $3 Rscript $scripts_dir"metag.R" $3
Ribotaper.sh (135-153行目, 167-177行目)
今回は、このスクリプトを分割して実行していますが、シェルスクリプトの中身は以下のように修正しています。
# $scripts_dir"/tracks_analysis.R" ccds $scripts_dir $n_of_cores Rscript $scripts_dir"/tracks_analysis.R" ccds $scripts_dir $n_of_cores echo "Running calculations exons_ccds..." # $scripts_dir"/tracks_analysis.R" exonsccds $scripts_dir $n_of_cores Rscript $scripts_dir"/tracks_analysis.R" exonsccds $scripts_dir $n_of_cores echo "Running calculations nonccds..." # $scripts_dir"/tracks_analysis.R" nonccds $scripts_dir $n_of_cores Rscript $scripts_dir"/tracks_analysis.R" nonccds $scripts_dir $n_of_cores # annotates the exons relative to ccds regions TO BE ADAPTED, CHECK WHICH FILES THEY NEED. echo "Annotate exons..." # $scripts_dir"/annotate_exons.R" $annot_dir $scripts_dir $n_of_cores Rscript $scripts_dir"/annotate_exons.R" $annot_dir $scripts_dir $n_of_cores echo "Making quality plots..." # $scripts_dir"/quality_check.R" $annot_dir Rscript $scripts_dir"/quality_check.R" $annot_dir
# $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores Rscript $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores echo "NONCCDS ORF finding..." # $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores Rscript $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores # Groups ORFs and creates BED files + protein fasta database echo "Grouping ORFs and creating protein fasta database..." # $scripts_dir"/create_protein_db.R" Rscript $scripts_dir"/create_protein_db.R"
今回実行に使うRiboTaperのシェルスクリプトは以下の通りに分割しておきます。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Ribo-seq_RiboTaper/scripts
- Ribo-seq-RiboTaper-3-define_ORF_STEP1.sh
- Ribo-seq-RiboTaper-3-define_ORF_STEP2-1.sh
- Ribo-seq-RiboTaper-3-define_ORF_STEP2-2.sh
- Ribo-seq-RiboTaper-3-define_ORF_STEP2-3.sh
- Ribo-seq-RiboTaper-3-define_ORF_STEP3.sh
- Ribo-seq-RiboTaper-3-define_ORF_STEP4.sh
上記が、実行するためのシェルスクリプトになります。これらを、ribotaperのscriptsディレクトリにあらかじめコピーしておきます。
実行時に注意すべき点は、RiboTaperはマルチコアで処理させる場合、使用するCPUのコア数によりますが、実行時に50-60GB程度のメモリを確保したほうがよいです（メモリが不足すると途中で処理が止まってしまいます）。
以下では、コマンドライン上で１つずつ実行したときの例を示します。
まず、最初に作成したribotaperの仮想環境に移行します。
source activate ribotaper
STEP1（所要時間: 2.5hr）
Ribotaper_step1.sh \ ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \ ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \ 27,28,29 \ 12,12,12 \ /home/akimitsu/software/RiboTaper_v1.3/scripts \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \ 8
- 1つ目の引数: Ribo-seqデータのBAMファイルを指定。
- 2つ目の引数: RNA-seqデータのBAMファイルを指定。
- 3つ目の引数: 先ほど作成したRiboTaper用のアノテーション情報のディレクトリを指定。
- 4つ目の引数: 使用するリード長（ノイズでなく、Ribosomeが結合していたリードのみを選抜）を指定（カンマ区切り）。先ほどのRibo-seqデータのチェック時に、緑色のバーに偏ってリードがマッピングされたリード長のみを選抜してくる。
- 5つ目の引数: P-sitesの計算に使用するOffset値を指定（カンマ区切り）。
- 6つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。
- 7つ目の引数: 使用するBedtoolsのパスを指定。
- 8つ目の引数: 使用するCPUのコア数。
実行後、以下のファイル群が得られるはずです。
P_sites_all Centered_RNA RIBO_unique_counts_ccds RIBO_best_counts_ccds RNA_unique_counts_ccds RNA_best_counts_ccds RIBO_unique_counts_exonsccds RIBO_best_counts_exonsccds RNA_unique_counts_exonsccds RNA_best_counts_exonsccds RIBO_unique_counts_nonccds RIBO_best_counts_nonccds RNA_unique_counts_nonccds RNA_best_counts_nonccds # data_tracks/ Centered_RNA_tracks_ccds Centered_RNA_tracks_exonsccds Centered_RNA_tracks_nonccds P_sites_all_tracks_ccds P_sites_all_tracks_exonsccds P_sites_all_tracks_nonccds Psit_Ribo_Rna_Cent_tracks_ccds Psit_Ribo_Rna_Cent_tracks_exonsccds Psit_Ribo_Rna_Cent_tracks_nonccds RIBO_tracks_ccds RIBO_tracks_exonsccds RIBO_tracks_nonccds RNA_tracks_ccds RNA_tracks_exonsccds RNA_tracks_nonccds index_tracks_ccds index_tracks_exonsccds index_tracks_nonccds
STEP2-1, 2-2, 2-3（所要時間: 2.5hr）
STEP2以降も同様に実行するだけです。STEP2-1, 2-2, 2-3は、並行して実行して構いません。Sun Grid Engineなどの並列計算システムを利用している場合、計算時間を短縮するために、並行してジョブを実行するようにしましょう。
以下には、STEP2-1の実行例を示しました。
Ribotaper_step2-1.sh \ ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \ ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \ 27,28,29 \ 12,12,12 \ /home/akimitsu/software/RiboTaper_v1.3/scripts \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \ 8
実行後、以下のファイル群が得られるはずです。
results_ccds results_exonsccds results_nonccds
STEP3（所要時間: 40min）
Ribotaper_step3.sh \ ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \ ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \ 27,28,29 \ 12,12,12 \ /home/akimitsu/software/RiboTaper_v1.3/scripts \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \ 8
実行後、以下のファイル群が得られるはずです。
results_nonccds_annot all_calculations_ccdsgenes_annot_new quality_check_plots.pdf
STEP4（所要時間: 8-9hr）
Ribotaper_step4.sh \ ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \ ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \ /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \ 27,28,29 \ 12,12,12 \ /home/akimitsu/software/RiboTaper_v1.3/scripts \ /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \ 8
実行後、以下のファイル群が得られるはずです。
# ORFs_CCDS/ # ORFs_NONCCDS/
ORF_maxORF_max_filt: filtered for multimapping and non-coding nonccds_coding_ORFs (see below) (recommended for further analysis)protein_db_max.fasta: (for proteomics search, not filtered for multimapping)translated_ORFs_sorted.bed: corresponding to ORF_maxtranslated_ORFs_filtered_sorted.bed: corresponding to ORF_max_filt filtered for multimapping and non-coding nonccds_coding_ORFs (see below) (recommended for further analysis)ORFs_genes_found: tab-separated values for number of ORFs found and their corresponding genesFinal_ORF_results.pdf: pdf file with length/coverage distributions for different ORF categories.
太文字で書かれたファイルが予測されたORFになります。ファイルの中身を確認すると以下の通り。
詳細については言及しませんが、Categoryを確認すると、主に、
- ORFs_ccds
- uORF
- ncORFS
があります。
ORFs_ccdsはmRNA由来のCanonical ORFを、uORFはmRNAの上流5'UTR中に存在するUpstream ORFを、ncORFsはNoncoding RNA中に予測されるORFを指しています。
