蛋白互作网络分析（STRING/Cytoscape/MCODE/Hub基因筛选）¶

BiocManager::install("STRINGdb")
library(STRINGdb)

# 初始化STRING对象（9606 = 人类）
string_db <- STRINGdb$new(
  version = "12.0",          # R包STRINGdb目前支持到12.0，网页版已更新到12.5
  species = 9606,
  score_threshold = 400,      # medium confidence
  network_type = "full",      # full = 所有证据; physical = 仅物理互作
  input_directory = ""
)

# 准备基因列表（需要有gene symbol列）
deg_list <- data.frame(gene = c("TP53", "BRCA1", "MYC", "EGFR", "CDK2",
                                 "RB1", "PTEN", "AKT1", "MTOR", "PIK3CA"))

# 映射到STRING ID
mapped <- string_db$map(deg_list, "gene", removeUnmappedRows = TRUE)
head(mapped)

# 获取互作网络
interactions <- string_db$get_interactions(mapped$STRING_id)
head(interactions)

# 绘制网络
string_db$plot_network(mapped$STRING_id)

# 富集分析（STRING内置）
enrichment <- string_db$get_enrichment(mapped$STRING_id)
head(enrichment[enrichment$category == "Process", ])

# 导出网络到本地文件
write.table(interactions, "PPI_network.tsv", sep = "\t", quote = FALSE, row.names = FALSE)

# 转换STRING ID回gene symbol用于Cytoscape
id_map <- mapped[, c("STRING_id", "gene")]
interactions_named <- merge(interactions, id_map, by.x = "from", by.y = "STRING_id")
colnames(interactions_named)[ncol(interactions_named)] <- "gene_from"
interactions_named <- merge(interactions_named, id_map, by.x = "to", by.y = "STRING_id")
colnames(interactions_named)[ncol(interactions_named)] <- "gene_to"
write.table(interactions_named[, c("gene_from", "gene_to", "combined_score")],
            "PPI_for_cytoscape.tsv", sep = "\t", quote = FALSE, row.names = FALSE)

1. 安装StringApp插件：
   Apps → App Manager → 搜索"stringApp" → Install

2. 直接从STRING导入（推荐）：
   File → Import → Network from Public Databases
   → 选择STRING → 输入基因列表 → 设置species和confidence

3. 或导入本地文件：
   File → Import → Network from File
   → 选择TSV文件 → 指定Source/Target/Interaction列

4. 调整布局：
   Layout → yFiles Organic Layout（推荐）
   或 Layout → Prefuse Force Directed Layout

5. 视觉映射（Visual Mapping）：
   - 节点大小 → 映射到 Degree
   - 节点颜色 → 映射到 logFC 或 P-value
   - 边粗细 → 映射到 Combined Score

BiocManager::install("RCy3")
library(RCy3)

# 确保Cytoscape已打开并运行
cytoscapePing()  # 测试连接

# 创建网络
nodes <- data.frame(
  id = c("TP53", "BRCA1", "MYC", "EGFR", "CDK2", "RB1", "PTEN"),
  logFC = c(-2.1, -1.5, 3.2, 2.8, 1.9, -1.2, -2.5),
  pvalue = c(0.001, 0.003, 0.0001, 0.002, 0.01, 0.005, 0.0002)
)

edges <- data.frame(
  source = c("TP53", "TP53", "BRCA1", "MYC", "EGFR", "PTEN"),
  target = c("BRCA1", "RB1", "TP53", "CDK2", "MYC", "AKT1"),
  score = c(0.99, 0.95, 0.99, 0.85, 0.78, 0.92)
)

createNetworkFromDataFrames(nodes, edges, title = "PPI Network")

# 设置视觉样式
setNodeSizeMapping("id", mapping.type = "d")
setNodeColorMapping("logFC", c(-3, 0, 3), c("#0000FF", "#FFFFFF", "#FF0000"))

# 设置布局
layoutNetwork("force-directed")

# 导出图片
exportImage("PPI_network.png", type = "PNG", resolution = 300)

Tools → Analyze Network → 选择 Undirected
结果面板会显示：
- 节点数、边数
- 平均degree
- 网络直径
- 聚类系数
- 各种分布图（degree distribution等）

library(igraph)

# 从边列表创建网络
ppi_data <- read.table("PPI_for_cytoscape.tsv", header = TRUE, sep = "\t")
g <- graph_from_data_frame(ppi_data[, c("gene_from", "gene_to")], directed = FALSE)

# 计算拓扑参数
V(g)$degree <- degree(g)
V(g)$betweenness <- betweenness(g)
V(g)$closeness <- closeness(g)
V(g)$eigenvector <- eigen_centrality(g)$vector
V(g)$clustering <- transitivity(g, type = "local")

# 汇总为数据框
topo_stats <- data.frame(
  gene = V(g)$name,
  degree = V(g)$degree,
  betweenness = V(g)$betweenness,
  closeness = V(g)$closeness,
  eigenvector = V(g)$eigenvector,
  clustering = V(g)$clustering
)

# 排序找Hub基因
topo_stats <- topo_stats[order(-topo_stats$degree), ]
head(topo_stats, 10)  # Top 10 hub genes

