教師なし学習（次元縮約）

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)

PCAで得られた軸が、もともとのデータのばらつきをどの程度説明できているかは、寄与率（それを累積した累積寄与率）を用います。今回は、PC1とPC2で32.6%のばらつきを説明できているということです。15121個の遺伝子をたった2個の軸にまとめていることを考えると、３割程度も説明できていれば、見方によっては十分な気がします。

(cumprop <- result_pca@R2cum)

    PC1     PC2 
0.19701 0.32669 

どれくらいの数の軸をとれば良いかは、累積寄与率をみて判断します。nPcs=10にすれば、10個の軸に要約できます。累積寄与率を見てみると、PC5までで、60%のばらつきを説明できていることになります。

result_pca2 <- pca(data_example,nPcs=10)
(cumprop2 <- result_pca2@R2cum)

    PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8     PC9    PC10 
0.19701 0.32669 0.41836 0.50605 0.57372 0.63551 0.67875 0.71249 0.74496 0.76931 

PCAの原理はとてもシンプルで、軸の回転です。一つ目の軸を選ぶには、まず一番ばらつきの大きい方向を15121個の軸からなる空間の中で見つけ、その方向を一つ目の軸(PC1)とします。二つ目の軸は、PC1に直角となるという制約のもとで、次にばらつきの大きい方向を見つけて、二つ目の軸(PC2)とします。これをすべての繰り返していけば最大で15121個の新しい軸を構成できます。このようにとった軸は、データのばらつきを説明できている順番に並べることができので、そのうちの上位を主成分として用いるわけです。

軸を回転させただけですから元々のデータ全体のばらつきと各軸のばらつきの和は同じです。

\[\sum_{i=1}^{15121}\mbox{var}(\mbox{gene}_i)=\sum_{i=1}^{15121}\mbox{var}(\mbox{PC}_i)\]

しかも、元々のデータの相関構造も壊していないことが理解できると思います。 ちなみに、先ほどの寄与率とは、全体のばらつきに対する各軸のばらつきの割合です。

\[\mbox{Proportion} = \frac{\mbox{var}(\mbox{PC}_i)}{\sum_{i=1}^{15121}\mbox{var}(\mbox{PC}_i)}\]

MDSについて説明します。PCAは多数の変数を、データのばらつきを最も説明する方向へ軸を回転させることで新しい変数を作り、寄与率の高い上位の軸を選ぶ方法でした。MDSでは、元々のデータにおける個体(サンプル)間の距離をできるだけ保持するようなユークリッド空間を構成するアプローチで新しい軸を構成します。例えて言えば、都市間の距離のみわかっている状況から地図を構成するような方法です。データのばらつきは、言い換えれば個体間の距離なので、結果はPCAと似たものとなります。

MDSを実行するには、cmdscale関数を用います。

d <- dist(data_example)
result_mds <- cmdscale(d,k = 2)

plot(result_mds,xlab="x1",ylab="x2",cex=0.1)
text(result_mds,labels = rownames(data_example),col=rep(c("black","blue"),each = ncol(data_tumor)))

先ほどのPCAと上下が逆になっており、また座標値も異なりますが、配置はほぼ同じだということがわかります。 PCAもMDSもほぼ同じ結果を与えます。違う点は、MDSは個体間の距離さえわかっていれば、座標空間を構成できるという点です。単純な表形式のデータだといいのですが、より複雑な形式のデータである場合には、明示的に表形式にはならないが、個体間の非類似性（距離）は計算できる場合には有効です。

このような可視化は、色々な場面で役に立ちます。例えば、クラスター分析の結果を確かめる場合に便利です。実際、MDSの場合は、クラスター分析で用いた非類似度に基づいて、座標空間を構成するので、クラスター分析の結果と点の配置が対応する傾向があります。また、TumorとNormalのサンプルの近さを見たように、明らかにおかしなサンプルをはずれ値として視覚的に見つけることも可能でしょう。さらに別の使い方としては、縮約した数次元の軸は、データのばらつきを最もよく表す軸であると考えられます。それを新たな特徴量を持ったデータと考えて、クラスター分析をしたりするなどの使い方も考えられます。