R/plotcopynumber.R

In ndbrown6/CNu: Copy Number Utilities

'plotcopynumber' <- function(path=NULL, platform="exome", 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 (!dir.exists(paste0(path, "/resu/mad"))) {
		stop("directory 'mad' absent")
	}
	if (!dir.exists(paste0(path, "/resu/copy"))) {
		dir.create(paste0(path, "/resu/copy"))
	}
	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
		}
	}
	if (platform=="exome") {
		pch = CNu::.CnuEnv$pch_bE
		cex = CNu::.CnuEnv$cex_bE
	} else if (platform=="targeted") {
		pch = CNu::.CnuEnv$pch_bT
		cex = CNu::.CnuEnv$cex_bT
	}
	registerDoMC(ncores)
	sampleNames = gsub(pattern=".RData", replacement="", x=dir(path=paste0(path, "/resu/mad/"), pattern=".RData", full.names=FALSE), fixed=TRUE)
	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(paste0(path, "/resu/mad/", sampleNames[i], ".RData"))
		sCN[,"N"] = cumsum(sCN[,"N"])
		col = rep("grey80", nrow(tCN))
		col[(tCN[,"Chromosome"] %% 2)==0] = "grey70"
		pdf(file=paste0(path, "/resu/copy/", sampleNames[i], ".pdf"), height=10, width=20)
		par(mar=c(5, 5, 4, 2)+.1)
		plot(tCN[,"Log2Ratio"], type="p", pch=pch, cex=cex, col=col, axes=FALSE, frame=TRUE, xlab="", ylab="", main="", ylim=c(-2,2))
		for (j in 2:nrow(sCN)) {
			lines(x=c(sCN[j-1,"N"], sCN[j,"N"]), y=rep(sCN[j,"SegmentedLog2Ratio"],2), lty=1, lwd=5, col="red")
		}
		lines(x=c(1, sCN[1,"N"]), y=rep(sCN[1,"SegmentedLog2Ratio"],2), lty=1, lwd=5, col="red")
		axis(2, at = NULL, cex.axis = 1.5, las = 1)
		mtext(side = 1, text = "Chromosome", line = 3, cex = 1.5)
		mtext(side = 2, text = expression(Log[2]~"Ratio"), line = 3, cex = 1.5)
		abline(v=1, col="goldenrod3")
		abline(h=0, col="red")
		for (k in 1:23) {
			v = max(which(tCN[,"Chromosome"]==k & tCN[,"Position"]>1))
			abline(v=v, col="goldenrod3")
		}
		start = NULL
		end = NULL
		for (k in 1:23) {
			start = c(start, min(which(tCN[,"Chromosome"]==k & tCN[,"Position"]>1)))
			end = c(end, max(which(tCN[,"Chromosome"]==k & tCN[,"Position"]>1)))
		}
		axis(1, at = .5*(start+end), labels=c(1:22,"X"), cex.axis = 1.5, las = 1)
		dipP = log2((rho*2 + (1-rho)*2)/(rho*psi + (1-rho)*2))
		dip0 = mean(sCN[sCN[,"IntegerCopies"]==2,"SegmentedLog2Ratio"])
		for (k in 1:10) {
			if (dipP>dip0) {
				abline(h=log2((rho*k + (1-rho)*2)/(rho*psi + (1-rho)*2))-(dipP-dip0), col="darkorange", lty=3)
				mtext(text=k, side=4, line=.5, at=log2((rho*k + (1-rho)*2)/(rho*psi + (1-rho)*2))-(dipP-dip0), las=2, cex=.75, col="orange")
			} else {
				abline(h=log2((rho*k + (1-rho)*2)/(rho*psi + (1-rho)*2))+(dip0-dipP), col="darkorange", lty=3)
				mtext(text=k, side=4, line=.5, at=log2((rho*k + (1-rho)*2)/(rho*psi + (1-rho)*2))+(dip0-dipP), las=2, cex=.75, col="orange")
			}
		}
		dev.off()
		return(1)
	}
	close(pb)
	if (sum(unlist(res))==length(sampleNames)) {
		cat("\n")
	}
	return(invisible(as.logical(unlist(res))))
}