kallisto や Salmon を利用して定量したデータを使って、edgeR や DESeq2 などで発現量の群間比較を行うことができる。この際に、Bioconductor の tximport パッケージ（Soneson et al, 2015）を利用することで、簡単に kallisto/Salmon の定量結果を edgeR/DESEq2 に渡すことができる。
kallisto
kallisto の定量結果は abundance.tsv に保存されている。また、kallisto で定量を行うときにブートストラップによる不確実性推測を行ったとき、abundance.tsv に加えて abundance.h5 も出力される。これらの定量結果ファイルは、kallisto 実行時に -o で指定されたフォルダに保存される。例えば、この解析手順で解析した場合は、6 サンプルの結果はそれぞれ次の 6 つのフォルダに保存される。
ls -l
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:04 SRR7508939_exp_kallisto
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:14 SRR7508940_exp_kallisto
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:25 SRR7508941_exp_kallisto
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:35 SRR7508942_exp_kallisto
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:45 SRR7508943_exp_kallisto
## drwxr-xr-x 2 _____ ________ 4096 Aug 18 12:55 SRR7508944_exp_kallistoabundance.tsv
abundance.tsv ファイルを読み込んで発現量データを取得する場合は、まず abundance.tsv ファイルへのパスを変数に保存する。次にその変数を tximport 関数に与えて、発現量データの読み込みを行う。
library(tximport)
kallisto.files <- file.path(list.files('.', pattern = 'exp_kallisto'), 'abundance.tsv')
names(kallisto.files) <- c('WT1', 'WT2', 'WT3', 'NSR1', 'NSR2', 'NSR3')
tx.exp <- tximport(kallisto.files, type = "kallisto", txOut = TRUE)
head(tx.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G19090.1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## AT1G18320.1 12.2794 34.0661 27.4029 12.2095 15.9555 12.2916
## AT5G11100.1 154.0000 183.0000 170.0000 162.0000 180.0000 168.0000
## AT4G35335.1 455.0000 509.0000 514.0000 400.0000 493.0000 467.0000
## AT1G60930.1 158.0000 200.0000 176.0000 122.0000 174.0000 157.0000
## AT1G17000.1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000転写産物レベルの発現量の差を比較したい場合は、tx.exp$counts に保存されている発現量データを整数にまるめて edgeR や DESeq2 に与えればよい。
abundance.h5
abundance.h5 ファイルを読み込んで発現量データを取得する場合は、tximport パッケージに加えて rhdf5 パッケージも利用する。
library(tximport)
library(rhdf5)
kallisto.files <- file.path(list.files('.', pattern = 'exp_kallisto'), 'abundance.h5')
names(kallisto.files) <- c('WT1', 'WT2', 'WT3', 'NSR1', 'NSR2', 'NSR3')
tx.exp <- tximport(kallisto.files, type = "kallisto", txOut = TRUE)
head(tx.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G19090.1 0.00000 0.00000 0.00000 0.00000 0.0000 0.00000
## AT1G18320.1 12.27942 34.06605 27.40287 12.20948 15.9555 12.29165
## AT5G11100.1 154.00000 183.00000 170.00000 162.00000 180.0000 168.00000
## AT4G35335.1 455.00000 509.00000 514.00000 400.00000 493.0000 467.00000
## AT1G60930.1 158.00000 200.00000 176.00000 122.00000 174.0000 157.00000
## AT1G17000.1 0.00000 0.00000 0.00000 0.00000 0.0000 0.00000遺伝子発現量
kallisto が出力した発現量データは転写産物レベルでの発現量である。遺伝子レベルでの解析を行う場合は、これら転写産物レベルでの発現量を遺伝子レベルの発現量に換算する必要がある。この換算は tximport パッケージの summarizeToGene 関数で実行できる。その際に、転写産物の名前と遺伝子の名前の対応表をデータフレーム型で作成する必要がある。
今回のデータは TAIR のデータで、1 つの遺伝子に由来する転写産物が 1 つ以上存在する場合、転写産物の名前は AT1G19090.1、AT1G19090.2 などのように遺伝子 ID の後ろに .1 や .2 などつけて命名されている。そのため、転写物名からピリオド以降の文字を消せば遺伝子名となるので、転写産物の名前と遺伝子の名前の対応は次のようにして簡単に作れる。ただし、ヒトやマウスなどのデータについては、遺伝子名が ENSG00000254647、転写物名が ENST00000381330.4 のようになっているので、遺伝子 ID と転写物 ID が必ずしも同じ番号となっていない。このとき、Ensembl BioMart などを使って遺伝子 ID と転写物 ID の対応表を作るとよい。
tx2gene <- data.frame(
TXNAME = rownames(tx.exp$counts),
GENEID = sapply(strsplit(rownames(tx.exp$counts), '\\.'), '[', 1)
)
head(tx2gene)
## TXNAME GENEID
## 1 AT1G19090.1 AT1G19090
## 2 AT1G18320.1 AT1G18320
## 3 AT5G11100.1 AT5G11100
## 4 AT4G35335.1 AT4G35335
## 5 AT1G60930.1 AT1G60930
## 6 AT1G17000.1 AT1G17000summarizeToGene 関数を使って、転写産物レベルの発現量を遺伝子レベルの発現量に変更する。
gene.exp <- summarizeToGene(tx.exp, tx2gene)
head(gene.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G01010 222.000 253 223.000 223 294.000 280.000
## AT1G01020 423.001 447 439.001 392 406.000 404.999
## AT1G01030 121.000 123 104.000 89 89.000 85.000
## AT1G01040 1404.000 1694 1609.000 1389 1542.000 1519.000
## AT1G01050 1599.000 1860 1794.000 1571 1863.000 1867.000
## AT1G01060 150.000 209 179.999 136 162.999 155.001ただし、この方法で取得した発現量データをそのまま使って edgeR や DESeq2 に入力するのは推奨されていない。次のように summarizeToGene を実行するときに countsFromAbundance を指定してカウントデータを取得する必要がある。countsFromAbundance に指定できるスケーリング方法は scaledTPM と lengthScaledTPM の 2 種類があるが、どちらを指定してもいい。
gene.exp <- summarizeToGene(tx.exp, tx2gene, countsFromAbundance = "scaledTPM")
head(gene.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G01010 188.86270 214.61478 189.75976 190.13068 249.63802 238.97528
## AT1G01020 497.72941 524.26361 518.76680 454.08255 481.00256 477.12631
## AT1G01030 90.70910 91.24827 77.58019 66.38409 66.10422 63.84052
## AT1G01040 303.46163 362.61675 342.59580 300.09212 330.97301 328.64488
## AT1G01050 2506.22423 2915.33620 2822.49119 2474.33500 2919.97415 2936.40219
## AT1G01060 79.22783 117.31857 96.77064 77.70996 91.45248 88.53382
library(edgeR)
d <- DGEList(counts = gene.exp$counts, group = c("WT", "WT", "WT", "NSR", "NSR", "NSR"))
d <- calcNormFactors(d)
# ...Salmon
Salmon の定量結果は quant.sf に保存されている。これらの定量結果ファイルは、quant 実行時に -o で指定されたフォルダに保存される。例えば、この解析手順で解析した場合は、6 サンプルの結果はそれぞれ次の 6 つのフォルダに保存される。
ls -l
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:04 SRR7508939_exp_salmon
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:05 SRR7508940_exp_salmon
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:06 SRR7508941_exp_salmon
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:07 SRR7508942_exp_salmon
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:09 SRR7508943_exp_salmon
## drwxr-xr-x 5 _____ ________ 4096 Aug 18 12:10 SRR7508944_exp_salmontximport パッケージで Salmon の定量結果を読み込んで発現量行列を取得する場合は、次のようにする。このとき tximport パッケージのほかに、jsonlite パッケージ、readr パッケージも必要である。
library(tximport)
library(jsonlite)
library(readr)
salmon.files <- file.path(list.files('.', pattern = 'exp_salmon'), 'quant.sf')
names(salmon.files) <- c('WT1', 'WT2', 'WT3', 'NSR1', 'NSR2', 'NSR3')
tx.exp <- tximport(salmon.files, type = "salmon", txOut = TRUE)
head(tx.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G19090.1 0.000 0.000 0.000 0.000 0.000 0.000
## AT1G18320.1 12.232 34.056 27.443 10.766 15.972 12.213
## AT5G11100.1 154.000 183.000 170.000 162.000 180.000 167.000
## AT4G35335.1 455.000 509.000 514.000 399.000 492.000 467.000
## AT1G60930.1 158.000 200.000 176.000 120.000 171.000 157.000
## AT1G17000.1 0.000 0.000 0.000 0.000 0.000 0.000kallisto が出力した発現量データは転写産物レベルでの発現量である。遺伝子レベルでの解析を行う場合は、これら転写産物レベルでの発現量を遺伝子レベルの発現量に換算する必要がある。この換算は tximport パッケージの summarizeToGene 関数で実行できる。その際に、転写産物の名前と遺伝子の名前の対応表をデータフレーム型で作成する必要がある。
tx2gene <- data.frame(
TXNAME = rownames(tx.exp$counts),
GENEID = sapply(strsplit(rownames(tx.exp$counts), '\\.'), '[', 1)
)
head(tx2gene)
## TXNAME GENEID
## 1 AT1G19090.1 AT1G19090
## 2 AT1G18320.1 AT1G18320
## 3 AT5G11100.1 AT5G11100
## 4 AT4G35335.1 AT4G35335
## 5 AT1G60930.1 AT1G60930
## 6 AT1G17000.1 AT1G17000
gene.exp <- summarizeToGene(tx.exp, tx2gene)
head(gene.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G01010 222.000 253 223.000 223 294.000 280.000
## AT1G01020 423.001 447 439.001 392 406.000 404.999
## AT1G01030 121.000 123 104.000 89 89.000 85.000
## AT1G01040 1404.000 1694 1609.000 1389 1542.000 1519.000
## AT1G01050 1599.000 1860 1794.000 1571 1863.000 1867.000
## AT1G01060 150.000 209 179.999 136 162.999 155.001
gene.exp <- summarizeToGene(tx.exp, tx2gene, countsFromAbundance = "scaledTPM")
head(gene.exp$counts)
## WT1 WT2 WT3 NSR1 NSR2 NSR3
## AT1G01010 189.14087 214.94639 189.94012 190.28102 249.90102 239.17130
## AT1G01020 497.27802 521.11795 518.73172 454.95068 480.78651 476.41325
## AT1G01030 90.73407 91.37981 77.56397 66.42176 66.15988 63.82483
## AT1G01040 303.53048 362.81946 342.18813 299.85386 330.76783 328.38150
## AT1G01050 2514.29196 2921.89373 2830.07796 2480.22428 2927.59919 2942.01994
## AT1G01060 77.95323 117.72696 96.64158 78.16146 91.81200 88.78277遺伝子レベルの発現量比較を行う場合は countsFromAbundance を使用した summarizeToGene 関数で取得したカウントデータを edgeR や DESeq2 などに与えればよい。countsFromAbundance に指定できるスケーリング方法は scaledTPM と lengthScaledTPM の 2 種類があるが、どちらを指定してもいい。
References
- Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2015, 7(3):562-78. DOI: 10.12688/f1000research.7563.1