TestAllPatternsの出力のうち、対数尤度、被捜索者リスト、候補者リストを並べて出力する
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- printTAPout
- タイトル
printTAPout
- 説明
TestAllPatternsの出力のうち、対数尤度、被捜索者リスト、候補者リストを並べて出力する
- 使用例
printTAPout(tap)
- ソース
printTAPout<-function(tap){ n<-length(tap$LogLike) for(i in 1:n){ tmp<-paste("LogLike",tap$LogLike[[i]][2],"SearchedID",c(tap$SearchedID[[i]]),"CandidateID",c(tap$CandidateID[[i]])) print(tmp) } }
- Rdファイル
\name{printTAPout} \alias{printTAPout} %- Also NEED an '\alias' for EACH other topic documented here. \title{ printTAPout } \description{ TestAllPatternsの出力のうち、対数尤度、被捜索者リスト、候補者リストを並べて出力する } \usage{ printTAPout(tap) } \arguments{ TestAllPatternsの出力オブジェクト } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ 標準出力 } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ p<-matrix( c(1:8, 0,0,1,0,3,3,3,3, 0,0,2,0,4,4,4,4, 1,0,1,0,1,1,0,0, 1,1,2,1,1,1,2,1), ncol=5) NL<-15 AandP<-MakeAlleleProb(NL=NL) Alleles<-AandP$Alleles Probs<-AandP$Probs G<-RandomGenotypeFamily(p,NL,Alleles,Probs) tap<-TestAllPatterns(p=p,G=G,Gpool=Gpool,candidates=candidates,Alleles=Alleles,Probs=Probs) printTAPout(tap) }
ある家系に1人以上の被捜索者がいるときにすべての捜索パターンをテストする
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- TestAllPatterns
- タイトル
TestAllPatterns
- 説明
ある家系に属する1人以上のメンバーが被捜索者となっていて、そのうちの0人以上、全員以下が候補者プールに含まれうる。被捜索者は候補者プールのうちの一部でもありえるし、全員でもあり得る。そのうえで、被捜索者は、1人から全員の全通りについて、また、候補者については、すべての場合について、尤度を計算する
- 使用例
TestAllPatterns(p=p,G=G2,Gpool=Gpool2,candidates=candidates,Alleles=Alleles,Probs=Probs)
- ソース
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 #print(tmpG2[,,1]) 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),,] } #print(tmptmpG2) #print("---") #print(tmpG2) #print(tobeEvaluated) #print(tmpp) #print(tmptmpG2) ret[[counter]]<-CalcLogLikelihoodFamily(tmpp,tmptmpG2,Alleles,Probs) ret2[[counter]]<-list(ToTestCurrent) ret3[[counter]]<-list() #print(ret[[counter]]) 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])){ #print("to be tested IDs in pedigree") #print(ToTestCurrent) #print("to be tested sampleIDs in pool") #print(tmpTestPattern[k,]) 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),,] } #print("---") #print(tobeEvaluated) #print(tmpp) #print(tmptmpG2) ret[[counter]]<-CalcLogLikelihoodFamily(tmpp,tmptmpG2,Alleles,Probs) #ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs) ret2[[counter]]<-list(ToTestCurrent) ret3[[counter]]<-list(tmpTestPattern[k,]) #print(ret[[counter]]) counter<-counter+1 } } } list(LogLike=ret,SearchedID=ret2,CandidateID=ret3,familyinfo=p,pedigree=ptemp,tobeSearched=tobeSearched,candidates=candidates,offeredGenotype=G,poolGenotype=Gpool,Alleles=Alleles,Probs=Probs) }
- Rdファイル
\name{TestAllPatterns} \alias{TestAllPatterns} %- Also NEED an '\alias' for EACH other topic documented here. \title{ TestAllPatterns } \description{ ある家系に属する1人以上のメンバーが被捜索者となっていて、そのうちの0人以上、全員以下が候補者プールに含まれうる。被捜索者は候補者プールのうちの一部でもありえるし、全員でもあり得る。そのうえで、被捜索者は、1人から全員の全通りについて、また、候補者については、すべての場合について、尤度を計算する } \usage{ TestAllPatterns(p=p,G=G,Gpool=Gpool,candidates=candidates,Alleles=Alleles,Probs=Probs) } \arguments{ \item{p}{家系関係行列} \item{G}{ジェノタイプ} \item{Gpool}{候補者を含むプールサンプルのジェノタイプ} \item{candidates}{被捜索者ごとのプール内候補ベクトルを要素とするリスト} \item{Alleles}{アレル名ベクトルのリスト} \item{Probs}{アレル頻度ベクトルのリスト} } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ \item{LogLike}{全捜索パターンの対数尤度} \item{SearchedID}{全捜索パターンの被捜索対象者をベクトルとするリスト} \item{CandidateID}{全捜索パターンの被捜索対象者ごとの候補者ID} \item{familyinfo}{家系情報行列} \item{pedigree}{家系pedigreeオブジェクト} \item{tobeSearched}{家系内の被捜索者ベクトル} \item{candidates}{被捜索者ごとの候補者IDベクトルを要素とするリスト} \item{offeredGenotype}{提供された家族メンバーのジェノタイプの行列} \item{poolGenotype}{候補者プールのジェノタイプ行列} \item{Alleles}{アレル名ベクトルのリスト} \item{Probs}{アレル頻度のベクトルのリスト} } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ p<-matrix( c(1:8, 0,0,1,0,3,3,3,3, 0,0,2,0,4,4,4,4, 1,0,1,0,1,1,0,0, 1,1,2,1,1,1,2,1), ncol=5) NL<-15 AandP<-MakeAlleleProb(NL=NL) Alleles<-AandP$Alleles Probs<-AandP$Probs G<-RandomGenotypeFamily(p,NL,Alleles,Probs) tap<-TestAllPatterns(p=p,G=G,Gpool=Gpool,candidates=candidates,Alleles=Alleles,Probs=Probs) printTAPout(tap) }
ある家系図であるマーカーのジェノタイプが観測される確率を計算する
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- CalcLogLikelihoodFamily
- タイトル
CalcLogLikelihoodFamily
- 説明
家系図のすべての親子トリオについての、その両親からその子が生まれる尤度を計算し、その総積をとる
- 使用例
CalcLogLikelihoodFamily(p,G,Alleles,Probs)
- ソース
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)) }
- Rdファイル
\name{CalcLogLikelihoodFamily} \alias{CalcLogLikelihoodFamily} %- Also NEED an '\alias' for EACH other topic documented here. \title{ CalcLogLikelihoodFamily } \description{ 尤度をローカスごとに計算し、その総積をとった上で、すべてのローカスについて総積をとる } \usage{ CalcLogLikelihoodFamily(p,G,Alleles,Probs) } \arguments{ \item{p}{家系関係行列} \item{G}{ジェノタイプ} \item{Alleles}{アレル名ベクトルのリスト} \item{Probs}{アレル頻度ベクトルのリスト} } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ \item{loglikelist}{各ローカスの尤度の対数} \item{loglikesum}{全ローカスの尤度の総積の対数} } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ p<-matrix( c(1:8, 0,0,1,0,3,3,3,3, 0,0,2,0,4,4,4,4, 1,0,1,0,1,1,0,0, 1,1,2,1,1,1,2,1), ncol=5) NL<-15 AandP<-MakeAlleleProb(NL=NL) Alleles<-AandP$Alleles Probs<-AandP$Probs G<-RandomGenotypeFamily(p,NL,Alleles,Probs) CalcLogLikelihoodFamilyLocus(p,G,Alleles,Probs) }
ある家系図であるマーカーのジェノタイプが観測される確率を計算する
書きかけ
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- CalcLogLikelihoodFamilyLocus2
- タイトル
CalcLogLikelihoodFamilyLocus2
- 説明
- 使用例
CalcLogLikelihoodFamilyLocus2(p,g,A,P)
- ソース
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) } CalcLogLikelihoodFamilyLocus2<-function(p,g,A,P){ unknown<-"0" Unfixed<-which(g[,1]==unknown) Nunfixed<-length(Unfixed) # 不明ジェノタイプ人数 ns<-length(p[,1]) na<-length(A) # すべての未確定ジェノタイプについて場合分けして計算し、その和をとる # 集団ジェノタイプ頻度の影響を受けるのは親が与えられていない人のみ # それ以外の個人で不明ジェノタイプの確率は「組み込まれる」 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-1,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)) #given<-rep(0,ns) for(i in 1:ns){ #tmp<-0 #if(g[i,1]!=unknown){ gvector[i,which(A==currentg[i,1]),1]<-1 #tmp<-tmp+1 #} #if(g[i,2]!=unknown){ gvector[i,which(A==currentg[i,2]),2]<-1 #tmp<-tmp+1 #} #if(tmp==2)given[i]<-1 } # すべての親子トリオについて計算 for(i in 1:ns){ tmp1<-p[i,2] tmp2<-p[i,3] if(tmp1*tmp2>0){ P1<-(gvector[tmp1,,1]+gvector[tmp1,,2])/2 P2<-(gvector[tmp2,,1]+gvector[tmp2,,2])/2 FromParents<-OffspringGenotypeProb(P1,P2) Offspring<-OffspringGenotypeProb(gvector[i,,1],gvector[i,,2]) #ret[i]<-log(sum(FromParents*Offspring)) # 常用対数 tmpret[i]<-sum(FromParents*Offspring) }else{ if(g[i,1]==unknown){ #print(gvector[i,,]) h1<-which(gvector[i,,1]==1) h2<-which(gvector[i,,2]==1) #print(h1) #print(h2) ttt<-P[h1]+P[h2] if(h1!=h2)ttt<-ttt*2 noParentProb<-noParentProb*ttt } } } tmprob<-prod(tmpret)*noParentProb ret<-ret+tmprob print(ret) #print(currentg) if(sum(counters)==0){ loop<-FALSE } } list(loglikesum=log(ret,10)) }
- Rdファイル
\name{CalcLogLikelihoodFamilyLocus2} \alias{CalcLogLikelihoodFamilyLocus2} %- Also NEED an '\alias' for EACH other topic documented here. \title{ CalcLogLikelihoodFamilyLocus2 } \description{ 家系図のすべての親子トリオについての、その両親からその子が生まれる尤度を計算し、その総積をとる } \usage{ CalcLogLikelihoodFamilyLocus2(p,g,A,P) } \arguments{ \item{p}{家系関係行列} \item{g}{ジェノタイプ} \item{A}{アレル名ベクトル} \item{P}{アレル頻度ベクトル} } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ \item{loglikelist}{各トリオの尤度の対数} \item{loglikesum}{全トリオの尤度の総積の対数} } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ p<-matrix( c(1:8, 0,0,1,0,3,3,3,3, 0,0,2,0,4,4,4,4, 1,0,1,0,1,1,0,0, 1,1,2,1,1,1,2,1), ncol=5) NL<-15 AandP<-MakeAlleleProb(NL=NL) Alleles<-AandP$Alleles Probs<-AandP$Probs g<-RandomGenotypeFamily(p,NL,Alleles,Probs) i<-1 CalcLogLikelihoodFamilyLocus2(p,g[,,i],Alleles[[i]],Probs[[i]]) }
メイティングして生じるディプロタイプ頻度を計算する
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- OffspringGenotypeProb
- タイトル
OffspringGenotypeProb
- 説明
任意のアレル数のマーカーについて、父方のアレル頻度ベクトルと母方のアレル頻度ベクトルとから生まれる子のディプロタイプの頻度を上三角行列で出す。 ある集団があって、HWEを想定するなら、2ベクトルに同じものを与えると、その集団のディプロタイプ頻度が出る。 ある個人の父由来のアレルと母由来のアレルをそれぞれ確率1としたベクトルを与えれば、その個人のディプロタイプ確率が1であるような行列が出る。 父と母のディプロタイプがわかっているときには、それぞれのアレル確率(2本の染色体を合わせた確率)を与えると、子のディプロタイプ確率が出る
- 使用例
OffspringGenotypeProb(p1,p2)
- ソース
OffspringGenotypeProb<-function(p1,p2){ V<-p1%*%t(p2) U<-V+t(V) U[lower.tri(U)]<-0 diag(U)<-diag(U)/2 U }
- Rdファイル
\name{OffspringGenotypeProb} \alias{OffspringGenotypeProb} %- Also NEED an '\alias' for EACH other topic documented here. \title{ OffspringGenotypeProb } \description{ 任意のアレル数のマーカーについて、父方のアレル頻度ベクトルと母方のアレル頻度ベクトルとから生まれる子のディプロタイプの頻度を上三角行列で出す。 ある集団があって、HWEを想定するなら、2ベクトルに同じものを与えると、その集団のディプロタイプ頻度が出る。 ある個人の父由来のアレルと母由来のアレルをそれぞれ確率1としたベクトルを与えれば、その個人のディプロタイプ確率が1であるような行列が出る。 父と母のディプロタイプがわかっているときには、それぞれのアレル確率(2本の染色体を合わせた確率)を与えると、子のディプロタイプ確率が出る } \usage{ OffspringGenotypeProb(p1,p2) } \arguments{ 自然数ベクトル } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ リスト。リストの各要素は行列。 } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ popProb<-c(0.2,0.3,0.5) OffspringGenotypeProb(popProb,popProb) dad<-c(0.5,0.5,0) mom<-c(0,0,1) OffspringGenotypeProb(dad,mom) }
複数の探索者を同じ候補者プールで探すとき
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- MakeTestPattern
- タイトル
MakeTestPattern
- 説明
探索対象者は複数。それぞれの探索対象の候補者も複数。探索対象の候補者には重複がある。でも、ある候補者が複数の探索対象の解になることは許されない。そんなときの、探索パターンを網羅的に作成
- 使用例
MakeTestPattern(candidates)
- ソース
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),]) }
- Rdファイル
\name{MakeTestPattern} \alias{MakeTestPattern} %- Also NEED an '\alias' for EACH other topic documented here. \title{ MakeTestPattern } \description{ 探索対象者は複数。それぞれの探索対象の候補者も複数。探索対象の候補者には重複がある。でも、ある候補者が複数の探索対象の解になることは許されない。そんなときの、探索パターンを網羅的に作成 } \usage{ MakeTestPattern(candidates) } \arguments{ 自然数ベクトル } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ リスト。リストの各要素は行列。 } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ candidates<-NULL tobeSearched<-c(2,4) for(i in 1:length(tobeSearched)){ candidates[[i]]<-1:Npool } MakeTestPattern(candidates) }
複数のメンバーが探索対象のときには、その全員もしくは一部を探索対象とすることが必要。そのパターンを網羅する
- 予定パッケージ名
- HigashiNopponLikelihood
- 関数名
- MakeSearchPattern
- タイトル
MakeSearchPattern
- 説明
探索対象k人について、i人探索のパターン数(k,i)として、(k,1)+(k,2)+...+(k,k)通りの探索パターンを作成する
- 使用例
MakeSearchPattern(v)
- ソース
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]]<-matrix(tmp2,nrow=length(tmp[,1])) X[[i]]<-tmp2 x<-sum(x,length(as.matrix(X[[i]])[,1])) } list(X=X,x=x) }
- Rdファイル
\name{MakeSearchPattern} \alias{MakeSearchPattern} %- Also NEED an '\alias' for EACH other topic documented here. \title{ MakeSearchPattern } \description{ 探索対象k人について、i人探索のパターン数(k,i)として、(k,1)+(k,2)+...+(k,k)通りの探索パターンを作成する } \usage{ MakeSearchPattern(v) } \arguments{ 自然数ベクトル } %- maybe also 'usage' for other objects documented here. \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ リスト。リストの各要素は行列。 } \references{ %% ~put references to the literature/web site here ~ } \author{ %% ~~who you are~~ } \note{ %% ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ } \examples{ # needs kinship package v<-c(3,5) MakeSearchPattern(v) }