Ошибка в if (min (gf)> = 0) {: пропущенное значение, где TRUE / FALSE требуется при оценке максимального правдоподобия кода функции частоты гаплотипа - PullRequest
Новая
       25

Ошибка в if (min (gf)> = 0) {: пропущенное значение, где TRUE / FALSE требуется при оценке максимального правдоподобия кода функции частоты гаплотипа

0 голосов
/ 08 ноября 2019

В 64-й строке кода выдается ошибка, и я понятия не имею, как избавиться от этого (и сделать дальнейший код исполняемым):

Ошибка в if (min (gf)> = 0) {: пропущенное значение там, где необходимо ИСТИНА / ЛОЖЬ

Код оценивает оценку максимального правдоподобия частоты гаплотипа. Это работает правильно, пока не попадет в проблемный цикл, тогда выполнение кода резко останавливается. В общем, я должен рассчитать логарифмическую вероятность.

Проблемный код: 

PA <- 0.7
PB <- 0.6
d<- 0.1
r <- d /sqrt(PA*Pa*PB*Pb)
make.hap.freq<-function(PA,PB,r){
    Pa <- 1-PA
    Pb <- 1-PB
    ret <- c(PA*PB,PA*Pb,Pa*PB,Pa*Pb)
    ret <- ret + c(d,-d,-d,d)
    ret
}

make.gen.freq <- function(PA,PB,r){
    h <- make.hap.freq(PA,PB,r)
    g <- matrix(0,3,3)
    g[1,1] <- h[1]^2
    g[1,3] <- h[2]^2
    g[3,1] <- h[3]^2
    g[3,3] <- h[4]^2
    g[1,2] <- 2*h[1] * h[2]
    g[2,1] <- 2*h[1] * h[3]
    g[2,3] <- 2*h[2] * h[4]
    g[3,2] <- 2*h[3] * h[4]
    g[2,2] <- 2*(h[1] * h[4] + h[2] * h[3])
    g
}


make.hap.freq(PA,PB,r)

make.gen.freq(PA,PB,r)


#########
## 7 Diplotype data simulation for 2 SNVs
#########

gf <- make.gen.freq(PA,PB,r)
N <- 1000
G <- sample(1:9, N, replace=TRUE,prob =gf)
tabulate(G)
matrix(tabulate(G), 3,3)

#########
## 8 Maximum likelihood estimation of haplotype frequency
#########
# Assume you have 3x3 diplotype counts data only.
# Based on the table, you can estimate haplotype frequency.


Gmat= matrix(c(20,24,5,14,22,6,2,2,5))
Gmat.rsum <- apply(Gmat,1,sum)
Gmat.csum <- apply(Gmat,2,sum)

PA <- (2*Gmat.rsum[1]+Gmat.rsum[2])/(2*sum(Gmat.rsum))
PB <- (2*Gmat.csum[1]+Gmat.csum[2])/(2*sum(Gmat.csum))


rs <- seq(from=-1,to=1,length=100)
logLike <- rep(NA,length(rs))
for(i in 1:length(rs)){
    tmp.r <- rs[i]
    gf <- make.gen.freq(PA,PB,tmp.r)


############
######This is the problematic loop:


    if(min(gf)>=0){
        logLike[i] <- dmultinom(c(Gmat),prob=c(gf),log=TRUE)
    }
}
plot(rs,logLike)

ml.r <- rs[which(logLike==max(logLike,na.rm=TRUE))]

abline(v=ml.r)
abline(h=max(logLike,na.rm=TRUE))

ml.hf <- make.hap.freq(PA,PB,ml.r)

Я исправил код. Тем не менее, я до сих пор не пришел к правильному ответу.

Используя команду 'dmultinom (c (Gmat), prob = c (gf), log = TRUE)', рассчитайте вероятность записи в журнал. Существует способ найти значение d или набор значений частоты гаплотипа, который максимизирует эту вероятность, и его результат называется оценкой максимальной вероятности частоты гаплотипа для 2 SNV. 

PA=0.7
PB=0.6
Pa=0.3
Pb=0.4
d=0.1
r <- d /sqrt(PA*Pa*PB*Pb)
make.hap.freq<-function(PA,PB,r){
    ret1 = c(PA*PB,PA*Pb,Pa*PB,Pa*Pb)
    ret = ret1 + c(d,-d,-d,d)
    ret
}

make.gen.freq <- function(PA,PB,r){
    h <- make.hap.freq(PA,PB,r)
    g <- matrix(0,3,3)
    g[1,1] <- h[1]^2
    g[1,3] <- h[2]^2
    g[3,1] <- h[3]^2
    g[3,3] <- h[4]^2
    g[1,2] <- 2*h[1] * h[2]
    g[2,1] <- 2*h[1] * h[3]
    g[2,3] <- 2*h[2] * h[4]
    g[3,2] <- 2*h[3] * h[4]
    g[2,2] <- 2*(h[1] * h[4] + h[2] * h[3])
    g
}


make.hap.freq(PA,PB,r)

make.gen.freq(PA,PB,r)


#########
## 7 Diplotype data simulation for 2 SNVs
#########

gf <- make.gen.freq(PA,PB,r)

#########
## 8 Maximum likelihood estimation of haplotype frequency
#########
# Assume you have 3x3 diplotype counts data only.
# Based on the table, you can estimate haplotype frequency.

N <- 1000
G <- sample(1:9, N, replace=TRUE,prob =gf)
tabulate(G)
matrix(tabulate(G), 3,3)

tabulate(G)
Gmat <- matrix(tabulate(G), 3,3)
Gmat.rsum <- apply(Gmat,1,sum)
Gmat.csum <- apply(Gmat,2,sum)

PA <- (2*Gmat.rsum[1]+Gmat.rsum[2])/(2*sum(Gmat.rsum))
PB <- (2*Gmat.csum[1]+Gmat.csum[2])/(2*sum(Gmat.csum))


rs <- seq(from=-1,to=1,length=100)
logLike <- rep(NA,length(rs))
for(i in 1:length(rs)){
    tmp.r <- rs[i]
    gf <- make.gen.freq(PA,PB,tmp.r)


############
######This is the problematic loop:


    if(min(gf)>=0){
        logLike[i] <- dmultinom(c(Gmat),prob=c(gf),log=TRUE)
    }
}
plot(rs,logLike)

ml.r <- rs[which(logLike==max(logLike,na.rm=TRUE))]

abline(v=ml.r)
abline(h=max(logLike,na.rm=TRUE))

ml.hf <- make.hap.freq(PA,PB,ml.r)
...