罹患同胞対解析のIBDのずれ

• 罹患同胞対を用いた遺伝因子解析では、ノンパラメトリック連鎖解析手法が用いられます
• ノンパラメトリック連鎖解析とは、パラメトリック連鎖解析との対比での呼称です
  • パラメトリック連鎖解析では、遺伝形式(優性・劣性など)や浸透率などを変数として与えた上で、家系データ(家系のつながりとジェノタイプ・フェノタイプ)について、尤度を最大にするゲノム上の位置を探します
• ノンパラでは、原因座位でのジェノタイプ別のリスクを変数化します
• その上で、アレルの伝達パターンそのものではなくて、アレルの伝達パターンをIBD数でグループ分けし、IBD数の３パターンについて、尤度を最大にする位置を探します
• 罹患同胞対があったとき、そして、親のジェノタイプがわからないときに、IBD数の割合が、帰無仮説での比である、0.25,0.5,0.25からどのようにずれるかを考えます
• 両親がリスクアレルを何本持っているかで変わります
  • 両親がリスクアレルを２本持っているときには、それが片親に偏っているか、両親が１本ずつ持つかで変わります
• AA,Aa,aaの３ジェノタイプでの発病リスクをとします
• 両親のリスクアレルが0本の場合と4本の場合は、子の遺伝子型はaaのみかAAのみなので、IBD数の割合は0.25,0.5,0.25です
• 両親のリスクアレルが0本で、片親がAA,もう片親がaaのとき、子はつねにAaで、この場合も0.25,0.5,0.25です
• 両親がA,a,a,aの場合
  • IBD数１の尤度は0.5
  • IBD数２の尤度は
  • IBD数０の尤度は
• 両親がa,A,A,Aの場合
  • IBD数１の尤度は0.5
  • IBD数２の尤度は
  • IBD数０の尤度は
• 両親のリスクアレルが0本で、両親ともAaの場合
  • IBD数１の尤度は
  • IBD数２の尤度は
  • IBD数０の尤度は
• 以上のようになるので、子どものジェノタイプにバリエーションが出るにしたがって、表現型がそろった同胞対を集めると、IBDが大きい方が、尤度が高くなります。
• 確かめのRソース

##########使用関数
#Parentsのアレル
#伝達行列も
ParentsGen<-function(n=2){
 N<-n^4
 d<-matrix(0,N,4)
 for(i in 2:N){
  d[i,]<-d[i-1,]
  d[i,1]<-d[i,1]+1
  for(j in 1:(4-1)){
   if(d[i,j]==n){
    d[i,j]<-0
    d[i,j+1]<-d[i,j+1]+1
   }else{
    j<-n
   }
  }
 }
 d
}



