0

我正在 R 中创建曼哈顿图，我的目标是根据 SNP 相关的特征突出显示 SNP。

这是我的主要数据框（gwas_plot）的尾部：

              SNP CHR        BP         P 
58008  rs77855001   1 249213034 8.665e-01
58009  rs74322942   1 249218520 7.437e-01
58010 rs114152373   1 249222547 6.504e-02
58011  rs11205304   1 149906423 4.000e-23
58012  rs11205305   1 149906433 2.000e-16
58013   rs1061956   1 149914327 4.000e-06

我有另一个数据集（gwascat_topsnp），其中包含 gwas_plot 中的 SNP 子集，并包含描述。它看起来像这样：

         SNP CHR        BP     P            Description
1 rs11205301   1 149906412 4e-23                 Height
2 rs11205301   1 149906412 2e-16                 Height
3  rs1061954   1 149914353 4e-06 Obesity-related Traits

我在允许您突出显示特定 SNP 的特定部分更改了曼哈顿图的源代码（原始源代码在这里：https ://github.com/stephenturner/qqman/tree/master/R）。我对其进行了更改，以便 pch 对于每个“描述”类别（在上述数据框中）都是唯一的。

这是我添加的部分：

  keyword2=gwascat_topsnp$Description
# Highlight2: highlight snps from a second character vector
  if (!is.null(highlight2)) {
    if (any(!(highlight2 %in% d$SNP))) warning("You're trying to highlight SNPs that don't exist in your results.")
    d.highlight2=d[which(d$SNP %in% highlight2), ]
    with(d.highlight2, points(pos, logp, col="red1", pch=1:24[as.factor(keyword2)], ...)) 
  }

所以这是我用来绘制曼哈顿图的代码：

par(xpd=T, mar=par()$mar+c(0,0,3,0))
manhattan_KB(gwas_plot,col=c("grey"),highlight=gwas_main[rownum,1],highlight2=gwascat_topsnp$SNP)
    par(xpd=TRUE)
    legend(0,40,legend=levels(as.factor(keyword2)),
           col=c("red1"),pch=1:24[as.factor(keyword2)],
           bty="n",cex=0.8,ncol=1,horiz=F)
    par(mar=c(5,4,4,2) + 0.1)

然而，在这里看到的输出中，描述为“高度”的两个 SNP 应该具有相同的 pch，但它们没有：其中一个有一个圆圈，一个有一个三角形，另一个特征 remians 在传奇。

有人可以帮我解决我哪里出错了吗？如果您需要更多信息，我将很乐意提供。

谢谢！

4

0 回答 0