高次の因子

• 昨日、LD領域における、SNP選択の話があった(こちら:まだ、記事が掲載されていないけれど…)。Lossoとベイズとを比較して、ベイズの方が、モデル適合性が高そうだ、という話だった(あくまでも紹介論文でのシミュレーションデータの話)
• LDの具合や、集団構造化の具合とか、入り組んだ、高次LDの影響を受けた結果のように思えるので、そのあたり、色々なLDを比較的思い通りに作りたい
• LDの評価はペアワイズLDで行うことが多い
• ペアワイズLDが共通でありながら、より高次のマーカーの組み合わせでのLDが異なるようなLD状態を作ろう
• 多次元分割表で、尺度をSNPとすれば、SNP１個に関する周辺度数を同じくするハプロタイプ頻度は、単SNPアレル頻度が共通
• SNPペアに関する周辺度数を同じくするハプロタイプ頻度は、SNPペアのハプロタイプ頻度が共通
• SNP n 個の組み合わせに関する周辺度数を同じくするハプロタイプ頻度は、SNP n 個組み合わせのハプロタイプ頻度が共通

CategoryVector<-function (d = 3) 
{
    if(d == 1){
     return(matrix(1,ncol=1,nrow=1))
    }
    df <- d - 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}
Simplex.Vector <- function(n){
	if(n==1){
		return(matrix(1,ncol=1,nrow=1))
	}
	X <- CategoryVector(n)
	X <- cbind(rep(1/sqrt(n),n),sqrt((n-1)/n)*X)
	X
}
Make.X <- function(rs){
	ret <- matrix(1,ncol = 1, nrow =1)
	for(i in 1:length(rs)){
		ret <- Simplex.Vector(rs[i]) %x% ret
	}
	t(ret)
}

make.snp.cmb <- function(k){
	snp.cmb <- matrix(0,nrow=1,ncol=k)

	for(i in 1:k){
		tmp <- rep(0,k)
		tmp[i] <- 1
		tobe.added <- t(t(snp.cmb) + tmp)
		snp.cmb <- rbind(snp.cmb, tobe.added)
	}
	snp.cmb
	
}

r.hap.freq <- function(hf,k){
	n <- log(length(hf),2)
	rs <- rep(2,n)
	X <- Make.X(rs)
	snp.cmb <- make.snp.cmb(n)
	rot.hf <- X %*% hf

	zeros <- which(apply(snp.cmb,1,sum) <=k)
	non.zeros <- which(apply(snp.cmb,1,sum) >=k+1)

	new.rot.hf <- rot.hf
	#jit.sign <- sample(c(-1,1),length(non.zeros),replace=TRUE)
	jit.value <- c(rdirichlet(1,rep(1,length(non.zeros))))
	#jit <- jit.sign * jit.value
	loop <- TRUE
	cnt <- 0
	while(loop){
		jit <- jit.value/(2^cnt) * sample(c(-1,1),length(non.zeros),replace=TRUE)
		new.rot.hf[non.zeros] <- rot.hf[non.zeros] + jit/(2^cnt)
		cnt <- cnt+1
		new.hf <- t(X) %*% new.rot.hf
		if(min(new.hf)>=0){
			loop <- FALSE
		}
	}
	new.hf

}

k <- 10
n.zero.h <- 10
n.non.zero.h <- 2^k-n.zero.h
f <- 0.999
tmp.hf.zero <- rdirichlet(1,rep(1,n.zero.h))*(1-f)
tmp.hf.non.zero <- rdirichlet(1,rep(1,n.non.zero.h))*f
hf <- sample(c(tmp.hf.zero,tmp.hf.non.zero))
#hf[c(1,sample(2:(length(hf)-1),n.non.zero.h-2),2^k)] <- tmp.hf
#hf

n.iter <- 10

hf.array <- array(hf,rep(2,k))
for(j in 1:n.iter){
	new.hf <- r.hap.freq(hf,3)
	new.hf.array <- array(new.hf,rep(2,k))
	print(snp.cmb[1,])
	print("-old---")
	print(sum(hf.array))
	print("-new---")
	print(sum(new.hf.array))

	for(i in 2:length(snp.cmb[,1])){
		print("=====")
		print(snp.cmb[i,])
		print("-old---")
		print(apply(hf.array,which(snp.cmb[i,]==1),sum))
		print("-new---")
		print(apply(new.hf.array,which(snp.cmb[i,]==1),sum))
		print("-----")
	}

}