2010-02-22

ジェノタイプデータのＰＣＡ補正

R 集団構造化 Eigenstrat PCA

昨日の続き
PCAにより、ＧＷＡＳのジェノタイプデータでいくつかの軸情報で個人に「位置情報」が与えられる
個人の位置情報に応じて、個人のフェノタイプ(ケース・コントロール)とＳＮＰのジェノタイプの値を補正する
補正したフェノタイプと補正したジェノタイプの間の相関係数x(サンプル数-(軸数+1))を検定統計量とする
これは0/1の形質と0/1/2のジェノタイプに関するトレンド検定を補正したもの
すべてのサンプルが同一地点にあると仮定したのが、構造化補正をしないときのトレンド検定
サンプル数から(軸数+1)で補正するのは(軸数+1)の亜集団の存在を仮定した場合に想定するので、その分の自由度補正をしている(ようなもの)か。
以下のＲのソースは不明コールがないことを想定して計算している。実際のEigenstratでは、その部分のケアをしているために、ソースが少し長いか。
構造化集団のデータを作り、ケース・コントロールの割り当てに偏りを入れた上で、実施。
確かに、構造化補正できていて、Ｐ値が均一分布からのそれになっている模様。
Eigenstratとの答えの一致が確認されました。

#構造化集団をシミュレート
Nm<-1000 #マーカー数
Npop<-4 #亜集団数
Ns<-c(100,200,200,200) #集団別人数
M<-NULL #全ジェノタイプデータを納める行列
#亜集団別にアレル頻度を振ってシミュレーション
for(j in 1:Npop){
 tmpM<-matrix(rep(0,Nm*Ns[j]),nrow=Nm)

 for(i in 1:Nm){
  af<-runif(1)*0.8+0.1
  f<-rnorm(1,sd=0.01)
  if(abs(f)>1) f=0
  df<-c(af^2,2*af*(1-af),(1-af)^2)
  df[1]<-df[1]+f/2*df[2]
  df[3]<-df[3]+f/2*df[2]
  df[2]<-1-df[1]-df[3]
  tmpM[i,]<-sample(c(0,1,2),Ns[j],replace=TRUE,prob=df)
 }
 #全データ行列に格納
 M<-cbind(M,tmpM)
}


##PCA (Eigenstratの方式)##
#マーカー別平均
mu<-apply(M,1,mean)
#マーカー平均で標準化
M<-M-mu
#分散で標準化
M<-M/sqrt(mu/2*(1-mu/2))
#個人間の分散・共分散行列
X<-1/Nm*t(M)%*%M
#固有値分解
eiout<-eigen(X)
#固有値をプロット
plot(eiout$values)
#意味のある固有値の数は、亜集団の数-1
#PLoS "Population Structure and Eigenanalysis by Nick Patterson"にも記述がある通り
#plotしてみる
eivect<-as.data.frame(eiout$vectors)
eilist<-1:(Npop+1)
plot(eivect[,eilist])

########################

#偏らせて形質を与える
phenotype<-c(sample(c(0,1),sum(Ns)/2,replace=TRUE,prob=c(0.45,0.55)),sample(c(0,1),sum(Ns)/2,replace=TRUE,prob=c(0.55,0.45)))
L<-3 # 考慮する軸数
#検定統計量とｐ値の格納用
Chisq<-rep(0,Nm)
CorrChisq<-rep(0,Nm)
Ps<-rep(0,Nm)
CorrPs<-rep(0,Nm)

#考慮する軸数分をeigenvector行列から抜き出し
Emat<-eiout$vectors[,1:L]
#抜けデータがなければ、これはすべて値が1
Esqs<-apply(Emat*Emat,2,sum)
#phenotypeの平均値補正
phenotype<-phenotype-mean(phenotype)
#補正項
Gamma<-apply(Emat*phenotype,2,sum)/Esqs
#補正形質
corrphenotype<-phenotype-Emat%*%Gamma

#マーカーごとにループ
for(i in 1:Nm){
 genotype<-M[i,]
#マーカーことの補正項
 Gamma<-apply(Emat*genotype,2,sum)/Esqs
#マーカー別の補正ジェノタイプ
 corrgenotype<-genotype-Emat%*%Gamma
#補正なしのトレンド検定統計量は、ジェノタイプ-フェノタイプのcor^2のサンプル数倍
 Chisq[i]<-(sum(Ns))*cor(genotype,phenotype)^2
#補正ありのそれは、補正後ジェノタイプ-補正後フェノタイプのcor^2の(サンプル−亜集団数)倍
 CorrChisq[i]<-(sum(Ns)-(L+1))*cor(corrgenotype,corrphenotype)^2
}
#自由度１でｐ値算出
Ps<-pchisq(Chisq,1,lower.tail=FALSE)
CorrPs<-pchisq(CorrChisq,1,lower.tail=FALSE)
#plot
ylim<-c(0,1)
plot(ppoints(length(Ps),a=0),sort(Ps),ylim=ylim)
par(new=T)
plot(ppoints(length(Ps),a=0),sort(CorrPs),ylim=ylim,col="red")


plot(Ps,CorrPs)