R/computecallsandcopy.R

In ndbrown6/CNu: Copy Number Utilities

'computecallsandcopy' <- function(path=NULL, algorithm="ascat", ncores=24)
{
	if (is.null(path)) {
		stop("path to facets output not supplied")
	}
	if (!dir.exists(paste0(path, "/resu/"))) {
		stop("directory 'resu' absent")
	}
	if (algorithm=="ascat") {
		if (!dir.exists(paste0(path, "/resu/ASCAT"))) {
			stop("directory 'ASCAT' absent")
		} else {
			demo = read.csv(file=paste0(path, "/resu/ASCAT/ASCAT_Params/ASCAT_Params.csv"), header=TRUE, sep=",", row.names=1, stringsAsFactors=FALSE)
			purity = demo[,"purity"]
			ploidy = demo[,"ploidy"]
		}
	} else if (algorithm=="gap") {
		if (!dir.exists(paste0(path, "/resu/GAP"))) {
			stop("directory 'GAP' absent")
		} else {
			demo = read.csv(file=paste0(path, "/resu/GAP/GAP_Params/GAP_Params.csv"), header=TRUE, sep=",", row.names=1, stringsAsFactors=FALSE)
			purity = 1-demo[,"p_BAF"]
			ploidy = demo[,"DNA_Index"]*2
		}
	}
	registerDoMC(ncores)
	sampleNames = rownames(demo)
	pb = txtProgressBar(min=1, max=length(sampleNames), style=3)
	res = foreach (i=1:length(sampleNames)) %dopar% {
		if (length(sampleNames)<=ncores) {
			if (i==length(sampleNames)) {
				setTxtProgressBar(pb, i)
			}
		} else {
			if ((i %% ncores)==0) {
				setTxtProgressBar(pb, i)
			}
		}
		rho = purity[i]
		psi = ploidy[i]
		load(file=paste0(path, "/resu/mad/", sampleNames[i], ".RData"))
		nonroundedcopy = .AbsoluteCopyNumber(rho=rho, psi=psi, gamma=1, x=sCN[,"SegmentedLog2Ratio"])
		integercopy = round(nonroundedcopy)
		calls3cat = rep(0, nrow(sCN))
		calls3cat[integercopy>round(psi)] = 1
		calls3cat[integercopy<round(psi)] = -1
		calls5cat = calls3cat
		calls5cat[integercopy>(round(psi)+2)] = 2
		calls5cat[integercopy<(round(psi)-1)] = -2
		if ("AbsoluteCopies" %in% colnames(sCN)) {
			sCN[,"AbsoluteCopies"] = nonroundedcopy
		} else {
			sCN = cbind(sCN, "AbsoluteCopies"=nonroundedcopy)
		}
		if ("IntegerCopies" %in% colnames(sCN)) {
			sCN[,"IntegerCopies"] = integercopy
		} else {
			sCN = cbind(sCN, "IntegerCopies"=integercopy)
		}
		if ("CalledData3Cat" %in% colnames(sCN)) {
			sCN[,"CalledData3Cat"] = calls3cat
		} else {
			sCN = cbind(sCN, "CalledData3Cat"=calls3cat)
		}
		if ("CalledData5Cat" %in% colnames(sCN)) {
			sCN[,"CalledData5Cat"] = calls5cat
		} else {
			sCN = cbind(sCN, "CalledData5Cat"=calls5cat)
		}
		save(sCN, tCN, file=paste0(path, "/resu/mad/", sampleNames[i], ".RData"))
		return(1)
	}
	close(pb)
	if (sum(unlist(res))==length(sampleNames)) {
		cat("\n")
	}
	return(invisible(as.logical(unlist(res))))
}