ある家系図であるマーカーのジェノタイプが観測される確率を計算する

書きかけ

  • 予定パッケージ名
    • 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]])
}