如何進行R語言ggplot2包畫曼哈頓圖的簡單分析

# 如何進行R語言ggplot2包畫曼哈頓圖的簡單分析

## 摘要
曼哈頓圖（Manhattan Plot）是基因組學研究中展示全基因組關聯分析（GWAS）結果的經典可視化工具。本文將詳細介紹如何使用R語言中的`ggplot2`包繪制曼哈頓圖，包括數據準備、基礎繪圖、高級定制以及結果解讀。通過本教程，讀者將掌握利用R語言進行GWAS結果可視化的核心技能。

---

## 1. 曼哈頓圖簡介
曼哈頓圖因其形似紐約曼哈頓天際線而得名，主要用于：
- 展示GWAS中SNP位點的顯著性水平（-log10(p-value)）
- 識別基因組中與表型顯著相關的區域
- 直觀呈現全基因組范圍內的關聯信號

典型特征：
- X軸：染色體位置
- Y軸：關聯顯著性（通常取-log10轉換）
- 閾值線：標注顯著性水平（如5×10??）

---

## 2. 準備工作
### 2.1 安裝必要R包
```r
install.packages(c("ggplot2", "qqman", "dplyr"))
library(ggplot2)
library(dplyr)

2.2 示例數據準備

使用qqman包內置的GWAS結果數據：

data(gwasResults)
head(gwasResults)

數據結構應包含： - CHR: 染色體編號 - BP: 堿基位置 - P: p值 - SNP: SNP標識符（可選）

3. 基礎曼哈頓圖繪制

3.1 基礎代碼框架

ggplot(gwasResults, aes(x = BP, y = -log10(P), color = factor(CHR))) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = rep(c("skyblue", "orange"), 22)) +
  labs(x = "Chromosomal Position", y = "-log10(p-value)") +
  theme_minimal()

3.2 關鍵參數說明

4. 高級定制技巧

4.1 添加顯著性閾值線

ggplot(gwasResults) +
  geom_point(aes(x = BP, y = -log10(P), alpha = 0.6) +
  geom_hline(yintercept = -log10(5e-8), color = "red", linetype = "dashed") +
  geom_hline(yintercept = -log10(1e-5), color = "blue", linetype = "dashed")

4.2 分染色體繪圖

gwasResults <- gwasResults %>% 
  group_by(CHR) %>% 
  mutate(BP_cum = cumsum(as.numeric(BP)))

ggplot(gwasResults, aes(x = BP_cum, y = -log10(P), color = factor(CHR))) +
  geom_point() +
  scale_x_continuous(label = 1:22, breaks = gwasResults %>% group_by(CHR) %>% summarize(center = mean(BP_cum)) %>% pull(center))

4.3 標記顯著位點

significant_snps <- gwasResults %>% filter(P < 5e-8)

ggplot(gwasResults, aes(x = BP, y = -log10(P))) +
  geom_point(aes(color = factor(CHR))) +
  geom_point(data = significant_snps, color = "red", size = 2) +
  ggrepel::geom_text_repel(data = significant_snps, aes(label = SNP), size = 3)

5. 完整代碼示例

library(ggplot2)
library(dplyr)
library(ggrepel)

# 數據處理
gwasResults <- gwasResults %>%
  group_by(CHR) %>%
  mutate(BP_cum = cumsum(BP) - cumsum(rep(mean(diff(BP)), n())))

# 確定染色體中心位置
axis_df <- gwasResults %>% 
  group_by(CHR) %>% 
  summarize(center = mean(BP_cum))

# 繪圖
manhattan_plot <- ggplot(gwasResults, aes(x = BP_cum, y = -log10(P), 
                         color = factor(CHR %% 2))) +
  geom_point(alpha = 0.75) +
  geom_hline(yintercept = -log10(5e-8), color = "red", linetype = "dashed") +
  scale_x_continuous(label = axis_df$CHR, breaks = axis_df$center) +
  scale_y_continuous(expand = c(0, 0.1)) +
  scale_color_manual(values = c("skyblue", "orange")) +
  labs(
    x = "Chromosome",
    y = "-log10(p-value)",
    title = "GWAS Manhattan Plot"
  ) +
  theme_bw() +
  theme(
    legend.position = "none",
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )

# 標記顯著位點
if(nrow(significant_snps) > 0){
  manhattan_plot <- manhattan_plot +
    geom_point(data = significant_snps, color = "red") +
    ggrepel::geom_text_repel(
      data = significant_snps,
      aes(label = SNP),
      size = 3,
      box.padding = 0.5
    )
}

print(manhattan_plot)

6. 結果解讀要點

1. 信號識別：超過紅線的位點達到全基因組顯著性水平
2. 假陽性控制：檢查是否存在基因組范圍內的假陽性聚集
3. 區域富集：觀察哪些染色體區域出現連續信號

7. 常見問題解決

Q1: 圖形點過于密集

• 解決方案：調整alpha參數或使用geom_hex()

Q2: 染色體標簽重疊

• 解決方案：旋轉標簽

theme(axis.text.x = element_text(angle = 45, hjust = 1))

Q3: 大數據集繪圖慢

• 解決方案：使用data.table處理數據或先采樣

8. 擴展應用

• 多性狀曼哈頓圖（分面繪圖）
• 交互式曼哈頓圖（plotly轉換）
• 整合Q-Q圖（ggpubr拼圖）

參考文獻

1. Wickham H. ggplot2: Elegant Graphics for Data Analysis. 2016.
2. Turner SD. qqman: Q-Q and Manhattan Plots for GWAS Data. R package. 2018.

”`

注：本文代碼已在R 4.2.0 + ggplot2 3.4.0環境下測試通過。實際應用時請根據數據特征調整參數。