Rスクリプト（PCA, クラスタリング） - yoshida1983/protc2026 GitHub Wiki

【重要】Mascotの Score について Mascot の Score は、ペプチドの同定の確からしさをもとにした値で、タンパク質の量を表す定量値ではありません。ここでは、Rでの一連の解析フロー（ファイルの結合 → Venn図 → PCA → 階層的クラスタリング）を体験することを目的に、Score を便宜的に強度の代わりに使います。結果の生物学的な解釈には使わないでください。また、Score が 0 のタンパク質は「そのサンプルでは同定されなかった」ことを意味します。

Mascot結果の準備

表計算ソフトは使わず、Rだけで進めます。

RStudio で新しいプロジェクトを、デスクトップに作成
- File → New Project → New Directory → New Project
- Directory name：適宜（半角英数字がおすすめ。例：protc2026）
- Create project as subdirectory of：Desktop
1つ前のコマでダウンロードした csvファイル（8サンプル分：A1, A2, B3, B4, C5, C6, D7, D8） を、作成したプロジェクトのフォルダにコピー＆ペースト
- A〜D はグループ、ファイル名がそのままサンプル名になります
RStudio の【右下】「Files」に、csvファイルが表示されていることを確認
- 表示されない場合は、Files の更新ボタン（丸い矢印）を押す
File → New File → R Script でスクリプトを開き、以下のコードを貼り付けて、上から順に実行（Ctrl + Enter）

プロジェクトを作成すると、そのフォルダが作業ディレクトリになります。setwd() は不要です。

Rの準備

パッケージの読み込み

# パッケージのインストール（初回のみ、先頭の「#」を消して実行）
# install.packages(c("tidyverse", "VennDiagram", "pheatmap"))

library(tidyverse)     # dplyr, readr, purrr, ggplot2 などをまとめて読み込む
library(VennDiagram)
library(pheatmap)

作業ディレクトリの確認（プロジェクトのフォルダになっているはず）

getwd()

csvファイルの読み込みと結合

Mascotのcsvから「Accession」と「Score」だけを取り出す関数を定義（実行するだけでOK）
- csvの先頭には検索条件などのヘッダーがあるので、「Accession」「Score」の見出し行を探して、そこから読み込みます

read_mascot <- function(path) {
  lines <- readLines(path, warn = FALSE, encoding = "UTF-8")

  # 「Accession」「Score」の見出し行を探す（それより前は検索条件などのヘッダー）
  start <- which(grepl("Accession", lines) & grepl("Score", lines))[1]
  if (is.na(start)) stop("見出し行が見つかりません: ", path)

  # 見出し行から次の空行まで（＝タンパク質の表の部分）だけを使う
  block <- lines[start:length(lines)]
  blank <- which(trimws(block) == "")
  if (length(blank) > 0) block <- block[seq_len(blank[1] - 1)]

  df <- read.csv(text = paste(block, collapse = "\n"), check.names = FALSE)

  tibble(Accession = as.character(df$Accession),
         Score     = as.numeric(df$Score)) |>
    group_by(Accession) |>
    summarise(Score = max(Score, na.rm = TRUE))   # 同じAccessionが重複していれば最大値
}

csvファイルの一覧を確認（A1〜D8 の8ファイル分が表示されればOK）

files <- list.files(pattern = "\\.csv$")

# 後で保存する結合結果（merged_proteome.csv）は入力ではないので除く
files <- setdiff(files, "merged_proteome.csv")
files

サンプル名（ファイル名から「.csv」を除いたもの）、グループ（A〜D）、サンプル番号（1〜8）
- A1 (1) のような名前があれば、同じファイルを2回ダウンロードしたものです。余分なファイルを削除してから、もう一度実行してください

sample_names <- tools::file_path_sans_ext(files)
sample_group <- sub("[0-9]+$", "", sample_names)                  # グループ（A1, A2 → A）
sample_no    <- as.numeric(sub("^[A-Za-z]+", "", sample_names))   # サンプル番号（A1 → 1）