# 网络全局特征
cat("Nodes:", vcount(g), "\n")
cat("Edges:", ecount(g), "\n")
cat("Average degree:", mean(degree(g)), "\n")
cat("Network density:", graph.density(g), "\n")
cat("Average path length:", mean_distance(g), "\n")
cat("Clustering coefficient:", transitivity(g, type = "global"), "\n")

# Scale-free性质检验
degree_dist <- degree_distribution(g)
# 如果degree分布近似power-law，则网络是scale-free的

1. 安装MCODE：Apps → App Manager → MCODE → Install
2. 运行：Apps → MCODE → 选择网络 → 参数设置：
   - Degree Cutoff: 2（默认）
   - Node Score Cutoff: 0.2（默认）
   - K-Core: 2（默认）
   - Max Depth: 100（默认）
3. 查看结果：每个模块会高亮显示，按score排序
4. 对每个模块做功能富集验证其生物学意义

# 使用igraph的社区检测算法（类似MCODE思想）
library(igraph)

# 方法1：基于边介数的Girvan-Newman算法
communities_gn <- cluster_edge_betweenness(g)
# 方法2：Louvain算法（大网络推荐）
communities_louvain <- cluster_louvain(g)
# 方法3：Fast greedy算法
communities_fg <- cluster_fast_greedy(g)

# 查看模块信息
membership(communities_louvain)
sizes(communities_louvain)
modularity(communities_louvain)

# 提取每个模块的基因
modules <- split(V(g)$name, membership(communities_louvain))
for (i in seq_along(modules)) {
  cat("Module", i, "(", length(modules[[i]]), "genes):", 
      paste(modules[[i]], collapse = ", "), "\n")
}

# 对每个模块做富集分析
library(clusterProfiler)
library(org.Hs.eg.db)

for (i in seq_along(modules)) {
  if (length(modules[[i]]) >= 5) {
    ego <- enrichGO(
      gene = modules[[i]],
      OrgDb = org.Hs.eg.db,
      keyType = "SYMBOL",
      ont = "BP",
      pAdjustMethod = "BH",
      qvalueCutoff = 0.05
    )
    if (nrow(ego) > 0) {
      cat("\n=== Module", i, "Top GO terms ===\n")
      print(head(ego@result[, c("Description", "p.adjust", "Count")], 5))
    }
  }
}

1. 安装：Apps → App Manager → cytoHubba → Install
2. 运行：Apps → cytoHubba → 选择算法 → Top N nodes
3. 推荐策略：
   - 分别用Degree, MCC, Betweenness, Closeness, EPC计算Top 20
   - 取交集作为最终Hub基因
4. 可视化：Hub基因用不同颜色/大小标注

# 用多种中心性指标筛选Hub基因
library(igraph)

# 计算多种中心性
hub_scores <- data.frame(
  gene = V(g)$name,
  degree = degree(g),
  betweenness = betweenness(g, normalized = TRUE),
  closeness = closeness(g, normalized = TRUE),
  eigenvector = eigen_centrality(g)$vector,
  page_rank = page_rank(g)$vector
)

# 每种方法取Top20
top_n <- 20
top_degree <- hub_scores$gene[order(-hub_scores$degree)][1:top_n]
top_between <- hub_scores$gene[order(-hub_scores$betweenness)][1:top_n]
top_close <- hub_scores$gene[order(-hub_scores$closeness)][1:top_n]
top_eigen <- hub_scores$gene[order(-hub_scores$eigenvector)][1:top_n]

# 取交集（至少在3种方法中出现）
all_hubs <- c(top_degree, top_between, top_close, top_eigen)
hub_freq <- table(all_hubs)
final_hubs <- names(hub_freq[hub_freq >= 3])
cat("Final Hub Genes:", paste(final_hubs, collapse = ", "), "\n")

# 绘制Venn图
library(VennDiagram)
venn.diagram(
  x = list(Degree = top_degree, Betweenness = top_between,
           Closeness = top_close, Eigenvector = top_eigen),
  filename = "hub_venn.png",
  fill = c("red", "blue", "green", "yellow"),
  alpha = 0.5
)

# ===== ggraph高级可视化 =====
library(ggraph)
library(tidygraph)

# 转换为tidygraph对象
tg <- as_tbl_graph(g) %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    community = as.factor(group_louvain()),
    is_hub = name %in% final_hubs
  )

# 绘制网络图
ggraph(tg, layout = "fr") +  # Fruchterman-Reingold布局
  geom_edge_link(aes(alpha = 0.3), show.legend = FALSE) +
  geom_node_point(aes(size = degree, color = community)) +
  geom_node_text(aes(label = ifelse(is_hub, name, "")),
                 repel = TRUE, size = 3) +
  scale_size_continuous(range = c(2, 10)) +
  theme_graph() +
  labs(title = "PPI Network with Hub Genes Highlighted")

# ===== 整合表达信息的可视化 =====
# 假设有差异表达结果
deg_data <- data.frame(
  gene = c("TP53", "BRCA1", "MYC", "EGFR"),
  logFC = c(-2.1, -1.5, 3.2, 2.8)
)

