Rスクリプト（GO） - yoshida1983/protc2026 GitHub Wiki

PCA・クラスタリングの続きとして、GO（Gene Ontology）のエンリッチメント解析を行います。表計算ソフトは使わず、Rだけで進めます。

目標：次の2つのセットで、サンプルごとにエンリッチした GO term が、どれだけ異なるかを比べます。

• 奇数番号のサンプル（A1, B3, C5, D7）：各班で異なるタグ由来のサンプル
• 偶数番号のサンプル（A2, B4, C6, D8）：同じライン由来のサンプル

進め方

• 同じプロジェクト、同じRの画面で続けます。PCA・クラスタリングで作った merged（GFPを除いたもの、Scoreが0のものは0で置換済み）と、sample_names、sample_no、even_samples を使います
• GO解析では、Score の値は使いません。検出されたかどうか（Score > 0）だけを使います

GO解析の準備

• 1つ前のコマ（Mascotと同じコマ）で、UniProtのブラウザで開いて、Mascotの結果といっしょにダウンロードした対応表を、プロジェクトのフォルダにコピー＆ペースト
  • シロイヌナズナの reviewed エントリの一覧です（ファイル名：uniprotkb_organism_id_3702_AND_reviewed_（ダウンロードした日）.tsv）
• RStudio の【右下】「Files」に、csvファイルといっしょに、表示されていることを確認

対応表の読み込み（Accession → AGIコード）

Mascotの Accession（RCA_ARATH など）を、GO解析に使うAGIコード（AT2G39730 など）に変換します。

• 対応表のファイルを探して、読み込み
  • ファイル名の日付は、ダウンロードした日によって変わるので、名前の最初の部分（uniprotkb_）で探します

uni_file <- list.files(pattern = "^uniprotkb_.*\\.tsv(\\.gz)?$")
uni_file

if (length(uni_file) != 1) stop("対応表のファイルが1つではありません。Filesを確認してください")

uni <- read_tsv(uni_file, show_col_types = FALSE)
names(uni)     # 列名を確認

• 「Entry Name」の列と、AGIコードの列（Gene Names (ordered locus)）から、対応表を作る
  • UniProtのAGIコードは At2g39730 のような形なので、大文字（AT2G39730）にそろえます
  • 1つのエントリに複数のAGIコードがある場合（At1g67370/At1g67380 など）は、すべて使います（別の行になります）

id_col  <- grep("Entry Name", names(uni), value = TRUE)[1]
oln_col <- grep("ordered locus", names(uni), value = TRUE)[1]
if (is.na(id_col) || is.na(oln_col)) stop("「Entry Name」または「Gene Names (ordered locus)」の列がありません")

id2agi <- tibble(
  Accession = uni[id_col](/yoshida1983/protc2026/wiki/id_col),
  AGI       = str_extract_all(toupper(uni[oln_col](/yoshida1983/protc2026/wiki/oln_col)), "AT[1-5CM]G[0-9]{5}")
) |>
  unnest_longer(AGI)

id2agi <- id2agi[!is.na(id2agi$AGI), ]     # AGIコードがないエントリを除く
head(id2agi)

• 検出されたタンパク質に、AGIコードを付ける

go_df <- left_join(merged, id2agi, by = "Accession")

# AGIコードに変換できなかったタンパク質
go_df$Accession[is.na(go_df$AGI)]

go_df <- go_df[!is.na(go_df$AGI), ]
dim(go_df)     # 複数のAGIコードをもつタンパク質は、複数行になります

サンプルごとのタンパク質リスト

• 奇数番号のサンプルを取り出す（偶数番号は、PCA・クラスタリングで作った even_samples）

odd_samples <- sample_names[sample_no %% 2 == 1]
odd_samples
even_samples

• サンプルごとに、検出された（Score > 0）タンパク質のAGIコードのリストを作る

gene_lists <- lapply(sample_names, function(s) unique(go_df$AGI[go_df[s](/yoshida1983/protc2026/wiki/s) > 0]))
names(gene_lists) <- sample_names