# 各グループが2サンプルずつになっているか確認
if (any(table(sample_group) != 2)) stop("各グループ2サンプルになっていません。ファイル名を確認してください")

# 偶数番号のサンプル（同じライン由来）
even_samples <- sample_names[sample_no %% 2 == 0]

data.frame(sample_names, sample_group, sample_no)
even_samples

各ファイルを読み込み、Score列の名前をサンプル名に変更

dfs <- lapply(seq_along(files), function(i) {
  d <- read_mascot(files[i])
  colnames(d) <- c("Accession", sample_names[i])
  d
})

「Accession」列で結合

merged <- reduce(dfs, full_join, by = "Accession")

空欄（NA）を0で置換

merged[is.na(merged)] <- 0

中身を確認

head(merged)
dim(merged)     # タンパク質数 × (1 + サンプル数)

GFP（ベイト）を確認して除く

Mascotの検索は、シロイヌナズナとGFPのデータベースに対して行っています。GFPはベイト（釣り餌）なので、精製がうまくいっていれば、どのサンプルでも高いScoreで検出されているはずです。まず確認してから、以降の解析では除きます。

Accessionの形を確認（先頭20個）
- シロイヌナズナのタンパク質は RCA_ARATH のような形です

head(merged$Accession, 20)

GFPのScoreを確認
- 全サンプルで高い値なら、精製がうまくいった証拠です

is_gfp <- str_detect(merged$Accession, regex("GFP", ignore_case = TRUE))

merged[is_gfp, ]

GFPを除く
- ベイトはどのサンプルにも共通して入っているので、サンプル間の違いを見る解析では除きます

if (sum(is_gfp) == 0) warning("GFPが見つかりません。Accessionの形を確認してください")

merged <- merged[!is_gfp, ]
dim(merged)     # 除いた後のタンパク質数 × (1 + サンプル数)

結合結果の保存（ファイル名はこのまま。変えると、次に実行したときにcsvの一覧に混ざります）

write_csv(merged, "merged_proteome.csv")

Venn図でサンプル間比較

検出されたタンパク質（Score > 0）で、次の2つを比較します。

各班の2サンプル（A1-A2, B3-B4, C5-C6, D7-D8）：同じ班で調製した2サンプルで、どれだけ同じタンパク質が検出されたか
偶数番号のサンプル（A2, B4, C6, D8）：同じライン由来のサンプルで、A〜D班の間でどれだけ再現性があるか

サンプルごとに、検出されたタンパク質のリストを取り出す関数を定義（実行するだけでOK）

detected <- function(s) merged$Accession[merged[[s]] > 0]

1) 各班の2サンプルの比較

円の大きさは、検出されたタンパク質の数に比例します。片方のサンプルで検出されたものが、すべてもう片方でも検出されている場合は、円が入れ子になります

# VennDiagramがログファイルを作らないようにする設定
invisible(futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger"))

groups <- unique(sample_group)     # A B C D

grid.newpage()
pushViewport(viewport(layout = grid.layout(ceiling(length(groups) / 2), 2)))   # 2列に並べる

for (i in seq_along(groups)) {
  cols <- sample_names[sample_group == groups[i]]        # 例：A1, A2

  venn.plot <- venn.diagram(
    x = setNames(lapply(cols, detected), cols),
    filename = NULL,      # ファイルに保存せず、メモリ上に保持
    fill = c("skyblue", "pink1"),
    alpha = 0.5,
    cex = 1.2,
    cat.cex = 1.2,
    main = paste("Group", groups[i])
  )

  pushViewport(viewport(layout.pos.row = (i - 1) %/% 2 + 1,
                        layout.pos.col = (i - 1) %% 2 + 1))
  grid.draw(venn.plot)
  popViewport()
}

popViewport()

2) 偶数番号のサンプルの比較

