GWAS解决时间批次效应

1. 识别批次效应

通过主成分分析（PCA）识别批次效应

原理：将高维基因型数据降维到低维空间（如前 2-3 个主成分），若样本在散点图中按批次聚集（同一批次样本紧密聚类，不同批次明显分离），则提示存在批次效应。
方法：基于基因型数据进行 PCA（用 PLINK 的--pca），提取前 10-20 个主成分（PCs）；
分析：通过散点图（PC1 vs PC2）观察不同时间段的样本是否在主成分上重叠（理想状态），或显著分离（提示种群分层）；
统计检验：用方差分析（ANOVA）检验 “时间标签” 是否与前 5 个主成分显著相关（p<0.05），若相关则需校正。
示例：用 PLINK 的--pca命令生成主成分，再用 R 的ggplot2绘制散点图，以批次为颜色分组，观察聚类模式。

2. 不同批次数据合并

2.1 不同批次进行基因组数据质量控制（QC)

使用plink软件对不同批次的数据在合并前进行QC

2.2 合并文件

检查重复样本：

1 2	`# 生成样本ID统计 plink --bfile batch1 --list-duplicate-vars ids-only suppress-first > duplicates.txt`

创建合并文件列表：

1
2
3

# 创建merge_list.txt，包含所有待合并文件的前缀
echo "batch1_qc" > merge_list.txt
echo "batch2_qc" >> merge_list.txt

合并文件：

# 尝试合并（自动检测并跳过冲突SNP）
plink --merge-list merge_list.txt --make-bed --out merged_data

# 如果存在冲突，生成报告
plink --bfile batch1_qc --bmerge batch2_qc.bed batch2_qc.bim batch2_qc.fam --recode --out temp

处理冲突：

# 手动移除冲突SNP后重新合并
plink --bfile batch1_qc --exclude temp-merge.missnp --make-bed --out batch1_clean
plink --bfile batch2_qc --exclude temp-merge.missnp --make-bed --out batch2_clean
plink --bfile batch1_clean --bmerge batch2_clean.bed batch2_clean.bim batch2_clean.fam --make-bed --out final_merged

2.3 关于合并冲突

在使用 PLINK 合并不同批次基因数据时，冲突是常见问题，主要源于 SNP 或样本信息的不一致。以下详细讲解冲突的产生原因、检测方法及解决流程：

一、冲突的类型及产生原因

合并冲突通常分为两类，需针对性处理：

1. SNP 冲突

定义：不同批次文件中存在同名 SNP（rsID 相同），但核心信息不一致。
常见原因：
- 等位基因分型错误（如 A/T vs T/A 颠倒，可能因检测平台链方向不同）。
- 参考基因组版本差异（如 hg19 vs hg38 导致 SNP 位置偏移）。
- 探针设计误差（同一 rsID 对应不同位点）。

2. 样本冲突

定义：不同批次文件中存在相同样本 ID，但表型或基因型数据不一致。
常见原因：
- 样本重复录入（同一受试者在不同批次被测序）。
- 样本 ID 命名错误（如拼写重复）。

二、冲突的检测方法

使用 PLINK 的–bmerge命令合并时，会自动检测冲突并生成报告文件，步骤如下：

1. 尝试合并并生成冲突报告

# 以两个批次为例（batch1_qc和batch2_qc为质控后的文件）
plink --bfile batch1_qc \
      --bmerge batch2_qc.bed batch2_qc.bim batch2_qc.fam \
      --recode \
      --out temp_merge

关键输出文件：
- temp_merge-merge.missnp：记录所有SNP 冲突（必看文件）。
- temp_merge-merge.missfam：记录样本 ID 冲突（若存在重复样本）。
- 日志信息（终端输出）：提示冲突数量及类型（如 “100 SNPs with non-matching alleles”）。

2. 解读冲突报告

temp_merge-merge.missnp 内容示例：

1 2	`rs12345 # 冲突SNP的rsID rs67890`

需进一步检查这些 SNP 在两个批次中的具体差异（通过查看.bim文件）：

# 查看batch1中冲突SNP的信息
grep "rs12345" batch1_qc.bim
# 输出格式：染色体  rsID  遗传距离  物理位置  等位基因1  等位基因2

# 对比batch2中的信息
grep "rs12345" batch2_qc.bim

