発展演習（Hard） - yoshida1983/protc2026 GitHub Wiki
実習では Mascot の Score を使いましたが、ここでは定量プロテオームデータで、PCA・クラスタリング・GO解析を行います。 表計算ソフトは使わず、Rだけで進めます。
- 測定法が違うため（ラベルフリー → TMT）、前処理の考え方から変わります
- 実習用スクリプトは、そのままでは使えません。どこを変える必要があるかを自分で判断してください
- 末尾に回答例があります。先に自分で組み立ててから確認してください
→ 先に 演習問題（Easy）、演習問題（Medium） に取り組むことをお勧めします
データ
- シロイヌナズナ葉の定量プロテオームデータ
- 出典：Yoshida et al. Plant J (2024)
- 論文のSupporting Information, Data S1のシートDを、本演習用に抜粋
- コントロール（Col; 野生型）、R6過剰発現体、R10過剰発現体の比較
- 反復数が違うことに注意（Col, 3; R6, 4; R10, 2）
- 値は正規化済みの log2 強度（再度の対数変換は不要）
- 10,497 タンパク質、うち ANOVA で有意に変化したもの 562
このデータは TMT で測定されています
- TMT10plex：10サンプルをそれぞれ異なる同位体タグで標識し、混ぜて1回の LC-MS/MS でまとめて測定する方法
- MS/MS で出るレポーターイオンの強度が、各サンプル由来の量に対応します
- mix チャネル：全サンプルを混ぜた参照で、正規化と品質確認に使います。群間の比較からは除きます
- 欠損値がほとんどありません：同じ測定の中で10チャネルを同時に読むためです。ラベルフリー（実習の Mascot 結果）で未検出が多数出るのとは対照的です
データの読み込み
ファイルを読むときの注意が3つあります。
- 見出しは3行目です（1行目はタイトル、2行目はTMTチャネル名）。2行読み飛ばします
- 列名が重複しています（
Col_Wが3つなど）。列名ではなく位置で取り出します AGI codeに欠損や重複がある行があります
library(tidyverse)
library(readxl)
# 3行目が見出しなので、2行読み飛ばす
d <- read_excel("tpj16765-sup-0001-DataS1_sheetD.xlsx", skip = 2)
dim(d) # 10497 × 30
names(d)
# 21〜30列目が10チャネル分の定量値。列名が重複しているので位置で取り出す
quant <- as.matrix(d[, 21:30])
rownames(quant) <- make.unique(d$`AGI code`)
colnames(quant) <- c(paste0("Col_", 1:3), paste0("R6_", 1:4),
paste0("R10_", 1:2), "mix")
head(quant)
ここから先は、実習用スクリプトを流用しながら自分で組み立ててください。
課題1「PCAで全体の傾向をつかむ」
目的：定量データのPCAを、実習との前処理の違いを意識して実行する。
課題
mixチャネルを除いた9サンプルで、PCAを実行する- スコアプロットを描き、3系統が分かれるか、同じ系統の反復がまとまるかを確認する
- PC1・PC2 の寄与率を軸ラベルに表示する
- ヒント：実習用スクリプト「PCA解析」の節をそのまま使えます。ただしlog2変換の行は不要です
scale. = TRUEとFALSEの両方を試し、結果を比較する- 考察：なぜ実習では
TRUE、ここではFALSEが既定なのか
- 考察：なぜ実習では
課題2「有意に変化したタンパク質でクラスタリング」
目的：統計的に有意な変動に絞り、系統間の違いを可視化する。
課題
ANOVA Significant列が+の行だけを取り出す（検定は済んでいるので、自分で計算する必要はありません）- 絞り込んだタンパク質でヒートマップを描き、3系統が分かれるかを確認する
- ヒント：実習用スクリプト「階層的クラスタリング」の節を使えます
- 距離尺度やクラスタリング手法を変えて、結果がどう変わるかを比較する
考察
- R6 と R10 では、変動の向きは同じですか、逆ですか
- R10 は反復が2つしかありません。結果を読むときに、どのような注意が必要ですか
課題3「GO解析で特徴を抽出する」
目的：背景遺伝子集合の設定を含めて、濃縮解析を行う。
課題
ANOVA Significantが+のタンパク質の AGIコードで、GO濃縮解析を行う- ヒント：
AGI code列があるので、UniProtの対応表は不要です
- ヒント：
- 背景（
universe）に、このデータで定量された全タンパク質を指定する- 全遺伝子を背景にすると、MSで検出されやすいタンパク質の term ばかりが上位に来ます（講義5）
universeを指定した場合としない場合で、結果を比較するontを"BP""CC""MF"と変えて、見えるものがどう変わるかを確認する
発展
- R6 と R10 で、変動の向きが揃うタンパク質／逆になるタンパク質を抽出してみる
-Log ANOVA p valueの列を使って、変化の大きさと有意性の関係を散布図にしてみる- 論文（Yoshida et al. 2024）の図と、自分の結果を見比べてみる
回答例
先に自分で組み立ててから、確認用に見てください。 あくまで一例であり、パラメーター設定やアルゴリズムの選択によって、結果の見え方が変わることに留意してください。
課題1：PCA
# mix を除いた9サンプル
mat <- quant[, 1:9]
group <- factor(c(rep("Col", 3), rep("R6", 4), rep("R10", 2)),
levels = c("Col", "R6", "R10"))
# すでに log2 済みなので、対数変換はしない
# 全サンプルで値が変わらないタンパク質は標準化できないので除く
mat <- mat[apply(mat, 1, sd) > 0, ]
dim(mat)
# スコアプロットを描く関数（scale. の比較に使う）
plot_pca <- function(p, title) {
df <- as.data.frame(p$x)
df$Sample <- rownames(p$x)
df$Group <- group
pct <- round(summary(p)$importance[2, 1:2] * 100, 1)
ggplot(df, aes(PC1, PC2, color = Group, label = Sample)) +
geom_point(size = 4) +
geom_text(vjust = -1.2, size = 4, show.legend = FALSE) +
theme_minimal() +
labs(title = title,
x = paste0("PC1 (", pct[1], "%)"),
y = paste0("PC2 (", pct[2], "%)"))
}
pca <- prcomp(t(mat), center = TRUE, scale. = FALSE)
summary(pca)
plot_pca(pca, "scale. = FALSE")
# スケーリングした場合と比較
plot_pca(prcomp(t(mat), center = TRUE, scale. = TRUE), "scale. = TRUE")
課題2：有意に変化したタンパク質でクラスタリング
library(pheatmap)
# ANOVA で有意なタンパク質（空欄が NA になるので which で囲む）
sig_rows <- which(d$`ANOVA Significant` == "+")
length(sig_rows) # 562
sig_id <- rownames(quant)[sig_rows]
mat_sig <- mat[rownames(mat) %in% sig_id, ]
dim(mat_sig)
ann <- data.frame(Group = group, row.names = colnames(mat_sig))
pheatmap(mat_sig,
scale = "row",
clustering_distance_rows = "correlation",
clustering_distance_cols = "correlation",
clustering_method = "average",
show_rownames = FALSE,
annotation_col = ann,
main = "correlation + average")
# 距離尺度と手法を変えて比較
pheatmap(mat_sig,
scale = "row",
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
show_rownames = FALSE,
annotation_col = ann,
main = "euclidean + ward.D2")
- R6 と R10 で、変動の向きを見る
# Col の平均からの差（log2）
mean_col <- rowMeans(mat_sig[, 1:3])
d_r6 <- rowMeans(mat_sig[, 4:7]) - mean_col
d_r10 <- rowMeans(mat_sig[, 8:9]) - mean_col
plot(d_r6, d_r10, pch = 16, col = "#00000060",
xlab = "R6 - Col (log2)", ylab = "R10 - Col (log2)")
abline(h = 0, v = 0, lty = 2)
課題3：GO解析
AGI code 列には、; で区切られた複数のAGIコードや、デコイ（REV__）が含まれる行があります。
取り出す前に整理します。
library(clusterProfiler)
library(org.At.tair.db)
# AGIコードを取り出す（デコイを除き、複数ある場合はすべて使う）
extract_agi <- function(x) {
x |>
str_remove_all("(REV__|CON__)[^;]*") |>
toupper() |>
str_extract_all("AT[1-5CM]G[0-9]{5}") |>
unlist() |>
unique()
}
# 背景：このデータで定量された全タンパク質
universe_agi <- extract_agi(d$`AGI code`)
length(universe_agi)
sig_genes <- extract_agi(d$`AGI code`[which(d$`ANOVA Significant` == "+")])
length(sig_genes)
ego <- enrichGO(
gene = sig_genes,
universe = universe_agi, # 背景を指定する
OrgDb = org.At.tair.db,
keyType = "TAIR",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
head(as.data.frame(ego))
dotplot(ego, showCategory = 20)
- 背景を指定した場合としない場合の比較
ego_nobg <- enrichGO(
gene = sig_genes,
OrgDb = org.At.tair.db, # universe を指定しない＝全遺伝子が背景
keyType = "TAIR",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
data.frame(
background = c("with_universe", "without_universe"),
n_terms = c(nrow(as.data.frame(ego)), nrow(as.data.frame(ego_nobg)))
)
- カテゴリーを変えて比較
for (o in c("BP", "CC", "MF")) {
e <- enrichGO(gene = sig_genes, universe = universe_agi,
OrgDb = org.At.tair.db, keyType = "TAIR", ont = o,
pAdjustMethod = "BH",
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
cat(o, ":", nrow(as.data.frame(e)), "terms\n")
}