case3 - nibb-gitc/gitc2020jun-rnaseq GitHub Wiki
Case study 3: Transcript-based RNA-Seq pipeline (alignment-free)
シロイヌナズナ(Arabidopsis thaliana)のRNA-seqを行った。ライブラリは2D sample (2days dark conditionで生育させた黄色芽生え)と2D2L sample (その後さらに2days light conditionで生育させた緑化芽生え）,4D sample (2days dark condition後引き続き 2days dark conditionで生育させたもの) の３条件、それぞれbiological replicates を３つ用意した。シーケンスはpaired-end（インサートの両端を読む）で76bpシークエンスした(MiSeq)。
これらのシーケンスファイルは、DDBJ/DRAに登録されている。公的データベースに登録されているNGSデータを取得して解析するパイプラインを学ぶために、sratools を使ってデータを取得するステップから実習する。また、この演習では、alignment-free のTranscript-base のパイプラインで解析を進める。３条件以上のRNAseq解析の例として、MDS plot やクラスタリング、ヒートマップ描画などの多変量解析も学ぶ。
基本パイプライン：Salmon -> edgeR
多変量解析は R
Data
(ファイルは、~/data/EX/case3 以下 にある)
Reference
- transcript sequence: Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa （シロイヌナズナのトランスクリプト。１遺伝子１トランスクリプトの代表配列。詳細は後述。）
Reads
- DRP004823_SRR_Acc_List.txt に書かれてあるaccession番号に基づいて、自分でダウンロードする。
Software
salmon (installed)
edgeR (installed)
Setup
case3 ディレクトリをつくり、以下の解析はその下で作業しよう。
$ mkdir case3
$ cd case3
Download Illumina sequences
DRP004823_SRR_Acc_List.txt に書かれてあるaccession番号のNGSデータをNCBI SRA から取得したい。
[DRP004823_SRR_Acc_List.txt]
DRR168786
DRR168787
DRR168788
DRR168789
DRR168790
DRR168791
DRR168792
DRR168793
DRR168794
まず、NCBI SRAのウェブサイトでどのような情報が登録されているか、ウェブ上で確認してみよう。
https://www.ncbi.nlm.nih.gov/sra/
にアクセスし、DRR168786 などのアクセション番号を検索してみよう。 情報をまとめると下記の通り。
Run Library_Name treatment
----------------------------------------
DRR168786 2D_rep1 2D
DRR168787 2D_rep2 2D
DRR168788 2D_rep3 2D
DRR168789 4D_rep1 4D
DRR168790 4D_rep2 4D
DRR168791 4D_rep3 4D
DRR168792 2D2L_rep1 2D2L
DRR168793 2D2L_rep2 2D2L
DRR168794 2D2L_rep3 2D2L
それでは、sratoolkit を使ってこれらのNGSデータをダウンロードしよう。
prefetch DRR168786
一括ダウンロードも可能
prefetch DRR168787 DRR168788 DRR168789 DRR168790
prefetch に -p オプションをつけると進行の度合いが表示されるのでオススメです。
sraフォーマットのファイルがダウンロードされる。ダウンロードされたファイルは、sratoolkitのデフォルト設定では、~/ncbi/public/sra に保存される。
9 file全てダウンロードしよう。ダウンロードが完了したらカレントディレクトリにコピーしよう。
cp ~/ncbi/public/sra/DRR1687*.sra ./
次に、sraフォーマットをfastqに変換しよう。
fastq-dump --split-files --origfmt DRR168786.sra
今回はpaird-endライブラリなので、それぞれRead1, Read2のファイルが生成される。
- DRR168786_1.fastq
- DRR168786_2.fastq
9 file全てについて、fastq変換をしよう。上記コマンドをDRR168786.sraの部分だけ変えて８回繰り返しても良いが、シェルスクリプトを使うと便利。
for f in DRR1687*sra
do
fastq-dump --split-files --origfmt $f
done
seqkit の stats サブコマンドを使って生成されたfastqのシーケンス数などの基礎情報を取得する。
seqkit stats DRR1687*fastq
file format type num_seqs sum_len min_len avg_len max_len
DRR168786_1.fastq FASTQ DNA 477,502 36,290,152 76 76 76
DRR168786_2.fastq FASTQ DNA 477,502 36,290,152 76 76 76
DRR168787_1.fastq FASTQ DNA 414,943 31,535,668 76 76 76
DRR168787_2.fastq FASTQ DNA 414,943 31,535,668 76 76 76
DRR168788_1.fastq FASTQ DNA 604,501 45,942,076 76 76 76
DRR168788_2.fastq FASTQ DNA 604,501 45,942,076 76 76 76
DRR168789_1.fastq FASTQ DNA 463,620 35,235,120 76 76 76
DRR168789_2.fastq FASTQ DNA 463,620 35,235,120 76 76 76
DRR168790_1.fastq FASTQ DNA 616,252 46,835,152 76 76 76
DRR168790_2.fastq FASTQ DNA 616,252 46,835,152 76 76 76
DRR168791_1.fastq FASTQ DNA 501,444 38,109,744 76 76 76
DRR168791_2.fastq FASTQ DNA 501,444 38,109,744 76 76 76
DRR168792_1.fastq FASTQ DNA 856,831 65,119,156 76 76 76
DRR168792_2.fastq FASTQ DNA 856,831 65,119,156 76 76 76
DRR168793_1.fastq FASTQ DNA 564,622 42,911,272 76 76 76
DRR168793_2.fastq FASTQ DNA 564,622 42,911,272 76 76 76
DRR168794_1.fastq FASTQ DNA 460,199 34,975,124 76 76 76
DRR168794_2.fastq FASTQ DNA 460,199 34,975,124 76 76 76
Prepare reference transcripts
tarsncrtip-base RNA-seq パイプラインで、referenceに何を使うかはとても重要です。ここでは、JGIのPhytozome (https://phytozome.jgi.doe.gov/pz/portal.html) が整理した、Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz をリファレンスとして使うことにします。１遺伝子につき代表的なスプライシングバリアントが１つ選ばれているので、得られた結果の解釈やその後の解析がgene単位でやりやすいためです。
cp ~/data/EX/case3/JGI/Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz ./
gunzip Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz
これでreferenceとreadsの両方が準備完了です。
Salmon: alignment-free quantification
build reference index
salmon index -t Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa -i Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.salmon.idx -k 31
この作業は一度やればよい。case1でも同じリファレンスを使ったので、コピーしても良い。
Mapping and quantification
DRR168786の場合、
salmon quant \
-i Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.salmon.idx \
-l A \
-1 DRR168786_1.fastq -2 DRR168786_2.fastq \
-o salmon_out_DRR168786
case1 と異なり、今回はpaired-end readsなので、オプションがいくつか異なる点に注意。
9 reads 全て同様に解析する。
salmon quant -i Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.salmon.idx -l A -1 DRR168786_1.fastq -2 DRR168786_2.fastq -o salmon_out_DRR168786
salmon quant -i Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.salmon.idx -l A -1 DRR168787_1.fastq -2 DRR168787_2.fastq -o salmon_out_DRR168787
...
(option, advanced) 上記ように、手動で一つずつコマンドを実行しても良いが、以下のようなシェルスクリプトでまとめて実行すると便利だし、ファイル名のミスなどのヒューマンエラーも防げる
for n in DRR168786 DRR168787 DRR168788 DRR168789 DRR168790 DRR168791 DRR168792 DRR168793 DRR168794
do
salmon quant -i Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.salmon.idx -l A -1 ${n}_1.fastq -2 ${n}_2.fastq -o salmon_out_${n} > run_salmon_${n}.log 2>&1
done
Results =>
salmon_out_DRR168786
salmon_out_DRR168787
salmon_out_DRR168788
salmon_out_DRR168789
salmon_out_DRR168790
salmon_out_DRR168791
salmon_out_DRR168792
salmon_out_DRR168793
salmon_out_DRR168794
の９つのディレクトリが生成されそれらの下にsalmonの結果が保存されたと仮定して、以降進める。
Expression analysis using edgeR
Prepare count matrix
このあとの解析がしやすいように、サンプルごとに別々のファイルに記録されているsalmonのカウントデータを、ひとつのファイルにまとめる。edgeRは、TPMでなくカウントデータを入力としなければいけない。それぞれの、quant.sf ファイルの、NumReads カラムを抜き出す。この作業にはやや煩雑なテキストデータ処理を要するので、筆者が用意したRubyスクリプト(case1でも利用)、merge_salmon_quant.rbを使って欲しい。スクリプトは、~/data/EX/case1/ にある。
(使い方)
$ruby merge_salmon_quant.rb dir1 dir2 dir3 ...
(dir1, dir2 はsalmonの出力ディレクトリ。それぞれ、quant.sf ファイルが含まれていることを前提としている)
(今回のケースでの実行例)
ruby merge_salmon_quant.rb salmon_out_DRR168786 salmon_out_DRR168787 salmon_out_DRR168788 salmon_out_DRR168789 salmon_out_DRR168790 salmon_out_DRR168791 salmon_out_DRR168792 salmon_out_DRR168793 salmon_out_DRR168794 > salmon.quant.matrix.txt
output: salmon.quant.matrix.txt
(ヘッダ省略）
id salmon_out_DRR168786 salmon_out_DRR168787 salmon_out_DRR168788 salmon_out_DRR168789 salmon_out_DRR168790 salmon_out_DRR168791 salmon_out_DRR168792 salmon_out_DRR168793 salmon_out_DRR168794
AT1G01010.1 5.000 5.000 5.000 6.000 7.000 11.000 0.000 0.000 4.000
AT1G01020.1 2.000 3.000 9.000 6.000 10.000 5.000 6.000 4.000 2.000
AT1G01030.1 2.000 2.000 0.000 2.000 2.000 0.000 0.000 0.000 0.000
...
Analysis with edgeR
準備
ここからRでの解析。カレントディレクトリに移動する必要がある。GUIの場合、メニューの、その他＞作業ディレクトリの変更... setwdコマンドも使える。 [R]
setwd("~/case3")
(作業ディレクトリが上記と異なる場合は各自の環境に合わせて入力すること)
edgeRとgplotsのライブラリを読み込んでおく。
libary(edgeR)
library(gplots)
Normalization
[R]
library(edgeR)
tab <- read.delim("salmon.quant.matrix.txt", comment.char="#", row.names=1)
colnames(tab) <- c("DRR168786", "DRR168787", "DRR168788", "DRR168789", "DRR168790", "DRR168791", "DRR168792", "DRR168793", "DRR168794")
treat <-factor(c(rep("2D", 3), rep("4D", 3), rep("2D2L", 3)))
treat <- relevel(treat, ref="2D")
dat <- DGEList(tab, group=treat)
dat <- calcNormFactors(dat, method="TMM")
dat$samples
group lib.size norm.factors
DRR168786 2D 238144 1.0602971
DRR168787 2D 217460 1.0614850
DRR168788 2D 244953 1.0300861
DRR168789 4D 306532 1.0297365
DRR168790 4D 369546 1.0407085
DRR168791 4D 295835 1.0318780
DRR168792 2D2L 478921 0.9354238
DRR168793 2D2L 367523 0.9051940
DRR168794 2D2L 321147 0.9211946
MDS plot
MDS plot はサンプル間の関係性を可視化する。これは unsupervised clustering（教師なしクラスタリング）の一種である。replicatesどうしが近くにクラスタリングされることが期待される。MDS plot は実験の妥当性の評価にも有用である。edgeRでは plotMDSコマンドでMDS plotを簡単に描くことができる。
plotMDS(dat)
DEG analysis
GLMを用いて、ANOVA-likeな解析を行う。３つの実験条件でいずれかの条件で発現変動のある遺伝子を検出する。
design <- model.matrix(~treat)
dat <- estimateDisp(dat, design)
fit <- glmFit(dat, design)
lrt <- glmLRT(fit, coef=2:3)
計算結果を確認しよう。
topTags(lrt)
DEG遺伝子はいくつあるだろうか？
lrt.fdr.sorted <- topTags(lrt, n=nrow(dat))$table
sum(lrt.fdr.sorted$FDR < 0.05)
sum(lrt.fdr.sorted$FDR < 0.01)
多変量解析
sampleの階層的クラスタリング
MDS plot でサンプル間の関係性を可視化することができたが、ここでは階層的クラスタリングでサンプル間の関係性を捉えてみよう。
ここではDEG遺伝子のみを使って、サンプルの階層的クラスタリングを行う。
deg <- subset(lrt.fdr.sorted, lrt.fdr.sorted$FDR<0.01)
logcpm <- cpm(dat, prior.count=1, log=TRUE)
logcpm.deg <- logcpm[rownames(deg), ]
d2 <- as.data.frame(t(logcpm.deg))
sampleTree <- hclust(dist(d2), method="ward.D2")
plot(sampleTree, hang=0.05)
MDS plot と同様、replicate間は非常に似ていることが確認できた。さらに、2Dと4Dがクラスターを形成しているが、全く光を当てない２条件が、光をあてる条件よりお互い似ているのはリーズナブルな結果といえよう。
heatmap
ここではDEG上位200位までの遺伝子を使ってヒートマップを描いてみよう。
top200 <- topTags(lrt, n=200)
logcpm <- cpm(dat, prior.count=1, log=TRUE)
logcpm.top200 <- logcpm[rownames(top200$table), ]
heatmapColors <- colorpanel(40, low="blue", high="yellow")
heatmap.2(as.matrix(logcpm.top200), col=heatmapColors, scale="row", trace="none")
heatmap.2のdefaultではクラスタリングにユークリッド距離が使われる。遺伝子のクラスタリングに相関係数 (Pearson correlation coefficient) を使って、ヒートマップ・クラスタリングをやり直してみよう。
hcl.genes <- hclust(as.dist(1-cor(t(logcpm.top200))))
heatmap.2(as.matrix(logcpm.top200), col=heatmapColors, scale="row", trace="none", Rowv=as.dendrogram(hcl.genes) )
最初に書いたヒートマップよりクラスタリングがよりリーズナブルになっているように見える。
さて、遺伝子の発現パターンは大まかにいくつに分類できるだろうか。ヒートマップとクラスタリングの樹形図を見ると３〜５がリーズナブルなグループ数に見える。ここでは４グループとして、遺伝子を発現パターンで分類してみよう。cutree関数を使うと階層的クラスタリングの結果に基づいてグルーピングしてくれる。
> mycl <- cutree(hcl.genes, k=4)
> head(mycl)
AT2G34430.1 AT1G29920.1 AT5G54770.1 AT3G54890.1 AT1G29910.1 AT2G34420.1
1 1 2 2 2 2
> table(mycl)
mycl
1 2 3 4
36 82 57 25
このグルーピング結果も合わせてヒートマップ図に描画する。
mycolhc <- rainbow(length(unique(mycl)), start=0.1, end=0.9)
mycolhc <- mycolhc[as.vector(mycl)]
heatmap.2(as.matrix(logcpm.top200), col=heatmapColors, scale="row", trace="none", Rowv=as.dendrogram(hcl.genes), RowSideColors=mycolhc )