case1 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki
Case study 1: Transcript-based RNA-seq pipeline
last updated: August 1, 2024
Copyright: Shuji Shigenobu [email protected] (NIBB)
Pipeline overview
- 受講生: 本練習問題は、RCCSのLinux環境で解析する。
- 聴講生: 解析用ソフトウェアを自身のPCにインストールする方法 を参考に、自身の端末で作業を行う。
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
RCCSのLinux環境で作業する。RCCSに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
(発展) 通常は上記の通り、transcriptのインデクシングのみで十分だが、より正確性を高めるためには、decoy-aware なindexingが推奨されている。具体的な方法は、こちら => case1-supp
Mapping and quantification
例として、D1_R1.fastq の処理の手順を示す。
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D1 -p 2 -r D1_R1.fastq
注）引数の「D1_R1.fastq」はファイルの存在するパスを指定するか、当該ファイルをカレントディレクトリにコピー（もしくはシンボリックリンク）する必要がある。
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 D1_R1.fastq
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D2 -p 2 -r D2_R1.fastq
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_D3 -p 2 -r D3_R1.fastq
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L1 -p 2 -r L1_R1.fastq
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L2 -p 2 -r L2_R1.fastq
$ salmon quant -i At.reference_transcripts.salmon.idx -l A -o salmon_out_L3 -p 2 -r L3_R1.fastq
Differential expression analysis using edgeR
Prepare count matrix
サンプルごとに別々のファイル(quant.sf)に記録されているsalmonのカウントデータを、ひとつのファイルにまとめる。edgeRは、TPMでなくカウントデータを入力としなければいけない。それぞれの、quant.sf ファイルの、NumReads カラムを抜き出す。salmon のquantmergeサブコマンドを使う。salmon_quant.txt という名前のファイルに出力することにする。
salmon quantmerge の使い方
salmon quantmerge --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はローカル環境で実行するので、RCCS上で得られた、salmon_quant.txt ファイルをscpでローカルにコピーする必要がある。
[from your computer]
UNIX環境
$ scp [email protected]:{PATH_TO_CASE1}/salmon_quant.txt ./
R環境
> library(edgeR)
> dat <- read.delim("salmon_quant.txt", 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")
edgeR Data 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")
> D$samples
group lib.size norm.factors
D1 D 232200 1.0893989
D2 D 211877 1.0862693
D3 D 237463 1.0644724
L1 L 362196 0.9139904
L2 L 316769 0.9179276
L3 L 328381 0.9462179
# 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)
注）edgeR ver.4以上では、calcNormFactorsの代わりにnormLibSizesを使う。
Differential expression analysis
> D <- estimateDisp(D) # estimate common dispersion and tagwise dispersion
> D$common.dispersion
[1] 0.02764591
> summary(D$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.02070 0.03587 0.11010 0.10972 0.18205 0.29842
> de <- exactTest(D, pair=c("D", "L")) # significance test to find differentially expressed genes
> topTags(de) # view the most significant genes
Comparison of groups: L-D
logFC logCPM PValue FDR
AT1G29920.1 8.145632 11.100253 1.362386e-151 3.719586e-147
AT3G21720.1 -6.078798 11.918235 9.810673e-147 1.339255e-142
AT2G34430.1 10.908939 10.819305 1.224629e-139 1.114494e-135
AT1G29910.1 6.309143 11.393867 2.486223e-134 1.696972e-130
AT2G38530.1 -5.603057 11.855005 4.809343e-131 2.626094e-127
AT5G54770.1 6.694113 10.885557 2.535386e-123 1.153685e-119
AT5G54190.1 -10.105614 9.702486 2.362122e-121 9.212949e-118
AT2G34420.1 5.315564 12.663228 4.755532e-120 1.622944e-116
AT3G54890.1 5.145765 11.691668 4.676092e-116 1.418518e-112
AT4G14130.1 -6.105307 9.702112 4.758250e-103 1.299098e-99
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も便利。(新しいedgeR ver.4以降では、decideTests)
使い方例
> summary(decideTestsDGE(de, p.value=0.01, lfc=1))
L-D
Down 504
NotSig 25891
Up 907
(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
2024-7-31
- 動作確認 on bias5 (salmon: v1.10.0; edgeR: v4.2.1)
- text minor updates
2024-2-26
- decoy-aware indexについての記述と、その解析方法についてcase1-suppを追加。
2023-2-25
- 動作確認 on bias5 (salmon: v1.8.0; edgeR: 3.32.1)
2022-8-28
- 動作確認＋minor update. salmon: v1.8.0; edgeR: 3.32.1
2022-2-26
- minor update. salmon v1.7.0
2021-9-11
- minor updates with the latest salmon (v1.4.0) / edgeR ()
2021-3-11
- minor changes. Salmon 1.4, R 4.0.4 今回のトレーニングコースの環境での実行結果を反映。
2021-3-10
- count table作成に、salmon quantmerge を使う方法へ変更。以前のrubyスクリプトの方法も記載は残してある。