BRIC-seqデータ解析（２群間比較, Unstranded, Single-end）

BRIC-seq BrideR2 NGS RNA decay

概要

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

使用するソフトウェア

解析ワークフロー

BRIC-seqデータからRNA半減期を算出し、さらに2群間の半減期を比較する方法を紹介したいと思います。 PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2016-09-06 16.00.52.png (152.5 kB)

サンプルコード

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半減期の算出、および２群間比較を行うスクリプト。

$ ./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変数に２群のラベル名を代入。 9行目: hour変数にタイムコースを代入。 35行目: comparisonFile変数に２群のラベル名を代入。 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に変更。

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

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

各ステップの説明

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

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

imamachi-n.hatenablog.com

$ conda install fastqc
$ conda install fastx_toolkit
$ conda install bowtie=1.12
$ conda install tophat
$ conda install samtools
$ conda install 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以上）に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。

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

Per tile sequence quality

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

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

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

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

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

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

Per sequence quality scores

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

Per base sequence content

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

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

Per sequence GC content

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

Per base N content

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

Sequence Length Distribution

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

Sequence Duplication Levels

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

Overrepresented sequences

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

Adapter Content

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

Kmer Content

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

3. 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というファイルが出力されます。ファイルの中身を見るとこんな感じ。

image.png (73.1 kB)

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ボタンをクリックすると、データをみることができます。

20170205212758.png (1.7 MB)

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の半減期の計算で出力されるファイル）を指定します。

表示させると、以下のようになります。

C__Users_Naoto_Downloads - Shiny 2017-02-06 19.41.22.png (115.4 kB)

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を開き、新しいプロジェクトを作成する。

無題 - ペイント 2017-02-06 19.46.39.png (175.0 kB)

保存したフォルダを開き、server.Rとui.Rを消します。代わりにapp.Rというファイルを作成してください（ここは必ずapp.Rという名前にしてください。別名にするとshinyapp.ioにデプロイできません）。

また、10. RNA半減期の計算の最後に生成された./BridgR_result_PUM1_study_mRNA/BridgeR_6_PUM1KD_halflife_pvalue_evaluation.txtのファイルをコピーして同じフォルダに入れます。

こんな感じになります。

BRIC-seq_shinyapp 2017-02-07 10.56.39.png (14.2 kB)

上図のように、保存先のフォルダに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になっているので、この関数１つで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をクリックして内容をコピーします。

shinyapps.io - Google Chrome 2017-02-07 11.20.19.png (248.3 kB)

[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と出たら、デプロイ成功です。

C__Users_Naoto_OneDrive_Shiny_app_BRIC-seq_shinyapp - RStudio 2017-02-06 20.05.00.png (31.4 kB)

表示されているURLにアクセスすると、今までRStudio上で見ていたRNAの分解曲線のデータをWeb上でも閲覧することが可能になります。

また、指定されたURLを他人に公開することで、RStudioをインストールしていないユーザーに対しても、Webブラウザを介してデータを共有することができます。

shinyapps.ioの自身のDashboardからもShinyアプリのURLを確認することができます。

Recent ApplicationsのところにアップロードしたShinyアプリが表示されており、それをクリックすると、URLなどの情報が記載されているページに飛ぶことができます。

shinyapps.io - Google Chrome 2017-02-07 14.59.04.png (249.6 kB)

内容は以上になります。