若差异为等位基因颠倒（如 batch1 为 A/T，batch2 为 T/A），属于链方向问题（可校正）
若差异为非互补等位基因（如 batch1 为 A/T，batch2 为 C/G），可能是分型错误（需排除）
temp_merge-merge.missfam 内容示例：

1	`sample123 # 重复的样本ID`

说明同一受试者出现在多个批次中，需确认是否为重复测序（保留一个）或 ID 录入错误（修正 ID）

三、解决冲突的具体操作

根据冲突类型选择以下方法，核心原则是保证合并后数据的一致性：

1. 解决 SNP 冲突

（1）排除冲突 SNP（最简单直接）

若冲突 SNP 数量少或无法校正，直接从所有批次中排除：

# 从batch1中排除冲突SNP
plink --bfile batch1_qc \
      --exclude temp_merge-merge.missnp \  # 导入冲突SNP列表
      --make-bed \
      --out batch1_clean

# 从batch2中排除相同的冲突SNP
plink --bfile batch2_qc \
      --exclude temp_merge-merge.missnp \
      --make-bed \
      --out batch2_clean

# 重新合并（此时无SNP冲突）
plink --bfile batch1_clean \
      --bmerge batch2_clean.bed batch2_clean.bim batch2_clean.fam \
      --make-bed \
      --out merged_no_conflict

（2）校正链方向冲突（保留 SNP，推荐）

若冲突源于等位基因链方向颠倒（如 A/T vs T/A），可通过–flip命令校正：

# 生成需翻转的SNP列表（即冲突SNP）
cp temp_merge-merge.missnp flip_list.txt

# 翻转batch2中冲突SNP的等位基因
plink --bfile batch2_qc \
      --flip flip_list.txt \  # 翻转指定SNP的等位基因
      --make-bed \
      --out batch2_flipped
      
# 使用批次1作为参考，翻转批次2中的等位基因（该条代码有待考证）
plink --bfile batch2_qc \
      --reference-allele batch1_qc.bim 4 2 \  # 关键参数：统一等位基因
      --extract common_snps.txt \
      --make-bed \
      --out batch2_aligned

# 重新合并（此时SNP等位基因一致）
plink --bfile batch1_qc \
      --bmerge batch2_flipped.bed batch2_flipped.bim batch2_flipped.fam \
      --make-bed \
      --out merged_flipped

注意：翻转后需再次检查等位基因是否一致，避免错误翻转。

（3）基于参考基因组校正（精准但复杂）

若因参考基因组版本不同（如 hg19 vs hg38）导致 SNP 位置或等位基因差异：

使用liftOver工具将 SNP 坐标转换为统一版本。
比对参考基因组的等位基因（如从 dbSNP 下载标准信息），修正批次中的错误分型。

2. 解决样本冲突

（1）移除重复样本

保留一个批次的重复样本（优先保留表型完整或测序质量高的批次）：

# 创建重复样本列表（从missfam文件提取）
cut -f1 temp_merge-merge.missfam > duplicate_samples.txt

# 从batch2中移除重复样本（假设保留batch1的样本）
plink --bfile batch2_qc \
      --remove duplicate_samples.txt \  # 移除指定样本
      --make-bed \
      --out batch2_no_dups

# 重新合并
plink --bfile batch1_qc \
      --bmerge batch2_no_dups.bed batch2_no_dups.bim batch2_no_dups.fam \
      --make-bed \
      --out merged_no_dups

（2）修正样本 ID（若为命名错误）

若样本 ID 重复是因命名规则不同（如 batch1 为 “ID123”，batch2 为 “123”），可批量修改 ID：

# 创建ID映射文件（old_id new_id）
echo -e "123\tID123" > id_map.txt  # 示例：将batch2的"123"改为"ID123"

# 替换batch2中的样本ID
plink --bfile batch2_qc \
      --update-ids id_map.txt \  # 更新样本ID
      --make-bed \
      --out batch2_renamed

# 重新合并
plink --bfile batch1_qc \
      --bmerge batch2_renamed.bed batch2_renamed.bim batch2_renamed.fam \
      --make-bed \
      --out merged_renamed

四、合并后验证

冲突解决后，需验证合并数据的完整性和一致性：

检查样本和 SNP 数量：

1
2
3

# 查看合并后的数据量
wc -l merged_final.fam  # 样本数（行数）
wc -l merged_final.bim  # SNP数（行数）

确保数量符合预期（总样本数 = 各批次样本数 - 重复样本数；总 SNP 数 = 共同 SNP 数 - 冲突排除数）。

