home <- [配布したフォルダをコピーしたパス]
setwd(home)

教師なし学習でクラスタリングと並んで重要なのが、次元縮約という方法論です。以下のように194サンプル(Tumor:97, Normal:97)の15121遺伝子の発現プロファイルを考えます。次元縮約の目的は、たくさんの変数(遺伝子)からなるデータを情報を可能な限り失わずに、少ない変数に要約し、解釈しやすくすることです。例えば、サンプル間の近さを可視化しようとすると、15121個の軸からなる座標空間を表示する必要がありますがそれは無理です。それを次元縮小で2,3個の軸に要約できれば、可視化が可能です。次元縮約の方法はたくさんありますが、今回は主成分分析(Principal component analysis;PCA)と多次元尺度構成法(Multidimensional scaling;MDS)を例に説明します。

input_file <- "data/data_expr_tumor_filtered.tsv"
data_tumor <- as.matrix(fread(input_file,header = TRUE,drop = 1))
colnames(data_tumor) <- paste0("Tumor",1:ncol(data_tumor))
input_file <- "data/data_expr_normal_filtered.tsv"
data_normal <- as.matrix(fread(input_file,header = TRUE,drop = 1))
colnames(data_normal) <- paste0("Normal",1:ncol(data_normal))
data_example <- cbind(data_tumor,data_normal)
rownames(data_example) <- paste0("Gene",1:nrow(data_example))
data_example <- t(data_example)
kable(data_example[1:10,1:10])
Gene1 Gene2 Gene3 Gene4 Gene5 Gene6 Gene7 Gene8 Gene9 Gene10
Tumor1 1.4159230 3.377278 0.7519159 98.35123 4.889642 15.20150 16.749664 0.3881095 1.4582108 139.74486
Tumor2 4.7858274 3.716744 0.2152176 100.90701 2.858270 16.31025 10.808428 0.3653588 0.3872724 24.27507
Tumor3 2.4368570 3.240557 0.3437818 40.82916 2.053996 16.05518 10.382871 0.0820921 0.2464288 35.18589
Tumor4 0.1958113 2.795901 0.8464482 46.42132 1.439745 13.99235 24.138840 0.1061655 0.5089498 37.40367
Tumor5 1.5079382 3.978028 0.2604621 234.46532 2.865083 11.50181 7.954990 0.7363401 0.5994274 30.36481
Tumor6 1.2348674 3.659223 0.4216911 327.46827 4.418752 13.27451 14.629319 0.3825851 0.9646250 19.29234
Tumor7 2.5003410 2.544354 0.8019765 180.69782 7.673741 10.38882 11.732746 0.2632887 0.5305637 19.26877
Tumor8 2.1449790 4.806979 0.2630553 246.70265 10.056010 14.63006 7.041516 0.2052525 2.3645559 24.56007
Tumor9 3.0588154 3.339891 6.4699673 89.54150 9.167577 10.79484 9.543649 0.2530093 2.4072224 17.29611
Tumor10 2.1071854 3.114608 0.5098996 109.82388 3.102324 23.12806 13.790673 0.5555666 0.6074557 23.08223

PCAから始めます。PCAでは、情報=データのばらつきと考えます。以下の図は、最初の20個の遺伝子のみについての対散布図です。遺伝子と遺伝子の間には、いくつかの相関関係が見られます。主成分分析では、このような構造も含めて、次元縮約を行います。さっそく15121個の遺伝子を2個の軸に要約してみましょう。今回はpcaMethodというパッケージ中のpca関数を用います。

pairs(data_example[,1:20])

library(pcaMethods)
(result_pca <- pca(data_example,nPcs = 2))
svd calculated PCA
Importance of component(s):
                PC1    PC2
R2            0.197 0.1297
Cumulative R2 0.197 0.3267
15121   Variables
194     Samples
0   NAs ( 0 %)
2   Calculated component(s)
Data was mean centered before running PCA 
Data was NOT scaled before running PCA 
Scores structure:
[1] 194   2
Loadings structure:
[1] 15121     2

要約した軸における値をスコアと呼びます。スコアを使って、2次元空間上でサンプルを表示して色分けしてみると、TumorとNormalの間で異なる傾向が見て取れます。

score <- result_pca@scores
score[1:10,]
                PC1         PC2
Tumor1   -2827.4559 -3197.35023
Tumor2     400.9473  8453.70884
Tumor3   -1571.8821  -655.33791
Tumor4   -4941.7396 -3255.54703
Tumor5     836.2170  2264.18164
Tumor6   -1399.2882   408.77359
Tumor7   -5014.9903    94.66862
Tumor8  -16852.6414 12812.95465
Tumor9   -3023.3209 -3246.54041
Tumor10  -3964.8193 -2034.51527
plot(score,cex=0.1)
text(score,labels=rownames(data_example),col=rep(c("black","blue"),each=ncol(data_tumor)))

得られた軸が何を表すかは必要に応じて、解析者が解釈を与えます。もともとの遺伝子の情報が軸に対して、どの程度、含まれているかは因子負荷量(Loadings)と呼ばれます。絶対値が大きい因子負荷量の遺伝子ほど、その軸に寄与しているという解釈ができます。例えば、Gene5はGene1よりもPC1に寄与しており、Gene4はGene5は符号が逆なので、何らかの逆の効果を持ってPC1に寄与しています。下図は各軸に対するloadingsを表示したものです。

loadings <- result_pca@loadings
loadings[1:10,]
                 PC1           PC2
Gene1  -4.510943e-05  9.483366e-05
Gene2   4.760426e-06 -1.528371e-05
Gene3  -1.171204e-04 -2.070240e-04
Gene4   5.027964e-03 -2.585365e-03
Gene5  -6.502928e-05 -3.937341e-05
Gene6  -4.661339e-05  2.453150e-04
Gene7   2.105485e-04 -8.920673e-06
Gene8   1.587441e-04 -6.743233e-05
Gene9  -2.211212e-05 -1.060542e-04
Gene10 -4.194461e-04  4.531030e-04
plot(loadings,cex=0.1)
text(loadings,labels=colnames(data_example),cex=0.5)