ex6: Clustering and PCA

ex6-1

データセットSato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt(61遺伝子x8遺伝子型）を使って、ユークリッド距離を使った場合とコサイン係数を使った距離でクラスタリングした時の違いを調べなさい（クラスタリング結果の違いの可視化はdendextendライブラリーを使うのが便利である）。このデータはシロイヌナズナ変異体にバクテリアを感染させた際の発現プロファイルを取り、野生型との比（log2)をとったものである。PR-1(At2g14610)は_P.syringae_に対する防御応答の主要なホルモンであるサリチル酸を介したシグナル伝達経路のマーカー遺伝子である。

コサイン係数での距離、クラスタリングは下記のカスタム関数を用いてよい。また、heatmapおよびheatmap.2(gplotsライブラリー）では引数に Rowv=as.dendrogram(“行のクラスタリング結果”), Colv=as.dendrogram(“列のクラスタリング結果”)と指定することで任意のクラスタリング結果でヒートマップを描くことができる。

準備

library(colorspace)
library(dendextend)
library(gplots)

source("~/data/MS/multivariate_analysis_source.R")

inputMatrix <- read.delim("~/data/MS/Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt", 
               header=TRUE, row.name=1)
heatmapColors <- colorpanel(10, low="blue", mid="white", high="orange")

実行

1. まずはユークリッド距離を使ってヒートマップを描いてみる。

heatmap.2(as.matrix(inputMatrix), 
          scale="none",           # 発現量比のスケーリング無し
          trace="none",           # heatmap.2デフォルトのトレースをキャンセル
          # ヒートマップのマス目設定
          sepcolor="black", colsep=0:ncol(inputMatrix), rowsep=0:nrow(inputMatrix),
          sepwidth=c(0.01, 0.01), 
          density.info="none",    # ヒストグラム
          col=heatmapColors,
          cexRow=(0.2 + 1/log10(nrow(inputMatrix)))/3*2,
          RowSideColors=ifelse(rownames(inputMatrix)=="At2g14610", "magenta", "grey")
)

2. 次にコサイン係数（ベクトルの角度）を距離尺度としたヒートマップを描いてみる。コサイン係数は関数が無いので自作する。

multivariate_analysis_source.R の中で定義してある。中身を確認してみよう。

上で、source()コマンドで読み込み済み。

3. 上記関数とas.dist関数を使ってコサイン係数の距離行列を求め、hclust関数でクラスタリングを実行する。

# 行のクラスタリング
rowClusters <- hclust(as.dist(cosine.table(as.matrix((t(inputMatrix))))))
# 列のクラスタリング
colClusters <- hclust(as.dist(cosine.table(as.matrix((inputMatrix)))))

4. ヒートマップを描く。Rowv, Colv引数にはas.dendrogram関数を介して上記クラスタリング結果を渡す。

heatmap.2(as.matrix(inputMatrix), Rowv=as.dendrogram(rowClusters), Colv=as.dendrogram(colClusters), 
scale="none", trace="none", sepcolor="black", colsep=0:ncol(inputMatrix), rowsep=0:nrow(inputMatrix), 
sepwidth=c(0.01, 0.01), density.info="none", col=heatmapColors, 
cexRow=(0.2 + 1/log10(nrow(inputMatrix)))/3*2, 
RowSideColors=ifelse(rownames(inputMatrix)=="At2g14610", "magenta", "grey"))

5. ユークリッド距離、コサイン係数を使ったクラスタリング結果の違いをdendextendパッケージに含まれる関数を使って可視化する。

rowClusters1 <- as.dendrogram(hclust(as.dist(dist(as.matrix((inputMatrix))))))
rowClusters2 <- as.dendrogram(hclust(as.dist(cosine.table(as.matrix(t(inputMatrix))))))
rowDendrogramList <- dendlist(rowClusters1, rowClusters2)
png("tanglegram_row.png", width=480*4, height=480*2, res=200)
tanglegram(rowDendrogramList, common_subtrees_color_branches = TRUE, 
columns_width= c(10,3,10), lab.cex=0.6, lwd=2, main_left="Euclidian", main_right="Cosine")
dev.off()

ex6-2

同じデータセットを主成分分析を用いて解析しなさい。この演習ではPR-1(At2g14610)の主成分得点、npr1-1, sid2-2の負荷量などサリチル酸シグナル伝達経路に注目して主成分分析の結果を評価しなさい。

1. 主成分分析を行うpca関数を定義する。

（出典：青木繁伸先生のwebsiteより、http://aoki2.si.gunma-u.ac.jp/R/pca.html )

multivariate_analysis_source.R の中で定義してある。中身を確認してみよう。

上で、source()コマンドで読み込み済み。

2. 主成分分析を実行する。

PCAresults <- pca(inputMatrix)

3. サリチル酸シグナル伝達経路のマーカー遺伝子(PR-1: At2g14610), 変異体（npr1-1, sid2-2)を手がかりに。結果を主成分得点、因子負荷量を用いて評価しなさい。

# 主成分得点（scores）から評価
png("ex6-2_PCA_scores.png", width=480*7, height=480*4, res=200)
par(mfrow=c(4,7))
for (kPC1 in 1:(ncol(inputMatrix)-1)){
    for (kPC2 in (kPC1+1):ncol(inputMatrix)){
        plot.pca(PCAresults, which="scores", pc.no=c(kPC1, kPC2), cex=0.5, markers="At2g14610")
    }
}

dev.off()

# 因子負荷量（loadings）から評価
png("ex6-2_PCA_loadings.png", width=480*7, height=480*4, res=200)
par(mfrow=c(4,7))
for (kPC1 in 1:(ncol(inputMatrix)-1)){
    for (kPC2 in (kPC1+1):ncol(inputMatrix)){
        plot.pca(PCAresults, which="loadings", pc.no=c(kPC1, kPC2), cex=0.5, 
                 markers=c("npr1.1", "sid2.2"))
    }
}

dev.off()

ex6-3

Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_full.txt（484遺伝子x22遺伝子型の実データセット）を使って同じ解析をし、より複雑なデータセットを使った場合のクラスタリング結果の違いを検討しなさい。

ex6-4

（発展問題: k-means法はトレーニングコースでは解説していません） coi1, dde2, jar1, jin1はジャスモン酸シグナル伝達経路、ein2-1, ein3はエチレンシグナル伝達経路、npr1-1, sid2-2はサリチル酸シグナル伝達経路に関わる遺伝子の変異体である。k-means法はデータをいくつのクラスターに分ければよいか見当が付いている場合によいクラスリング方法である。これらの変異体遺伝子発現プロファイルをk-means法を使って解析し、クラスター数の妥当性を考察せよ。kの数は3から始めなさい。同じ処理

例　

kmeans(t(inputMatrix), centers=3)$cluster

を繰り返し、結果の安定性やクラスタリングのされ方を指標に結果を評価しなさい。

ヒント: centers引数、iter.max引数を調整することにより、結果が変化します。