tg2 <- tg %>%
  activate(nodes) %>%
  left_join(deg_data, by = c("name" = "gene"))

ggraph(tg2, layout = "fr") +
  geom_edge_link(alpha = 0.2) +
  geom_node_point(aes(size = degree, color = logFC)) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +
  geom_node_text(aes(label = ifelse(degree > 5, name, "")),
                 repel = TRUE, size = 2.5) +
  theme_graph()

# ===== 生存分析验证Hub基因 =====
library(survival)
library(survminer)

# 假设有TCGA临床数据
for (hub in final_hubs) {
  # 按中位数分高低表达组
  median_exp <- median(expr_data[hub, ])
  group <- ifelse(expr_data[hub, ] > median_exp, "High", "Low")

  surv_data <- data.frame(
    time = clinical$OS_time,
    status = clinical$OS_status,
    group = group
  )

  fit <- survfit(Surv(time, status) ~ group, data = surv_data)
  p <- surv_pvalue(fit)$pval

  if (p < 0.05) {
    cat(hub, "is prognostic, P =", p, "\n")
  }
}

# ===== ROC分析验证诊断价值 =====
library(pROC)
for (hub in final_hubs) {
  roc_obj <- roc(group_label, as.numeric(expr_data[hub, ]))
  cat(hub, "AUC =", auc(roc_obj), "\n")
}

# ===== 独立数据集验证 =====
# 在GEO独立数据集中验证hub基因表达模式

# ========================================
# PPI网络分析完整流程
# ========================================

# 1. 安装和加载包
packages <- c("STRINGdb", "igraph", "ggraph", "tidygraph",
              "clusterProfiler", "org.Hs.eg.db", "ggplot2")
BiocManager::install(packages)
lapply(packages, library, character.only = TRUE)

# 2. 准备差异基因列表
deg_genes <- read.csv("DEG_results.csv")
gene_list <- deg_genes$gene[deg_genes$adj.P.Val < 0.05 & abs(deg_genes$logFC) > 1]

# 3. STRING映射
string_db <- STRINGdb$new(version = "12.0",          # R包STRINGdb目前支持到12.0，网页版已更新到12.5 species = 9606, score_threshold = 400)
mapped <- string_db$map(data.frame(gene = gene_list), "gene", removeUnmappedRows = TRUE)
interactions <- string_db$get_interactions(mapped$STRING_id)

# 4. 构建igraph网络
# 转换ID
id_map <- setNames(mapped$gene, mapped$STRING_id)
edges <- data.frame(
  from = id_map[interactions$from],
  to = id_map[interactions$to],
  score = interactions$combined_score
)
edges <- edges[complete.cases(edges), ]
g <- graph_from_data_frame(edges[, 1:2], directed = FALSE)
g <- simplify(g)  # 去除自环和多重边

# 5. 拓扑分析
topo <- data.frame(
  gene = V(g)$name,
  degree = degree(g),
  betweenness = betweenness(g, normalized = TRUE),
  closeness = closeness(g, normalized = TRUE),
  clustering = transitivity(g, type = "local")
)
topo <- topo[order(-topo$degree), ]
write.csv(topo, "network_topology.csv", row.names = FALSE)

# 6. Hub基因（多算法交集）
top20_degree <- topo$gene[1:20]
top20_between <- topo$gene[order(-topo$betweenness)][1:20]
top20_close <- topo$gene[order(-topo$closeness)][1:20]
all_top <- c(top20_degree, top20_between, top20_close)
hub_genes <- names(which(table(all_top) >= 2))
cat("Hub genes:", paste(hub_genes, collapse = ", "), "\n")

# 7. 模块检测
comm <- cluster_louvain(g)
modules <- split(V(g)$name, membership(comm))

# 8. 可视化
tg <- as_tbl_graph(g) %>%
  activate(nodes) %>%
  mutate(degree = centrality_degree(),
         community = as.factor(group_louvain()),
         is_hub = name %in% hub_genes)

p <- ggraph(tg, layout = "fr") +
  geom_edge_link(alpha = 0.15) +
  geom_node_point(aes(size = degree, color = community)) +
  geom_node_text(aes(label = ifelse(is_hub, name, "")), repel = TRUE, size = 3) +
  scale_size(range = c(1, 8)) +
  theme_graph() +
  labs(title = "PPI Network Analysis")
ggsave("PPI_network.pdf", p, width = 12, height = 10)

import networkx as nx
import pandas as pd

# 读取互作数据
edges = pd.read_csv("PPI_edges.tsv", sep="\t")
G = nx.from_pandas_edgelist(edges, "gene_from", "gene_to")

# 计算中心性
degree = nx.degree_centrality(G)
betweenness = nx.betweenness_centrality(G)
closeness = nx.closeness_centrality(G)

# 社区检测
from networkx.algorithms.community import louvain_communities
communities = louvain_communities(G)

# Hub基因
sorted_degree = sorted(degree.items(), key=lambda x: x[1], reverse=True)
top_hubs = [x[0] for x in sorted_degree[:10]]