シミュレーションデータ

RNA-seq のリードカウントデータから発現変動遺伝子の検出法を開発するとき、一般的に、先にシミュレーションデータを利用して手法開発し、性能テストを行う。開発した手法がシミュレーションデータにおいて、十分に精度が確認された後、実際の RNA-seq のリードカウントデータに適応し、動作確認を行う。

シミュレーションデータは一般に負の二項分布に従うようにランダム生成され、さらに外れ値を加えたものが多い。以下に、シミュレーションデータを生成する方法を紹介する。二群間比較のシミュレーションデータを生成する関数は TCC、edgeR や compcodeR などのパッケージに実装されている。また、多群間比較あるいは多因子、時系列データのようなシミュレーションデータを生成する関数は TCC パッケージに実装されている。

二群間比較のシミュレーションデータ

TCC パッケージによるシミュレーションデータの生成

library(TCC)
packageVersion("TCC")
## [1] ‘1.8.0’

Ngene          <- 10000          # 遺伝子数
replicates     <- c(3, 3)        # 各実験群の biological replicate 数
PDEG           <- 0.1            # データ全体のおける発現変動遺伝子の割合
DEG.assign     <- c(0.8, 0.2)    # 発現変動遺伝子の分布
DEG.foldchange <- c(4, 4)        # 発現変動遺伝子の発現量の fold change

set.seed(2015)
tcc <- simulateReadCounts(
    Ngene = Ngene, replicates = replicates, PDEG = PDEG,
    DEG.assign = DEG.assign, DEG.foldchange = DEG.foldchange
)

カウントデータは tcc$count に保存される。

head(tcc$count)
##        G1_rep1 G1_rep2 G1_rep3 G2_rep1 G2_rep2 G2_rep3
## gene_1       5       0      11       5       2       0
## gene_2     356     478     783     249     103     120
## gene_3      31      32      32       4      10      11
## gene_4      73     160     111      45      61      37
## gene_5       3       2       2       0       0       0
## gene_6     313     321     552     128      76     107

どの遺伝子が発現変動遺伝子なのかを示す情報は tcc$simulation に保存されている。TCC パッケージが生成するカウントデータは遺伝子のデータセットの最初の方にまとめられている。

head(tcc$simulation$trueDEG)
## gene_1 gene_2 gene_3 gene_4 gene_5 gene_6
##      1      1      1      1      1      1
tail(tcc$simulation$trueDEG)
## gene_9995  gene_9996  gene_9997  gene_9998  gene_9999 gene_10000
##         0          0          0          0          0          0

発現変動遺伝子の割合を 10% としているので、実際の発現変動遺伝子の個数は 1000 個である。

sum(tcc$simulation$trueDEG == 1)
## [1] 1000
sum(tcc$simulation$trueDEG == 0)
## [1] 9000

生成されたシミュレーションデータの MA プロットは以下のようになる。赤でプロットされた遺伝子は発現変動遺伝子である。発現変動遺伝子の fold change を 4 に設定してあるので、赤のプロットは ± log24 に集中しているのが確認できる。

edgeR パッケージによるシミュレーションデータの生成

edgeR-robust の手法開発にシミュレーションデータが使われた。そのようなシミュレーションデータを生成する関数は Robinson らが公開している R コードに書かれている。

source("http://130.60.190.4/robinson_lab/edgeR_robust/robust_simulation.R")
library(tweeDEseqCountData)
library(edgeR)
packageVersion("edgeR")
## [1] ‘3.10.0’
packageVersion("tweeDEseqCountData")
## [1] ‘1.6.0’

dataset     <- pickrell       # サンプリングする母集団
nTags       <- 10000          # 遺伝子数
group       <- factor(c(1, 1, 1, 2, 2, 2)) # 実験群の識別ラベル
pDiff       <- 0.1            # データ全体のおける発現変動遺伝子の割合
pUp         <- 0.8            # 発現変動遺伝子の分布
foldDiff    <- 4              # 発現変動遺伝子の発現量の fold change
pOutlier    <- 0.1            # 外れ値の割合
outlierMech <- "S"            # 外れ値の生成アルゴリズム("S", "R", "M")

data(pickrell)
pickrell <- as.matrix(exprs(pickrell.eset))

set.seed(2015)
simdata <- NBsim(
   foldDiff = foldDiff, dataset = dataset, nTags = nTags,
   group = group, pDiff = pDiff, pUp = pUp,
   pOutlier = pOutlier, outlierMech = outlierMech
)

