library(kinship)
library(MCMCpack)
library(gtools)
library(sets)
library(paramlink)
MakePedigreeFromFamilyInfo<-function(p){
ns<-length(p[,1])
affected<-status<-rep(1,ns)
affected[which(p[,5]==2)]<-0
affected[which(p[,5]==3)]<-0
status[which(p[,5]==1)]<-0
status[which(p[,5]==2)]<-0
ptemp<-pedigree(id=p[,1],dadid=p[,3],momid=p[,2],sex=p[,4],affected=affected,status=status)
if(sum(ptemp$affected)==0)ptemp$affected<-affected
ptemp
}
MakeAlleleProb<-function(NL=15){
Alleles<-NULL
Probs<-NULL
library(MCMCpack)
for(i in 1:NL){
na<-sample(2:5,1)
st<-sample(10:20,1)
tmp<-NULL
for(j in 1:na){
a<-paste("",st+j,sep="")
tmp<-c(tmp,a)
}
Alleles[[i]]<-tmp
Probs[[i]]<-c(rdirichlet(1,rep(1,na)))
}
list(Alleles=Alleles,Probs=Probs)
}
RandomGenotype<-function(N,n,Alleles,Probs){
ns<-N
g<-array(0,c(ns,2,n))
noparents<-1:ns
for(i in noparents){
for(j in 1:n){
g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
}
}
g
}
RandomGenotypeFamily<-function(d,n,Alleles,Probs,missing=FALSE){
ns<-length(d[,1])
g<-array(0,c(ns,2,n))
noparents<-which(d[,2]==0,d[,3]==0)
for(i in noparents){
for(j in 1:n){
g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
}
}
given<-rep(0,ns)
given[noparents]<-1
loop<-TRUE
if(sum(given)==ns)loop<-FALSE
while(loop){
yet<-which(given==0)
for(i in yet){
prt1<-d[i,2]
prt2<-d[i,3]
if(given[prt1]+given[prt2]==2){
for(j in 1:n){
g[i,1,j]<-sample(g[prt1,,j],1)
g[i,2,j]<-sample(g[prt2,,j],1)
}
given[i]<-1
}
}
if(sum(given)==ns)loop<-FALSE
}
if(missing){
Type<-p[,5]
G2<-array(0,c(length(p[,1]),2,NL))
G2[which(Type==1),,]<-g[which(Type==1),,]
g<-G2
}
g
}
MakeSearchPattern<-function(v){
library(gtools)
nt<-length(v)
x<-0
X<-NULL
for(i in 1:nt){
tmp<-combinations(nt,i)
tmp2<-tmp
tmp2[1:length(tmp2)]<-c(v[tmp])
X[[i]]<-tmp2
x<-sum(x,length(as.matrix(X[[i]])[,1]))
}
list(X=X,x=x)
}
MakeTestPattern<-function(candidates){
gr<-expand.grid(candidates)
n<-length(gr[,1])
ok<-rep(1,n)
for(i in 1:n){
tmp<-unlist(gr[i,])
tmp2<-outer(tmp,tmp,"-")
diag(tmp2)<-1
if(prod(tmp2)==0){
ok[i]<-0
}
}
as.matrix(gr[which(ok==1),])
}
OffspringGenotypeProb<-function(p1,p2){
V<-p1%*%t(p2)
U<-V+t(V)
U[lower.tri(U)]<-0
diag(U)<-diag(U)/2
U
}
LogLikelihoodTrio<-function(p1,p2,c){
V<-p1%*%t(p2)
U<-V+t(V)
U[lower.tri(U)]<-0
diag(U)<-diag(U)/2
}
CalcLogLikelihoodFamily<-function(p,g,Alleles,Probs){
nl<-length(g[1,1,])
ret<-rep(0,nl)
for(i in 1:length(g[1,1,])){
tmp<-CalcLogLikelihoodFamilyLocus(p,g[,,i],Alleles[[i]],Probs[[i]])
ret[i]<-tmp$loglikesum
}
list(loglikelist=ret,loglikesum=sum(ret))
}
LimitDiplotypeZ<-function(p,g,A){
unknown="0"
ns<-length(p[,1])
na<-length(A)
D<-matrix(0,ns,na)
H<-array(0,c(ns,2,na))
for(i in 1:ns){
if(g[i,1]!=unknown){
tmp<-rep(0,na)
tmp[which(A==g[i,1])]<-1
tmp[which(A==g[i,2])]<-1
D[i,which(tmp==1)]<-1
D[i,which(tmp==0)]<-(-1)
H[i,1,which(tmp==0)]<-(-1)
H[i,2,which(tmp==0)]<-(-1)
}
}
loop<-TRUE
while(loop){
tmpH<-H
tmpD<-D
for(i in 1:ns){
if(p[i,2]!=0){
H[i,1,which(D[p[i,2],]==-1)]<-(-1)
if(length(which(D[p[i,2],]>=0))==1){
H[i,1,which(D[p[i,2],]==1)]<-1
H[i,1,which(D[p[i,2],]!=1)]<-(-1)
}
}
if(length(which(H[i,1,]==-1))==(na-1))H[i,1,which(H[i,1,]>(-1))]<-1
if(length(which(H[i,1,]==1))==1)H[i,1,which(H[i,1,]!=1)]<--1
D[i,which(H[i,1,]==1)]<-1
if(p[i,3]!=0){
H[i,2,which(D[p[i,3],]==-1)]<-(-1)
if(length(which(D[p[i,3],]>=0))==1){
H[i,2,which(D[p[i,3],]==1)]<-1
H[i,2,which(D[p[i,3],]!=1)]<-(-1)
}
}
if(length(which(H[i,2,]==-1))==(na-1))H[i,2,which(H[i,2,]>(-1))]<-1
if(length(which(H[i,2,]==1))==1)H[i,2,which(H[i,2,]!=1)]<--1
D[i,which(H[i,2,]==1)]<-1
H[i,1,which(((D[i,]==1) * (H[i,2,]==-1))==1)]<-1
H[i,2,which(((D[i,]==1) * (H[i,1,]==-1))==1)]<-1
D[p[i,2],which(H[i,1,]==1)]<-1
D[p[i,3],which(H[i,2,]==1)]<-1
if(length(which(D[i,]==1))==2)D[i,which(D[i,]!=1)]<-(-1)
for(i in 1:ns){
D[i,which((H[i,1,]==-1) * (H[i,2,]==-1)==1)]<-(-1)
if(length(which(D[i,]==-1))>=(na-1))D[i,which(D[i,]>(-1))]<-1
}
}
if(isTRUE(all.equal(D,tmpD)) & isTRUE(all.equal(H,tmpH)))loop<-FALSE
}
dset<-dsetA<-list()
dtypeset<-dtypesetA<-list()
hset1<-hset2<-hset1A<-hset2A<-list()
for(i in 1:length(D[,1])){
tmpd1<-as.set(which(D[i,]==1))
tmpd0<-as.set(which(D[i,]==0))
tmpd1A<-as.set(A[which(D[i,]==1)])
tmpd0A<-as.set(A[which(D[i,]==0)])
if(set_is_empty(tmpd1)){
dset[[i]]<-set_combn(tmpd0,1) | set_combn(tmpd0,2)
dsetA[[i]]<-set_combn(tmpd0A,1) | set_combn(tmpd0A,2)
}else{
if(set_is_empty(tmpd0)){
dset[[i]]<-set(tmpd1)
dsetA[[i]]<-set(tmpd1A)
}else{
dset[[i]]<-dsetA[[i]]<-set()
tmpd01<-tmpd0 | tmpd1
for(j in tmpd01){
dset[[i]]<-dset[[i]] | set(tmpd1 | j)
}
tmpd01A<-tmpd0A | tmpd1A
for(j in tmpd01A){
dsetA[[i]]<-dsetA[[i]] | set(tmpd1A | j)
}
}
}
hset1[[i]]<-as.set(which(H[i,1,]>=0))
hset2[[i]]<-as.set(which(H[i,2,]>=0))
dtypeset[[i]]<-as.set(hset1[[i]]*hset2[[i]])
hset1A[[i]]<-as.set(A[which(H[i,1,]>=0)])
hset2A[[i]]<-as.set(A[which(H[i,2,]>=0)])
dtypesetA[[i]]<-as.set(hset1A[[i]]*hset2A[[i]])
}
for(i in 1:length(D[,1])){
tmpdset<-set()
for(j in dset[[i]]){
offsOK<-TRUE
oyajudge<-FALSE
oyanashi<-TRUE
for(k in 1:length(p[,1])){
if(p[k,2] == i | p[k,3] == i){
offsOKperchild<-FALSE
oyajudge<-TRUE
oyanashi<-FALSE
for(l in dset[[k]]){
if(!set_is_empty(j & l))offsOKperchild<-TRUE
}
if(!offsOKperchild)offsOK<-FALSE
}
}
hapsOK<-FALSE
if(!(set_is_empty(j & hset1[[i]])) & !(set_is_empty(j & hset2[[i]])) &
length(j & (hset1[[i]] | hset2[[i]]) )==length(j) )hapsOK<-TRUE
if(offsOK & hapsOK){
tmpdset<-tmpdset | set(j)
}
if(oyanashi & hapsOK){
oyajudge=TRUE
tmpdset<-tmpdset | set(j)
}
}
dset[[i]]<-tmpdset
}
for(i in 1:length(D[,1])){
tmpdset<-set()
for(j in dsetA[[i]]){
offsOK<-TRUE
oyajudge<-FALSE
oyanashi<-TRUE
for(k in 1:length(p[,1])){
if(p[k,2] == i | p[k,3] == i){
offsOKperchild<-FALSE
oyajudge<-TRUE
oyanashi<-FALSE
for(l in dsetA[[k]]){
if(!set_is_empty(j & l))offsOKperchild<-TRUE
}
if(!offsOKperchild)offsOK<-FALSE
}
}
hapsOK<-FALSE
if(!(set_is_empty(j & hset1A[[i]])) & !(set_is_empty(j & hset2A[[i]])) &
length(j & (hset1A[[i]] | hset2A[[i]]) )==length(j) )hapsOK<-TRUE
if(offsOK & hapsOK){
tmpdset<-tmpdset | set(j)
}
if(oyanashi & hapsOK){
oyajudge=TRUE
tmpdset<-tmpdset | set(j)
}
}
dsetA[[i]]<-tmpdset
}
dlist<-dlistA<-list()
for(i in 1:length(dset)){
cnt<-1
dlist[[i]]<-list()
for(j in dset[[i]]){
dlist[[i]][[cnt]]<-as.set(j)
cnt<-cnt+1
}
}
for(i in 1:length(dsetA)){
cnt<-1
dlistA[[i]]<-list()
for(j in dsetA[[i]]){
dlistA[[i]][[cnt]]<-as.set(j)
cnt<-cnt+1
}
}
list(D=D,H=H,dlist=dlist,dlistA=dlistA,dset=dset,dsetA=dsetA,dtypeset=dtypeset,dtypesetA=dtypesetA,
hset1=hset1,hset2=hset2,hset1A=hset1A,hset2A=hset2A)
}
TestAllPatterns<-function(p,G,Gpool,candidates,Alleles,Probs){
ptemp<-MakePedigreeFromFamilyInfo(p)
tobeSearched<-which(p[,5]==2)
SearchPattern<-MakeSearchPattern(tobeSearched)
ret<-NULL
ret2<-NULL
ret3<-NULL
counter<-1
for(i in 1:length(SearchPattern$X)){
spx<-as.matrix(SearchPattern$X[[i]])
for(j in 1:length(spx[,1])){
ToTestCurrent<-spx[j,]
tmpG2<-G
tobeEvaluated<-rep(0,length(p[,1]))
tobeEvaluated[which(p[,5]==1)]<-1
tobeEvaluated[which(p[,5]==3)]<-1
tobeEvaluated[ToTestCurrent]<-1
tmpp<-p[which(tobeEvaluated==1),]
tmptmpG2<-tmpG2[which(tobeEvaluated==1),,]
if(length(tmpG2[1,1,])==1){
tmptmpG2<-array(0,c(length(which(tobeEvaluated==1)),2,1))
tmptmpG2[,,1]<-tmpG2[which(tobeEvaluated==1),,]
}
ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list()
counter<-counter+1
numtested<-length(spx[j,])
tmpcandidates<-NULL
for(k in 1:numtested){
selectedToBeTested<-which(tobeSearched==spx[j,k])
tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
}
tmpTestPattern<-MakeTestPattern(tmpcandidates)
counter<-counter+1
numtested<-length(spx[j,])
tmpcandidates<-NULL
for(k in 1:numtested){
selectedToBeTested<-which(tobeSearched==spx[j,k])
tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
}
tmpTestPattern<-MakeTestPattern(tmpcandidates)
for(k in 1:length(tmpTestPattern[,1])){
tmpG2<-G
tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
tobeEvaluated<-rep(0,length(p[,1]))
tobeEvaluated[which(p[,5]==1)]<-1
tobeEvaluated[which(p[,5]==3)]<-1
tobeEvaluated[ToTestCurrent]<-1
tmpp<-p[which(tobeEvaluated==1),]
tmptmpG2<-tmpG2[which(tobeEvaluated==1),,]
if(length(tmpG2[1,1,])==1){
tmptmpG2<-array(0,c(length(which(tobeEvaluated==1)),2,1))
tmptmpG2[,,1]<-tmpG2[which(tobeEvaluated==1),,]
}
ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list(tmpTestPattern[k,])
counter<-counter+1
}
}
}
LRmatrix<-matrix(0,counter,counter)
for(i in 1:(counter-1)){
for(j in 1:(counter-1)){
LRmatrix[i,j]<-ret[[i]][[2]][1]-ret[[j]][[2]][1]
}
}
list(LogLike=ret,LRmatrix=LRmatrix,SearchedID=ret2,CandidateID=ret3,familyinfo=p,pedigree=ptemp,tobeSearched=tobeSearched,candidates=candidates,offeredGenotype=G,poolGenotype=Gpool,Alleles=Alleles,Probs=Probs)
}
counterplus<-function(x,t){
x[1]<-x[1]+1
for(i in 1:(length(x))){
if(x[i]==t[i]){
x[i]<-0
if(i<length(x)){
x[i+1]<-x[i+1]+1
}
}else{
return(x)
}
}
return(x)
}
counterplus2<-function(x,t){
x[1]<-x[1]+1
for(i in 1:(length(x))){
if(x[i]==(t[i]+1)){
x[i]<-1
if(i<length(x)){
x[i+1]<-x[i+1]+1
}
}else{
return(x)
}
}
return(x)
}
subnucs <- function(ped) {
parents <- unique(ped[, 2:3])
parents = parents[-match(0, parents[, 1]), , drop = FALSE]
list1 <- lapply(nrow(parents):1, function(i) {
par = parents[i, ]
c(fa = par[[1]], mo = par[[2]], offs = as.vector(ped[,
1])[which(ped[, 2] == par[[1]] & ped[, 3] ==
par[[2]], useNames = FALSE)])
})
res = list()
i = 1
k = 1
while (length(list1) > 1) {
if (i > length(list1)) {
warning("Loop detected, likelihood calculations will not work.")
return(NULL)
}
link = ((sub <- list1[[i]]) %in% unlist(list1[-i]))
if (sum(link) == 1) {
res[[k]] <- c(pivot = sub[[which(link)]], sub)
list1 <- list1[-i]
k <- k + 1
i <- 1
}
else i <- i + 1
}
res[[k]] <- c(pivot = 0, list1[[1]])
res
}
CalcLogLikelihoodFamilyLocus2<-function(p,g,A,P){
unknown<-"0"
Unfixed<-which(g[,1]==unknown)
Nunfixed<-length(Unfixed)
ns<-length(p[,1])
na<-length(A)
if(Nunfixed==0){
tmpret<-rep(1,ns)
currentg<-g
gvector<-array(0,c(ns,na,2))
ret<-0
for(i in 1:ns){
gvector[i,which(A==currentg[i,1]),1]<-1
gvector[i,which(A==currentg[i,2]),2]<-1
}
trioloop<-TRUE
triocnt<-1
while(trioloop){
P1<-gvector[triocnt,,1]
P2<-gvector[triocnt,,2]
tmp1<-p[triocnt,2]
tmp2<-p[triocnt,3]
if(tmp1!=0){
P1<-(gvector[tmp1,,1]+gvector[tmp1,,2])/2
}
if(tmp2!=0){
P2<-(gvector[tmp2,,1]+gvector[tmp2,,2])/2
}
FromParents<-OffspringGenotypeProb(P1,P2)
tmpallele1<-which(A==currentg[triocnt,1])
tmpallele2<-which(A==currentg[triocnt,2])
lessa<-min(tmpallele1,tmpallele2)[1]
morea<-max(tmpallele1,tmpallele2)[1]
tmptmp<-FromParents[lessa,morea]
tmpret[triocnt]<-tmptmp
if(tmptmp==0)trioloop<-FALSE
triocnt<-triocnt+1
if(triocnt==(ns+1))trioloop<-FALSE
}
tmprob<-prod(tmpret)
ret<-ret+tmprob
}else{
Ng<-na*(na+1)/2
diplotype<-matrix(0,Ng,2)
diplotypeProb<-rep(0,Ng)
cnt<-1
for(i in 1:na){
for(j in 1:i){
diplotype[cnt,1]<-A[i]
diplotype[cnt,2]<-A[j]
diplotypeProb[cnt]<-P[i]*P[j]
if(i!=j)diplotypeProb[cnt]<-2*diplotypeProb[cnt]
cnt<-cnt+1
}
}
counters<-rep(0,Nunfixed)
nchild<-rep(Ng,Nunfixed)
loop<-TRUE
ret<-0
noParentProb<-1
while(loop){
tmpret<-rep(1,ns)
currentgenotype<-counters+1
currentg<-g
counters<-counterplus(counters,nchild)
currentg[Unfixed,]<-diplotype[currentgenotype,]
gvector<-array(0,c(ns,na,2))
for(i in 1:ns){
gvector[i,which(A==currentg[i,1]),1]<-1
gvector[i,which(A==currentg[i,2]),2]<-1
}
trioloop<-TRUE
triocnt<-1
while(trioloop){
P1<-P2<-P
tmp1<-p[triocnt,2]
tmp2<-p[triocnt,3]
if(tmp1!=0){
P1<-(gvector[tmp1,,1]+gvector[tmp1,,2])/2
}
if(tmp2!=0){
P2<-(gvector[tmp2,,1]+gvector[tmp2,,2])/2
}
FromParents<-OffspringGenotypeProb(P1,P2)
Offspring<-OffspringGenotypeProb(gvector[triocnt,,1],gvector[triocnt,,2])
tmptmp<-sum(FromParents*Offspring)
tmpret[triocnt]<-tmptmp
if(tmptmp==0)trioloop<-FALSE
triocnt<-triocnt+1
if(triocnt==(ns+1))trioloop<-FALSE
}
tmprob<-prod(tmpret)
ret<-ret+tmprob
if(sum(counters)==0){
loop<-FALSE
}
}
}
list(loglikesum=log(ret,10),like=ret)
}
gameteFreq<-function(P){
(apply(P,1,sum)+apply(P,2,sum))/2
}
OffspringDiplotypeFreq<-function(m,f){
gameteM<-gameteFreq(m)
gameteF<-gameteFreq(f)
outer(gameteM,gameteF,FUN="*")
}
LimitDiplotype<-function(p,g,A){
unknown="0"
ns<-length(p[,1])
na<-length(A)
a<-array(0,c(ns,na,na))
for(i in 1:ns){
if(g[i,1]==unknown){
a[i,,]<-1/(na^2)
}else{
tmpa1<-which(A==g[i,1])
tmpa2<-which(A==g[i,2])
a[i,tmpa1,tmpa2]<-a[i,tmpa2,tmpa1]<-1
}
}
for(i in 1:ns){
tmpm<-p[i,2]
tmpf<-p[i,3]
if(tmpm!=0){
print(i)
print(tmpm)
print(tmpf)
print(a[tmpm,,])
print(a[tmpf,,])
a[i,,]<-a[i,,]*OffspringDiplotypeFreq(a[tmpm,,],a[tmpf,,])
}
a[i,,]<-a[i,,]+t(a[i,,])
}
sign(a)
}
MakeVectorFromDiplotypeMat<-function(m){
tmp<-m
tmp[,]<-1
tmp[upper.tri(tmp)]<-0
m[which(tmp==1,arr.ind=TRUE)]
}
DiplotypeListFromMat<-function(m){
which(MakeVectorFromDiplotypeMat(m)==1)
}
CalcLogLikelihoodFamily<-function(p,g,Alleles,Probs){
nl<-length(g[1,1,])
ret<-rep(0,nl)
for(i in 1:length(g[1,1,])){
tmp<-CalcLogLikelihoodFamilyLocus2(p,g[,,i],Alleles[[i]],Probs[[i]])
ret[i]<-tmp$loglikesum
}
list(loglikelist=ret,loglikesum=sum(ret))
}
printTAPout<-function(tap){
n<-length(tap$LogLike)
for(i in 1:n){
if(tap$LogLike[[i]][[2]][1]!=-Inf){
tmp<-paste("LogLike",tap$LogLike[[i]][[2]][1],"Prob",10^(tap$LogLike[[i]][[2]][1]),"SearchedID",c(tap$SearchedID[[i]]),"CandidateID",c(tap$CandidateID[[i]]))
print(tmp)
}
}
}
printSMPFout<-function(SMMFout){
n<-length(SMMFout$LogLike)
plot(SMMFout$pedigree,id=unlist(SMMFout$IndNames),main=SMMFout$FamilyName)
n<-length(SMMFout$LogLike)
for(i in 1:n){
if(SMMFout$LogLike[[i]]!=-Inf){
tmp<-paste("LogLike",SMMFout$LogLike[[i]],"Prob",10^(SMMFout$LogLike[[i]]),"SearchedID",c(SMMFout$SearchedID[[i]]),"CandidateID",c(SMMFout$CandidateID[[i]]))
print(tmp)
}
}
}
TestAllPatterns2<-function(p,G,Gpool,candidates,Alleles,Probs){
ptemp<-MakePedigreeFromFamilyInfo(p)
nucs<-subnucs(p)
tobeSearched<-which(p[,5]==2)
SearchPattern<-MakeSearchPattern(tobeSearched)
LikeFromGenPop<-rep(0,length(Gpool[,1,1]))
for(i in 1:length(Alleles)){
tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
for(j in 1:length(LikeFromGenPop)){
tmpallele1<-which(Alleles[[i]]==Gpool[j,1,i])
tmpallele2<-which(Alleles[[i]]==Gpool[j,2,i])
tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
}
}
ret<-NULL
ret2<-NULL
ret3<-NULL
counter<-1
for(i in 1:length(SearchPattern$X)){
spx<-as.matrix(SearchPattern$X[[i]])
for(j in 1:length(spx[,1])){
ToTestCurrent<-spx[j,]
tmpG2<-G
BaseLike<-CalcLogLikelihoodFamilyX(p,nucs,G,Alleles,Probs)$loglikesum
numtested<-length(spx[j,])
tmpcandidates<-NULL
for(k in 1:numtested){
selectedToBeTested<-which(tobeSearched==spx[j,k])
tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
}
print("done for none")
tmpTestPattern<-MakeTestPattern(tmpcandidates)
for(k in 1:length(tmpTestPattern[,1])){
tmpG2<-G
tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
LoutTmp<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list(tmpTestPattern[k,])
counter<-counter+1
}
}
}
LRmatrix<-matrix(0,counter,counter)
for(i in 1:(counter-1)){
for(j in 1:(counter-1)){
LRmatrix[i,j]<-ret[[i]]-ret[[j]]
}
}
list(LogLike=ret,LRmatrix=LRmatrix,SearchedID=ret2,CandidateID=ret3,familyinfo=p,pedigree=ptemp,tobeSearched=tobeSearched,candidates=candidates,offeredGenotype=G,poolGenotype=Gpool,Alleles=Alleles,Probs=Probs)
}
printTAPout2<-function(tap){
n<-length(tap$LogLike)
for(i in 1:n){
if(tap$LogLike[[i]]!=-Inf){
tmp<-paste("LogLike",tap$LogLike[[i]],"Prob",10^(tap$LogLike[[i]]),"SearchedID",c(tap$SearchedID[[i]]),"CandidateID",c(tap$CandidateID[[i]]))
print(tmp)
}
}
}
MakeInfo<-function(LDZout,gnew){
ret<-LDZout$dlistA
unknown="0"
for(i in 1:length(gnew[,1])){
if(gnew[i,1]!=unknown){
ret[[i]]<-list(as.set(gnew[i,]))
}
}
ret
}
MakeDiplotypePrior<-function(p,g,A,P,LDZout,info){
ret<-list()
for(i in 1:length(p[,1])){
ret[[i]]<-rep(1,length(LDZout$dlist[[i]]))
if(length(info[[i]])>1){
if(p[i,2]=="0"){
tmph1<-LDZout$hset1A[[i]]
tmph2<-LDZout$hset2A[[i]]
for(j in 1:length(LDZout$dlist[[i]])){
tmpp<-1
heterocheck<-TRUE
for(k in LDZout$dlistA[[i]][[j]]){
if(!(k %e% (tmph1 & tmph2)))heterocheck<-FALSE
tmpp<-tmpp*P[which(A==k)]
}
if(length(LDZout$dlistA[[i]][[j]])==2){
if(heterocheck){
tmpp<-tmpp*2
}
}
ret[[i]][j]<-tmpp
}
ret[[i]]<-ret[[i]]/sum(ret[[i]])
}
}else{
ret[[i]]<-rep(0,length(LDZout$dlist[[i]]))
for(j in 1:length(LDZout$dlist[[i]])){
if(info[[i]][[1]] == LDZout$dlistA[[i]][[j]]){
ret[[i]][j]=1
}
}
}
}
ret
}
CalcProbNucs<-function(p,g,A,P,nucs,LDZout){
rets<-list()
for(nn in 1:length(nucs)){
nuc<-nucs[[nn]]
dimvector<-nuc[2:length(nuc)]
fa<-nuc[2]
mo<-nuc[3]
offs<-nuc[4:length(nuc)]
pivot<-nuc[1]
if(pivot!=0){
who<-which(nuc[2:length(nuc)]==pivot)
dimvector<-nuc[2:length(nuc)]
pivotLoc<-which(dimvector==pivot)
dimvector<-c(dimvector[-pivotLoc],pivot)
}
genotypeNumVector<-rep(0,length(dimvector))
for(i in 1:length(dimvector)){
genotypeNumVector[i]<-length(LDZout$dlist[[dimvector[i]]])
}
genotypeNumVector
ret<-array(1,genotypeNumVector)
faid<-which(dimvector==fa)
moid<-which(dimvector==mo)
ofid<-which(dimvector!=fa & dimvector!=mo)
counter<-array(1:prod(genotypeNumVector),genotypeNumVector)
for(i in 1:length(ret)){
x<-which(counter==i,arr.ind=TRUE)
tmpprob<-1
faset<-LDZout$dlist[[dimvector[faid]]][[x[faid]]]
moset<-LDZout$dlist[[dimvector[moid]]][[x[moid]]]
for(j in ofid){
ofset<-LDZout$dlist[[dimvector[j]]][[x[j]]]
tmptmpprob<-LikeTrio(faset,moset,ofset)
tmpprob<-tmpprob*tmptmpprob
}
ret[x]<-tmpprob
}
rets[[nn]]<-list(ret,dimvector)
}
rets
}
LikeTrio<-function(setf,setm,seto){
couple<-setf*setm
cnt<-0
for(i in couple){
if(as.set(i) == seto)cnt<-cnt+1
}
cnt/length(couple)
}
LimitAlleles<-function(g,A,P){
tmp<-rep(0,length(A))
for(i in 1:length(g)){
tmp[which(A==g[i])]<-1
}
A2<-c(A[which(tmp==1)],"pool")
P2<-c(P[which(tmp==1)])
P2<-c(P2,1-sum(P2))
list(A=A2,P=P2)
}
LikeNucWithPrior<-function(cpnout,nucs,DiplotypePrior){
prob<-list()
for(nn in 1:length(nucs)){
nucDimVector<-cpnout[[nn]][[2]]
tmp<-DiplotypePrior[[nucDimVector[1]]]
for(i in 2:length(nucDimVector)){
tmp<-tmp%o%DiplotypePrior[[nucDimVector[i]]]
}
prob[[nn]]<-tmp * cpnout[[nn]][[1]]
}
prob
}
SumPivot<-function(cpnout,like,info){
ret<-NULL
pivotted<-rep(0,length(info))
cumulProb<-list()
for(nn in 1:length(cpnout)){
tmpdim<-dim(cpnout[[nn]][[1]])
tmp<-rep(1,tmpdim[1])
if(pivotted[[cpnout[[nn]][[2]][[1]]]]==1){
tmp<-cumulProb[[cpnout[[nn]][[2]][[1]]]]
}
for(i in 2:length(tmpdim)){
tmp2<-rep(1,tmpdim[i])
if(pivotted[[cpnout[[nn]][[2]][[i]]]]==1){
tmp2<-cumulProb[[cpnout[[nn]][[2]][[i]]]]
}
tmp<-tmp%o%tmp2
}
tmp3<-like[[nn]]*tmp
ret<-apply(tmp3,length(cpnout[[nn]][[2]]),sum)
pivotted[cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]]]<-1
cumulProb[[cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]]]]<-ret
}
prob<-sum(ret)
list(like=prob,loglike=log(prob,10))
}
CalcLikeZ<-function(p,G,nucs,Alleles,Probs){
tmpret<-0
for(na in 1:length(Alleles)){
A<-unlist(Alleles[[na]])
P<-Probs[[na]]
g<-G[,,na]
gpool<-G[,,na]
A2P2<-LimitAlleles(gpool,A,P)
A2<-A2P2[[1]]
P2<-A2P2[[2]]
LDZout2<-LimitDiplotypeZ(p,g,A2)
MCheck<-TRUE
for(i in 1:length(LDZout2$dlistA)){
if(length(LDZout2$dlistA[[i]])==0){
print("checkFALSE")
print(p)
print(g)
print(LDZout2$dlistA)
MCheck<-FALSE
}
}
if(!MCheck){
return(-Inf)
}
cpnout2<-CalcProbNucs(p,g,A2,P2,nucs,LDZout2)
info2<-LDZout2$dlistA
DiplotypePrior2<-MakeDiplotypePrior(p,g,A2,P2,LDZout2,info2)
like2<-LikeNucWithPrior(cpnout2,nucs,DiplotypePrior2)
tmpret<-tmpret+SumPivot(cpnout2,like2,info2)$loglike
}
tmpret
}
SearchMissingsMultiFamily<-function(pedigrees,genotypesFamily,Gpool,candidatesList,FamilyNames,IndNames,Alleles,Probs){
output<-list()
LikeFromGenPop<-rep(0,length(Gpool[,1,1]))
for(i in 1:length(Alleles)){
tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
for(j in 1:length(LikeFromGenPop)){
tmpallele1<-which(Alleles[[i]]==Gpool[j,1,i])
tmpallele2<-which(Alleles[[i]]==Gpool[j,2,i])
tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
}
}
for(ip in 1:length(pedigrees)){
output[[ip]]<-SearchMissingsPerFamily(pedigrees[[ip]],genotypesFamily[[ip]],
Gpool,candidatesList[[ip]],FamilyNames[[ip]],IndNames[[ip]],Alleles,Probs,LikeFromGenPop)
}
output
}
SearchMissingsPerFamily<-function(p,G,Gpool,candidatesL,FamilyNs,IndNs,Alleles,Probs,LikeFromGenPop){
ptemp<-MakePedigreeFromFamilyInfo(p)
tmpPed<-p[,1:5]
tmpPed<-data.frame(ID=tmpPed[,1],FID=tmpPed[,3],MID=tmpPed[,2],SEX=tmpPed[,4]+1,AFF=tmpPed[,5]-1)
tmpPed<-linkdat(tmpPed,model=1)
plot(ptemp,id=unlist(IndNs))
pdffile<-paste(FamilyNs,".pdf")
pdf(pdffile,family="Japan1")
plot(ptemp,id=unlist(IndNs))
dev.off()
nucs<-subnucs(p)
searched<-which(p[,5]==2)
candidates<-list()
for(i in 1:length(searched)){
candidates[[i]]<-candidatesL[[searched[[i]]]]
}
BaseMade<-FALSE
SearchPattern<-MakeSearchPattern(searched)
ret<-NULL
ret2<-NULL
ret3<-NULL
counter<-1
for(i in 1:length(SearchPattern$X)){
spx<-as.matrix(SearchPattern$X[[i]])
for(j in 1:length(spx[,1])){
ToTestCurrent<-spx[j,]
numtested<-length(ToTestCurrent)
tmpcandidates<-list()
if(numtested>0){
for(k in 1:numtested){
selectedToBeTested<-which(searched==ToTestCurrent[k])
tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
}
}
tmpTestPattern<-MakeTestPattern(tmpcandidates)
for(k in 1:length(tmpTestPattern[,1])){
tmpG2<-G
tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
tttmpG<-tmpG2[,,1]
for(i in 2:length(tmpG2[1,1,])){
tttmpG<-cbind(tttmpG,tmpG2[,,i])
}
tmpPed2<-setMarkers(tmpPed, m=tttmpG)
mCout<-mendelianCheck(tmpPed2)
if(length(mCout)!=0){
ret[[counter]]<--Inf
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list(tmpTestPattern[k,])
counter<-counter+1
}else{
if(!BaseMade){
BaseLike<-CalcLikeZ(p,G,nucs,Alleles,Probs)
BaseMade<-TRUE
}
LoutTmp<-CalcLikeZ(p,tmpG2,nucs,Alleles,Probs)
ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list(tmpTestPattern[k,])
counter<-counter+1
}
}
}
}
list(LogLike=ret,SearchedID=ret2,CandidateID=ret3,pedigree=ptemp,FamilyName=FamilyNs,IndNames=IndNs)
}
CalcLikeForCandidates<-function(p,G,candidates,Gcandidates,Alleles,Probs){
LikeFromGenPop<-rep(0,length(Gpool2[,1,1]))
for(i in 1:length(Alleles)){
tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
for(j in 1:length(LikeFromGenPop)){
tmpallele1<-which(Alleles[[i]]==Gpool2[j,1,i])
tmpallele2<-which(Alleles[[i]]==Gpool2[j,2,i])
tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
}
}
nucs<-subnucs(p)
BaseLike<-CalcLikeZ(p,G2,nucs,Alleles,Probs)
SearchPattern<-MakeSearchPattern(searched)
ret<-NULL
ret2<-NULL
ret3<-NULL
counter<-1
for(i in 1:length(SearchPattern$X)){
spx<-as.matrix(SearchPattern$X[[i]])
for(j in 1:length(spx[,1])){
ToTestCurrent<-spx[j,]
tmpG2<-G
numtested<-length(ToTestCurrent)
tmpcandidates<-list()
if(numtested>0){
for(k in 1:numtested){
selectedToBeTested<-which(searched==ToTestCurrent[k])
tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
}
}
tmpTestPattern<-MakeTestPattern(tmpcandidates)
for(k in 1:length(tmpTestPattern[,1])){
tmpG2<-G
tmpG2[ToTestCurrent,,]<-Gpool2[tmpTestPattern[k,],,]
LoutTmp<-CalcLikeZ(p,tmpG2,nucs,Alleles,Probs)
ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
ret2[[counter]]<-list(ToTestCurrent)
ret3[[counter]]<-list(tmpTestPattern[k,])
counter<-counter+1
}
}
}
list(LogLike=ret,SearchedID=ret2,CandidateID=ret3)
}
Identifiler15<-c("D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01","D13S317","D16S539","D2S1338","D19S433","vWA","TPROX","D18S51","D5S818","FGA")
AllelesId<-NULL
ProbsIDYoshidaN<-NULL
ProbsIDYoshidaFreq<-NULL
ProbsIDYoshidaFreq2<-NULL
AllelesId[[1]]<-c("7",9,10,11,12,13,14,15,16,17,18)
AllelesId[[2]]<-c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2)
AllelesId[[3]]<-c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15)
AllelesId[[4]]<-c("7",8,9,10,11,12,13,14,15,16)
AllelesId[[5]]<-c("12",13,14,15,16,17,18,19,21)
AllelesId[[6]]<-c("5",6,7,8,9,9.3,10)
AllelesId[[7]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[8]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[9]]<-c("16",17,18,19,20,21,22,23,24,25,26,27,28)
AllelesId[[10]]<-c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18)
AllelesId[[11]]<-c("13",14,15,16,17,18,19,20,21,22)
AllelesId[[12]]<-c("7",8,9,10,11,12,13,14)
AllelesId[[13]]<-c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27)
AllelesId[[14]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[15]]<-c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28)
ProbsIDYoshidaN[[1]]<-c(1,2,358,294,331,607,553,363,173,17,1)
ProbsIDYoshidaN[[2]]<-c(3,114,13,666,910,12,3,279,192,52,321,7,6,109,2,11)
ProbsIDYoshidaN[[3]]<-c(8,342,123,1,591,1,2,887,634,94,16,1)
ProbsIDYoshidaN[[4]]<-c(27,3,126,603,561,1130,186,48,14,2)
ProbsIDYoshidaN[[5]]<-c(5,2,78,1069,827,539,171,8,1)
ProbsIDYoshidaN[[6]]<-c(4,603,720,179,1076,94,24)
ProbsIDYoshidaN[[7]]<-c(3,721,348,310,597,546,138,35,2)
ProbsIDYoshidaN[[8]]<-c(1,5,955,531,505,482,196,22,3)
ProbsIDYoshidaN[[9]]<-c(23,263,429,564,285,40,136,396,291,166,78,23,6)
ProbsIDYoshidaN[[10]]<-c(1,3,10,2,109,14,776,81,944,238,137,310,14,52,8,1)
ProbsIDYoshidaN[[11]]<-c(1,524,72,497,763,609,199,27,7,1)
ProbsIDYoshidaN[[12]]<-c(1,1227,308,83,980,96,3,2)
ProbsIDYoshidaN[[13]]<-c(6,12,129,538,600,454,339,220,1,130,98,59,42,35,20,10,4,2,1)
ProbsIDYoshidaN[[14]]<-c(7,18,232,542,789,635,450,24,3)
ProbsIDYoshidaN[[15]]<-c(9,58,181,240,353,544,5,554,13,424,3,198,5,86,21,6)
for(i in 1:15){
ProbsIDYoshidaFreq[[i]]<-ProbsIDYoshidaN[[i]]/sum(ProbsIDYoshidaN[[i]])
ProbsIDYoshidaFreq2[[i]]<-(ProbsIDYoshidaN[[i]]+1)/(sum(ProbsIDYoshidaN[[i]])+1)
ProbsIDYoshidaFreq2[[i]][which(ProbsIDYoshidaN[[i]]<5)]<-5/sum(ProbsIDYoshidaN[[i]])
}