Добрый день, я выполняю геномное предсказание (GBLUP) у тетраплоидных видов.Я начал использовать зоммер в прошлом году до изменений.Тогда вы могли использовать функцию mmer с аргументом ZETA , чтобы указать ковариацию случайного члена.Теперь есть спецификация с random ~ vs (...)
Я сравнил это с функцией MEMMA , которая использует ZETA аргумент, и, насколько я понимаю, он используется внутри Mmer.Однако по какой-то причине я получаю разные ответы, например
library(devtools)
install_github("covaruber/sommer")
library(sommer)
install.packages("AGHmatrix")
library(AGHmatrix)
data(DT_polyploid)
# ####=========================================####
# ####### convert markers to numeric format
# ####=========================================####
numo <- atcg1234(data=GT, ploidy=4);
# ###=========================================####
# ###### plants with both genotypes and phenotypes
# ###=========================================####
common <- intersect(DT$Name,rownames(numo$M))
#
# ###=========================================####
# ### get the markers and phenotypes for such inds
# ###=========================================####
marks <- numo$M[common,]; marks[1:5,1:5]
DT2 <- DT[match(common,DT$Name),];
# ###=========================================####
G.A <- Gmatrix(marks, method="VanRaden", ploidy=4, missingValue=NA, impute.method = TRUE)
G.D <- Gmatrix(marks, method="Endelman", ploidy=4, missingValue=NA, impute.method = TRUE)
T.pheno <- DT2[,c(1,9)]
T.pheno$Name <- as.factor(T.pheno$Name)
T.pheno$Name2 <- T.pheno$Name
set.seed(1892)
rrn <- sample(1:187, 50, replace = F)
T.pheno$sucrose[rrn] <- NA
ans.A <- mmer(sucrose ~ 1,
random=~vs(Name, Gu=G.A),
rcov = ~units,
data=T.pheno)
cor(ans.A$U$`u:Name`$sucrose[which(is.na(T.pheno$sucrose))],
DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")
[1] 0,2205831
Затем применяю тот же анализ с использованием MEMMA
Z1 <- diag(length(T.pheno$sucrose))
Z2 <- diag(length(T.pheno$sucrose))
ETA.A <- list(list(Z=Z1, K=G.A) )
ETA.AD <- list(list(Z=Z1, K=G.A), list(Z=Z2, K=G.D) )
ans.A <- MEMMA(Y=T.pheno$sucrose, ZETA=ETA.A)
cor(ans.A$fitted.y[which(is.na(T.pheno$sucrose))],
DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")
[1] 0,2778689
Почему есть разница?Для добавления матриц доминирования и эпистазиса я нахожу формулировку в MEMMA более простой и менее подверженной ошибкам, и, конечно, я воодушевлен более высокой точностью прогнозирования.Во-вторых, функция MEMMA обеспечивает установленные значения, которые находятся в том же масштабе и, следовательно, сопоставимы с наблюдаемыми значениями.Однако MEMMA довольно медленная ... и все же должно быть объяснение разницы в точности.Ниже приведен код и результаты для mmer и MEMMA при включении доминирования.
ans.AD <- mmer(sucrose ~ 1,
random= ~ vs(Name, Gu=G.A) + vs(Name2, Gu=G.D),
rcov = ~ units,
data=T.pheno)
cor(ans.AD$U$`u:Name`$sucrose[which(is.na(T.pheno$sucrose))],
DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")
ans.AD <- MEMMA(Y=T.pheno$sucrose, ZETA=ETA.AD)
cor(ans.AD$fitted.y[which(is.na(T.pheno$sucrose))],
DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")
[1] 0.2357571 (mmer с доминированием)
[1] 0.2785493 (MEMMA с доминированием)
Спасибо, что уделили время