読者です 読者をやめる 読者になる 読者になる

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

使用するソフトウェア

解析ワークフロー

RNA-seqデータを2群間比較する方法として、Cuffdiffを用いた方法とfeatureCount-edgeRを用いた方法を紹介したいと思います。

TopHat-Cuffdiff workflow

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2016-09-02 11.33.15.png (62.3 kB)

TopHat-featureCounts-edgeR workflow

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2016-09-01 18.22.58.png (156.1 kB)

サンプルコード

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_listgencode.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行目: annoListgencode.v19.annotation_filtered_symbol_type_list.txt(後述)を指定。
11行目: saveDirに保存先のディレクトリ名を指定。

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

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

権限の問題

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

$ chmod 755 hoge.sh

改行コードの問題

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

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

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

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

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

各ステップの説明

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

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

imamachi-n.hatenablog.com

$ conda install fastqc
$ conda install fastx_toolkit
$ conda install bowtie=1.12
$ conda install tophat
$ conda install samtools
$ conda install rseqc
$ conda install bedtools
$ conda install 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という項目にリストアップしてあります。

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

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

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

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

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

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

SRR4081222.sraのURLをコピー。

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

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

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

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

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

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

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

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

2. Quality check

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

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

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

実行例

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

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

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

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

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

Basic Statistics

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

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

などなど

Per base sequence quality

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

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

Per tile sequence quality

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

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

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

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

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

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

Per sequence quality scores

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

Per base sequence content

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

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

Per sequence GC content

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

Per base N content

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

Sequence Length Distribution

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

Sequence Duplication Levels

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

Overrepresented sequences

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

Adapter Content

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

Kmer Content

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

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

image.png (73.1 kB)

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

20170205212758.png (1.7 MB)

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列目のfeatureexonである行を使用するように指定している。
  • -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つ目の引数で指定したサンプルに対して、ラベル付け(例えば、ControlKnockdown)を行う。データはカンマ区切りで指定する。
  • 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比較

2017-02-08_cuffdiff_vs_edgeR.xlsx - Excel 2017-02-08 20.03.59.png (23.9 kB)

CuffdiffとedgeRのFold-changeの比較

image.png (85.9 kB)

発現量が低い遺伝子については、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を使うようにしましょう。