2012年にnatureでtophatとcufflinksを使ったRNA-seq解析の手法が発表されている。
http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html
tophat、bowtie、cufflinksなどを使ったいわゆるタキシードプロトコルはRNA seqの1つのワークフローとして知られており、検索すると多くのプロトコルが見つかる(検索例)。
その1つが以下になる。
上記の手法にのっとり、RNA seq解析を行ってみよう。使えそうなRのプログラムも試してみることにする。
生き物はなんでもよいが、今回は高精度なgtfファイルが手に入るシロイヌナズナのRNA seqデータを使うことにする。SRAからアラビドプシスのRNA-seqデータを探すと、以下の2群比較のデータが見つかった。これを使ってみる。
Transcriptome analysis of arabidopsis kamchatica under normal and stress conditions
アラビドプシス属ではあるが、シロイヌナズナではない。Arabidopsis kamchatica(画像)を実験材料に使っている。このシーケンスデータをナズナゲノムに当てて定量していく。
シロイヌナズナのfastaファイルとgtfファイルはEnsembl genomeの
http://plants.ensembl.org/info/website/ftp/index.html
からダウンロード。
nature protocolに載ったワークフローは以下になる。
Prestep、Quality check
まずデータのクオリティを確認しておく。末端でクオリティが極端に悪くなる部位があればこの段階でクオリティトリムする。
ペアリードのR1を調べる。
fastqc -t 12 -o fastx_test_R1/ -f fastq DRR013381_1.fastq
- -t Specifies the number of files which can be processed simultaneously.
- -o Create all output files in the specified output directory.
- -f Bypasses the normal sequence file format detection and forces the program to use the specified format. Valid formats are bam,sam,bam_mapped,sam_mapped and fastq
- -k Specifies the length of Kmer to look for in the Kmer content module. Specified Kmer length must be between 2 and 10. Default length is 7 if not specified.
出力されたhtmlを開く。
トリミングは必要なさそうである。マッピングに移ろう。
step1、Mapping
cufflinksまでは、biological replicatesがあっても独立したサンプルとして扱っていく。今回はテストなので、chr1だけに限定してマッピングを行う(本当は良くない。アライメントの仕方次第でアーティファクトが出る可能性もゼロではない)。
bowtie2のidnex作成
bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa chr1
最後のchr1はユーザーが決めた名前である。この名前を使いtophatを走らせる。
tophat2でpaired-endリードをmapping
tophat2 -p 4 -I 50000 -i 20 --library-type fr-firststrand \
-o tophat_control1 -G Arabidopsis_thaliana.TAIR10.36.gtf chr1 \
DRR013381_1.fastq DRR013381_2.fastq
- -p num-threads
- -I --max-intron-length (default: 500000)
- -i --min-intron-length (default: 50)
- -G --GTF
- -o --output-dir (default: ./tophat_out)
(gtfをtophatの時点で指定してmappingする方法もある。その場合-Gをつける)。
アライメントファイルが出力される。今回はリファレンスとデータで種が違うので、アライメントは非常に困難になる。エラーは覚悟の上、--read-mismatches 4 --read-gap-length 4 --read-edit-dist 4も上のtohat2につけてランした(How to control the alignment of reads in terms of number of mismatches, gap length etc. ? を参考。edit distはwikiに説明がある)。ナズナはエコタイプの違いでもたくさん塩基置換が起きている。このパラメータでマッピング感度は上がるが、ランダムマッチも増える。どれくらい有効 or 無謀な方法なのかは分からない。または、bowtie2のvery-sensitiveオプションの方がいいかもしれない。tophatはbowtie2のオプションを --b2-のprefixをつけると認識できるようになっている。bowtie2には以下のような一括設定のフラグがある。
--very-fast Same as: -D 5 -R 1 -N 0 -L 22 -i S,0,2.50
--fast Same as: -D 10 -R 2 -N 0 -L 22 -i S,0,2.50
--sensitive Same as: -D 15 -R 2 -L 22 -i S,1,1.15 (default in --end-to-end mode)
--very-sensitive Same as: -D 20 -R 3 -N 0 -L 20 -i S,1,0.50
なのでtophatで--very-sensitive をつけるなら--b2-very-sensitive というフラグをつければいい(リンク)。
追記;マッピング率は40%ほどだった(default parameterで20%ほど)。
tophatの出力はaccepted_hits.bamという定型名のため、サンプルの区別がつかなくなる恐れがある。出力されるbamをサンプル名が分かる名前にリネームする。
mv tophat_output/accepted_hits.bam tophat_output/control1.bam
のちにIGVで可視化するために、bamに indexをつけておく。
samtools index tophat_output/control1.bam
gtfファイルを指定して解析した場合、次の2-4は飛ばしてstep5に進む。gtfを指定していない場合、cufflinksにアライメント状況からfeatureファイル(gtf)を作ってもらう必要がある。step2-4を実行する。(詳細については、topにリンクを貼ったnature protocolの論文box1. Cを参照)。
今回はcuffmergeしてアライメント領域を1から再定義した方が賢明と思われるが、時間があまりない。でナズナのgtfをそのまま使う。
2、Cufflinks
アライメント部位の情報に基づき、orfを予測してgtfファイルを出力する。
A、 mRNAを再構成して、gftファイルを出力する。
cufflinks --no-update-check --overlap-radius 1 --library-type fr-firststrand -o cufflinks control1.bam
- –overlap-radius <int> Transfrags that are separated by less than this distance get merged together, and the gap is filled. Default: 50 (in bp).
- –no-update-check Turns off the automatic routine that contacts the Cufflinks server to check for a more recent version.
B、出力されるgtfファイルをサンプル名がわかる名前にリネームする。
mv cufflinks/transcripts.gtf cufflinks/control1.gtf
これをサンプル全てに対して行う。サンプルが多いならシェルスクリプトを書いて自動処理を考える。
3-4、Cuffmerge
このstepもgtfファイルをアライメント時に指定している場合は必要ない。cufflinksに読み込んで全サンプルのgtfを統合し、1つのgtfファイルを出力する。最初に読み込むファイルのリストを作る(サンプルが少ないならコマンド実行時に明示してもよい)。
echo cufflinks/control1.gtf > assemblies.txt
echo cufflinks/control2.gtf >> assemblies.txt
echo cufflinks/control3.gtf >> assemblies.txt
echo cufflinks/stress1.gtf >> assemblies.txt
echo cufflinks/stress2.gtf >> assemblies.txt
echo cufflinks/stress3.gtf >> assemblies.txt
catで中身を確認しておく。
gtfのマージを実行する。
cuffmerge -s Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa assemblies.txt
merged_asm/merged.gtfができる。このmeggeしたgtfファイルを使って以下の計算を進める。
5、Cuffdiff
cuffdiffコマンドで サンプル間の発現量差を計算する。
使い方は
cuffdiff [options] <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> [... sampleN_hits.sam]
実際のランは以下のように行った。
cuffdiff -o cuffdiff_output -p 12 --library-type fr-firststrand \
-u Arabidopsis_thaliana.TAIR10.36.gtf -b \
Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa \
-N --compatible-hits-norm -L control,stress \
tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam \
tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam
- -o --output-dir
- -p --num-threads
- -b use bias correction - reference fasta required
- -u use 'rescue method' for multi-reads
- --library-type Library prep used for input reads
- -N --compatible-hits-norm count hits compatible with reference RNAs only
- -L comma-separated list of condition labels
-Lでつけた名前を今後使うことになる。区別がつくシンプルな名前が良い。 n > 2の実験でreplicatesがある場合、replicatesは","でつなぐ。異なるサンプルはスペースでつなぐ(manual)。上のデータはcontrolとstressのtriplicateのデータなので、上のboxの様に記述している。そのほか、末端にある""はコマンドを複数行にわたって記述する時のマークである。
出力のcuffdiff_outputにはたくさんのファイルができる。
発現量の差を調べたファイルの中を見てみる。
gene_exp.diff
最初の10行だけ表示。 value(FPKMではない)とp_valueが表示されている。これだけでも十分解析できそうである。
FPKMはgenes.fplm.trackingの中にある。
おまけ
igvでbamを開いて、データを確認してみる。igvがない人は
brew install igv
で入る。
fastaとgtf、bamファイルをコマンドから指定して開く。
igv -g Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa
Arabidopsis_thaliana.TAIR10.36.gtf,tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam,tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam
gtfファイル以降は ","でつないで記載する(gtf,bam1,bam2,bam3...)。
gtfがソートされてないと怒られるが、そのまま開く。(gtfのソートについてはバイオインフォ道場さんの記事が参考になります http://bioinfo-dojo.net/2017/04/01/igv_gtf_bed/)。
bedも選択した場合、gtfの下に追加表示される。
pop upが邪魔なら上のボタンからnever showをチェック。
6-18、CummeRbund
Rにデータを読み込む。ターミナルでRとタイプ
uesaka$ R
R version 3.3.0 (2016-05-03) -- "Supposedly Educational"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
R は、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。
R は多くの貢献者による共同プロジェクトです。
詳しくは 'contributors()' と入力してください。
また、R や R のパッケージを出版物で引用する際の形式については
'citation()' と入力してください。
'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。
>
ターミナル内でRを起動する。3.3が起動している。
library("cummeRbund")を打ってcummeRbundを読み込む。cummeRbundを入れてない人はBiocunducterのパッケージをインストールする。
source("http://bioconductor.org/biocLite.R") #時間がかかります
biocLite("cummeRbund")
biocLite("org.At.tair.db") #GO解析用
library("cummeRbund")
library("org.At.tair.db")
cuffdiff出力ディレクトリを読み込む。
cuff = readCufflinks('diff_output')
cummeRbundのグラフコマンドを一通り実行してみる。
csDensity(genes(cuff))
csScatter(genes(cuff), 'stress1', 'control1')
csScatterMatrix(genes(cuff))
csVolcanoMatrix(genes(cuff))
cufflinks2.2.1にバグがあるのか、cummeRbundで正しく読み込めない。というかsqliteのデータのファイルサイズが妙に小さいので、正しくデータベースが構築できていない気がする。以前cuffldiffで解析したデータは正常にcummeRbundに読み込めたため、cummeRbundは問題ない可能性が高い。正しくcummeRbundに読み込め次第追記します。
cufflinksを使わず、featurecountsでリードをカウントし、edgeRで正規化、検定してDEGを出す流れは以下の2記事を参考にしてください。
リードのカウント
edgeRで検定
全体の流れを再現するならこちらをどうぞ。2017年夏に少人数で勉強会を行った時の資料です。コピペで再現できるはずです。
引用
RNA seq workflow
https://en.m.wikibooks.org/wiki/Next_Generation_Sequencing_(NGS)/RNA