case1 - nibb-gitc/gitc2021mar-rnaseq GitHub Wiki
Case study 1: Transcript-based RNA-seq pipeline
updated: Mar 10, 2021
Copyright: Shuji Shigenobu [email protected] (NIBB)
Pipeline overview
transcript base のRNA-seq解析の基本的なパイプラインを学ぶ。発現量推定と発現比較解析までのプロトコルを具体的なスクリプトやコマンドを追って説明する。
実験デザイン：シロイヌナズナ Arabidopsis thaliana を明条件(L)と暗条件(D)で育て、それぞれ3個体からRNAを抽出して (three biological replicates)、illumina TruSeq kitでRNA-seqライブラリを作製した（ライブラリの平均長は200bp、SD=±50bpであった）。それらを、イルミナシーケンサーMiSeqでシーケンスした(片側76base シーケンス)。明暗 (L vs D) の条件の差で発現の異なる遺伝子を同定したい。
[Strategy]
1. mapping reads to transcriptome reference and abundance estimation
tool: salmon
2. differential expression analysis
tool: edgeR
Setup
bias5のLinux環境で作業する。bias5にsshでログインし、case1ディレクトリを作りその下で解析しよう。
$ mkdir case1
$ cd case1
Software
本解析に必要なソフトウェアは以下の通り。
- salmon, R, edgeR,
Data set
[~/gitc/data/EX/case1/]
data/
└── EX
└── case1
├── Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz
├── IlluminaReads
│ ├── D1_R1.fastq
│ ├── D2_R1.fastq
│ ├── D3_R1.fastq
│ ├── L1_R1.fastq
│ ├── L2_R1.fastq
│ └── L3_R1.fastq
└── merge_salmon_quant.rb
Prepare reference transcripts
transcript-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/case1/Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz ./
$ gunzip Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa.gz
解凍したreference transcript のファイルには幾つtranscriptが含まれているか確認しておきましょう。
# 方法１
$ grep "^>" Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa |wc
# 方法２
$ seqkit stats Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa
# => 結果
file format type num_seqs sum_len min_len avg_len max_len
Athaliana_....fa FASTA DNA 27,416 40,622,613 21 1,481.7 16,347
Prepare Illumina read files
Six Illumina read files are saved in ~/gitc/data/EX/case1/IlluminaReads
まずはlsコマンドでファイルサイズなどを確認
$ ls -lh ~/gitc/data/EX/case1/IlluminaReads/
total 486M
-rw-r--r-- 1 shige analyins 77M Mar 1 2018 D1_R1.fastq
-rw-r--r-- 1 shige analyins 68M Mar 1 2018 D2_R1.fastq
-rw-r--r-- 1 shige analyins 88M Mar 1 2018 D3_R1.fastq
-rw-r--r-- 1 shige analyins 97M Mar 1 2018 L1_R1.fastq
-rw-r--r-- 1 shige analyins 82M Mar 1 2018 L2_R1.fastq
-rw-r--r-- 1 shige analyins 77M Mar 1 2018 L3_R1.fastq
それぞれのファイルが何本のリードを含んでいるかを確認しておく。
$ seqkit stats ~/gitc/data/EX/case1/IlluminaReads/*fastq
# => output
file format type num_seqs sum_len min_len avg_len max_len
data/EX/case1/IlluminaReads/D1_R1.fastq FASTQ DNA 382,799 29,047,372 51 75.9 76
data/EX/case1/IlluminaReads/D2_R1.fastq FASTQ DNA 336,384 25,474,502 51 75.7 76
data/EX/case1/IlluminaReads/D3_R1.fastq FASTQ DNA 436,695 33,107,457 51 75.8 76
data/EX/case1/IlluminaReads/L1_R1.fastq FASTQ DNA 485,913 36,864,248 51 75.9 76
data/EX/case1/IlluminaReads/L2_R1.fastq FASTQ DNA 407,812 30,943,819 51 75.9 76
data/EX/case1/IlluminaReads/L3_R1.fastq FASTQ DNA 382,669 28,988,314 51 75.8 76
Mapping Illumina Reads to Reference and quantification
Salmon を使って、リードのマッピングとカウントを同時に行う。
Build salmon index
まず、reference をindexing。この作業は一度やればよい。
$ salmon index -t Athaliana_167_TAIR10.transcript_primaryTranscriptOnly.fa -i At.reference_transcripts.salmon.idx -k 31
Mapping and quantification
例として、D1_R1.fastq の処理の手順を示す。
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D1 -p 2 -r data/EX/case1/IlluminaReads/D1_R1.fastq --fldMean 200 --fldSD 50
outputs => salmon_out_D1/
quant.sf が結果ファイル。中身を確認しておきましょう。
D1_R1.fastq D2_R1.fastq D3_R1.fastq L1_R1.fastq L2_R1.fastq L3_R1.fastq ６つのシーケンスファイルすべてについて、同様にマッピングを行なう。
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D1 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/D1_R1.fastq --fldMean 200 --fldSD 50
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D2 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/D2_R1.fastq --fldMean 200 --fldSD 50
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D3 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/D3_R1.fastq --fldMean 200 --fldSD 50
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L1 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/L1_R1.fastq --fldMean 200 --fldSD 50
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L2 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/L2_R1.fastq --fldMean 200 --fldSD 50
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L3 -p 2 -r ~/gitc/data/EX/case1/IlluminaReads/L3_R1.fastq --fldMean 200 --fldSD 50
Differential expression analysis using edgeR
Prepare count matrix
サンプルごとに別々のファイル(quant.sf)に記録されているsalmonのカウントデータを、ひとつのファイルにまとめる。edgeRは、TPMでなくカウントデータを入力としなければいけない。それぞれの、quant.sf ファイルの、NumReads カラムを抜き出す。salmon のquantmergeサブコマンドを使う。salmon_quant.txt という名前のファイルに出力することにする。
salmon quantmerge の使い方
salmon quant merge --column numreads -o OUTPUT --quants DIR1 DIR2 ...
- --column : どのカラムのデータを抽出するか。今回 NumReads を抽出するので明示的に指定する必要がある。
- -o : 出力結果ファイルの名前
- --quants : 抽出したいquant.sf を含むディレクトリの名前 複数スペースで区切って列挙する。
$salmon quantmerge --column numreads -o salmon_quant.txt --quants salmon_out_D1 salmon_out_D2 salmon_out_D3 salmon_out_L1 salmon_out_L2 salmon_out_L3
出力結果のsalmon_quant.txt を確認しよう
$ head salmon_quant.txt
Name salmon_out_D1 salmon_out_D2 salmon_out_D3 salmon_out_L1 salmon_out_L2 salmon_out_L3
AT5G67620.1 0 0 0 0 0 1
AT5G67580.1 5.02 5 2.202 11 5.289 8
AT5G67550.1 0 0 0 0 0 0
AT5G67540.2 2 3 3 7 2 3
AT5G67520.1 2 2 0 0 3 0
AT5G67490.1 2 8 10 14 11 11
AT5G67460.1 0 0 1 0 0 0
AT5G67450.1 0 0 1 0 1 3
AT5G67440.1 3 2 9 15 10 5
...
〜〜〜参考スタート〜〜〜
注）以下は前回のコースまでに使ったrubyスクリプトを使ったカウントマトリックスの作成法です。今回は上記の通りsalmon quantmerge で同じことが実現できますので、rubyスクリプトは不要です。ただ、参考までに記録を残しておきます。（スクリプトは改造して応用が効くので、参考したい方のために）。
サンプルごとに別々のファイルに記録されているsalmonのカウントデータを、ひとつのファイルにまとめる。edgeRは、TPMでなくカウントデータを入力としなければいけない。それぞれの、quant.sf ファイルの、NumReads カラムを抜き出す。この作業にはやや煩雑なテキストデータ処理を要するので、筆者が用意したRubyスクリプトmerge_salmon_quant.rb`を使って欲しい。スクリプトは、~/gitc/data/EX/case1/ にある。
(使い方)
$ruby merge_salmon_quant.rb dir1 dir2 dir3 ...
(dir1, dir2 はsalmonの出力ディレクトリ。それぞれ、quant.sf ファイルが含まれていることを前提としている)
(今回のケースでの実行例)
ruby merge_salmon_quant.rb salmon_out_D1 salmon_out_D2 salmon_out_D3 salmon_out_L1 salmon_out_L2 salmon_out_L3 > salmon_quant.txt
=> salmon_quant.txt
〜〜〜参考終わり〜〜〜
Data check and scatter plot
複雑な統計計算で発現変動解析をあれこれ行なう前に、簡単なデータチェックをしておく。
以下、R環境で。Rはローカル環境で実行するので、bias5上で得られた、salmon_quant.txt ファイルをscpでローカルにコピーする必要がある。
[from your computer]
$ scp [email protected]:case1/salmon_quant.txt ./
> library(edgeR)
> dat <- read.delim("salmon_quant.txt", comment.char="#", row.name=1)
> head(dat)
salmon_out_D1 salmon_out_D2 salmon_out_D3 salmon_out_L1
AT1G01010.1 6 5 5 0
AT1G01020.1 2 3 10 4
AT1G01030.1 2 2 0 0
AT1G01040.2 11 22 9 8
AT1G01050.1 22 17 12 35
AT1G01060.1 0 1 2 1
salmon_out_L2 salmon_out_L3
AT1G01010.1 4 1
AT1G01020.1 2 8
AT1G01030.1 0 5
AT1G01040.2 10 10
AT1G01050.1 28 15
AT1G01060.1 0 1
カラム名の"salmon_out_"が目障りなので削除して、D1, D2などの名前に変更する。
> new.colnames <- sub("salmon_out_", "", colnames(dat))
> colnames(dat) <- new.colnames
> head(dat)
D1 D2 D3 L1 L2 L3
AT5G67620.1 0.00 0 0.000 0 0.000 1
AT5G67580.1 5.02 5 2.202 11 5.289 8
AT5G67550.1 0.00 0 0.000 0 0.000 0
AT5G67540.2 2.00 3 3.000 7 2.000 3
AT5G67520.1 2.00 2 0.000 0 3.000 0
AT5G67490.1 2.00 8 10.000 14 11.000 11
OK.
scatter plot を描いてみよう
(example of scatter plot)
> plot(dat$D1 + 1, dat$D2+1, log="xy")
(example of all-vs-all scatter plot)
> pairs(dat+1, log="xy")
(example of comparison between D1vsD2 and D1vsL1)
> par(mfrow=c(1,2))
> plot(dat$D1 + 1, dat$D2+1, log="xy")
> plot(dat$D1 + 1, dat$L1+1, log="xy")
edgeRData import
> library(edgeR)
> category <- c("D", "D", "D", "L", "L", "L")
> D <- DGEList(dat, group=category) # import table into edgeR
TMM normalization
> D <- calcNormFactors(D, method="TMM") # TMM normalization
> D$samples
group lib.size norm.factors
D1 D 232200 1.0902235
D2 D 211877 1.0865568
D3 D 237463 1.0633594
L1 L 362196 0.9137676
L2 L 316769 0.9182738
L3 L 328381 0.9461152
# dump normalized count data
> D.cpm.tmm <- cpm(D, normalized.lib.size=T)
> write.table(D.cpm.tmm, file="cpm.tmm.txt", sep="\t", quote=F)
Differential expression analysis
> D <- estimateCommonDisp(D) # estimate common dispersion
> D$common.dispersion
[1] 0.02656185
> D <- estimateTagwiseDisp(D) # estimate tagwise dispersion
> summary(D$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000415 0.011661 0.032213 0.095506 0.143547 1.299024
> de <- exactTest(D, pair=c("D", "L")) # significance test to find differentially expressed genes
> topTags(de) # view the most significant genes
# dump DE analysis result
> de.sorted <- topTags(de, n=nrow(de$table))
> write.table(de.sorted$table, "de.txt", sep="\t", quote=F)
Result evaluation
有意に発現変動している遺伝子はいくつあるのか。例えば、LでDより２倍以上高発現し (logFC > log2(2))、FDR < 0.01 の遺伝子の数は、以下で求められる。
> sum(de.sorted$table$FDR<0.01 & de.sorted$table$logFC > log2(2))
edgeRの command decideTestsDGEも便利。
使い方例
> summary(decideTestsDGE(de, p.value=0.01, lfc=1))
L-D
Down 513
NotSig 25867
Up 922
(edgeRのバージョンによって多少結果は異なる)。
plotSmear (edgeRに含まれるMA描画ツール) でMA plotを描いてみよう。
de.names <- row.names(D[decideTestsDGE(de, p.value=0.01) !=0, ])
plotSmear(D, de.tags=de.names)
MDS plotを行なうと、ライブラリ間の発現パターンの類似性をおおざっぱにとらえることができる。
plotMDS(D)
Revision History
2021-3-10
- count table作成に、salmon quantmerge を使う方法へ変更。以前のrubyスクリプトの方法も記載は残してある。
2021-3-11
- minor changes. Salmon 1.4, R 4.0.4 今回のトレーニングコースの環境での実行結果を反映。