Nothing
### This file contains codes modified from
### "tzeng-main.function.forRDscoreTest.public.r".
### Reference: J.-Y. Tzeng (2005) Genetic Epidemiology.
### J.-Y. Tzeng, C.-H. Wang, J.-T. K, and C.-H. K. Hsiao (2006) Am. J. Hum. Genet.
### Modified: Wei-Chen Chen 2009-09-17
### X: array[N, L], N haplotypes, L loci.
### ploidy: 1 or 2 => index of X = rep(1:(N/ploidy), each = ploidy)
### This function returns a fake object "haplo" to replace the output of haplo.em().
haplo.phase.known <- function(X, ploidy = 2){
X.names <- apply(X, 1, paste, collapse="")
unique.X <- unique(X)
unique.X.names <- apply(unique.X, 1, paste, collapse="")
X.summary <- as.data.frame(table(X.names))
X.summary$prob <- X.summary$Freq / sum(X.summary$Freq)
haplo <- NULL
haplo$haplotype <- unique.X
haplo$hap.prob <- X.summary$prob[match(unique.X.names, as.character(X.summary$X.names))]
haplo$post <- rep(1, nrow(X))
if(ploidy == 1){
haplo$hap1code <- match(X.names, unique.X.names)
haplo$hap2code <- haplo$hap1code
haplo$indx.subj <- 1:nrow(X)
} else if(ploidy == 2){
tmp <- matrix(match(X.names, unique.X.names), nrow = 2)
haplo$hap1code <- tmp[1,]
haplo$hap2code <- tmp[2,]
haplo$indx.subj <- rep(1:(nrow(X)/2), each = 2)
} else{
stop("ploidy = 1 or 2")
}
haplo
} # End of haplo.phase.know().
### This function returns posterior probabilities based on the fake "haplo" object.
haplo.post.prob <- function(X, ploidy = 2, skip.haplo = 1e-7, K = NULL){
### Add pos for arbitrarily choosing K. Use default getbut.fun() if K = NULL.
### Modified: Wei-Chen Chen 2009-12-05.
haplo <- haplo.phase.known(X, ploidy) # The fake "haplo" object.
nhap <- nrow(X)
n.loci <- ncol(X)
n.subj <- nrow(X)
digit <- 1
H.names<-apply(haplo$haplotype, 1, paste, collapse="")
Pcscn <-haplo$hap.prob; names(Pcscn)<-H.names
pcscn <-rev(sort(Pcscn))
if(is.null(K)){
#WCC pos <-getcut.fun(pcscn,nhap, plot=0)
pos <-getcut.fun(pcscn, 2*nhap, plot=0)
} else{
tl.K <- length(pcscn)
if(K < 1 || K > (tl.K - 1)){
stop(paste("K should be in [", 1, ", ", tl.K - 1, "]", sep = ""))
}
pos <- K
}
reserv <-names(pcscn)[1:pos]
##---- prepare to get X, X*, E(X*|G), E(X|G) ----
# n.subj <- length(y)
hap1 <- haplo$hap1code
hap2 <- haplo$hap2code
indx <- haplo$indx.subj
post <- haplo$post
uhap <- sort(unique(c(hap1, hap2)))
which.haplo<- haplo$hap.prob >= skip.haplo
p.cscn <- Pcscn[which.haplo]
h.names <- names(p.cscn)
uhap <- uhap[which.haplo]
##---- obtain bigMatB for Reducing Dim ----
BigMatB <- final.BigMatB.matC.fun(p.cscn, reserv, n.loci, digit) ##"BB.cscn""CC.cscn""BigMatB.cscn""subPI.cscn"
matB <- BigMatB$BigMatB.cscn
# CC <- BigMatB$CC.cscn
# BB <- BigMatB$BB.cscn
# subPI <- BigMatB$subPI.cscn
# p.cscn.FD <- p.cscn[match(rownames(matB), h.names)]
# p.cscn.RD <- p.cscn.FD %*% matB
# RR <- nrow(matB)
# RRstar <- ncol(matB)
##---- get E(X*|G) = xRD.post ----
x <- 0.5* (outer(hap1, uhap, "==") + outer(hap2, uhap, "==")); colnames(x)<-h.names
x <- x[,match(rownames(matB), h.names)]
## the "0.5" is because in Schaid's coding the rowsum of X vec is 2 while mine is 1
## xRD <- x %*% matB[,-1] ## modified on May 17, 2005
xRD <- x %*% matB
n.xRD <- ncol(xRD)
xRD.post <- matrix(rep(NA, n.subj * n.xRD), ncol = n.xRD);for(j in 1:n.xRD){ xRD.post[,j]<-tapply(xRD[,j]*post,indx,sum)}
##---- get E(X|G) = x.post ----
n.x <- ncol(x)
x.post <- matrix(rep(NA, n.subj * n.x), ncol = n.x);for(j in 1:n.x){ x.post[,j]<-tapply(x[,j]*post,indx,sum)};
list(haplo = haplo,
FD.id = match(rownames(matB), H.names),
RD.id = match(colnames(matB), H.names),
FD.post = x.post,
RD.post = xRD.post,
g.truncate = ncol(matB))
} # End of haplo.post.prob()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.