ジェノタイプデータのＰＣＡその２

昨日の続き

では、正方行列を作らずに、非正方行列のままsvd()をかけるとどうなるかもやってみます。
同じ構造を表す固有値と固有ベクトルが取れました。

#構造化集団をシミュレート
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)
#マーカー平均で標準化
M2<-M-mu
#分散で標準化
M2<-M2/sqrt(mu/2*(1-mu/2))
#個人間の分散・共分散行列
X<-1/Nm*t(M2)%*%M2
#固有値分解
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)
plot(eivect[,eilist])
#############################
##標準化したM2をsvdしてみる##
#############################
svdout2<-svd(M2)
#固有値をプロット
plot(svdout2$d)
#意味のある固有値の数は、亜集団の数-1

plot(eiout$values,svdout2$d)
#eigenの固有値と同じ

svdvect2<-as.data.frame(svdout2$v)

plot(svdvect2[,eilist])
################################
#2つのeigenevectorの関係を見る#
################################

comparison<-as.data.frame(cbind(eivect[,eilist],svdvect2[,eilist]))
plot(comparison)