ex601: Clustering and PCA

ex601-1: 階層クラスタリング

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

コサイン係数での距離、クラスタリングは下記のカスタム関数を用いてよい。また、heatmapおよびheatmap.2(gplotsライブラリー）では、引数に

• Rowv=as.dendrogram(行のクラスタリング結果オブジェクト)
• Colv=as.dendrogram(列のクラスタリング結果オブジェクト) と指定することで任意のクラスタリング結果でヒートマップを描くことができる。

準備

library(dendextend)
library(gplots)

source("nibb_gitc_rnaseq_20210310/gitc/data/MS/multivariate_analysis_source.R")

inputMatrix <- read.delim(
  "nibb_gitc_rnaseq_20210310/gitc/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. ユークリッド距離、コサイン係数を使ったクラスタリング。コサイン係数は自作関数cosine.tableを使って計算している。

rowClusters1 <- hclust(
  as.dist( dist( as.matrix(inputMatrix) ) ) 
)
rowClusters2 <- hclust(
  as.dist( cosine.table( as.matrix( t(inputMatrix) ) ) )
)

2. ユークリッド距離でのクラスタリング結果の可視化

heatmap.2(as.matrix(inputMatrix), 
          Rowv=as.dendrogram(rowClusters1),
          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")
)

3. 次にコサイン係数（ベクトルの角度）を距離尺度としたクラスタリング結果を可視化する。

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

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

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

# 行のクラスタリング: rowClusters2として計算済み
# 列のクラスタリング：
colClusters <- hclust(
  as.dist( cosine.table( as.matrix(inputMatrix) ) )
)

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

heatmap.2(
  as.matrix(inputMatrix), 
  Rowv=as.dendrogram(rowClusters2), 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パッケージに含まれる関数を使って可視化する。

rowDendrogramList <- dendlist(as.dendrogram(rowClusters1), as.dendrogram(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()

ex601-2

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

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

(pca関数はmultivariate_analysis_source.R の中で定義している。 ここまでの作業をしていれば、source()コマンドで読み込み済み。 http://aoki2.si.gunma-u.ac.jp/R/pca.html のpca関数を改変して使用している)

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

PCAresults <- pca(inputMatrix)

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

# 主成分得点（scores）から評価
png("ex601-2_PCA_scores.png", 
     width=480*14, height=480*8, 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("ex601-2_PCA_loadings.png", 
    width=480*14, height=480*8, 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()

ex601-3

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