二群間比較(DESeq)

DESeq を利用して発現変動遺伝子の検出する例を示す。二群間比較では nbinomTest と尤度比検定の 2 種類の検定法を利用できる。

R を起動して DESeq を呼び出す。呼び出しが完了したら DESeq のバージョンをチェックする。このページではバージョン 1.20.0 の解析結果を示している。バージョンが異なると解析結果が異なる場合がある。

library(DESeq)
packageVersion("DESeq")
##[1] ‘1.20.0’

サンプルデータ

サンプルデータは、カウントデータが負の二項分布に従うように乱数を生成して作成した。A 群と B 群の 2 群からなり、それぞれの実験群には 3 つの biological replicate が含まれる。また、全遺伝子数は 1,000 個であり、その発現パターンは以下のように設定した。

gene_1 ~ gene_160 B 群に比べ、A 群での発現量が高い
gene_161 ~ gene_200 A 群に比べ、B 群で発現量が高い
gene_201 ~ gene_1000 両群において差異なし

サンプルデータを読み込んでから行列型に変換し count 変数に代入する。

count <- read.table("http://bi.biopapyrus.net/data/count2g.txt", 
                    sep = "\t", header = T, row.names = 1)
count <- as.matrix(count)

dim(count)
## [1] 1000    6

head(count)
##         A1  A2  A3  B1 B2 B3
## gene_1   3   2  40   4  2  0
## gene_2 680 807 594 165 91 73
## gene_3  35  38  31   7  9 11
## gene_4 103 165 110  38 31 28
## gene_5   2   9   1   3  1  0
## gene_6 419 373 399 133 89 58

実験群の識別ラベルを作成し、group 変数に代入する。

group <- factor(c("A", "A", "A", "B", "B", "B"))

nbinomTest

カウントデータと実験群の識別ラベルを 1 つの変数に保存する。

cds <- newCountDataSet(countData = count, conditions = group)

データを正規化するためのサイズファクターを計算し、その次に分散を計算する。

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

最後に、nbinomTest により発現変動遺伝子を検出する。A 群と B 群の比較が目的であるから、以下のように順に引数を与える。

res <- nbinomTest(cds, "A", "B")

結果は res に保存される。

head(res)
##       id   baseMean  baseMeanA  baseMeanB foldChange log2FoldChange         pval        padj
## 1 gene_1   8.059985  14.035830   2.084139  0.1484871      -2.751591 2.996011e-01 1.000000000
## 2 gene_2 383.730021 653.208090 114.251952  0.1749090      -2.515324 3.284198e-07 0.000108707
## 3 gene_3  20.998313  32.633121   9.363505  0.2869326      -1.801216 3.421741e-02 0.263394460
## 4 gene_4  76.197165 118.734093  33.660237  0.2834926      -1.818617 1.861466e-03 0.024645816
## 5 gene_5   2.589502   3.788293   1.390712  0.3671078      -1.445724 5.757148e-01 1.000000000
## 6 gene_6 235.296812 373.396805  97.196818  0.2603044      -1.941729 1.001718e-04 0.002786959

解析結果をファイルに保存するには write.table を利用する。

write.table(res, file = "result.txt", col.names = T, row.names = T, sep = "\t")

biological replicates がない場合

A 群に 1 サンプル、B 群に 1 サンプルあるような biological replicate がない場合は estimateDispersions を実行するときエラーがでる場合がある。

nbr_count <- count[, c(1, 4)]
nbr_group <- group[c(1, 4)]
cds <- newCountDataSet(countData = nbr_count, conditions = nbr_group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
## Error in .local(object, ...) : 
##  None of your conditions is replicated. Use method='blind' 
##  to estimate across conditions, or 'pooled-CR', if you have 
##  crossed factors.

このとき、エラー文をみると method = "blind" あるいは method = "pooled-CR" と指定して、群間のサンプルを一時的に biological replicate とみなして分散を推測する。

cds <- estimateDispersions(cds, method = "blind")
res <- nbinomTest(cds, "A", "B")
## Warning message:
## In nbinomTest(cds, "A", "B") :
##   You have used 'method="blind"' in estimateDispersion without
##   also setting 'sharingMode="fit-only"'. This will not yield 
##   useful results.

結果は res に保存されるが、警告メッセージをみると estimateDispersions 関数で method = "blind" を利用するとき、sharingMode = "fit-only" も合わせて利用した方がいいと書いてある。そこで、再実行する。

cds <- newCountDataSet(countData = nbr_count, conditions = nbr_group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only")
res <- nbinomTest(cds, "A", "B")
head(res)
##       id   baseMean  baseMeanA  baseMeanB foldChange log2FoldChange      pval padj
## 1 gene_1   3.515138   2.918317   4.111958  1.4090168      0.4946888 0.9893173    1
## 2 gene_2 415.551790 661.485298 169.618282  0.2564203     -1.9634174 0.3403113    1
## 3 gene_3  20.621482  34.047037   7.195927  0.2113525     -2.2422768 0.4534654    1
## 4 gene_4  69.629586 100.195567  39.063604  0.3898736     -1.3589217 0.5486856    1
## 5 gene_5   2.514757   1.945545   3.083969  1.5851439      0.6646138 0.9914388    1
## 6 gene_6 272.157146 407.591677 136.722615  0.3354402     -1.5758727 0.4448613    1

尤度比検定(nbinomGLMTest)

尤度比検定の詳細に関して説明を省略する。newCountDatSet を作成するとき、実験群の識別ラベルは conditions 引数にベクトルの形で与えている。このとき、full model は ~ condition(最後に s が付かない!)とすれば、conditions の情報をもとにモデルが作成される。reduced model はすべての情報を省略するという意味で ~ 1 で指定する。

cds <- newCountDataSet(countData = count, conditions = group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
f.model <- fitNbinomGLMs(cds, count ~ condition)
r.model <- fitNbinomGLMs(cds, count ~ 1)
res <- nbinomGLMTest(f.model, r.model)
head(res)
## [1] 2.251768e-01 2.665830e-07 2.826232e-02 1.581978e-03 4.651769e-01 8.747931e-05

biological replicate がないとき、estimateDispersions 関数に与える引数を上述のようにオプションを与えて実行すればよい。

cds <- newCountDataSet(countData = nbr_count, conditions = nbr_group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only")
f.model <- fitNbinomGLMs(cds, count ~ condition)
r.model <- fitNbinomGLMs(cds, count ~ 1)
res <- nbinomGLMTest(f.model, r.model)
head(res)
## [1] 0.9040478 0.2945206 0.3625542 0.4958876 0.8882824 0.3987974

解析結果を整形して、result.txt ファイルに保存する。

table <- data.frame(
  gene.name = rownames(count),
  p.value = res,
  q.value = p.adjust(res, method = "BH")
)
write.table(table, file = "result.txt", col.names = T, row.names = F, sep = "\t")

References

  • Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11(10):R106. PubMed Abstract