R- Как построить правильные круговые диаграммы в сетях гаплотипов haploNet {pegas} {ape} {adegenet}

При использовании пакета haploNet для создания графиков в сети гаплотипов я использовал для этого скрипт, доступный в Интернете. Однако я думаю, что что-то не так. Скрипт доступен в виде примера Woodmouse. Код, который я использовал:

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, поэтому постарайтесь не быть слишком абстрактным, а вместо этого предоставьте код.

Если есть кто-то, кто очень хорошо знает этот пакет, пожалуйста, объясните также, почему за реальными диаграммами (с цифрами на них) есть эти странные лишние линии, поскольку они не были видны в примере с лесной мышью (может быть, это из-за того, что не так слишком?)

Спасибо заранее


r genetics phylogeny dna-sequence
person Arn De Grauwe    schedule 04.07.2015    source источник

Ответы (1)


arrow_upward
5
arrow_downward

Я боролся с той же проблемой, но считаю, что нашел решение.

Проблема в том, что шаг, делающий 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)

введите здесь описание изображения

И все теперь хорошо!

Дополнительно:

Хотя это не запрашивается ОП, я подумал оставить это здесь, если это представляет интерес для кого-либо еще. Иногда мне удобно маркировать гаплотипы по их частоте. Это можно сделать, изменив метки гаплотипов, чтобы они были равны их частотам:

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

введите здесь описание изображения

person hugot    schedule 13.10.2015
comment
Привет, есть ли способ изменить цвет круговых диаграмм? В функции plot есть параметр col, но он у меня не работает. - person GabrielMontenegro; 08.06.2017
comment
@ JavierM88, вы можете использовать опцию bg. См. справку по функции с ?plot.haploNet - person hugot; 21.06.2017