IBDsCalc<-function(parents,parentsProb,trvec,GRisk){ 
 PairsIBD<-matrix(0,3,3)
 PairsIBDs<-list(a0=PairsIBD,a1=PairsIBD,a2HoHo=PairsIBD,a2HeHe=PairsIBD,a3=PairsIBD,a4=PairsIBD,aAll=PairsIBD)

 for(i in 1:length(parents[,1])){
  for(j in 1:length(trvec[,1])){
   #S1とS2の父・母由来染色体が祖父由来か祖母由来
   father1<-trvec[j,1]
   mother1<-trvec[j,2]
   father2<-trvec[j,3]
   mother2<-trvec[j,4]
   #祖父・祖母由来に照らして、Aかaかを確認
   g1f<-parents[i,father1+1]
   g1m<-parents[i,mother1+3]
   g2f<-parents[i,father2+1]
   g2m<-parents[i,mother2+3]
   #IBDの数は、father1,father2の一致、mother1,mother2の一致、のその和
   ibd<-0
   if(father1==father2){
    ibd<-ibd+1
   }
   if(mother1==mother2){
    ibd<-ibd+1
   }
   #S1,S2のジェノタイプは伝達アレルの和
   gtype1<-g1f+g1m
   gtype2<-g2f+g2m
   #S1,S2のジェノタイプに照らして、2人が罹患対か、片方罹患か、2人とも非罹患かの確率を計算
   AfSib<-GRisk[gtype1+1]*GRisk[gtype2+1]
   NafSib<-(1-GRisk[gtype1+1])*(1-GRisk[gtype2+1])
   AfNafSib<-1-AfSib-NafSib
   numRiskAllele<-4-sum(parents[i,])
   if(numRiskAllele==0){
    PairsIBDs$a0[1,ibd+1]<-PairsIBDs$a0[1,ibd+1]+parentsProb[i]*AfSib/16
    PairsIBDs$a0[2,ibd+1]<-PairsIBDs$a0[2,ibd+1]+parentsProb[i]*AfNafSib/16
    PairsIBDs$a0[3,ibd+1]<-PairsIBDs$a0[3,ibd+1]+parentsProb[i]*NafSib/16
   }else if(numRiskAllele==1){
    PairsIBDs$a1[1,ibd+1]<-PairsIBDs$a1[1,ibd+1]+parentsProb[i]*AfSib/16
    PairsIBDs$a1[2,ibd+1]<-PairsIBDs$a1[2,ibd+1]+parentsProb[i]*AfNafSib/16
    PairsIBDs$a1[3,ibd+1]<-PairsIBDs$a1[3,ibd+1]+parentsProb[i]*NafSib/16
   }else if(numRiskAllele==2){
    if(parents[i,1]==parents[i,2]){ #父の２アレルが同じ：リスクアレルが片親に偏っている
     PairsIBDs$a2HoHo[1,ibd+1]<-PairsIBDs$a2HoHo[1,ibd+1]+parentsProb[i]*AfSib/16
     PairsIBDs$a2HoHo[2,ibd+1]<-PairsIBDs$a2HoHo[2,ibd+1]+parentsProb[i]*AfNafSib/16
     PairsIBDs$a2HoHo[3,ibd+1]<-PairsIBDs$a2HoHo[3,ibd+1]+parentsProb[i]*NafSib/16
    }else{# ２親ともヘテロ
     PairsIBDs$a2HeHe[1,ibd+1]<-PairsIBDs$a2HeHe[1,ibd+1]+parentsProb[i]*AfSib/16
     PairsIBDs$a2HeHe[2,ibd+1]<-PairsIBDs$a2HeHe[2,ibd+1]+parentsProb[i]*AfNafSib/16
     PairsIBDs$a2HeHe[3,ibd+1]<-PairsIBDs$a2HeHe[3,ibd+1]+parentsProb[i]*NafSib/16
    }

   }else if(numRiskAllele==3){
    PairsIBDs$a3[1,ibd+1]<-PairsIBDs$a3[1,ibd+1]+parentsProb[i]*AfSib/16
    PairsIBDs$a3[2,ibd+1]<-PairsIBDs$a3[2,ibd+1]+parentsProb[i]*AfNafSib/16
    PairsIBDs$a3[3,ibd+1]<-PairsIBDs$a3[3,ibd+1]+parentsProb[i]*NafSib/16
   }else if(numRiskAllele==4){
    PairsIBDs$a4[1,ibd+1]<-PairsIBDs$a4[1,ibd+1]+parentsProb[i]*AfSib/16
    PairsIBDs$a4[2,ibd+1]<-PairsIBDs$a4[2,ibd+1]+parentsProb[i]*AfNafSib/16
    PairsIBDs$a4[3,ibd+1]<-PairsIBDs$a4[3,ibd+1]+parentsProb[i]*NafSib/16
   }
  }
 }
 PairsIBDs$aAll<-PairsIBDs$a0+PairsIBDs$a1+PairsIBDs$a2HoHo+PairsIBDs$a2HeHe+PairsIBDs$a3+PairsIBDs$a4
 PairsIBDs$a0<-PairsIBDs$a0/apply(PairsIBDs$a0,1,sum)
 PairsIBDs$a1<-PairsIBDs$a1/apply(PairsIBDs$a1,1,sum)
 PairsIBDs$a2HoHo<-PairsIBDs$a2HoHo/apply(PairsIBDs$a2HoHo,1,sum)
 PairsIBDs$a2HeHe<-PairsIBDs$a2HeHe/apply(PairsIBDs$a2HeHe,1,sum)
 PairsIBDs$a3<-PairsIBDs$a3/apply(PairsIBDs$a3,1,sum)
 PairsIBDs$a4<-PairsIBDs$a4/apply(PairsIBDs$a4,1,sum)
 PairsIBDs$aAll<-PairsIBDs$aAll/apply(PairsIBDs$aAll,1,sum)

 PairsIBDs
}

makeSNPrisk<-function(r=2,dom=0.5){
 DRisk<-c(r,0,1)
 DRisk[2]<-dom*DRisk[1]+(1-dom)*DRisk[3]
 DRisk
}
#集団のジェノタイプ別人数(ケース・コントロール別)を作成

makePop<-function(prev,DFreq,DRisk){
 CaseMat<-DFreq*DRisk
 CaseMatRatio<-CaseMat/sum(CaseMat)
 CaseMat<-CaseMatRatio*prev

 ContMat<-DFreq-CaseMat
 ContMatRatio<-ContMat/sum(ContMat)
 return(list(Case=CaseMat,CaseRatio=CaseMatRatio,Cont=ContMat,ContRatio=ContMatRatio))
}

##########実行

#条件
n<-2
Risk<-2 #リスクホモのRR
dom<-0.3 #優性効果
prevalence<-0.01 #有病率
AFreq<-0.1

#ジェノタイプ別リスク
DRisk<-makeSNPrisk(Risk,dom)

#伝達行列を作成
trvec<-ParentsGen(n=n)

#両親のアレル保有パターンを作成
parents<-ParentsGen(n=n)

#両親のジェノタイプ別頻度(集団内での)を作成

parentsProb<-rep(1,length(parents[,1]))

for(i in 1:length(parentsProb)){
 for(j in 1:length(parents[i,])){
  tmp<-AFreq
  if(parents[i,j]==1){
   tmp<-1-AFreq
  }
  parentsProb[i]<-parentsProb[i]*tmp
 }
}

#ジェノタイプ別
DFreq<-c(AFreq^2,2*AFreq*(1-AFreq),(1-AFreq)^2)
#集団のジェノタイプ別人数(ケース・コントロール別)を作成

pop<-makePop(prevalence,DFreq,DRisk)

#集団のジェノタイプ別発病率
GRisk<-pop$Case/DFreq

#両親のアレル保有パターン別・同胞の表現型共有場合別、IBD数比率
IBDs<-IBDsCalc(parents,parentsProb,trvec,GRisk)

IBDs