怎样使用Python绘制GWAS分析中的曼哈顿图和QQ图

发布时间：2021-10-11 17:05:41 作者：柒染
来源：亿速云阅读：903

# 怎样使用Python绘制GWAS分析中的曼哈顿图和QQ图

## 摘要
全基因组关联分析（GWAS）是研究遗传变异与表型关联的重要方法，而曼哈顿图和QQ图是展示GWAS结果的两种关键可视化形式。本文详细介绍如何使用Python中的`matplotlib`、`seaborn`和`qqman`等库，从数据预处理到图形定制化绘制高质量的曼哈顿图和QQ图，并提供完整的代码示例和常见问题解决方案。

---

## 1. 引言

### 1.1 GWAS分析简介
全基因组关联分析（Genome-Wide Association Study, GWAS）通过检测全基因组范围内的遗传变异（通常是SNP）与表型之间的关联，来识别潜在的致病基因或位点。典型的GWAS分析会产生每个SNP的p值，用于衡量其与表型关联的显著性。

### 1.2 可视化需求
- **曼哈顿图**：展示全基因组范围内SNP的显著性水平（-log10(p值)），X轴为染色体位置，Y轴为p值的负对数。
- **QQ图**：评估p值分布是否偏离期望，用于检测群体分层或技术偏差。

Python因其丰富的科学计算库（如`pandas`、`numpy`）和可视化库（如`matplotlib`、`seaborn`），成为GWAS可视化的理想工具。

---

## 2. 准备工作

### 2.1 安装必要库
```bash
pip install numpy pandas matplotlib seaborn qqman

2.2 数据格式要求

曼哈顿图和QQ图的输入数据通常为包含以下列的DataFrame： - 曼哈顿图：CHR（染色体编号）、BP（碱基位置）、P（p值） - QQ图：仅需P列

示例数据前5行：

import pandas as pd

data = pd.DataFrame({
    'CHR': [1, 1, 2, 2, 3],
    'BP': [1000, 2000, 1500, 3000, 2500],
    'P': [1e-4, 5e-7, 0.01, 1e-6, 0.1]
})

3. 绘制曼哈顿图

3.1 基础曼哈顿图

使用qqman库快速绘制：

from qqman import manhattan

manhattan(data, chr='CHR', bp='BP', p='P', snp='SNP')
plt.show()

3.2 自定义曼哈顿图（matplotlib实现）

import matplotlib.pyplot as plt
import numpy as np

# 预处理：按染色体分组并计算累计位置
data['-logP'] = -np.log10(data['P'])
data = data.sort_values(['CHR', 'BP'])
data['pos'] = data.groupby('CHR')['BP'].transform(lambda x: x - x.min())
chromosomes = data['CHR'].unique()
colors = ['#1f77b4', '#ff7f0e']  # 交替颜色

fig, ax = plt.subplots(figsize=(12, 6))
for i, chr in enumerate(chromosomes):
    chr_data = data[data['CHR'] == chr]
    ax.scatter(
        x=chr_data['pos'] + i * 1e7,  # 染色体间留空
        y=chr_data['-logP'],
        color=colors[i % 2],
        s=10
    )

# 添加显著线（通常为5e-8）
ax.axhline(y=-np.log10(5e-8), color='r', linestyle='--')
ax.set_xticks([i * 1e7 for i in range(len(chromosomes))])
ax.set_xticklabels(chromosomes)
ax.set_xlabel('Chromosome')
ax.set_ylabel('-log10(p-value)')
plt.show()

3.3 高级定制

分染色体颜色：使用seaborn的hue参数
标注显著位点：

sig_snps = data[data['P'] < 5e-8]
for _, row in sig_snps.iterrows():
    ax.text(row['pos'], row['-logP'], row['SNP'], fontsize=8)

4. 绘制QQ图

4.1 基础QQ图

import statsmodels.api as sm

sm.qqplot(data['P'], line='45')
plt.title('QQ Plot of GWAS p-values')
plt.show()

4.2 自定义QQ图

observed = -np.log10(np.sort(data['P']))
expected = -np.log10(np.linspace(1/len(data), 1, len(data)))

plt.figure(figsize=(8, 8))
plt.scatter(expected, observed, s=5)
plt.plot([0, max(expected)], [0, max(expected)], 'r--')
plt.xlabel('Expected (-log10(p))')
plt.ylabel('Observed (-log10(p))')
plt.title('QQ Plot with λGC: {:.3f}'.format(np.median(observed)/np.median(expected)))

4.3 解释QQ图

理想情况：点沿对角线分布
膨胀（λGC > 1）：可能存在群体分层
截尾偏离：提示真实关联信号

5. 代码优化与常见问题

5.1 大数据处理技巧

使用numpy的memmap处理超大数据
分染色体绘制后合并：

fig, axes = plt.subplots(3, 4, figsize=(16, 12))  # 12染色体分块绘制

5.2 常见错误解决

p值为0或1：过滤或替换为极小值


data['P'] = data['P'].clip(lower=1e-300)

染色体顺序错乱：转换为分类变量


data['CHR'] = pd.Categorical(data['CHR'], categories=list(range(1, 23)) + ['X', 'Y'])

6. 完整代码示例

[此处可放置一个从数据加载到生成两图的完整脚本]

7. 总结

本文演示了利用Python绘制GWAS曼哈顿图和QQ图的完整流程，涵盖： - 使用qqman快速绘图 - 基于matplotlib的自定义实现 - 图形美学调整与生物学解释 - 大数据场景下的优化方案

通过灵活调整代码参数，研究者可以生成符合出版物要求的高质量图像。

参考文献

Turner, S.D. (2014). qqman: Q-Q and manhattan plots for GWAS data. R package.
Purcell et al. (2007). PLINK: A toolset for whole-genome association analysis. AJHG.

”`

注：实际文章需补充完整代码示例部分并调整细节至3850字左右。此框架已包含技术要点和逻辑结构，可根据需要扩展以下内容： - 不同GWAS软件（PLINK、GCTA）的输出格式处理 - 更多自定义选项（如3D曼哈顿图） - 与其他工具（如R的ggplot2）的对比 - 实际生物学案例分析