シミュレーションデータは simdata$counts に保存されている。

head(simdata$counts)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## ids1    3    1    0    1    1    1
## ids2    3    1    0    0    1    0
## ids3    8    5    1    9   11    6
## ids4    2    0    0    0    2    1
## ids5  950 1017  925  886 1457 1171
## ids6   65   26   50   35   50   58

simdata 変数の中にはカウントデータのみならず、その他の情報も保存されている。

names(simdata)
##  [1] "dataset"     "nTags"       "lib.size"    "nlibs"       "group"      
##  [6] "design"      "pDiff"       "pUp"         "foldDiff"    "outlierMech"
## [11] "min.factor"  "max.factor"  "name"        "Lambda"      "Dispersion" 
## [16] "indDE"       "indNonDE"    "mask_nonDE"  "mask_DEdown" "mask_DEup"  
## [21] "mask_DE"     "counts"    

どの遺伝子が発現変動遺伝子かどうかの情報は mask_DEmask_DEdownmask_DEup に保存されている。mask_DE には発現変動遺伝子かどうかの情報が直接に保存されている。シミュレーションデータの生成条件では nTags = 10000pDiff = 0.1 としているため、発現変動遺伝子の数は 1000 と期待される。実際にシミュレーション中の発現変動遺伝子の数を計算すると、1000 個であることがわかる。

is.deg <- apply(simdata$mask_DE, 1, all)
head(is.deg)
## [1] FALSE  TRUE FALSE FALSE FALSE FALSE
sum(is.deg)
## [1] 1000

どの実験群で高発現であるかの情報は mask_DEup に保存されている。発現変動遺伝子の高発現と低発現のパターンはを表示すると以下のようになっている。

head(simdata$mask_DEup[is.deg, ])
##      [,1] [,2] [,3]  [,4]  [,5]  [,6]
## [1,] TRUE TRUE TRUE FALSE FALSE FALSE
## [2,] TRUE TRUE TRUE FALSE FALSE FALSE
## [3,] TRUE TRUE TRUE FALSE FALSE FALSE
## [4,] TRUE TRUE TRUE FALSE FALSE FALSE
## [5,] TRUE TRUE TRUE FALSE FALSE FALSE
## [6,] TRUE TRUE TRUE FALSE FALSE FALSE

head(simdata$counts[is.deg, ])
##       [,1] [,2] [,3] [,4] [,5] [,6]
## ids2     3    1    0    0    1    0
## ids16  573  382  592  129  126  230
## ids39   21    8   18    0   22    3
## ids49    0    0    0    0    0    0
## ids50   29   17   43   13    4   11
## ids63   44   38   40    8    8   14

発現変動遺伝子がどのように分布されているかを調べる。発現変動遺伝子のうち、最初の実験群(左の 3 biological replicate)で高発現する遺伝子数は 812 個である。発現変動遺伝子の約 80% にあたる。これは pUp = 0.8 と設定しているからである。

colSums(simdata$mask_DEup[is.deg, ])
## [1] 812 812 812 188 188 188

生成されたシミュレーションデータの MA プロットは以下のようになる。赤でプロットされた遺伝子は発現変動遺伝子である。発現変動遺伝子の fold change を 4 に設定してあるので、赤のプロットは ± log24 に集中しているのが確認できる。

compcodeR パッケージによるシミュレーションデータの生成

compcodeR は発現変動遺伝子の検出法を比較するために作成されるパッケージである。パッケージの中ではシミュレーションデータを生成する関数が用意されている。

library(compcodeR)
packageVersion("compcodeR")
## [1] ‘1.4.0’

dataset          <- "simdata_seed_2015"    # データセットの名前
samples.per.cond <- 3             # 1 つの実験群のサンプル数
n.vars           <- 10000         # 遺伝子数
n.diffexp        <- 1000          # 発現変動遺伝子の数
effect.size      <- 3             # 発現変動遺伝子の発現量の fold change
random.outlier.high.prob <- 0.1   # 外れ値を入れる確率の最大値
random.outlier.low.prob  <- 0.0   # 外れ値を入れる確率の最小値
fraction.upregulated     <- 0.2   # 第 2 群で高発現する発現変動遺伝子の数