sapply(gene_lists, length)     # サンプルごとのタンパク質の数

GO解析（エンリッチメント解析）

• パッケージの読み込み（初回のみインストール）
  • Bioconductorのパッケージは、dplyr の関数（select() など）と名前がぶつかるので、データの準備を先に済ませてから読み込みます

# 初回のみ、先頭の「#」を消して実行（時間がかかります）
# if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
# BiocManager::install(c("clusterProfiler", "org.At.tair.db", "ggtangle"))

library(clusterProfiler)
library(org.At.tair.db)

• library() で「there is no package called '○○'」と出たら、そのパッケージを install.packages("○○") でインストールしてから、もう一度 library() を実行してください（例：ggtangle）

• サンプルごとに、enrichGOを実行

  • 背景（バックグラウンド）は、8サンプルのいずれかで検出された全タンパク質

universe_agi <- unique(go_df$AGI)
length(universe_agi)

ego_list <- lapply(gene_lists, function(g) {
  enrichGO(
    gene          = g,
    universe      = universe_agi,
    OrgDb         = org.At.tair.db,
    keyType       = "TAIR",            # AGIコードを使う
    ont           = "BP",              # BP, CC, MF
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.05
  )
})

# 1つ目のサンプルの結果を確認
head(as.data.frame(ego_list[1](/yoshida1983/protc2026/wiki/1)))

• エンリッチした GO term（GO ID）を、サンプルごとに取り出す

go_terms <- lapply(ego_list, function(e) {
  if (is.null(e)) character(0) else as.data.frame(e)$ID
})

n_terms <- sapply(go_terms, length)
n_terms                                # サンプルごとの、エンリッチした GO term の数

if (sum(n_terms) == 0) stop("エンリッチした GO term がありませんでした")

奇数・偶数のサンプルの比較

• 奇数番号（異なるタグ由来）と偶数番号（同じライン由来）のセットで、GO term の重なりを数える
  • 全サンプル共通：セット内のすべてのサンプルでエンリッチした GO term
  • 1サンプルのみ：セット内の1つのサンプルだけでエンリッチした GO term

odd_terms  <- go_terms[odd_samples]
even_terms <- go_terms[even_samples]

summarise_terms <- function(term_list) {
  all_terms <- unique(unlist(term_list))
  n_samples <- sapply(all_terms, function(t) sum(sapply(term_list, function(x) t %in% x)))
  data.frame(
    "全GO term数"   = length(all_terms),
    "全サンプル共通" = sum(n_samples == length(term_list)),
    "1サンプルのみ"  = sum(n_samples == 1),
    check.names = FALSE
  )
}

rbind(
  "奇数（異なるタグ由来）" = summarise_terms(odd_terms),
  "偶数（同じライン由来）" = summarise_terms(even_terms)
)

• Venn図で比較

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

venn_odd <- venn.diagram(
  x = odd_terms,
  filename = NULL,
  fill = c("skyblue", "pink1", "mediumorchid", "lightgreen"),
  alpha = 0.5,
  cex = 1.2,
  cat.cex = 1.2,
  main = "GO terms: odd-numbered samples"
)

venn_even <- venn.diagram(
  x = even_terms,
  filename = NULL,
  fill = c("skyblue", "pink1", "mediumorchid", "lightgreen"),
  alpha = 0.5,
  cex = 1.2,
  cat.cex = 1.2,
  main = "GO terms: even-numbered samples"
)

# 左右に並べて描画
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))

pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1))
grid.draw(venn_odd)
popViewport()

pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 2))
grid.draw(venn_even)
popViewport()

popViewport()

サンプルごとの AGIコードを書き出す

for (s in names(gene_lists)) {
  writeLines(gene_lists[s](/yoshida1983/protc2026/wiki/s), paste0("shinygo_", s, ".txt"))
}
writeLines(universe_agi, "shinygo_background.txt")   # 背景（Background）用

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

