ex601 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki
ex601: Clustering and PCA
ex601-0: データ可視化
データセットSato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt(61遺伝子x8遺伝子型）をRで読み込んで、image関数、heatmap関数で可視化し、image関数とheatmap関数の違いを確認する。
# getwd() で現在のディレクトリを確認すること。
dataset_dir <- "~/Desktop/gitc/data/MS/dataset"
data_small_path <- file.path(dataset_dir, "Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt")
inputMatrix <- read.delim(
data_small_path,
header=TRUE, row.name=1
)
str(inputMatrix) # 読み込んだデータの構造を確認する
image 関数で読み込んだ数値行列をカラーコードとして可視化する。
image(t(inputMatrix))
heatmap 関数で階層クラスタリングしてデータを並び替えて可視化する。
heatmap(as.matrix(inputMatrix))
ex601-1: 階層クラスタリング
データセットSato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt(61遺伝子x8遺伝子型）を使って、ユークリッド距離、コサイン係数でクラスタリングした時の違いを調べなさい。なお、クラスタリング結果の違いの可視化はdendextendライブラリーを使うのが便利である。このデータはシロイヌナズナ変異体に_P.syringae_を感染させた際の発現プロファイルを取り、野生型との対数比（底：2)である。PR-1(At2g14610)は_P.syringae_に対する防御応答の主要なホルモンであるサリチル酸を介したシグナル伝達経路のマーカー遺伝子である。
コサイン係数での距離、クラスタリングは下記のカスタム関数を用いてよい。また、heatmapおよびheatmap.2(gplotsパッケージ）では、引数に
Rowv=as.dendrogram(行のクラスタリング結果オブジェクト)Colv=as.dendrogram(列のクラスタリング結果オブジェクト) と指定することで任意のクラスタリング結果でヒートマップを描くことができる。
準備
library(coop)
library(dendextend)
library(gplots)
source("~/Desktop/gitc/data/MS/multivariate_analysis_source.R")
実行
1. ユークリッド距離、コサイン係数を使ってクラスタリングする。
rowClusters1 <- hclust(
as.dist( dist( as.matrix(inputMatrix) ) )
)
rowClusters2 <- hclust(
as.dist(1-tcosine(as.matrix(inputMatrix)))
)
それぞれのクラスタリング結果を可視化する。
layout(matrix(1:2, ncol=2))
plot(rowClusters1)
plot(rowClusters2)
[補足] コサイン係数の扱いを検討する。コサイン係数は-1から1までの範囲の有理数をとり、コサイン係数=1は２つのプロファイルの角が0°、コサイン係数=-1は180°、コサイン係数=0は90°（直角、独立）である。今回のデータは野生型と変異体の発現プロファイルの対数比であり、ある遺伝子Aから計算したコサイン係数=1の遺伝子Bとコサイン係数=-1の遺伝子Cがある場合には、遺伝子B、Cは対数比の+/-が反転した遺伝子である。
研究者がこれらA、B、Cをどのようにクラスリングしたいかは目的による。ある転写制御メカニズムによって同時に制御される状況（活性化・抑制）をまとめる形で解析したい場合、分けて解析したい場合などがそれにあたる。ここでは、それぞれを想定してコサイン係数の扱いを例示する。受講者には可視化による確認をしてもらいたい。
# 活性化・抑制を同じクラスターにまとめたい場合
rowClusters2.1 <- hclust(
as.dist(1-abs(tcosine(as.matrix(inputMatrix))))
)
layout(matrix(1:2, ncol=2))
plot(rowClusters2) # 活性化・抑制を分けてクラスターしたい場合
plot(rowClusters2.1)
2. ユークリッド距離でのクラスタリング結果の可視化
heatmapColors <- colorpanel(10, low="blue", mid="white", high="orange")
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. 次にコサイン係数（ベクトルの角度）を距離尺度としたクラスタリング結果を可視化する。
コサイン係数を距離尺度とした行データ（遺伝子）のクラスタリング結果は上記でrowClusters2として定義してある。
列データ（サンプル）も同様にコサイン係数の距離行列を求め、hclust関数でクラスタリングを実行する。
# 行のクラスタリング: rowClusters2として計算済み
# 列のクラスタリング：
colClusters <- hclust(
as.dist(1-cosine(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()コマンドで読み込み済み。)
source("~/Desktop/gitc/data/MS/multivariate_analysis_source.R")
2. 主成分分析の実行内容を理解できる低次元データで試してみる。
２次元で可視化できるようにするために２サンプル分だけのデータで主成分分析を行う。ここではサリチル酸シグナル伝達変異体のnpr1-1とsid2-2を用いる。
test_PCA <- PCA(inputMatrix[,7:8])
plot(inputMatrix[,7], inputMatrix[,8]) #元データ
abline(h=0, v=0)
plot.PCA(test_PCA,"scores") #主成分得点
主成分得点が元データ（npr1-1, sid2-2の遺伝子発現データ）を、主成分分析で回転させた座標に投射したデータになったことがわかる。 主成分分析はこのようにデータの分散を説明する軸を見つけて回転させるアルゴリズムである。
3. 主成分分析を実行する。
PCAresults <- PCA(inputMatrix)
4. サリチル酸シグナル伝達経路のマーカー遺伝子(PR-1: At2g14610), 変異体（npr1-1, sid2-2)を手がかりに。結果を主成分得点、因子負荷量を用いて評価しなさい。
# 因子負荷量（loadings）と主成分得点(scores)から評価
for (kPC1 in 1:(ncol(inputMatrix)-1)){
for (kPC2 in (kPC1+1):ncol(inputMatrix)){
png(paste0("ex601-2_PCA_loadings_PC", kPC1, "-", kPC2, ".png"),
width=210/2, height=210/2, unit="mm", res=200 # 適宜、図の大きさ・解像度はこれら３つの引数で調節する
)
plot.PCA(PCAresults, which="loadings", pc.no=c(kPC1, kPC2), cex=0.5, markers=c("npr1.1", "sid2.2"))
dev.off()
png(paste0("ex601-2_PCA_scores_PC", kPC1, "-", kPC2, ".png"),
width=210/2, height=210/2, unit="mm", res=200 # 適宜、図の大きさ・解像度はこれら３つの引数で調節する
)
plot.PCA(PCAresults, which="scores", pc.no=c(kPC1, kPC2), cex=0.5, markers=c("At2g14610"))
dev.off()
}
}
ex601-3
Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_full.txt（484遺伝子x22遺伝子型の実データセット）を使って同じ解析をし、より複雑なデータセットを使った場合のクラスタリング結果の違いを検討しなさい。
source関数で読み込んだ関数のコード内容
multivariate_analysis_source.R
PCA <- function(input_matrix, scale=FALSE, center=FALSE){
input_rownames <- rownames(input_matrix)
INPUT_NROW <- nrow(input_matrix) # sample size
INPUT_NCOL <- ncol(input_matrix) # No of variables
if (is.null(input_rownames)) input_rownames <- paste("#", 1:INPUT_NROW, sep="")
input_matrix <- subset(input_matrix, complete.cases(input_matrix)) # remove cases with NAs
if (is.null(colnames(input_matrix))) {
colnames(input_matrix) <- paste("X", 1:INPUT_NCOL, sep="")
}
variable_names <- colnames(input_matrix)
mean_values <- colMeans(input_matrix)
variances <- apply(input_matrix, 2, var)
SDs <- sqrt(variances)
r <- cor(input_matrix)
PCA_result <- prcomp(input_matrix, scale=scale, center=center)
eigenvalues <- sqrt(PCA_result$sdev)
eigenvector <- PCA_result$rotation # this is rotation matrix
contribution <- eigenvalues/INPUT_NCOL*100
cumulative_contribution <- cumsum(contribution)
PCA_loadings <- t(sqrt(eigenvalues)*t(eigenvector))
PCA_scores <-
scale(input_matrix) %*% eigenvector * sqrt(INPUT_NROW/(INPUT_NROW-1))
names(mean_values) <-
names(variances) <- names(SDs) <-
rownames(r) <- colnames(r) <- rownames(PCA_loadings) <-
colnames(input_matrix)
names(eigenvalues) <-
names(contribution) <- names(cumulative_contribution) <-
colnames(PCA_loadings) <- colnames(PCA_scores) <-
paste("PC", 1:INPUT_NCOL, sep="")
return(structure(
list(
mean=mean_values, variance=variances,
standard.deviation=SDs, r=r,
factor.loadings=PCA_loadings, eval=eigenvalues,
evec=eigenvector, nr=INPUT_NROW, # added for subsequent PCA projection
contribution=contribution,
cum.contribution=cumulative_contribution, fs=PCA_scores),
class="pca"))
}
# print メソッド
print.PCA <- function(
pca_returned_object,
N_PCs_to_visualize=NULL,
digits_to_visualize=3)
{
eval <- pca_returned_object$eval
nv <- length(eval)
if (is.null(N_PCs_to_visualize)) {
N_PCs_to_visualize <- sum(eval >= 1)
}
eval <- eval[1:N_PCs_to_visualize]
cont <- eval/nv
cumc <- cumsum(cont)
fl <- pca_returned_object$factor.loadings[, 1:N_PCs_to_visualize, drop=FALSE]
rcum <- rowSums(fl^2)
vname <- rownames(fl)
max.char <- max(nchar(vname), 12)
fmt1 <- sprintf("%%%is", max.char)
fmt2 <- sprintf("%%%is", digits_to_visualize+5)
fmt3 <- sprintf("%%%i.%if", digits_to_visualize+5, digits_to_visualize)
cat("\nResult of PCA\n\n")
cat(sprintf(fmt1, ""),
sprintf(fmt2, c(sprintf("PC%i", 1:N_PCs_to_visualize), " Contribution")), "\n", sep="", collapse="")
for (i in 1:nv) {
cat(sprintf(fmt1, vname[i]),
sprintf(fmt3, c(fl[i, 1:N_PCs_to_visualize], rcum[i])),
"\n", sep="", collapse="")
}
cat(sprintf(fmt1, "Eigenvalue"), sprintf(fmt3, eval[1:npca]), "\n", sep="", collapse="")
cat(sprintf(fmt1, "Contribution"), sprintf(fmt3, cont[1:npca]), "\n", sep="", collapse="")
cat(sprintf(fmt1, "Cum.contrib."), sprintf(fmt3, cumc[1:npca]), "\n", sep="", collapse="")
}
# summary メソッド
summary.PCA <- function(pca_returned_object, output_digits=5){
print.default(pca_returned_object, digits=output_digits)
}
# plot メソッド
plot.PCA <- function(
pca_returned_object,
which=c("loadings", "scores"),
pc.no=c(1,2),
draw_axes=TRUE, # 座標軸を描き込むかどうか
label.cex=0.6, # 主成分負荷量のプロットのラベルのフォントサイズ
markers=NULL,
draw.names=FALSE,
col="black",
...) # plot に引き渡す引数
{
which <- match.arg(which)
if (which == "loadings") {
d <- pca_returned_object$factor.loadings
} else {
d <- pca_returned_object$fs
}
label <- sprintf("PC%i", pc.no)
plot(d[, pc.no[1]], d[, pc.no[2]],
xlab=label[1], ylab=label[2], col=col, ...)
if (which == "loadings") {
labelPosition <- ifelse(d[, pc.no[1]] < 0, 4, 2)
if (is.null(markers) == FALSE){
for (marker in markers){
points(x=d[marker, pc.no[1]], y=d[marker, pc.no[2]],
col="magenta", pch=16)
}
}
if (draw.names == TRUE){
text(x=d[, pc.no[1]], y=d[, pc.no[2]],
labels=rownames(pca_returned_object$factor.loadings),
pos=labelPosition, cex=label.cex
)
}
} else {
if (which == "scores" && is.null(markers) == FALSE){
for (marker in markers){
points(x=d[marker, pc.no[1]], y=d[marker, pc.no[2]], col="magenta", pch=16)
#text(x=d[marker, pc.no[1]], y=d[marker, pc.no[2]],
# labels=marker, col="red")# At2g14610
}
}
}
}