set.seed(2015)
simdata <- generateSyntheticData(
    dataset = dataset, n.vars = n.vars, n.diffexp = n.diffexp,
    effect.size = effect.size, samples.per.cond = samples.per.cond, 
    random.outlier.high.prob = random.outlier.high.prob,
    random.outlier.low.prob = random.outlier.low.prob,
    fraction.upregulated = fraction.upregulated
)

simdata 変数にはシミュレーションデータのほかに、シミュレーション条件などが保存されている。

head(simdata@count.matrix)
##    sample1 sample2 sample3 sample4 sample5 sample6
## g1     382    2388    3745    5922    3999    4731
## g2     986     885     184    6899   14920    1192
## g3     288      70     154     295     421     610
## g4      14       7       9      43      39      62
## g5       4      19       7     378      35      19
## g6    1520    1193    1389    3748    5425    4913

その遺伝子が発現変動遺伝子であるかどうかの情報は simdata@variable.annotations の中に保存されている。

head(simdata@variable.annotations$differential.expression)
## [1] 1 1 1 1 1 1
tail(simdata@variable.annotations$differential.expression)
## [1] 0 0 0 0 0 0 
sum(simdata@variable.annotations$differential.expression)
## [1] 1000

どの実験群で高発現であるかどうかの情報も保存されている。fraction.upregulaed = 0.2 としているため、第 2 群で高発現する発現変動遺伝子の数は 200 個と期待される。

sum(simdata@variable.annotations$upregulation)
## [1] 200
sum(simdata@variable.annotations$downregulation)
## [1] 800

生成されたシミュレーションデータの MA プロットは以下のようになる。赤でプロットされた遺伝子は発現変動遺伝子である。

sSeq パッケージによるシミュレーションデータの生成

sSeq パッケージ中の sim 関数を利用したシミュレーションデータの生成方法。上述のいくつかの方法では遺伝子発現量の期待値を実データから自動的にサンプリングして作るが、sSeq の場合は自分で用意する必要がある。

library(sSeq)
packageVersion("sSeq")
## [1] ‘1.6.0’

ng                   <- 10000       # 遺伝子数 
mean_DE              <- 2           # 発現変動遺伝子の発現量の差
sd_DE                <- 1           # 発現変動遺伝子の発現量の差のばらつき(?)
true_isDE_proportion <- 0.10        # 発現変動遺伝子の割合
conds      <- c("A", "A", "A", "B", "B", "B")    # 各実験群の biological replicate 数
true_mean1 <- round(rexp(ng, rate=1/200))        # 遺伝子発現量の期待値を用意
alpha      <- function(m){1/(m+100)}             # 発現量のばらつき
s0         <- runif(length(conds), 0, 2)         # size factor (ライブラリーを調整する因子)

set.seed(2015)
simdata <- sim(ng = ng, conds = conds, mean_DE = mean_DE, sd_DE = sd_DE,
               true_mean1 = true_mean1, alpha = alpha, s0 = s0,
               true_isDE_proportion = true_isDE_proportion)

生成されたデータは simdata$countsTable に保存されている。

head(simdata$countsTable)
##         A1  A2  A3   B4  B5   B6
## 1_FALSE  21  27  21   58  13   23
## 2_FALSE 208 310 257  680  71  338
## 3_TRUE    1   2   0 4155 413 1953
## 4_FALSE 217 305 229  675  82  351
## 5_FALSE 303 453 351 1038 122  484
## 6_FALSE 110 131 113  329  47  152

その遺伝子が発現変動遺伝子なのか否かの情報は simdata$true_isDEG に保存されている。

head(simdata$true_isDE)
## [1] FALSE FALSE  TRUE FALSE FALSE FALSE

生成されたシミュレーションデータの MA プロットは以下のようになる。赤でプロットされた遺伝子は発現変動遺伝子である。

References

  1. Sun J, Nishiyama T, Shimizu K, Kadota K. TCC: an R package for comparing tag count data with robust normalization strategies. BMC Bioinformatics. 2013, 14:219. PubMed Abstract
  2. Zhou X, Lindsay H, Robinson MD. Robustly detecting differential expression in RNA sequencing data using observation weights. Nucleic Acids Research. 2014, 42(11):e91. PubMed Abstract
  3. Soneson C, Delorenzi M. A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics. 2013, 14:91. PubMed Abstract
  4. Yu D, Huber W, Vitek O. Shrinkage estimation of dispersion in Negative Binomial models for RNA-seq experiments with small sample size. Bioinformatics. 2013, 29(10):1275-82. PubMed Abstract