主成分分析

マイクロアレイあるいは RNA-Seq を用いた解析で、多数のサンプル(またはライブラリー)を対象としているとき、サンプル同士の類似度を調べたい場合がある。類似度を調べる方法としてはクラスタリング、主成分分析(PCA)、独立成分分析(ICA)、非負行列分解(NMF)などの方法が挙げられる。ここで乱数の生成により作成したデータを利用して PCA を行ってみる。

サンプルデータ

サンプルデータは条件を設定して負の二項分布に従うように乱数生成により作成した。A 群と B 群の 2 つの実験群からなる。A 群には薬剤 A を投与し、0 時間、3 時間、6 時間の 3 つの時刻においてサンプルを採取した。B 群には薬剤 B を投与し、0 時間、3 時間、6 時間の 3 つの時刻においてサンプルを採取した。各実験の各時刻において 2 つの biological replicate があり、実験全体で計 12 biological replicate を使用したとする。また、全遺伝子数は 1,000 であり、その発現パターンは以下のように設定した。

遺伝子 A 群の発現パターン B 群の発現パターン
gene_1 ~ gene_40 A6h > A3h > A0h, B6h > B3h > B0h
A0h = B0h, A3h = B3h, A0h = B0h
gene_41 ~ gene_80 A6h = A3h > A0h A0h = B0h = B3h = B6h
gene_81 ~ gene_120 A0h = A3h = A6h A0h = B6h = B3h > B0h
gene_121 ~ gene_160 A0h = A6h < A3h A0h = B0h = B3h = B6h
gene_161 ~ gene_200 A0h = A6h < A3h A0h = B0h = B3h < B6h
gene_201 ~ gene_1000 A0h = A3h = A6h = B0h = B3h = B6h

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

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

dim(count)
## [1] 1000    12

head(count)
##        A1_0h A2_0h A3_3h A4_3h A5_6h A6_6h B1_0h B2_0h B3_3h B4_3h B5_6h B6_6h
## gene_1     1     7    11     0    22    24     0     0     9     3    22    16
## gene_2   170   182   744   683   918   637   128   112   645   517  1248  1234
## gene_3     8    12    37    46   100    94     9     8    30    46    47    63
## gene_4    44    38   109   187   354   225    40    97   103   146   505   343
## gene_5     0     0     5     0     4    10     0     0     1     0     7     4
## gene_6    88   112   372   459  1010   611   129    80   486   360   678   767

ここで、分散のもっとも大きい遺伝子上位 200 個に対して、PCA を行う。

count.log <- log10(count + 1)
count.log.var <- apply(count.log, 1, var)
count.log.pca <- count.log[order(count.log.var, decreasing = TRUE), ][1:200, ]

res <- prcomp(count.log.pca)

PCA の結果が res に保存される。

summary(res)
## Importance of components:
##                           PC1     PC2     PC3     PC4    PC5     PC6     PC7     PC8     PC9   PC10    PC11    PC12
## Standard deviation     2.3707 0.60435 0.57019 0.51679 0.4882 0.47561 0.46070 0.43716 0.42221 0.4073 0.38261 0.34214
## Proportion of Variance 0.6979 0.04535 0.04037 0.03316 0.0296 0.02809 0.02635 0.02373 0.02214 0.0206 0.01818 0.01454
## Cumulative Proportion  0.6979 0.74325 0.78362 0.81678 0.8464 0.87447 0.90082 0.92455 0.94669 0.9673 0.98546 1.00000

各主成分の座標を取り出して散布図を描く。サンプルデータは A 群と B 群の 2 タイプがあり、各群にはさらに 0h、3h、6h のデータがあるので、これらが区別しやすいように散布図を描く。

col <- c("red", "red", "red", "red", "red", "red", "blue", "blue", "blue", "blue", "blue", "blue")
time <- c("0", "0", "3", "3", "6", "6", "0", "0", "3", "3", "6", "6")

pc1 <- res$rotation[, 1]
pc2 <- res$rotation[, 2]
pc3 <- res$rotation[, 3]

plot(pc1, pc2, type = "n", asp = TRUE)
text(pc1, pc2, labels = time, col = col)

plot(pc2, pc3, type = "n", asp = TRUE)
text(pc2, pc3, labels = time, col = col)

遺伝子の発現量について、A 群と B 群はそれほど異なっているわけではないことが、第 1 主成分から読み取られる。また、第 2 主成分および第 3 主成分からは、A 群と B 群の遺伝子の発現量は、0 時間において似ていることがわかる。また、3 時間後および 6 時間後では、A 群と B 群の遺伝子の発現量に大きな違いが見られるようになる。

PCA結果  PCA結果

各成分の寄与率は以下のようにプロットできる。

propvar <- summary(res)$importance[2, ]
propvar
##     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8     PC9    PC10    PC11    PC12
## 0.69789 0.04535 0.04037 0.03316 0.02960 0.02809 0.02635 0.02373 0.02214 0.02060 0.01818 0.01454

barplot(propvar)
PCA結果