多群間比較(DESeq)

DESeq を利用して発現変動遺伝子の検出する例を示す。R を起動して DESeq を呼び出す。呼び出しが完了したら DESeq のバージョンをチェックする。

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

サンプルデータ

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

遺伝子 A 群で高発現 B 群で高発現 C 群で高発現
gene_1 ~ gene_50 Yes No No
gene_51 ~ gene_100 No Yes No
gene_101 ~ gene_150 No No Yes
gene_151 ~ gene_200 Yes Yes No
gene_201 ~ gene_250 No Yes Yes
gene_251 ~ gene_300 Yes No Yes
gene_201 ~ gene_1000 No No No

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

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

dim(count)
## [1] 1000    9

head(count)
##         A1  A2   A3  B1  B2  B3  C1  C2  C3
## gene_1   3   2   19   0   0   0  13   0   0
## gene_2 680 807 1011 168 132 138 135 188 202
## gene_3  35  38   30  11   9  11   8   9   8
## gene_4 103 165  124  35  37  31  36  36  30
## gene_5   2   9    0   0   0   0   0   0   0
## gene_6 419 373  478 114 117 109  89 142 106

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

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

検定(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] 4.516739e-02 1.448991e-08 2.662722e-02 3.392796e-04 4.849386e-02 2.218799e-05

解析結果を整形して、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")

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

nbr_count <- count[, c(1, 4, 7)]
nbr_group <- group[c(1, 4, 7)]

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.1747847 0.2034519 0.4908590 0.5176419 0.6087631 0.2477292

References

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