3

当使用 haploNet 包在单倍型网络上绘制一些图时，我使用了 Internet 上可用的脚本来执行此操作。不过我觉得有问题。该脚本以木鼠示例的形式提供。我使用的代码是：

x <- read.dna(file="Masto.fasta",format="fasta")
h <- haplotype(x)
net <- haploNet(h)
plot(net)

plot(net, size = attr(net, "freq"), fast = TRUE)
plot(net, size = attr(net, "freq"))
plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.8

table(rownames(x))

ind.hap<-with(
    stack(setNames(attr(h, "index"), rownames(h))), 
    table(hap=ind, pop=rownames(x)[values])
)
ind.hap 

plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.8, pie=ind.hap)
legend(50,50, colnames(ind.hap), col=rainbow(ncol(ind.hap)), pch=20)

legend(x=7,y=10,c("Baeti ero","Felege weyni","Golgole naele","Hagare selam","Ruba feleg","Ziway"),c("red","yellow","green","turquoise","blue","magenta"))

但是，在绘制 ind.hap 时，您会注意到某些行不在正确的位置。你可以在这里看到这个：

      pop
hap    Baetiero ETH022 ETH742 Felegeweyni Golgolenaele Rubafeleg
  I           0      0      1           0            0         0
  II          0      1      0           0            0         0
  III         1      0      0           1            0         1
  IV          2      0      0           0            0         3
  IX          0      0      0           1            0         0
  V           4      0      0           0            2         0
  VI          4      0      0           1            0         4
  VII         2      0      0           1            0         0
  VIII        0      0      0           1            0         1
  X           3      0      0           0            1         0
  XI          0      0      0           0            1         1
  XII         0      0      0           1            0         0
  XIII        0      0      0           0            0         1

您可以看到第 IX 行不在正确的位置。这不会有太大问题，但程序需要第 9 行来绘制 IX 的饼图，即 VIII 的数据。结果是这样的：（我无法插入图像，因为我的声誉低于 10...，无论如何您都可以通过执行整个文件来获得图像）

您可以看到，对于 V 直到 IX，它并不是应有的状态（这些是交换的行）。例如：IX 中只有 1 个单倍型，但有 2 个单倍型的饼图（两者都占图表的 50%），它是使用 VIII 数据生成的。由于行是按字母顺序而不是升序排序的，但这是包固有的，我不知道该怎么做。我远不是 R 的大师，所以尽量不要太抽象，而是提供代码。

如果有人非常了解这个包，请解释为什么在真实图表后面有这些奇怪的额外线条（上面有数字），因为它们在木鼠示例中不可见（可能是因为出了什么问题也？）

提前感谢

4

1 回答 1

5

我一直在努力解决同样的问题，但相信我想出了一个解决方案。

问题在于，table对每个“种群”进行单倍型计数的步骤按字母顺序排列单倍型。因此，例如，单倍型“IX”出现在“V”之前。另一方面，该函数haplotype()按其“数字”顺序对单倍型进行排序。这就是在绘图时造成差异的原因。

这可以通过按“标签”对单倍型对象进行排序来解决，如?haplotype帮助中所述。

我将使用woodmouse示例数据来举例说明：

# Sample 9 distinct haplotypes
library(pegas)
data(woodmouse)
x <- woodmouse[sample(9, 100, replace = T), ]

为了简化，我创建了一个函数来创建单倍型计数表（基于这篇文章）：

countHap <- function(hap = h, dna = x){
    with(
        stack(setNames(attr(hap, "index"), rownames(hap))),
        table(hap = ind, pop = attr(dna, "dimnames")[[1]][values])
    )
}

现在，让我们看看没有对单倍型进行排序的结果：

h <- haplotype(x) # create haplotype object
net <- haploNet(h) # create haploNet object

plot(net, pie = countHap(), size = attr(net, "freq"), legend = T)

在此处输入图像描述

现在，让我们看看我们的计数表，以检查这些结果：

countHap(h, x)

      pop
hap    No0906S No0908S No0909S No0910S No0912S No0913S No304 No305 No306
  I          0       0       0       0       0       0     0     8     0
  II         0       0       0       0       0       0     9     0     0
  III        0       0       0       0       0       0     0     0    10
  IV        16       0       0       0       0       0     0     0     0
  IX         0       0       0       0       0       8     0     0     0
  V          0      12       0       0       0       0     0     0     0
  VI         0       0      10       0       0       0     0     0     0
  VII        0       0       0      13       0       0     0     0     0
  VIII       0       0       0       0      14       0     0     0     0

事情不匹配：例如，单倍型“V”应该出现在个体“No0908S”中，而是被着色为个体“No0913S”（这应该是单倍型“IX”的标签）。

现在，让我们对单倍型进行排序：

h <- haplotype(x)
h <- sort(h, what = "labels") # This is the extra step!!
net <- haploNet(h)

plot(net, pie = countHap(), size = attr(net, "freq"), legend = T)

在此处输入图像描述

现在一切都很好！

额外：

尽管 OP 没有要求这样做，但如果其他人感兴趣，我想把它留在这里。有时，我发现按频率标记单倍型很方便。这可以通过将单倍型标签更改为等于它们的频率来完成：

attr(h, "labels") <- attr(h, "freq")
plot(net, pie = countHap(), size = attr(net, "freq"), legend = T)

在此处输入图像描述

于 2015-10-13T15:55:04.090 回答