venn.plot <- venn.diagram(
  x = setNames(lapply(even_samples, detected), even_samples),
  filename = NULL,
  fill = c("skyblue", "pink1", "mediumorchid", "lightgreen"),
  alpha = 0.5,
  cex = 1.2,
  cat.cex = 1.2,
  main = "Even-numbered samples (same line)"
)

grid.newpage()
grid.draw(venn.plot)

【参考】Rにおける色指定

[統計グラフの色（三重大学名誉教授奥村晴彦先生）]

PCA解析

行列に変換 → Log2変換 → 全サンプルで同じ値のタンパク質を除外

# タンパク質 × サンプルの行列（行名 = Accession）
score_mat <- merged |>
  column_to_rownames("Accession") |>
  as.matrix()

# Log2変換（+1 は log(0) を避けるため）
log_mat <- log2(score_mat + 1)

# 全サンプルで同じ値のタンパク質は標準化できないので除く
log_mat <- log_mat[apply(log_mat, 1, sd) > 0, , drop = FALSE]
dim(log_mat)

行／列を入れ換え（行 = サンプル）、分散をそろえてPCAを実行

pca_res <- prcomp(t(log_mat), center = TRUE, scale. = TRUE)
summary(pca_res)

PCA結果の図示

pca_df <- as.data.frame(pca_res$x)
pca_df$Sample <- rownames(pca_res$x)
pca_df$Group  <- sample_group[match(pca_df$Sample, sample_names)]

pct <- round(summary(pca_res)$importance[2, 1:2] * 100, 1)   # PC1, PC2 の寄与率（%）

ggplot(pca_df, aes(x = PC1, y = PC2, color = Group, label = Sample)) +
  geom_point(size = 4) +
  geom_text(vjust = -1.2, size = 4) +
  theme_minimal() +
  labs(
    title = "PCA of Proteome Data (log2 + scaled)",
    x = paste0("PC1 (", pct[1], "%)"),
    y = paste0("PC2 (", pct[2], "%)")
  )

階層的クラスタリング

Venn図と同じ方針で、次の2つを見ます。

全8サンプル：各班の2サンプル（A1-A2, B3-B4, C5-C6, D7-D8）が、近くにまとまるか
- 2サンプルだけではクラスタリングの意味がないので、全サンプルのクラスタリングで、各班の2サンプルが隣り合うかを見ます
偶数番号のサンプル（A2, B4, C6, D8）：同じライン由来のサンプルが、A〜D班の間でどれだけまとまるか

サンプル（列）に付けるグループの情報

ann_col <- data.frame(Group = sample_group[match(colnames(log_mat), sample_names)],
                      row.names = colnames(log_mat))

1) 全8サンプル

階層的クラスタリングの図示

pheatmap(
  log_mat,
  scale = "row",                       # 行方向のz-score（パターン重視）
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "correlation",
  clustering_method = "average",
  show_rownames = FALSE,
  annotation_col = ann_col
)

階層的クラスタリングの図示（2サンプル以上で検出されたタンパク質のみ）

# 検出されたサンプル数（Score > 0）を数える
det_cnt <- rowSums(score_mat[rownames(log_mat), ] > 0)

# 2サンプル以上で検出されたタンパク質だけを残す
mat_2 <- log_mat[det_cnt >= 2, , drop = FALSE]

pheatmap(
  mat_2,
  scale = "row",
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "correlation",
  clustering_method = "average",
  show_rownames = FALSE,
  annotation_col = ann_col,
  annotation_row = data.frame(Detected = det_cnt[rownames(mat_2)])
)

2) 偶数番号のサンプル

# 偶数番号のサンプルだけを取り出す
mat_even <- log_mat[, even_samples, drop = FALSE]

# 取り出したサンプルで値が変わらない（すべて0など）タンパク質は、z-scoreを計算できないので除く
mat_even <- mat_even[apply(mat_even, 1, sd) > 0, , drop = FALSE]
dim(mat_even)

