DESeq 正規化

DESeq の正規化は size factor と呼ばれる係数を求めることによって行います。TMM 法などの normalization factor とほぼ同じように使う。ただし、後で説明するように両者は同一なものではない。

m 個のライブラリーが存在するとき、ライブラリー j の size factor は次のように求める。また、kij はライブラリー j における遺伝子 i のカウントデータを表す。

\[ s_{j} = median \frac{k_{ij}}{\left(\prod_{v=1}^{m}k_{iv}\right)^{\frac{1}{m}}} \]

このライブラリーサイズを利用して、データを正規化するとき、次の手順により正規化後のカウントデータを qij を計算する。

\[q_{i\rho} = \frac{k_{ij}}{s_{j}} \]

実際に R で実行した際のスクリプトは次のようになる。 サンプルデータは A 群と B 群の 2 群間を想定し、biological replicate が 2 とする。また、遺伝子数を 5 個とする。

x <- data.frame(
        A1 = c(20, 40, 30, 120, 80),
        A2 = c(25, 50, 30, 100, 90),
        B1 = c(25, 20, 30, 420, 90),
        B2 = c(20, 25, 30, 390, 80)
)

# size factors を求める
sizefactors <- apply(x / apply(x, 1, function(y){(prod(y))^(1/4)}), 2, median)
sizefactors
##       A1       A2       B1       B2 
## 0.942809 1.060660 1.060660 0.942809 



# 実際に DESeq パッケージで確かめてみる
library(DESeq)
d <- newCountDataSet(x, condition=c(1,1,2,2))
d <- estimateSizeFactors(d)
d@phenoData@data$sizeFactor
## [1] 0.942809 1.060660 1.060660 0.942809

normalization factors と size factors

発現変動遺伝子同定を行うパッケージとして edgeR と DESeq がよく知られている。両者では検定法が異なるほか、正規化法も異なっている。edgeR は TMM 法による正規化を行う。実際には、normalization factors とよばれる係数が計算され、正規化はこの係数とライブラリーサイズの積を利用して正規化を行う。一方、DESeq は上述の方法により正規化を行い、size factors とよばれる係数が計算され、正規化はこの係数を利用する。

normalization factors と size factors はともに正規化に関わる係数であるが、両者は異なる概念である。normalization factors はライブラリーサイズに依存するが、size factors は依存しない。従って、両者を直接に比較できない。もし、比較する場合は、normalization factors とライブラリーサイズの積を size factors と比較すべきである(参照:grokbase)。

library(edgeR)
library(DESeq)

# サンプルデータの読み込みと実験群のラベルの用意
sample_data <- read.table("https://bi.biopapyrus.jp/https://stats.biopapyrus.jp/glm/taylor-series.html", header = T)
group_label <- factor(c("A", "A", "A", "B", "B", "B"))
count_data <- sample_data[, -1]            # 2 列目以降のカウントデータを count_data に代入
rownames(count_data) <- sample_data[, 1]   # count_data に遺伝子名を付ける 
head(count_data)
##        A1  A2 A3 B1 B2 B3
## gene_1  4   2  2  0  0  0
## gene_2  1   1  0  9  3  1
## gene_3 33  54 75 17 13 12
## gene_4  0   0  4  0  1  0
## gene_5 24 146 54 23 15 22
## gene_6  3   2  6  0  0  0


# normalization factors を計算
d1 <- DGEList(count = count_data, group = group_label)
d1 <- calcNormFactors(d1)
norm.factors <- d1$samples$norm.factors
# normalization factors を比較するために調整する
adj_norm.factors <- norm.factors * colSums(count_data)
adj_norm.factors <- adj_norm.factors / mean(adj_norm.factors)   # 小数値に調整(係数なので、このように調整しても副効果はない)

# size factors を計算
d2 <- newCountDataSet(count_data, group_label)
d2 <- estimateSizeFactors(d2)
size.factors <- sizeFactors(d2)

# normlization factors と size factors の関係
plot(
  size.factors,
  adj_norm.factors,
  xlab = "size.factors",
  ylab = "norm.factors * lib.sizes", cex = 4, pch = 20,
  col = c("orange", "orange", "orange", "darkgrey", "darkgrey", "darkgrey")
)
size factors と normalization factors の関係

References

  • Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11(10):R106. PubMed Abstract
  • Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11(3):R25. PubMed Abstract