case3 - nibb-gitc/gitc2021mar-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
(ファイルは、~/gitc/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フォーマットのファイルがダウンロードされる。ダウンロードされたファイルは、カレントディレクトリにそれぞれのIDのディレクトリの下に保存される。
注）以前のバージョンのsratoolkitのデフォルト設定では、~/ncbi/public/sra に保存される。という挙動だった。sratoolkitのバージョンに注意。
9 file全てダウンロードしよう。
次に、sraフォーマットをfastqに変換しよう。
fastq-dump --split-files --origfmt DRR168786/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 ~/gitc/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のカウントデータを、ひとつのファイルにまとめる。salmon quantmerge サブコマンドを使う。salmon.quant.matrix.txt ファイルに書き出すことにする。
$salmon quantmerge --column numreads -o salmon.quant.matrix.txt --quants 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
output: salmon.quant.matrix.txt
$head salmon.quant.matrix.txt
Name 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
AT5G67620.1 0 0 0 0 0 0 0 0 0
AT5G67580.1 5 5 2.185 9 11.809 12.082 13 11.865 5
AT5G67550.1 0 0 0 0 0 1 0 0 0
AT5G67540.2 2 2 3 12 5 6 9 7 2
AT5G67520.1 2 2 0 2 3 0 4 0 3
AT5G67490.1 2 8 10 9 10 6 12 14 11
AT5G67460.1 0 0 1 0 0 1 4 0 0
AT5G67450.1 0 0 1 3 0 1 0 0 1
AT5G67440.1 3 2 9 8 9 10 16 15 9
Analysis with edgeR
準備
ここからRでの解析。RをローカルPCで解析する場合、salmon.quant.matrix.txt をサーバからコピーしておく。scpを使用。
Rを起動し、scpのあるディレクトリで作業をする。setwdコマンドなどを使ってワーキングディレクトリを設定する。
[R] （例）
setwd("~/Documents／case3")
(作業ディレクトリが上記と異なる場合は各自の環境に合わせて入力すること)
edgeRとgplotsのライブラリを読み込んでおく。
libary(edgeR)
library(gplots)
Normalization
[R]
library(edgeR)
tab <- read.delim("salmon.quant.matrix.txt", comment.char="#", row.names=1)
new.colnames <- sub("salmon_out_", "", colnames(tab))
colnames(tab) <- new.colnames
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 233143 1.0753172
DRR168787 2D 212285 1.0847719
DRR168788 2D 238443 1.0507711
DRR168789 4D 300049 1.0496378
DRR168790 4D 361092 1.0118684
DRR168791 4D 288957 1.0429869
DRR168792 2D2L 473029 0.9035653
DRR168793 2D2L 362166 0.9036133
DRR168794 2D2L 316781 0.9020538
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 1 2 1 1
> table(mycl)
mycl
1 2 3 4
78 43 55 24
このグルーピング結果も合わせてヒートマップ図に描画する。
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 )
Revision History
2021-3-11
- minor updates. 今回のコースの環境に合わせてアップデート。