pheatmap(
  mat_even,
  scale = "row",
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "correlation",
  clustering_method = "average",
  show_rownames = FALSE,
  annotation_col = ann_col[even_samples, , drop = FALSE]
)

うまく動かないときは

エラーが出たり、図が表示されなかったりして、原因が分からないときは、次の順で試してください。

ほうきマークで、消す（RStudio の画面には、ほうきのマークが3か所にあります）
- 【右上】Environment のほうきマーク：これまでに作ったデータ（merged など）を、すべて消します。前回の実行の残りが、原因になっていることがあります
- 【右下】Plots のほうきマーク：表示されている図を、すべて消します。図が白紙のときなどに試します
- 【左下】Console のほうきマーク：Console の表示を消します（見やすくするだけで、動作には影響しません）
Rを再起動する
- メニューの Session → Restart R（ショートカット：Ctrl + Shift + F10、Mac は Cmd + Shift + F10）
- 読み込んだパッケージ（library()）も、作ったデータも、すべて消えます。プロジェクトのフォルダにあるファイル（csvなど）は、消えません
スクリプトの上から、順に実行し直す
- 再起動すると、library() の読み込みも消えているので、「Rの準備」から、もう一度実行します

【動作確認用】予備実験サンプル（pre1〜pre5）で行う場合

予備実験のcsv（pre1.csv〜pre5.csv）で、Venn図とクラスタリングまでを、まとめて確認するためのコードです。

「Rの準備」（パッケージの読み込み）と、read_mascot() の定義を実行したあとに、実行してください
A〜D班のcsvと同じフォルダに置くと、上の list.files() に pre のファイルが混ざってしまうので、別のプロジェクトで実行してください
結合結果の保存はしません（保存したcsvが、次の実行で入力に混ざるのを避けるため）

# ---- csvの読み込みと結合 ----
files_pre <- list.files(pattern = "^pre[0-9]+\\.csv$")
files_pre                                     # pre1〜pre5 の5ファイルが表示されればOK
names_pre <- tools::file_path_sans_ext(files_pre)
 
dfs_pre <- lapply(seq_along(files_pre), function(i) {
  d <- read_mascot(files_pre[i])
  colnames(d) <- c("Accession", names_pre[i])
  d
})
 
merged_pre <- reduce(dfs_pre, full_join, by = "Accession")
merged_pre[is.na(merged_pre)] <- 0
 
# ---- GFP（ベイト）を確認して除く ----
is_gfp_pre <- str_detect(merged_pre$Accession, regex("GFP", ignore_case = TRUE))
merged_pre[is_gfp_pre, ]                      # GFPのScore（全サンプルで高いはず）

merged_pre <- merged_pre[!is_gfp_pre, ]
dim(merged_pre)                               # 除いた後のタンパク質数 × (1 + サンプル数)
 
# ---- Venn図（5サンプル）----
invisible(futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger"))
 
venn_pre <- venn.diagram(
  x = setNames(lapply(names_pre, function(s) merged_pre$Accession[merged_pre[[s]] > 0]), names_pre),
  filename = NULL,
  fill = c("skyblue", "pink1", "mediumorchid", "lightgreen", "orange"),
  alpha = 0.5,
  cex = 1.2,
  cat.cex = 1.2,
  main = "Pre-experiment samples"
)
 
grid.newpage()
grid.draw(venn_pre)
 
# ---- 階層的クラスタリング ----
mat_pre <- merged_pre |>
  column_to_rownames("Accession") |>
  as.matrix()
 
log_pre <- log2(mat_pre + 1)
 
# 全サンプルで同じ値のタンパク質は、z-scoreを計算できないので除く
log_pre <- log_pre[apply(log_pre, 1, sd) > 0, , drop = FALSE]
dim(log_pre)
 
pheatmap(
  log_pre,
  scale = "row",
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "correlation",
  clustering_method = "average",
  show_rownames = FALSE
)