ジェノタイプデータのＰＣＡ eigen()

集団構造化があるときに、PCAして、プロットすることがある。
その情報を使って、形質マッピング検定に用いる前座のようなもの。
この論文がEigenstratのそれですが。
ここでやっている、PCA部分をRでなぞってみます

#構造化集団をシミュレート
Nm<-1000 #マーカー数
Npop<-4 #亜集団数
Ns<-c(100,150,200,250) #集団別人数
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])