予備実験のcsv（pre1.csv〜pre5.csv）で、GO解析（サンプルごとのエンリッチメント解析と、GO term の比較）までを、まとめて確認するためのコードです。

• PCA・クラスタリングの【動作確認用】のコードを、同じプロジェクトで実行して、merged_pre と names_pre を作ったあとに、実行してください
• 対応表（uniprotkb_...tsv）も、同じプロジェクトのフォルダに置いてください
• pre1〜pre5 の5サンプルを、まとめて比べます（奇数・偶数には分けません）

# ---- 対応表の読み込み（Accession → AGIコード）----
uni_file <- list.files(pattern = "^uniprotkb_.*\\.tsv(\\.gz)?$")
if (length(uni_file) != 1) stop("対応表のファイルが1つではありません。Filesを確認してください")

uni <- read_tsv(uni_file, show_col_types = FALSE)

id_col  <- grep("Entry Name", names(uni), value = TRUE)[1]
oln_col <- grep("ordered locus", names(uni), value = TRUE)[1]
if (is.na(id_col) || is.na(oln_col)) stop("「Entry Name」または「Gene Names (ordered locus)」の列がありません")

id2agi <- tibble(
  Accession = uni[id_col](/yoshida1983/protc2026/wiki/id_col),
  AGI       = str_extract_all(toupper(uni[oln_col](/yoshida1983/protc2026/wiki/oln_col)), "AT[1-5CM]G[0-9]{5}")
) |>
  unnest_longer(AGI)
id2agi <- id2agi[!is.na(id2agi$AGI), ]

go_pre <- left_join(merged_pre, id2agi, by = "Accession")
go_pre$Accession[is.na(go_pre$AGI)]          # AGIコードに変換できなかったタンパク質
go_pre <- go_pre[!is.na(go_pre$AGI), ]
dim(go_pre)

# ---- サンプルごとのタンパク質リスト ----
gene_pre <- lapply(names_pre, function(s) unique(go_pre$AGI[go_pre[s](/yoshida1983/protc2026/wiki/s) > 0]))
names(gene_pre) <- names_pre
sapply(gene_pre, length)                     # サンプルごとのタンパク質の数

# ---- GO解析（エンリッチメント解析）----
library(clusterProfiler)
library(org.At.tair.db)

# 背景：pre1〜pre5 のいずれかで検出された全タンパク質
universe_pre <- unique(go_pre$AGI)
length(universe_pre)

ego_pre <- lapply(gene_pre, function(g) {
  enrichGO(
    gene          = g,
    universe      = universe_pre,
    OrgDb         = org.At.tair.db,
    keyType       = "TAIR",
    ont           = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.05
  )
})

head(as.data.frame(ego_pre[1](/yoshida1983/protc2026/wiki/1)))            # 1つ目のサンプルの結果を確認

terms_pre <- lapply(ego_pre, function(e) {
  if (is.null(e)) character(0) else as.data.frame(e)$ID
})

n_pre <- sapply(terms_pre, length)
n_pre                                        # サンプルごとの、エンリッチした GO term の数

if (sum(n_pre) == 0) stop("エンリッチした GO term がありませんでした")

# ---- GO term の重なり ----
all_terms <- unique(unlist(terms_pre))
n_samples <- sapply(all_terms, function(t) sum(sapply(terms_pre, function(x) t %in% x)))

data.frame(
  "全GO term数"   = length(all_terms),
  "全サンプル共通" = sum(n_samples == length(terms_pre)),
  "1サンプルのみ"  = sum(n_samples == 1),
  check.names = FALSE
)

# Venn図（5サンプル）
invisible(futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger"))

venn_pre_go <- venn.diagram(
  x = terms_pre,
  filename = NULL,
  fill = c("skyblue", "pink1", "mediumorchid", "lightgreen", "orange"),
  alpha = 0.5,
  cex = 1.2,
  cat.cex = 1.2,
  main = "GO terms: pre samples"
)

grid.newpage()
grid.draw(venn_pre_go)