基于Limma交互模型分城市污染物CpG位点关联分析方法咨询

你用的log2Pb*City交互模型完全可行，既能识别不同城市中Pb对CpG甲基化影响存在差异的位点，也能单独分析每个城市内Pb与CpG的关联。不过当前代码在结果提取环节有问题，以下是具体解释和调整建议：

模型逻辑说明 #

你的模型公式~log2Pb*City + sex + Age + CD8T + CD4T + NK + Bcell + Mono涵盖了关键要素：

• log2Pb主效应：所有城市整体上Pb暴露对甲基化的平均影响
• City主效应：不同城市间甲基化水平的基线差异
• log2Pb:City交互效应：不同城市中Pb影响甲基化的程度差异
• 人口学（性别、年龄）和细胞比例协变量：有效控制混杂因素，这部分设置很规范

代码调整方案 #

你当前用coef=2提取的是模型中第二个变量（大概率是log2Pb的主效应），但要分析每个城市单独的Pb关联，需要针对对应城市的交互项提取结果：

# 构建模型矩阵
var <- model.matrix(~log2Pb*City + sex + Age + CD8T + CD4T + NK + Bcell + Mono, data=targets2)

# 查看所有变量名，确认每个城市对应的交互项位置
colnames(var)

# 拟合模型
fit <- lmFit(mval, var)
fit2 <- eBayes(fit, trend=TRUE, robust=TRUE)

# 示例：提取各城市的Pb关联结果（需替换为实际交互项名称）
# 参考城市的Pb效应直接取log2Pb主效应
probe_ref_city <- topTable(fit2, adjust="BH", coef="log2Pb", num=Inf)
# 提取City2中Pb的关联结果（交互项名称需匹配colnames输出）
probe_city2 <- topTable(fit2, adjust="BH", coef="log2Pb:CityCity2", num=Inf)
# 提取City3中Pb的关联结果
probe_city3 <- topTable(fit2, adjust="BH", coef="log2Pb:CityCity3", num=Inf)
# 提取City4中Pb的关联结果
probe_city4 <- topTable(fit2, adjust="BH", coef="log2Pb:CityCity4", num=Inf)

关键注意事项 #

• 必须先通过colnames(var)确认交互项的准确名称：因为City作为因子变量时，模型矩阵会自动生成对比项（比如城市名为Beijing、Shanghai，会生成CityShanghai这类变量，交互项就是log2Pb:CityShanghai）
• 若要同时检验所有城市的Pb效应差异，可通过正则匹配提取所有交互项：topTable(fit2, coef=grep("log2Pb:City", colnames(var)), num=Inf)
• 4个城市的交互项会有3个（参考城市是因子的第一个水平），其他城市的绝对Pb效应需要将log2Pb主效应与对应交互项系数相加，才能得到该城市中每单位log2Pb变化对应的甲基化改变量

内容的提问来源于stack exchange，提问作者Yogesh Gupta