随机抽查 SNP 一致性：

# 随机选择10个SNP，检查在合并前后的等位基因是否一致
shuf merged_final.bim | head -n 10 | cut -f2 > random_snps.txt
grep -f random_snps.txt batch1_qc.bim
grep -f random_snps.txt batch2_clean.bim  # 或处理后的batch2文件

五、关键参数总结

命令 / 参数	作用	适用场景
–bmerge	合并二进制文件并检测冲突	初始合并步骤
–exclude	排除指定 SNP	解决 SNP 冲突
–flip	翻转指定 SNP 的等位基因	校正链方向导致的 SNP 冲突
–remove	移除指定样本	解决样本重复冲突
–update-ids	批量修改样本 ID	修正样本 ID 命名错误
–make-bed	输出二进制文件（.bed/.bim/.fam）	所有步骤的输出格式指定

通过以上步骤，可系统解决 PLINK 合并中的冲突问题。实际操作中，建议优先通过校正（如–flip）保留数据，减少信息损失；仅在无法校正时才选择排除。若冲突数量过大，需回溯检查原始数据的质控流程（如芯片型号、参考基因组版本是否一致）。

3. 校正批次效应

3.1 将批次作为分类协变量

原理：将 “时间批次” 作为分类协变量直接纳入 GWAS 的回归模型中，通过统计模型控制批次对表型的独立影响，从而分离 SNP 与表型的真实关联。
适用场景：时间批次划分明确（如 “2020 年组”“2021 年组”“2022 年组”），且批次效应可通过线性关系校正。

方法：

先进行QC，再将批次作为分类协变量，同时将PCs作为连续型协变量加入模型，并使用GWAS中的gemma软件，模型使用MLM

(1) 使用plink进行QC

(2) 将批次作为分类协变量添加到协变量文件中

将分类协变量转换为哑变量

1	`plink - file b - covar cov1.txt - write-covar - dummy-coding`

将哑变量添加到协变量文件中

(3) 利用主成分分析 (PCA) 校正

基于所有样本的基因型数据计算主成分 (PCs)。
选择与批次效应显著相关的前K个PCs（通常通过可视化或统计检验确定K）。
在关联分析模型中，将这些选出的PCs作为连续型协变量加入模型。

(4) 计算遗传关系矩阵

1	`gemma -bfile your_data -gk 1 -o relatedness_matrix`

(3) 运行MLM模型（带批次协变量）

gemma -bfile [prefix] \
      -k [grm_file] \             # 上一步生成的关系矩阵文件
      -c [covariates.txt] \       # 包含批次信息的协变量文件
      -p [phenotype.txt] \        # 表型文件
      -lmm 4 \                    # 使用高效的 LMM 算法（LMM-LOGISTIC）
      -o [output_prefix]

3.2 通过R语言的Combat算法

在 GWAS 中，可使用 R 软件中的相关包，如 sva、limma 等，结合 Combat 算法、混合线性模型等来去除不同时间批次样本产生的误差
ComBat是一种基于经验贝叶斯的批次校正方法，直接对基因型数据进行 “批次标准化”（而非仅在模型中加协变量），消除不同时间批次间的系统性差异（如不同年份样本的 SNP 检出率差异），校正后的数据可直接用于常规 GWAS 分析。
适用场景:
- 时间批次间存在显著的系统性偏差（如 2020 年样本的平均测序深度显著低于 2022 年）；
- 希望直接修正数据而非通过模型控制

4. 校正后验证

检查重复样本：批次间重复样本基因型一致性应接近批次内水平。

QQ图/Lambda值

1
2
3

# 计算校正后λ值
awk 'NR>1 {print $15}' gwas_corrected.assoc.txt | grep -v "NA" > pvals.txt
Rscript -e "median(qchisq(scan('pvals.txt'), df=1, lower.tail=F)/0.456"

λ ≈ 1 表示批次效应已控制（理想值1.0-1.05）。

PCA重新可视化：校正后PCA中批次聚类现象应消失。

Bioinformatics

#Bioinformatics #GWAS

GWAS解决时间批次效应

https://oldstory.cn/2025/08/13/gwas_jie_jue_shi_jian_pi_ci_xiao_ying/

作者

Ricardo

发布于

2025年8月13日

许可协议

GWAS流程上一篇

Gemma使用下一篇