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つ目の引数: 上記のスクリプトから得られた出力ファイルを指定します。
以上になります。