R/functions-crlmm.R
In benilton/oligo-release: Preprocessing tools for oligonucleotide arrays.
### WARNING: IF ROBUST TO DIMNAMES,
### REMOVE AND SAVE LOTS
### OF MEMORY... ADD DIMNAMES
### AT THE END
rowEntropy <- function(p) rowMeans(rowSums(log2(p^p), dims=2))
getSnpFragmentLength <- function(object, enzyme){
## returns data.frame with 3 columns:
## SNP | Enzyme | Length
conn <- db(get(annotation(object)))
snps <- paste("(", paste("'", featureNames(object), "'", collapse=", ", sep=""), ")")
sql <- paste("SELECT man_fsetid, enzyme, length",
"FROM fragmentLength",
"INNER JOIN featureSet USING(fsetid)",
"WHERE man_fsetid IN", snps)
if (!missing(enzyme)){
stopifnot(length(enzyme) == 1)
sql <- paste(sql, "AND enzyme =", enzyme)
}
dbGetQuery(conn, sql)
}
snpGenderCall <- function(object){
XIndex <- getChrXIndex(object)
theA <- getA(object)
if (!is.ff(theA)){
tmp <- median(theA, na.rm=TRUE)
}else{
nc <- ncol(object)
tmp <- rep(NA, nc)
for (i in 1:nc) tmp[i] <- median(getA(object[,i]), na.rm=TRUE)
tmp <- median(tmp)
}
if (bothStrands(object)){
a <- apply(theA[XIndex,,], 2, median, na.rm=TRUE)
}else{
a <- apply(theA[XIndex,], 2, median, na.rm=TRUE)
}
kfit=kmeans(a, c(min(a),tmp))
return(factor(c("female","male")[as.numeric(kfit$cluster==1)+1]))
}
getChrXIndex <- function(object){
sql <- "SELECT man_fsetid, chrom, physical_pos FROM featureSet WHERE man_fsetid LIKE 'SNP%'"
chrs <- dbGetQuery(db(get(annotation(object))), sql)
chrs <- chrs[order(chrs[["chrom"]], chrs[["physical_pos"]], chrs[["man_fsetid"]], na.last=FALSE),]
return(which(chrs[["chrom"]] == "X"))
}
##gender in pData keeps male female
fitAffySnpMixture <- function(object, df1=3, df2=5,
probs=rep(1/3,3), eps=50,
subSampleSize=10^4, verbose=TRUE){
if(is.null(object$gender)){
maleIndex <- snpGenderCall(object)=="male"
}else{
maleIndex <- object$gender=="male"
}
XIndex=getChrXIndex(object)
I <- dim(object)[1]
J <- dim(object)[2]
set.seed(1)
## tmp <- c( (1:I)[-XIndex],((I+1):(2*I))[-XIndex])
## tmp <- c( (1:I),((I+1):(2*I)))
if (I > subSampleSize){
idx <- sort(sample(I, subSampleSize))
}else{
idx <- 1:I
}
## rm(tmp)
if (bothStrands(object)){
pis <- array(0,dim=c(I,J,3,2))
fs <- array(0,dim=c(I,J,2))
snr <- array(0,dim=J)
}else{
pis <- array(0,dim=c(I,J,3))
fs <- array(0,dim=c(I,J))
snr <- array(0,dim=J)
}
##### dimnames(fs)<-list(featureNames(object),
##### sampleNames(object),
##### c("antisense","sense"))
##### dimnames(pis)<-list(featureNames(object),
##### sampleNames(object),
##### c("AA","AB","BB"),
##### c("antisense","sense"))
##### names(snr) <- sampleNames(object)
if(verbose) cat("Fitting mixture model to ", J, " arrays. Epsilon must reach ", eps, ".\n",sep="")
## getSnpFragmentLength now returns a table SNP | enzyme | length
## and it is possible that there are 2+ lengths by SNPxEnzyme combo
## 21 Oct 2011
L <- getSnpFragmentLength(object)
L <- aggregate(L$length, by=list(SNP=L$man_fsetid), mean)
L <- L[match(featureNames(object), L$SNP), 2]
fix <- which(is.na(L))
L[fix] <- median(L, na.rm=T)
rm(fix)
if (bothStrands(object)){
L <- cbind(L,L)
l <- as.numeric(L[idx, ])
L <- as.numeric(L)
}else{
l <- L[idx]
}
L <- L-mean(L)
l <- l-mean(l)
matL <- ns(l, df1)
for(j in 1:J){
Y <- getM(object[,j])
A <- getA(object[,j])
if (bothStrands(object)){
Y <- Y[,1,]
A <- A[,1,]
}
fix <- which(is.na(Y))
Y[fix] <- median(Y, na.rm=T)
A[fix] <- median(A, na.rm=T)
rm(fix)
mus <- quantile(Y,c(1,3,5)/6);mus[2]=0
sigmas <- rep(mad(c(Y[Y<mus[1]]-mus[1],Y[Y>mus[3]]-mus[3])),3)
sigmas[2] <- sigmas[2]/2
if (bothStrands(object)){
a <- as.numeric(A[idx, ])
y <- as.numeric(Y[idx, ])
}else{
a <- A[idx]
y <- Y[idx]
}
A <- A-mean(A)
a <- a-mean(a)
weights <- apply(cbind(mus, sigmas), 1, function(p) dnorm(y, p[1], p[2]))
PreviousLogLik <- -Inf
change <- eps+1
itmax <- 0
matA <- ns(a,df2)
while (change > eps & itmax < 100){
itmax <- itmax+1
## E
z <- sweep(weights, 2, probs, "*")
LogLik <- rowSums(z)
z <- sweep(z, 1, LogLik, "/")
LogLik <- sum(log(LogLik))
change <- abs(LogLik-PreviousLogLik)
if (verbose){
if (itmax > 1 | j > 1) cat(del)
message <- paste("Array ",j,": epsilon = ", signif(change,2), " ", sep="")
del <- paste(rep("\b", nchar(message)), collapse="")
cat(message)
}
PreviousLogLik <- LogLik
probs <- colMeans(z)
## M
fit1 <- lm(y~matL+matA,weights=z[,1])
fit2 <- sum(z[,2]*y)/sum(z[,2])
fit3 <- lm(y~matL+matA,weights=z[,3])
sigmas[1] <- sqrt(sum(z[,1]*residuals(fit1)^2)/sum(z[,1]))
sigmas[2] <- sqrt(sum(z[,2]*(y-fit2)^2)/sum(z[,2]))
sigmas[3] <- sqrt(sum(z[,3]*residuals(fit3)^2)/sum(z[,3]))
weights[,1] <- dnorm(y,fitted(fit1),sigmas[1])
weights[,2] <- dnorm(y,fit2,sigmas[2])
weights[,3] <- dnorm(y,fitted(fit3),sigmas[3])
weights[y >= 0, 1] <- 0
weights[y <= 0, 3] <- 0
}
## gc()
bigX <- cbind(1, ns(L, knots=as.numeric(attr(matL, "knots")), Boundary.knots=attr(matL, "Boundary.knots")),
ns(A, knots=as.numeric(attr(matA, "knots")), Boundary.knots=attr(matA, "Boundary.knots")))
rm(matA); ## gc()
pred1 <- bigX%*%coef(fit1)
pred2 <- rep(fit2,length(Y))
pred3 <- bigX%*%coef(fit3)
rm(bigX); ## gc()
weights <- matrix(0,length(Y),3)
weights[,1] <- dnorm(Y,pred1,sigmas[1])
weights[,2] <- dnorm(Y,pred2,sigmas[2])
weights[,3] <- dnorm(Y,pred3,sigmas[3])
weights[Y >= pred2, 1] <- 0
## if (maleIndex[j]) weights[XIndex, 2] <- 0
weights[Y <= pred2, 3] <- 0
z <- sweep(weights, 2, probs, "*")
LogLik <- rowSums(z)
z <- sweep(z, 1, LogLik, "/")
if (bothStrands(object)){
fs[,j,] <- matrix((pred3-pred1)/2,ncol=2)
for(k in 1:3){
pis[,j,k,] <- matrix(z[,(4-k)],ncol=2) ##4-k cause 3is1,2is2 and 1is3
}
snr[j] <- median(fs[,j,])^2/(sigmas[1]^2+sigmas[2]^2)
}else{
fs[, j] <- (pred3-pred1)/2
for(k in 1:3){
pis[, j, k] <- z[,(4-k)] ##4-k cause 3is1,2is2 and 1is3
}
snr[j] <- median(fs[,j])^2/(sigmas[1]^2+sigmas[2]^2)
}
}
fs[fs < 0] <- 0
if(verbose) cat("Done.\n")
return(list(f0=median(fs),fs=fs, pis=pis, snr=snr))
}
getInitialAffySnpCalls <- function(object,subset=NULL,
concordanceCutoff=0.0001,
cutoffs=c(0.7,0.5,0.7),
returnProbs=FALSE,
verbose=FALSE,
sqsClass = "SnpQSet"){
if(is.null(subset)) subset <- 1:(dim(object$pis)[1])
if (sqsClass == "SnpQSet"){
pi1 <- object$pis[subset,,,1]
pi2 <- object$pis[subset,,,2]; rm(object); ## gc()
## fixing SNPs that have probes only on one strand:
idx <- is.na(pi1[,1,1])
pi1[idx,,] <- pi2[idx,,]
idx <- is.na(pi2[,1,1])
pi2[idx,,] <- pi1[idx,,]
rm(idx); ## gc()
if(verbose) cat("Picking good starting value: ")
if(verbose) cat("Computing entropy, ")
E1 <- rowEntropy(pi1)
E2 <- rowEntropy(pi2); ## gc()
tmpN <- dim(pi1)
tmpcall1 <- tmpcall2 <- matrix(NA, nrow=tmpN[1], ncol=tmpN[2])
##### rownames(tmpcall1) <- rownames(tmpcall2) <- dimnames(pi1)[[1]]
##### colnames(tmpcall1) <- colnames(tmpcall2) <- dimnames(pi1)[[2]]
## gc()
if(verbose) cat("calculating calls, ")
for (i in 1:tmpN[1]){
for (j in 1:tmpN[2]){
tmpcall1[i, j] <- which.max(pi1[i,j,])
tmpcall2[i, j] <- which.max(pi2[i,j,])
}
}
if(verbose) cat("determining non-concordant calls, ")
concordance <- rowIndepChiSqTest(tmpcall1,tmpcall2); ## gc()
if(verbose) cat("deciding which strand(s) to use")
tc1 <- tmpcall1 == 1
tc2 <- tmpcall1 == 2
tc3 <- tmpcall1 == 3
noABIndex1 <- which((rowSums(tc2)<3)*(rowSums(tc3)>0)*(rowSums(tc1)>0) == 1)
tc1 <- tmpcall2 == 1
tc2 <- tmpcall2 == 2
tc3 <- tmpcall2 == 3
noABIndex2 <- which((rowSums(tc2)<3)*(rowSums(tc3)>0)*(rowSums(tc1)>0) == 1)
rm(tc1, tc2, tc3); ## gc()
E1[noABIndex1] <- -Inf
E2[noABIndex2] <- -Inf
jointprobs <- (pi1+pi2)/2
## gc()
noInfoIndex <- intersect(noABIndex1, noABIndex2)
rm(noABIndex1, noABIndex2)
jointprobs[noInfoIndex,,] <- 1/3
rm(noInfoIndex)
notBoth <- concordance < concordanceCutoff
E1bE2 <- E1 > E2
rm(E1, E2)
i1 <- notBoth & E1bE2
i2 <- notBoth & !E1bE2
rm(notBoth, E1bE2)
jointprobs[i1,,] <- pi1[i1,,]
rm(i1, pi1)
jointprobs[i2,,] <- pi2[i2,,]
rm(i2, pi2); ## gc()
tmpN <- dim(jointprobs)
tmpcall <- tmpmax <- matrix(NA, nrow=tmpN[1], ncol=tmpN[2])
##### rownames(tmpcall) <- rownames(tmpmax) <- dimnames(jointprobs)[[1]]
##### colnames(tmpcall) <- colnames(tmpmax) <- dimnames(jointprobs)[[2]]
## gc()
if(verbose) cat(" finalizing"); ## gc()
for (i in 1:tmpN[1]){
for (j in 1:tmpN[2]){
tmpcall[i, j] <- which.max(jointprobs[i, j, ])
tmpmax[i, j] <- jointprobs[i, j, tmpcall[i,j]]
}
}
}else{
tmpcall <- apply(object$pis[subset,, ], 1:2, which.max)
tmpmax <- apply(object$pis[subset,, ], 1:2, max)
}
for(i in 1:3)
tmpcall[tmpcall==i & tmpmax<cutoffs[i]] <- NA
if(verbose) cat("\nCompleted!\n")
if(sqsClass == "SnpQSet"){
if(returnProbs) return(list(calls=tmpcall,probs=jointprobs)) else return(tmpcall)
}else{
return(tmpcall)
}
}
rowIndepChiSqTest <- function(call1,call2){
tmp=vector("numeric",nrow(call1))
for(i in 1:nrow(call1)){
tmpt= cbind(table(factor(call1[i,],levels=1:3)),
table(factor(call2[i,],levels=1:3)))
tmpt=tmpt[!rowSums(tmpt)==0,,drop=FALSE]
if(nrow(tmpt)>1){
rowtot <- rowSums(tmpt)
coltot <- colSums(tmpt)
e=outer(rowtot,coltot)/sum(tmpt)
stat = sum((tmpt-e)^2/e)
tmp[i] = 1 - pchisq(stat, (nrow(tmpt) - 1) * (ncol(tmpt) - 1))
}
else{
tmp[i]=0
}
}
return(tmp)
}
getGenotypeRegionParams <- function(M, initialcalls, f=0, verbose=TRUE, subset){
if (!missing(subset)){
M <- M[subset,, drop=FALSE]
initialcalls <- initialcalls[subset,, drop=FALSE]
f <- f[subset,, drop=FALSE]
}
if(verbose) cat("Computing centers and scales for 3 genotypes")
## tmp <- .Call("R_HuberMatrixRows2", M+(initialcalls-2)*f,
## as.integer(initialcalls), 1.5)
tmp <- .Call("R_trimmed_stats", M+(initialcalls-2)*f,
as.integer(initialcalls), 0.025)
##### centers <- scales <- N <- array(NA, dim=c(nrow(M),3))
##### scales <- N <- array(NA, dim=c(nrow(M),3))
##### dimnames(centers) <- dimnames(scales) <- dimnames(N) <- list(rownames(M), c("AA","AB","BB"))
##### centers[,] <- tmp[[1]]
tmpN <- tmp[[3]]; tmp[[3]] <- NULL
tmpN[tmpN <= 1] <- NA
tmpS <- tmp[[2]]; tmp[[2]] <- NULL
tmpS[is.na(tmpN)] <- NA
tmpV <- sqrt(rowSums((tmpN[,-2]-1)*tmpS[,-2]^2, na.rm=T)/rowSums(tmpN[,-2]-1, na.rm=T))
tmpV[!is.finite(tmpV)] <- NA
tmpS[,-2] <- tmpV
rm(tmpV)
tmpC <- tmp[[1]]; rm(tmp)
tmpC[is.na(tmpN)] <- NA
tmpN[is.na(tmpN)] <- 0
if(verbose) cat(" Done\n")
return(list(centers=tmpC, scales=tmpS, N=tmpN))
}
getAffySnpGenotypeRegionParams<-function(object,initialcalls,f=NULL,
subset=1:(dim(object)[1]),
verbose=FALSE, sqsClass = "SnpQSet"){
if(is.null(f)) f=fitAffySnpMixture(object,verbose=verbose)$fs
if(verbose) cat("Computing centers and scales:\n")
if (sqsClass == "SnpQSet"){
N <- scales <- centers <- array(NA,dim=c(length(subset),3,2))
##### dimnames(N) <- dimnames(scales) <- dimnames(centers) <- list(featureNames(object)[subset],
##### c("AA","AB","BB"), c("antisense","sense"))
for(s in 1:2){
## gc()
if(verbose) cat(c("\tantisense","\tsense....")[s],": ",sep="")
tmp <- getGenotypeRegionParams(getM(object)[subset,,s],
initialcalls[subset,],
f[subset,,s],
verbose=verbose)
centers[,,s] <- tmp$centers
scales[,,s] <- tmp$scales
N[,,s] <- tmp$N
}
N <- apply(N, 1:2, max, na.rm=TRUE) ## returns -Inf if c(NA, NA)
}else{
tmp <- getGenotypeRegionParams(getM(object)[subset,],
initialcalls[subset,],
f[subset,],
verbose=verbose)
centers <- tmp$centers
scales <- tmp$scales
N <- tmp$N
}
return(list(centers=centers,scales=scales,N=N,f0=median(f, na.rm=TRUE)))
}
getAffySnpPriors <- function(object,minN=20,subset=1:(dim(object$centers)[1]),
verbose=TRUE, sqsClass="SnpQSet"){
if(verbose) cat("Computing priors.\n")
if (sqsClass == "SnpQSet"){
N <- cbind(object$N,object$N)
Index <- subset[which(rowMeans(N[subset,]>minN)==1)]
N <- N[Index,]
mus <- cbind(object$centers[Index,,1],object$centers[Index,,2])
sigmas <- cbind(object$scales[Index,,1],object$scales[Index,,2])
sigmas[is.na(mus)] <- NA
maxsigmas <- c(quantile(sigmas[,c(1,3,4,6)],.99, na.rm=TRUE),
quantile(sigmas[,c(2,5)],.99, na.rm=TRUE))
}else{
N <- object$N
Index <- subset[which(rowMeans(N[subset,] > minN) == 1)]
N <- N[Index,]
mus <- object$centers[Index,]
sigmas <- object$scales[Index,]
maxsigmas <- c(quantile(sigmas[,c(1, 3)],.99, na.rm=TRUE),
quantile(sigmas[, 2],.99, na.rm=TRUE))
}
##variance for prior for mu
##need to make robust
V <- cov(mus,use="pairwise.complete")
##prior for sigma
zgs <- log(sigmas^2); rm(sigmas); ## gc()
dgs <- N-1
egs <- zgs-digamma(dgs/2)+log(dgs/2); rm(zgs); ## gc()
n <- length(Index)
d0s <- 2*trigammaInverse(colMeans(n*(egs-colMeans(egs, na.rm=TRUE))^2/(n-1)-trigamma(dgs/2), na.rm=TRUE))
s20 <- exp(colMeans(egs, na.rm=TRUE)+digamma(d0s/2)-log(d0s/2))
return(list(V=V,d0s=d0s,s20=s20,maxsigmas=maxsigmas))
}
##minN- minimum number of points in cluster required for use
##max$Sigma
updateAffySnpParams <- function(object, priors, missingStrandIndex, minN=3,
maxHomoSigma=priors$maxsigma[1],
maxHeteSigma=priors$maxsigma[2],
subset=1:(dim(object$centers)[1]),
d0s=80, verbose=FALSE){
object$centers <- object$centers[subset,,]
object$scales <- object$scales[subset,,]
object$N <- object$N[subset,]
missingStrandIndex <- missingStrandIndex[subset]
if(verbose) cat("Updating centers and scales")
##First variances
for(i in 1:2){
for(j in 1:2){ ##1 and 3 are the same
if(j==2) N <- object$N[,2] else N <- rowSums(object$N[,-2],na.rm=TRUE)
s <- object$scales[,j,i]
if (is.null(d0s))
d0s <- priors$d0s[3*(i-1)+j]
s20 <- priors$s20[3*(i-1)+j] ##notice the ad-hoc choice of 3
Index <- N>minN & !is.na(s)
N <- N[Index]; s <- s[Index]
object$scales[Index,j,i] <- sqrt ( ( (N-1)*s^2 + d0s*s20 ) / (d0s+N-1) )
object$scales[!Index & missingStrandIndex != i,j,i] <- sqrt(s20)
}
}
object$scales[,3,] <- object$scales[,1,] ##AA=BB
object$scales[,2,][object$scales[,2,]>maxHeteSigma] <- maxHeteSigma
object$scales[,-2,][object$scales[,-2,]>maxHomoSigma] <- maxHomoSigma
if(verbose) cat(".")
##Means
updateMean <- function(strandIndex, strand){
## strandIndex is a vector where
## 0 - none strands are missing on the array
## 1 - antisense strand is missing
## 2 - sense strand is missing
## strand argument is what do you want to be updated
if (strand == "both"){
snps <- strandIndex == 0; idx <- 1:6
N <- cbind(object$N[snps,],object$N[snps,])
mu <- cbind(object$centers[snps,,1],object$centers[snps,,2])
}else if(strand == "antisense"){
snps <- strandIndex == 2; idx <- 1:3
N <- object$N[snps,]
mu <- object$centers[snps,,1]
}else if(strand == "sense"){
snps <- strandIndex == 1; idx <- 4:6
N <- object$N[snps,]
mu <- object$centers[snps,,2]
}
if (sum(snps) == 1){
N <- matrix(N, nrow=1)
mu <- matrix(mu, nrow=1)
}
Vinv <- solve(priors$V[idx, idx])
NSinv <- t(N)/priors$s20[idx]
tmp <- t(sapply(1:nrow(mu),function(i){
if(verbose & i%%5000==0) cat(".")
mus=mu[i,]; Ns=N[i,]
mus[Ns<minN] <- 0
mus[is.na(mus)] <- 0
return(solve(Vinv+diag(NSinv[,i]))%*%(NSinv[,i]*mus))
}))
return(tmp)
}
if (any(missingStrandIndex == 0)){
tmp <- updateMean(missingStrandIndex, "both")
object$centers[missingStrandIndex == 0,,1] <- tmp[,1:3]
object$centers[missingStrandIndex == 0,,2] <- tmp[,4:6]; rm(tmp);
}
if (any(missingStrandIndex == 2))
object$centers[missingStrandIndex == 2,,1] <- updateMean(missingStrandIndex, "antisense")
if (any(missingStrandIndex == 1))
object$centers[missingStrandIndex == 1,,2] <- updateMean(missingStrandIndex, "sense")
if(verbose) cat("Done.\n")
return(object)
}
getAffySnpDistance <- function(object,params,f=0,subset=1:(dim(object)[1]),
w=NULL,verbose=FALSE){
x=getM(object[subset,])
rm(object)
Dist <- array(NA,dim=c(dim(x)[1],dim(x)[2],3,2))
if(verbose) cat("Calculating likelihood-based distances")
for(i in 1:2){
if(verbose) cat(".")
for(j in 1:3){
tmp <- x[,,i]+(j-2)*f[subset,,i]
Dist[,,j,i] <- 2*log(params$scales[subset,j,i]) +
((tmp-params$centers[subset,j,i])/params$scales[subset,j,i])^2
if(!is.null(w)) Dist[,,j,i] <- Dist[,,j,i] - 2*log(w[subset,,j,i])
}
}
##### dimnames(Dist) <- list(dimnames(x)[[1]],
##### dimnames(x)[[2]], 1:3,
##### dimnames(x)[[3]])
if(verbose) cat("Done.\n")
return(Dist)
}
getAffySnpCalls <- function(Dist,XIndex,maleIndex,subset=1:(dim(Dist)[1]),
verbose=FALSE, sqsClass="SnpQSet"){
if (sqsClass == "SnpQSet"){
Dist <- Dist[subset,,,,drop=FALSE]
}else{
Dist <- Dist[subset,,, drop=FALSE]
}
#### XIndex <- which(subset%in%XIndex)
## gc()
## Here we can even put default value to 3 for the new code. /HB
res <- array(as.integer(-1),dim=dim(Dist)[c(1,2)]); ## gc()
##### dimnames(res) <- dimnames(Dist)[1:2]
if (sqsClass == "SnpQSet")
Dist <- rowSums(Dist, na.rm=TRUE, dims=3)
## Dist[XIndex, maleIndex, 2] <- Inf
##the following is slow!
if(verbose) cat("Making calls for ", ncol(res), " arrays");
for(j in 1:ncol(res)){
if(verbose) cat(".");
D1 <- Dist[,j,1];
D2 <- Dist[,j,2];
D3 <- Dist[,j,3];
d12 <- (D1 < D2);
d23 <- (D2 < D3); rm(D2);
d13 <- (D1 < D3); rm(D3);
d <- rep(as.integer(3), length(D1)); rm(D1)
d[( d12 & d13)] <- as.integer(1); rm(d13)
d[(!d12 & d23)] <- as.integer(2); rm(d12, d23)
res[,j] <- d;
rm(d);
}
if(verbose) cat("Done\n")
return(res)
}
getAffySnpConfidence <- function(Dist, Calls, XIndex, maleIndex,
subset=1:nrow(Calls), verbose=TRUE, sqsClass = "SnpQSet"){
if (sqsClass == "SnpQSet"){
Dist <- rowSums(Dist[subset,,,,drop=FALSE], dims=3, na.rm=T)
}else{
Dist <- Dist[subset,,, drop=FALSE]
}
Calls <- Calls[subset,,drop=FALSE]
#### XIndex <- which(subset%in%XIndex)
res <- array(NA,dim=dim(Dist)[c(1,2)])
##### dimnames(res) <- list(dimnames(Dist)[[1]],dimnames(Dist)[[2]])
##### Dist <- rowSums(Dist, dims=3, na.rm=T)
if (verbose) cat("Computing confidence for calls on ",ncol(res)," arrays")
##apply is faster apply but takes too much memory
#### Index <- 1:nrow(Calls)
Index2 <- 1:nrow(Calls)
for(j in 1:ncol(res)){
## if(maleIndex[j]){
## Index2 <- Index[-XIndex]
## }else{
## Index2 <- Index
## }
## Index2 <- Index
if (verbose) cat(".")
tmpdist <- cbind(abs(Dist[,j,1]-Dist[,j,2]),abs(Dist[,j,2]-Dist[,j,3]))
tmpIndex <- split(Index2, factor(Calls[Index2,j], levels=1:3), drop=FALSE)
if (length(tmpIndex[[1]])>0) res[tmpIndex[[1]],j] <- tmpdist[tmpIndex[[1]],1]
if (length(tmpIndex[[3]])>0) res[tmpIndex[[3]],j] <- tmpdist[tmpIndex[[3]],2]
if (length(tmpIndex[[2]])>0) res[tmpIndex[[2]],j] <- pmin(tmpdist[tmpIndex[[2]], 1],
tmpdist[tmpIndex[[2]], 2])
rm(tmpIndex, tmpdist); ## gc()
## if(maleIndex[j]){
## Index2 <- Index[XIndex]
## res[Index2,j] <- abs(Dist[Index2,j,1]-Dist[Index2,j,3])
## }
}
if (verbose) cat("Done\n")
return(res)
}
replaceAffySnpParams <- function(object,value,subset){
object$centers[subset,,] <- value$centers
object$scales[subset,,] <- value$scales
object$N[subset,] <- value$N
return(object)
}
crlmm.old <- function(object, correction=NULL, recalibrate=TRUE,
minLLRforCalls=c(5, 1, 5),
verbose=TRUE, correctionFile=NULL, prefix="tmp.crlmm.", balance=1.5){
library(annotation(object), character.only=TRUE)
if(is.null(correctionFile))
stop("Provide correctionFile.\nIf the correctionFile is not found, it will be created and it will contain the EM results.")
if(file.exists(correctionFile)){
if (verbose) cat("File with correction information - from EM - was found.", "Loading...", sep="\n")
load(correctionFile)
}else{
if(verbose) cat("M correction not provided. Calculating. Will take several minutes.\n")
correction <- fitAffySnpMixture(object,verbose=verbose)
save(correction, file=correctionFile)
}
if (verbose) cat("Done.\n")
snr <- correction$snr
if(is.null(object$gender)){
if(verbose) cat("Gender not provided... using data to predict.\n")
warning("Gender not provided... using data to predict.")
maleIndex <- snpGenderCall(object)=="male"
}else{
maleIndex <- object$gender=="male"
}
annotname <- annotation(object)
load(system.file(paste("extdata/", annotname, "CrlmmInfo.rda", sep=""), package=annotname))
myenv <- get(paste(annotname,"Crlmm",sep="")); rm(list=paste(annotname,"Crlmm",sep=""))
thePriors <- get("priors", myenv)
## Index <- which(!get("hapmapCallIndex",myenv) | get("badCallIndex",myenv) | get("badRegions", myenv))
Index <- which(!get("hapmapCallIndex",myenv))
myCalls <- matrix(NA,dim(object)[1],dim(object)[2])
myCalls[Index,] <- getInitialAffySnpCalls(correction,Index,verbose=verbose, sqsClass=class(object))
fs <- correction$fs; rm(correction)
rparams <- getAffySnpGenotypeRegionParams(object, myCalls, fs,
subset=Index,verbose=verbose, sqsClass=class(object))
rm(myCalls); ## gc()
if (bothStrands(object)){
oneStrand <- apply(is.na(getM(object[,1])[,1,]), 1,
function(v) ifelse(length(ll <- which(v))==0, 0, ll))
rparams <- updateAffySnpParams(rparams, thePriors, oneStrand, verbose=verbose)
params <- replaceAffySnpParams(get("params",myenv), rparams, Index)
}else{
rparams <- updateAffySnpParamsSingle(rparams, thePriors, verbose=verbose)
params <- replaceAffySnpParamsSingle(get("params",myenv), rparams, Index)
}
## save(rparams, params, file=paste(prefix, "ParamsBeforeRec.rda", sep=""))
rm(myenv, Index)
if (bothStrands(object)){
myDist <- getAffySnpDistance(object, params, fs)
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myDist[,,-2,] <- balance*myDist[,,-2,]
}else{
myDist <- getAffySnpDistanceSingle(object, params, fs)
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myDist[,,-2] <- balance*myDist[,,-2]
}
if (!recalibrate)
## rm(params)
rm(fs); ## gc()
XIndex <- getChrXIndex(object)
myCalls <- getAffySnpCalls(myDist,XIndex,maleIndex,verbose=verbose, sqsClass=class(object))
LLR <- getAffySnpConfidence(myDist,myCalls,XIndex,maleIndex,verbose=verbose, sqsClass=class(object))
rm(myDist)
load(correctionFile)
fs <- correction$fs; rm(correction)
if(recalibrate){
if(verbose) cat("Recalibrating.")
for(k in 1:3)
myCalls[myCalls == k & LLR < minLLRforCalls[k]] <- NA
rm(LLR)
myCalls[, snr < 3.675] <- NA
rparams <- getAffySnpGenotypeRegionParams(object, myCalls,
fs, verbose=verbose, sqsClass=class(object))
rm(myCalls)
## gc()
if (bothStrands(object)){
rparams <- updateAffySnpParams(rparams, thePriors, oneStrand, verbose=verbose)
myDist <- getAffySnpDistance(object, rparams, fs, verbose=verbose)
myDist[,,-2,] <- balance*myDist[,,-2,]
}else{
rparams <- updateAffySnpParamsSingle(rparams, thePriors, verbose=verbose)
myDist <- getAffySnpDistanceSingle(object, rparams, fs, verbose=verbose)
myDist[,,-2] <- balance*myDist[,,-2]
}
## save(rparams, file=paste(prefix, "ParamsAfterRec.rda", sep=""))
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myCalls <- getAffySnpCalls(myDist,XIndex, maleIndex, verbose=verbose, sqsClass = class(object))
LLR <- getAffySnpConfidence(myDist,myCalls,XIndex,maleIndex,verbose=verbose, sqsClass = class(object))
rm(fs, myDist)
pacc <- LLR2conf(myCalls, LLR, snr, annotation(object))
}
## correction$snr
## maleIndex
gender <- rep("female", length(maleIndex))
gender[maleIndex] <- "male"
if (is.null(object$gender)){
addPhenoData <- new("AnnotatedDataFrame",
data=cbind(pData(object),
data.frame(crlmmSNR=as.numeric(snr),
gender=gender,
row.names=sampleNames(object))),
varMetadata= rbind(varMetadata(object),
data.frame(labelDescription=c("crlmmSNR", "gender"),
row.names=c("crlmmSNR", "gender"))))
}else{
addPhenoData <- new("AnnotatedDataFrame",
data=cbind(pData(object),
data.frame(crlmmSNR=as.numeric(snr),
row.names=sampleNames(object))),
varMetadata= rbind(varMetadata(object),
data.frame(labelDescription=c("crlmmSNR"),
row.names=c("crlmmSNR"))))
}
return(new("SnpSet",
phenoData=addPhenoData,
experimentData=experimentData(object),
annotation=annotation(object),
call=myCalls,
callProbability=pacc,
LLR=LLR))
}
###this one just for us
### addRegions <- function(i,params,...){
### require(ellipse)
### ADD <- params$f0*c(1,0,-1)
### idx <- which(!apply(is.na(params$centers[i,,]), 2, all))
### if (length(idx) == 1) idx <- rep(idx, 2)
### for(k in 1:3){
### points(t(params$centers[i,k,idx])+ADD[k],pch="+",col=k)
### lines(ellipse(diag(2),scale=params$scales[i,k,idx],centre=params$centers[i,k,idx]+ADD[k]),col=k,...)
### }
### }
LLR2conf <-function(theCalls, LLR, SNR, annot){
thePath <- system.file("extdata", package=annot)
theFile <- paste(annot, ".spline.params.rda", sep="")
load(file.path(thePath, theFile))
Het <- as.vector(theCalls==2)
## dst <- rep(Dst, ncol(theCalls))
LLR <- as.vector(sqrt(LLR))
conf <- vector("numeric", length(LLR))
## geting rid of NOTES:
HmzK3 <- get("HmzK3")
lm1 <- get("lm1")
HmzK2 <- get("HmzK2")
HtzK3 <- get("HtzK3")
lm2 <- get("lm2")
HtzK2 <- get("HtzK2")
SNRK <- get("SNRK")
SNRlm <- get("SNRlm")
if (any(!Het)){
tmp <- pmin(LLR[!Het], HmzK3)
conf[!Het] <- lm1$coef[1]+lm1$coef[2]*tmp+lm1$coef[3]*(tmp-HmzK2)*I(tmp>HmzK2)
}
if (any(Het)){
tmp <- pmin(LLR[Het], HtzK3)
conf[Het] <- lm2$coef[1]+lm2$coef[2]*tmp+lm2$coef[3]*(tmp-HtzK2)*I(tmp>HtzK2)
}
conf <- matrix(conf, ncol=ncol(theCalls))
X <- pmin(log(SNR), SNRK)
SNRfix <- SNRlm$coef[1]+SNRlm$coef[2]*X
conf <- sweep(conf, 2, SNRfix, FUN="+")
conf <- 1/(1+exp(-conf))
conf[,SNR<=3] <- 1/3
conf[conf<1/3] <- 1/3
return(conf)
}
benilton/oligo-release documentation built on May 12, 2019, 10:59 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### WARNING: IF ROBUST TO DIMNAMES,
### REMOVE AND SAVE LOTS
### OF MEMORY... ADD DIMNAMES
### AT THE END
rowEntropy <- function(p) rowMeans(rowSums(log2(p^p), dims=2))
getSnpFragmentLength <- function(object, enzyme){
## returns data.frame with 3 columns:
## SNP | Enzyme | Length
conn <- db(get(annotation(object)))
snps <- paste("(", paste("'", featureNames(object), "'", collapse=", ", sep=""), ")")
sql <- paste("SELECT man_fsetid, enzyme, length",
"FROM fragmentLength",
"INNER JOIN featureSet USING(fsetid)",
"WHERE man_fsetid IN", snps)
if (!missing(enzyme)){
stopifnot(length(enzyme) == 1)
sql <- paste(sql, "AND enzyme =", enzyme)
}
dbGetQuery(conn, sql)
}
snpGenderCall <- function(object){
XIndex <- getChrXIndex(object)
theA <- getA(object)
if (!is.ff(theA)){
tmp <- median(theA, na.rm=TRUE)
}else{
nc <- ncol(object)
tmp <- rep(NA, nc)
for (i in 1:nc) tmp[i] <- median(getA(object[,i]), na.rm=TRUE)
tmp <- median(tmp)
}
if (bothStrands(object)){
a <- apply(theA[XIndex,,], 2, median, na.rm=TRUE)
}else{
a <- apply(theA[XIndex,], 2, median, na.rm=TRUE)
}
kfit=kmeans(a, c(min(a),tmp))
return(factor(c("female","male")[as.numeric(kfit$cluster==1)+1]))
}
getChrXIndex <- function(object){
sql <- "SELECT man_fsetid, chrom, physical_pos FROM featureSet WHERE man_fsetid LIKE 'SNP%'"
chrs <- dbGetQuery(db(get(annotation(object))), sql)
chrs <- chrs[order(chrs[["chrom"]], chrs[["physical_pos"]], chrs[["man_fsetid"]], na.last=FALSE),]
return(which(chrs[["chrom"]] == "X"))
}
##gender in pData keeps male female
fitAffySnpMixture <- function(object, df1=3, df2=5,
probs=rep(1/3,3), eps=50,
subSampleSize=10^4, verbose=TRUE){
if(is.null(object$gender)){
maleIndex <- snpGenderCall(object)=="male"
}else{
maleIndex <- object$gender=="male"
}
XIndex=getChrXIndex(object)
I <- dim(object)[1]
J <- dim(object)[2]
set.seed(1)
## tmp <- c( (1:I)[-XIndex],((I+1):(2*I))[-XIndex])
## tmp <- c( (1:I),((I+1):(2*I)))
if (I > subSampleSize){
idx <- sort(sample(I, subSampleSize))
}else{
idx <- 1:I
}
## rm(tmp)
if (bothStrands(object)){
pis <- array(0,dim=c(I,J,3,2))
fs <- array(0,dim=c(I,J,2))
snr <- array(0,dim=J)
}else{
pis <- array(0,dim=c(I,J,3))
fs <- array(0,dim=c(I,J))
snr <- array(0,dim=J)
}
##### dimnames(fs)<-list(featureNames(object),
##### sampleNames(object),
##### c("antisense","sense"))
##### dimnames(pis)<-list(featureNames(object),
##### sampleNames(object),
##### c("AA","AB","BB"),
##### c("antisense","sense"))
##### names(snr) <- sampleNames(object)
if(verbose) cat("Fitting mixture model to ", J, " arrays. Epsilon must reach ", eps, ".\n",sep="")
## getSnpFragmentLength now returns a table SNP | enzyme | length
## and it is possible that there are 2+ lengths by SNPxEnzyme combo
## 21 Oct 2011
L <- getSnpFragmentLength(object)
L <- aggregate(L$length, by=list(SNP=L$man_fsetid), mean)
L <- L[match(featureNames(object), L$SNP), 2]
fix <- which(is.na(L))
L[fix] <- median(L, na.rm=T)
rm(fix)
if (bothStrands(object)){
L <- cbind(L,L)
l <- as.numeric(L[idx, ])
L <- as.numeric(L)
}else{
l <- L[idx]
}
L <- L-mean(L)
l <- l-mean(l)
matL <- ns(l, df1)
for(j in 1:J){
Y <- getM(object[,j])
A <- getA(object[,j])
if (bothStrands(object)){
Y <- Y[,1,]
A <- A[,1,]
}
fix <- which(is.na(Y))
Y[fix] <- median(Y, na.rm=T)
A[fix] <- median(A, na.rm=T)
rm(fix)
mus <- quantile(Y,c(1,3,5)/6);mus[2]=0
sigmas <- rep(mad(c(Y[Y<mus[1]]-mus[1],Y[Y>mus[3]]-mus[3])),3)
sigmas[2] <- sigmas[2]/2
if (bothStrands(object)){
a <- as.numeric(A[idx, ])
y <- as.numeric(Y[idx, ])
}else{
a <- A[idx]
y <- Y[idx]
}
A <- A-mean(A)
a <- a-mean(a)
weights <- apply(cbind(mus, sigmas), 1, function(p) dnorm(y, p[1], p[2]))
PreviousLogLik <- -Inf
change <- eps+1
itmax <- 0
matA <- ns(a,df2)
while (change > eps & itmax < 100){
itmax <- itmax+1
## E
z <- sweep(weights, 2, probs, "*")
LogLik <- rowSums(z)
z <- sweep(z, 1, LogLik, "/")
LogLik <- sum(log(LogLik))
change <- abs(LogLik-PreviousLogLik)
if (verbose){
if (itmax > 1 | j > 1) cat(del)
message <- paste("Array ",j,": epsilon = ", signif(change,2), " ", sep="")
del <- paste(rep("\b", nchar(message)), collapse="")
cat(message)
}
PreviousLogLik <- LogLik
probs <- colMeans(z)
## M
fit1 <- lm(y~matL+matA,weights=z[,1])
fit2 <- sum(z[,2]*y)/sum(z[,2])
fit3 <- lm(y~matL+matA,weights=z[,3])
sigmas[1] <- sqrt(sum(z[,1]*residuals(fit1)^2)/sum(z[,1]))
sigmas[2] <- sqrt(sum(z[,2]*(y-fit2)^2)/sum(z[,2]))
sigmas[3] <- sqrt(sum(z[,3]*residuals(fit3)^2)/sum(z[,3]))
weights[,1] <- dnorm(y,fitted(fit1),sigmas[1])
weights[,2] <- dnorm(y,fit2,sigmas[2])
weights[,3] <- dnorm(y,fitted(fit3),sigmas[3])
weights[y >= 0, 1] <- 0
weights[y <= 0, 3] <- 0
}
## gc()
bigX <- cbind(1, ns(L, knots=as.numeric(attr(matL, "knots")), Boundary.knots=attr(matL, "Boundary.knots")),
ns(A, knots=as.numeric(attr(matA, "knots")), Boundary.knots=attr(matA, "Boundary.knots")))
rm(matA); ## gc()
pred1 <- bigX%*%coef(fit1)
pred2 <- rep(fit2,length(Y))
pred3 <- bigX%*%coef(fit3)
rm(bigX); ## gc()
weights <- matrix(0,length(Y),3)
weights[,1] <- dnorm(Y,pred1,sigmas[1])
weights[,2] <- dnorm(Y,pred2,sigmas[2])
weights[,3] <- dnorm(Y,pred3,sigmas[3])
weights[Y >= pred2, 1] <- 0
## if (maleIndex[j]) weights[XIndex, 2] <- 0
weights[Y <= pred2, 3] <- 0
z <- sweep(weights, 2, probs, "*")
LogLik <- rowSums(z)
z <- sweep(z, 1, LogLik, "/")
if (bothStrands(object)){
fs[,j,] <- matrix((pred3-pred1)/2,ncol=2)
for(k in 1:3){
pis[,j,k,] <- matrix(z[,(4-k)],ncol=2) ##4-k cause 3is1,2is2 and 1is3
}
snr[j] <- median(fs[,j,])^2/(sigmas[1]^2+sigmas[2]^2)
}else{
fs[, j] <- (pred3-pred1)/2
for(k in 1:3){
pis[, j, k] <- z[,(4-k)] ##4-k cause 3is1,2is2 and 1is3
}
snr[j] <- median(fs[,j])^2/(sigmas[1]^2+sigmas[2]^2)
}
}
fs[fs < 0] <- 0
if(verbose) cat("Done.\n")
return(list(f0=median(fs),fs=fs, pis=pis, snr=snr))
}
getInitialAffySnpCalls <- function(object,subset=NULL,
concordanceCutoff=0.0001,
cutoffs=c(0.7,0.5,0.7),
returnProbs=FALSE,
verbose=FALSE,
sqsClass = "SnpQSet"){
if(is.null(subset)) subset <- 1:(dim(object$pis)[1])
if (sqsClass == "SnpQSet"){
pi1 <- object$pis[subset,,,1]
pi2 <- object$pis[subset,,,2]; rm(object); ## gc()
## fixing SNPs that have probes only on one strand:
idx <- is.na(pi1[,1,1])
pi1[idx,,] <- pi2[idx,,]
idx <- is.na(pi2[,1,1])
pi2[idx,,] <- pi1[idx,,]
rm(idx); ## gc()
if(verbose) cat("Picking good starting value: ")
if(verbose) cat("Computing entropy, ")
E1 <- rowEntropy(pi1)
E2 <- rowEntropy(pi2); ## gc()
tmpN <- dim(pi1)
tmpcall1 <- tmpcall2 <- matrix(NA, nrow=tmpN[1], ncol=tmpN[2])
##### rownames(tmpcall1) <- rownames(tmpcall2) <- dimnames(pi1)[[1]]
##### colnames(tmpcall1) <- colnames(tmpcall2) <- dimnames(pi1)[[2]]
## gc()
if(verbose) cat("calculating calls, ")
for (i in 1:tmpN[1]){
for (j in 1:tmpN[2]){
tmpcall1[i, j] <- which.max(pi1[i,j,])
tmpcall2[i, j] <- which.max(pi2[i,j,])
}
}
if(verbose) cat("determining non-concordant calls, ")
concordance <- rowIndepChiSqTest(tmpcall1,tmpcall2); ## gc()
if(verbose) cat("deciding which strand(s) to use")
tc1 <- tmpcall1 == 1
tc2 <- tmpcall1 == 2
tc3 <- tmpcall1 == 3
noABIndex1 <- which((rowSums(tc2)<3)*(rowSums(tc3)>0)*(rowSums(tc1)>0) == 1)
tc1 <- tmpcall2 == 1
tc2 <- tmpcall2 == 2
tc3 <- tmpcall2 == 3
noABIndex2 <- which((rowSums(tc2)<3)*(rowSums(tc3)>0)*(rowSums(tc1)>0) == 1)
rm(tc1, tc2, tc3); ## gc()
E1[noABIndex1] <- -Inf
E2[noABIndex2] <- -Inf
jointprobs <- (pi1+pi2)/2
## gc()
noInfoIndex <- intersect(noABIndex1, noABIndex2)
rm(noABIndex1, noABIndex2)
jointprobs[noInfoIndex,,] <- 1/3
rm(noInfoIndex)
notBoth <- concordance < concordanceCutoff
E1bE2 <- E1 > E2
rm(E1, E2)
i1 <- notBoth & E1bE2
i2 <- notBoth & !E1bE2
rm(notBoth, E1bE2)
jointprobs[i1,,] <- pi1[i1,,]
rm(i1, pi1)
jointprobs[i2,,] <- pi2[i2,,]
rm(i2, pi2); ## gc()
tmpN <- dim(jointprobs)
tmpcall <- tmpmax <- matrix(NA, nrow=tmpN[1], ncol=tmpN[2])
##### rownames(tmpcall) <- rownames(tmpmax) <- dimnames(jointprobs)[[1]]
##### colnames(tmpcall) <- colnames(tmpmax) <- dimnames(jointprobs)[[2]]
## gc()
if(verbose) cat(" finalizing"); ## gc()
for (i in 1:tmpN[1]){
for (j in 1:tmpN[2]){
tmpcall[i, j] <- which.max(jointprobs[i, j, ])
tmpmax[i, j] <- jointprobs[i, j, tmpcall[i,j]]
}
}
}else{
tmpcall <- apply(object$pis[subset,, ], 1:2, which.max)
tmpmax <- apply(object$pis[subset,, ], 1:2, max)
}
for(i in 1:3)
tmpcall[tmpcall==i & tmpmax<cutoffs[i]] <- NA
if(verbose) cat("\nCompleted!\n")
if(sqsClass == "SnpQSet"){
if(returnProbs) return(list(calls=tmpcall,probs=jointprobs)) else return(tmpcall)
}else{
return(tmpcall)
}
}
rowIndepChiSqTest <- function(call1,call2){
tmp=vector("numeric",nrow(call1))
for(i in 1:nrow(call1)){
tmpt= cbind(table(factor(call1[i,],levels=1:3)),
table(factor(call2[i,],levels=1:3)))
tmpt=tmpt[!rowSums(tmpt)==0,,drop=FALSE]
if(nrow(tmpt)>1){
rowtot <- rowSums(tmpt)
coltot <- colSums(tmpt)
e=outer(rowtot,coltot)/sum(tmpt)
stat = sum((tmpt-e)^2/e)
tmp[i] = 1 - pchisq(stat, (nrow(tmpt) - 1) * (ncol(tmpt) - 1))
}
else{
tmp[i]=0
}
}
return(tmp)
}
getGenotypeRegionParams <- function(M, initialcalls, f=0, verbose=TRUE, subset){
if (!missing(subset)){
M <- M[subset,, drop=FALSE]
initialcalls <- initialcalls[subset,, drop=FALSE]
f <- f[subset,, drop=FALSE]
}
if(verbose) cat("Computing centers and scales for 3 genotypes")
## tmp <- .Call("R_HuberMatrixRows2", M+(initialcalls-2)*f,
## as.integer(initialcalls), 1.5)
tmp <- .Call("R_trimmed_stats", M+(initialcalls-2)*f,
as.integer(initialcalls), 0.025)
##### centers <- scales <- N <- array(NA, dim=c(nrow(M),3))
##### scales <- N <- array(NA, dim=c(nrow(M),3))
##### dimnames(centers) <- dimnames(scales) <- dimnames(N) <- list(rownames(M), c("AA","AB","BB"))
##### centers[,] <- tmp[[1]]
tmpN <- tmp[[3]]; tmp[[3]] <- NULL
tmpN[tmpN <= 1] <- NA
tmpS <- tmp[[2]]; tmp[[2]] <- NULL
tmpS[is.na(tmpN)] <- NA
tmpV <- sqrt(rowSums((tmpN[,-2]-1)*tmpS[,-2]^2, na.rm=T)/rowSums(tmpN[,-2]-1, na.rm=T))
tmpV[!is.finite(tmpV)] <- NA
tmpS[,-2] <- tmpV
rm(tmpV)
tmpC <- tmp[[1]]; rm(tmp)
tmpC[is.na(tmpN)] <- NA
tmpN[is.na(tmpN)] <- 0
if(verbose) cat(" Done\n")
return(list(centers=tmpC, scales=tmpS, N=tmpN))
}
getAffySnpGenotypeRegionParams<-function(object,initialcalls,f=NULL,
subset=1:(dim(object)[1]),
verbose=FALSE, sqsClass = "SnpQSet"){
if(is.null(f)) f=fitAffySnpMixture(object,verbose=verbose)$fs
if(verbose) cat("Computing centers and scales:\n")
if (sqsClass == "SnpQSet"){
N <- scales <- centers <- array(NA,dim=c(length(subset),3,2))
##### dimnames(N) <- dimnames(scales) <- dimnames(centers) <- list(featureNames(object)[subset],
##### c("AA","AB","BB"), c("antisense","sense"))
for(s in 1:2){
## gc()
if(verbose) cat(c("\tantisense","\tsense....")[s],": ",sep="")
tmp <- getGenotypeRegionParams(getM(object)[subset,,s],
initialcalls[subset,],
f[subset,,s],
verbose=verbose)
centers[,,s] <- tmp$centers
scales[,,s] <- tmp$scales
N[,,s] <- tmp$N
}
N <- apply(N, 1:2, max, na.rm=TRUE) ## returns -Inf if c(NA, NA)
}else{
tmp <- getGenotypeRegionParams(getM(object)[subset,],
initialcalls[subset,],
f[subset,],
verbose=verbose)
centers <- tmp$centers
scales <- tmp$scales
N <- tmp$N
}
return(list(centers=centers,scales=scales,N=N,f0=median(f, na.rm=TRUE)))
}
getAffySnpPriors <- function(object,minN=20,subset=1:(dim(object$centers)[1]),
verbose=TRUE, sqsClass="SnpQSet"){
if(verbose) cat("Computing priors.\n")
if (sqsClass == "SnpQSet"){
N <- cbind(object$N,object$N)
Index <- subset[which(rowMeans(N[subset,]>minN)==1)]
N <- N[Index,]
mus <- cbind(object$centers[Index,,1],object$centers[Index,,2])
sigmas <- cbind(object$scales[Index,,1],object$scales[Index,,2])
sigmas[is.na(mus)] <- NA
maxsigmas <- c(quantile(sigmas[,c(1,3,4,6)],.99, na.rm=TRUE),
quantile(sigmas[,c(2,5)],.99, na.rm=TRUE))
}else{
N <- object$N
Index <- subset[which(rowMeans(N[subset,] > minN) == 1)]
N <- N[Index,]
mus <- object$centers[Index,]
sigmas <- object$scales[Index,]
maxsigmas <- c(quantile(sigmas[,c(1, 3)],.99, na.rm=TRUE),
quantile(sigmas[, 2],.99, na.rm=TRUE))
}
##variance for prior for mu
##need to make robust
V <- cov(mus,use="pairwise.complete")
##prior for sigma
zgs <- log(sigmas^2); rm(sigmas); ## gc()
dgs <- N-1
egs <- zgs-digamma(dgs/2)+log(dgs/2); rm(zgs); ## gc()
n <- length(Index)
d0s <- 2*trigammaInverse(colMeans(n*(egs-colMeans(egs, na.rm=TRUE))^2/(n-1)-trigamma(dgs/2), na.rm=TRUE))
s20 <- exp(colMeans(egs, na.rm=TRUE)+digamma(d0s/2)-log(d0s/2))
return(list(V=V,d0s=d0s,s20=s20,maxsigmas=maxsigmas))
}
##minN- minimum number of points in cluster required for use
##max$Sigma
updateAffySnpParams <- function(object, priors, missingStrandIndex, minN=3,
maxHomoSigma=priors$maxsigma[1],
maxHeteSigma=priors$maxsigma[2],
subset=1:(dim(object$centers)[1]),
d0s=80, verbose=FALSE){
object$centers <- object$centers[subset,,]
object$scales <- object$scales[subset,,]
object$N <- object$N[subset,]
missingStrandIndex <- missingStrandIndex[subset]
if(verbose) cat("Updating centers and scales")
##First variances
for(i in 1:2){
for(j in 1:2){ ##1 and 3 are the same
if(j==2) N <- object$N[,2] else N <- rowSums(object$N[,-2],na.rm=TRUE)
s <- object$scales[,j,i]
if (is.null(d0s))
d0s <- priors$d0s[3*(i-1)+j]
s20 <- priors$s20[3*(i-1)+j] ##notice the ad-hoc choice of 3
Index <- N>minN & !is.na(s)
N <- N[Index]; s <- s[Index]
object$scales[Index,j,i] <- sqrt ( ( (N-1)*s^2 + d0s*s20 ) / (d0s+N-1) )
object$scales[!Index & missingStrandIndex != i,j,i] <- sqrt(s20)
}
}
object$scales[,3,] <- object$scales[,1,] ##AA=BB
object$scales[,2,][object$scales[,2,]>maxHeteSigma] <- maxHeteSigma
object$scales[,-2,][object$scales[,-2,]>maxHomoSigma] <- maxHomoSigma
if(verbose) cat(".")
##Means
updateMean <- function(strandIndex, strand){
## strandIndex is a vector where
## 0 - none strands are missing on the array
## 1 - antisense strand is missing
## 2 - sense strand is missing
## strand argument is what do you want to be updated
if (strand == "both"){
snps <- strandIndex == 0; idx <- 1:6
N <- cbind(object$N[snps,],object$N[snps,])
mu <- cbind(object$centers[snps,,1],object$centers[snps,,2])
}else if(strand == "antisense"){
snps <- strandIndex == 2; idx <- 1:3
N <- object$N[snps,]
mu <- object$centers[snps,,1]
}else if(strand == "sense"){
snps <- strandIndex == 1; idx <- 4:6
N <- object$N[snps,]
mu <- object$centers[snps,,2]
}
if (sum(snps) == 1){
N <- matrix(N, nrow=1)
mu <- matrix(mu, nrow=1)
}
Vinv <- solve(priors$V[idx, idx])
NSinv <- t(N)/priors$s20[idx]
tmp <- t(sapply(1:nrow(mu),function(i){
if(verbose & i%%5000==0) cat(".")
mus=mu[i,]; Ns=N[i,]
mus[Ns<minN] <- 0
mus[is.na(mus)] <- 0
return(solve(Vinv+diag(NSinv[,i]))%*%(NSinv[,i]*mus))
}))
return(tmp)
}
if (any(missingStrandIndex == 0)){
tmp <- updateMean(missingStrandIndex, "both")
object$centers[missingStrandIndex == 0,,1] <- tmp[,1:3]
object$centers[missingStrandIndex == 0,,2] <- tmp[,4:6]; rm(tmp);
}
if (any(missingStrandIndex == 2))
object$centers[missingStrandIndex == 2,,1] <- updateMean(missingStrandIndex, "antisense")
if (any(missingStrandIndex == 1))
object$centers[missingStrandIndex == 1,,2] <- updateMean(missingStrandIndex, "sense")
if(verbose) cat("Done.\n")
return(object)
}
getAffySnpDistance <- function(object,params,f=0,subset=1:(dim(object)[1]),
w=NULL,verbose=FALSE){
x=getM(object[subset,])
rm(object)
Dist <- array(NA,dim=c(dim(x)[1],dim(x)[2],3,2))
if(verbose) cat("Calculating likelihood-based distances")
for(i in 1:2){
if(verbose) cat(".")
for(j in 1:3){
tmp <- x[,,i]+(j-2)*f[subset,,i]
Dist[,,j,i] <- 2*log(params$scales[subset,j,i]) +
((tmp-params$centers[subset,j,i])/params$scales[subset,j,i])^2
if(!is.null(w)) Dist[,,j,i] <- Dist[,,j,i] - 2*log(w[subset,,j,i])
}
}
##### dimnames(Dist) <- list(dimnames(x)[[1]],
##### dimnames(x)[[2]], 1:3,
##### dimnames(x)[[3]])
if(verbose) cat("Done.\n")
return(Dist)
}
getAffySnpCalls <- function(Dist,XIndex,maleIndex,subset=1:(dim(Dist)[1]),
verbose=FALSE, sqsClass="SnpQSet"){
if (sqsClass == "SnpQSet"){
Dist <- Dist[subset,,,,drop=FALSE]
}else{
Dist <- Dist[subset,,, drop=FALSE]
}
#### XIndex <- which(subset%in%XIndex)
## gc()
## Here we can even put default value to 3 for the new code. /HB
res <- array(as.integer(-1),dim=dim(Dist)[c(1,2)]); ## gc()
##### dimnames(res) <- dimnames(Dist)[1:2]
if (sqsClass == "SnpQSet")
Dist <- rowSums(Dist, na.rm=TRUE, dims=3)
## Dist[XIndex, maleIndex, 2] <- Inf
##the following is slow!
if(verbose) cat("Making calls for ", ncol(res), " arrays");
for(j in 1:ncol(res)){
if(verbose) cat(".");
D1 <- Dist[,j,1];
D2 <- Dist[,j,2];
D3 <- Dist[,j,3];
d12 <- (D1 < D2);
d23 <- (D2 < D3); rm(D2);
d13 <- (D1 < D3); rm(D3);
d <- rep(as.integer(3), length(D1)); rm(D1)
d[( d12 & d13)] <- as.integer(1); rm(d13)
d[(!d12 & d23)] <- as.integer(2); rm(d12, d23)
res[,j] <- d;
rm(d);
}
if(verbose) cat("Done\n")
return(res)
}
getAffySnpConfidence <- function(Dist, Calls, XIndex, maleIndex,
subset=1:nrow(Calls), verbose=TRUE, sqsClass = "SnpQSet"){
if (sqsClass == "SnpQSet"){
Dist <- rowSums(Dist[subset,,,,drop=FALSE], dims=3, na.rm=T)
}else{
Dist <- Dist[subset,,, drop=FALSE]
}
Calls <- Calls[subset,,drop=FALSE]
#### XIndex <- which(subset%in%XIndex)
res <- array(NA,dim=dim(Dist)[c(1,2)])
##### dimnames(res) <- list(dimnames(Dist)[[1]],dimnames(Dist)[[2]])
##### Dist <- rowSums(Dist, dims=3, na.rm=T)
if (verbose) cat("Computing confidence for calls on ",ncol(res)," arrays")
##apply is faster apply but takes too much memory
#### Index <- 1:nrow(Calls)
Index2 <- 1:nrow(Calls)
for(j in 1:ncol(res)){
## if(maleIndex[j]){
## Index2 <- Index[-XIndex]
## }else{
## Index2 <- Index
## }
## Index2 <- Index
if (verbose) cat(".")
tmpdist <- cbind(abs(Dist[,j,1]-Dist[,j,2]),abs(Dist[,j,2]-Dist[,j,3]))
tmpIndex <- split(Index2, factor(Calls[Index2,j], levels=1:3), drop=FALSE)
if (length(tmpIndex[[1]])>0) res[tmpIndex[[1]],j] <- tmpdist[tmpIndex[[1]],1]
if (length(tmpIndex[[3]])>0) res[tmpIndex[[3]],j] <- tmpdist[tmpIndex[[3]],2]
if (length(tmpIndex[[2]])>0) res[tmpIndex[[2]],j] <- pmin(tmpdist[tmpIndex[[2]], 1],
tmpdist[tmpIndex[[2]], 2])
rm(tmpIndex, tmpdist); ## gc()
## if(maleIndex[j]){
## Index2 <- Index[XIndex]
## res[Index2,j] <- abs(Dist[Index2,j,1]-Dist[Index2,j,3])
## }
}
if (verbose) cat("Done\n")
return(res)
}
replaceAffySnpParams <- function(object,value,subset){
object$centers[subset,,] <- value$centers
object$scales[subset,,] <- value$scales
object$N[subset,] <- value$N
return(object)
}
crlmm.old <- function(object, correction=NULL, recalibrate=TRUE,
minLLRforCalls=c(5, 1, 5),
verbose=TRUE, correctionFile=NULL, prefix="tmp.crlmm.", balance=1.5){
library(annotation(object), character.only=TRUE)
if(is.null(correctionFile))
stop("Provide correctionFile.\nIf the correctionFile is not found, it will be created and it will contain the EM results.")
if(file.exists(correctionFile)){
if (verbose) cat("File with correction information - from EM - was found.", "Loading...", sep="\n")
load(correctionFile)
}else{
if(verbose) cat("M correction not provided. Calculating. Will take several minutes.\n")
correction <- fitAffySnpMixture(object,verbose=verbose)
save(correction, file=correctionFile)
}
if (verbose) cat("Done.\n")
snr <- correction$snr
if(is.null(object$gender)){
if(verbose) cat("Gender not provided... using data to predict.\n")
warning("Gender not provided... using data to predict.")
maleIndex <- snpGenderCall(object)=="male"
}else{
maleIndex <- object$gender=="male"
}
annotname <- annotation(object)
load(system.file(paste("extdata/", annotname, "CrlmmInfo.rda", sep=""), package=annotname))
myenv <- get(paste(annotname,"Crlmm",sep="")); rm(list=paste(annotname,"Crlmm",sep=""))
thePriors <- get("priors", myenv)
## Index <- which(!get("hapmapCallIndex",myenv) | get("badCallIndex",myenv) | get("badRegions", myenv))
Index <- which(!get("hapmapCallIndex",myenv))
myCalls <- matrix(NA,dim(object)[1],dim(object)[2])
myCalls[Index,] <- getInitialAffySnpCalls(correction,Index,verbose=verbose, sqsClass=class(object))
fs <- correction$fs; rm(correction)
rparams <- getAffySnpGenotypeRegionParams(object, myCalls, fs,
subset=Index,verbose=verbose, sqsClass=class(object))
rm(myCalls); ## gc()
if (bothStrands(object)){
oneStrand <- apply(is.na(getM(object[,1])[,1,]), 1,
function(v) ifelse(length(ll <- which(v))==0, 0, ll))
rparams <- updateAffySnpParams(rparams, thePriors, oneStrand, verbose=verbose)
params <- replaceAffySnpParams(get("params",myenv), rparams, Index)
}else{
rparams <- updateAffySnpParamsSingle(rparams, thePriors, verbose=verbose)
params <- replaceAffySnpParamsSingle(get("params",myenv), rparams, Index)
}
## save(rparams, params, file=paste(prefix, "ParamsBeforeRec.rda", sep=""))
rm(myenv, Index)
if (bothStrands(object)){
myDist <- getAffySnpDistance(object, params, fs)
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myDist[,,-2,] <- balance*myDist[,,-2,]
}else{
myDist <- getAffySnpDistanceSingle(object, params, fs)
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myDist[,,-2] <- balance*myDist[,,-2]
}
if (!recalibrate)
## rm(params)
rm(fs); ## gc()
XIndex <- getChrXIndex(object)
myCalls <- getAffySnpCalls(myDist,XIndex,maleIndex,verbose=verbose, sqsClass=class(object))
LLR <- getAffySnpConfidence(myDist,myCalls,XIndex,maleIndex,verbose=verbose, sqsClass=class(object))
rm(myDist)
load(correctionFile)
fs <- correction$fs; rm(correction)
if(recalibrate){
if(verbose) cat("Recalibrating.")
for(k in 1:3)
myCalls[myCalls == k & LLR < minLLRforCalls[k]] <- NA
rm(LLR)
myCalls[, snr < 3.675] <- NA
rparams <- getAffySnpGenotypeRegionParams(object, myCalls,
fs, verbose=verbose, sqsClass=class(object))
rm(myCalls)
## gc()
if (bothStrands(object)){
rparams <- updateAffySnpParams(rparams, thePriors, oneStrand, verbose=verbose)
myDist <- getAffySnpDistance(object, rparams, fs, verbose=verbose)
myDist[,,-2,] <- balance*myDist[,,-2,]
}else{
rparams <- updateAffySnpParamsSingle(rparams, thePriors, verbose=verbose)
myDist <- getAffySnpDistanceSingle(object, rparams, fs, verbose=verbose)
myDist[,,-2] <- balance*myDist[,,-2]
}
## save(rparams, file=paste(prefix, "ParamsAfterRec.rda", sep=""))
save(myDist, file=paste(prefix, "distFile.rda", sep=""))
myCalls <- getAffySnpCalls(myDist,XIndex, maleIndex, verbose=verbose, sqsClass = class(object))
LLR <- getAffySnpConfidence(myDist,myCalls,XIndex,maleIndex,verbose=verbose, sqsClass = class(object))
rm(fs, myDist)
pacc <- LLR2conf(myCalls, LLR, snr, annotation(object))
}
## correction$snr
## maleIndex
gender <- rep("female", length(maleIndex))
gender[maleIndex] <- "male"
if (is.null(object$gender)){
addPhenoData <- new("AnnotatedDataFrame",
data=cbind(pData(object),
data.frame(crlmmSNR=as.numeric(snr),
gender=gender,
row.names=sampleNames(object))),
varMetadata= rbind(varMetadata(object),
data.frame(labelDescription=c("crlmmSNR", "gender"),
row.names=c("crlmmSNR", "gender"))))
}else{
addPhenoData <- new("AnnotatedDataFrame",
data=cbind(pData(object),
data.frame(crlmmSNR=as.numeric(snr),
row.names=sampleNames(object))),
varMetadata= rbind(varMetadata(object),
data.frame(labelDescription=c("crlmmSNR"),
row.names=c("crlmmSNR"))))
}
return(new("SnpSet",
phenoData=addPhenoData,
experimentData=experimentData(object),
annotation=annotation(object),
call=myCalls,
callProbability=pacc,
LLR=LLR))
}
###this one just for us
### addRegions <- function(i,params,...){
### require(ellipse)
### ADD <- params$f0*c(1,0,-1)
### idx <- which(!apply(is.na(params$centers[i,,]), 2, all))
### if (length(idx) == 1) idx <- rep(idx, 2)
### for(k in 1:3){
### points(t(params$centers[i,k,idx])+ADD[k],pch="+",col=k)
### lines(ellipse(diag(2),scale=params$scales[i,k,idx],centre=params$centers[i,k,idx]+ADD[k]),col=k,...)
### }
### }
LLR2conf <-function(theCalls, LLR, SNR, annot){
thePath <- system.file("extdata", package=annot)
theFile <- paste(annot, ".spline.params.rda", sep="")
load(file.path(thePath, theFile))
Het <- as.vector(theCalls==2)
## dst <- rep(Dst, ncol(theCalls))
LLR <- as.vector(sqrt(LLR))
conf <- vector("numeric", length(LLR))
## geting rid of NOTES:
HmzK3 <- get("HmzK3")
lm1 <- get("lm1")
HmzK2 <- get("HmzK2")
HtzK3 <- get("HtzK3")
lm2 <- get("lm2")
HtzK2 <- get("HtzK2")
SNRK <- get("SNRK")
SNRlm <- get("SNRlm")
if (any(!Het)){
tmp <- pmin(LLR[!Het], HmzK3)
conf[!Het] <- lm1$coef[1]+lm1$coef[2]*tmp+lm1$coef[3]*(tmp-HmzK2)*I(tmp>HmzK2)
}
if (any(Het)){
tmp <- pmin(LLR[Het], HtzK3)
conf[Het] <- lm2$coef[1]+lm2$coef[2]*tmp+lm2$coef[3]*(tmp-HtzK2)*I(tmp>HtzK2)
}
conf <- matrix(conf, ncol=ncol(theCalls))
X <- pmin(log(SNR), SNRK)
SNRfix <- SNRlm$coef[1]+SNRlm$coef[2]*X
conf <- sweep(conf, 2, SNRfix, FUN="+")
conf <- 1/(1+exp(-conf))
conf[,SNR<=3] <- 1/3
conf[conf<1/3] <- 1/3
return(conf)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.