Nothing
cmap <- function(...)
UseMethod("cmap")
cmap.default <- function(...)
return(cmap.xyz(...))
cmap.xyz <-
function(xyz, grpby=NULL, dcut=4, scut=3, pcut=1, binary=TRUE, mask.lower = TRUE,
collapse=TRUE, gc.first = FALSE, ncore=1, nseg.scale=1, ...) {
# Parallelized by parallel package (Mon Apr 22 16:32:19 EDT 2013)
ncore <- setup.ncore(ncore, bigmem = TRUE)
if(ncore > 1) {
# Issue of serialization problem
# Maximal number of cells of a double-precision matrix
# that each core can serialize: (2^31-1-61)/8
R_NCELL_LIMIT_CORE = 2.68435448e8
R_NCELL_LIMIT = ncore * R_NCELL_LIMIT_CORE
if(nseg.scale < 1) {
warning("nseg.scale should be 1 or a larger integer\n")
nseg.scale=1
}
}
if (!(is.numeric(pcut) && pcut >= 0 && pcut <= 1)) {
stop("Input 'pcut' should a number between 0 and 1")
}
xyz=as.xyz(xyz)
nxyz =nrow(xyz)
if(nrow(xyz)>1) {
pb <- txtProgressBar(min=0, max=nrow(xyz), style=3)
if(is.null(grpby)) {
nres <- ncol(xyz)/3
} else {
inds <- bounds(grpby, dup.inds = TRUE)
nres <- nrow(inds)
}
if(ncore > 1) {
## shared memory to follow progress bar
iipb <- bigmemory::big.matrix(1, nrow(xyz), init=NA)
ni = nrow(xyz)
RLIMIT = floor(R_NCELL_LIMIT/(0.5*nres*(nres+1)))
nDataSeg = floor((ni-1)/RLIMIT)+1
nDataSeg = floor(nDataSeg * nseg.scale)
lenSeg = floor(ni/nDataSeg)
cmap.list <- NULL
for(i in 1:nDataSeg) {
istart = (i-1)*lenSeg + 1
iend = if(i<nDataSeg) i*lenSeg else ni
cmap.list <- c(cmap.list, mclapply(istart:iend, function(j) {
if(gc.first) gc()
dmat <- dm.xyz(xyz[j,], grpby, scut, mask.lower=TRUE)
## progress bar
iipb <- get("iipb")
iipb[1,j] <- 1
k <- length(which(!is.na(bigmemory::as.matrix(iipb))))
setTxtProgressBar(pb, k)
return(as.numeric(dmat[!lower.tri(dmat)] < dcut))
}) )
}
} else {
cmap.list <- lapply(1:nrow(xyz), function(j) {
dmat <- dm.xyz(xyz[j,], grpby, scut, mask.lower=TRUE)
setTxtProgressBar(pb, j)
return(as.numeric(dmat[!lower.tri(dmat)] < dcut))
})
}
rm(xyz)
if(gc.first) gc()
if(collapse) {
cmap.t <- rowMeans(do.call(cbind, cmap.list), na.rm=TRUE)
cmap.t[!is.finite(cmap.t)] <- NA
rm(cmap.list)
if(gc.first) gc()
if(binary) cmap.t <- as.numeric(cmap.t >= pcut )
cont.map <- matrix(NA, nrow=nres, ncol=nres)
cont.map[!lower.tri(cont.map)] <- cmap.t
if(!mask.lower)
cont.map[lower.tri(cont.map)] <- t(cont.map)[lower.tri(cont.map)]
}
else {
cont.map <- array(NA, dim=c(nres, nres, nxyz))
cmap.t <- matrix(NA, nres, nres)
for(i in 1:nxyz) {
cmap.t[!lower.tri(cmap.t)] <- cmap.list[[i]]
if(!mask.lower)
cmap.t[lower.tri(cmap.t)] <- t(cmap.t)[lower.tri(cmap.t)]
cont.map[,,i] <- cmap.t
}
rm(cmap.list)
if(gc.first) gc()
}
rm(cmap.t)
if(gc.first) gc()
close(pb)
} else {
## Distance matrix (all-atom)
dmat <- dm.xyz( xyz, grpby, scut, mask.lower = mask.lower, ncore=ncore)
## Contact map
return(matrix(as.numeric(dmat < dcut),
ncol = ncol(dmat),
nrow = nrow(dmat)))
}
return (cont.map)
}
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.