R/plotcopybychrom.R

In ndbrown6/CNu: Copy Number Utilities

'plotcopybychrom' <- 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/bychrom"))) {
		dir.create(paste0(path, "/resu/bychrom"))
	}
	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% {
		data(CytoBand)
		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]
		if (!file.exists(paste0(path, "/resu/bychrom/", sampleNames[i]))) {
			dir.create(paste0(path, "/resu/bychrom/", sampleNames[i]))
		}
		load(paste0(path, "/resu/mad/", sampleNames[i], ".RData"))
		sCN[,"N"] = cumsum(sCN[,"N"])
		for (j in 1:23) {
			tmp2 = subset(tCN, tCN[,"Chromosome"]==j)
			tmp3 = subset(sCN, sCN[,"Chromosome"]==j)
			pdf(file=paste0(path, "/resu/bychrom/", sampleNames[i], "/chrom_", j, ".pdf"), width=7, height=7)
			par(mar=c(5, 5, 4, 2)+.1)
			figs = split.screen(figs=matrix(c(0,1,0.1,1, 0,1,0,.3), nrow=2, ncol=4, byrow=TRUE))
			screen(figs[1])
			plot(tmp2[,"Position"], tmp2[,"Log2Ratio"], type="p", pch=pch, cex=cex, col="grey85", axes=FALSE, frame=TRUE, xlab="", ylab="", main="", ylim=c(-2,2), xlim=c(0,max(CytoBand[CytoBand[,"Chromosome"]==j,"End"])))
			for (k in 1:nrow(tmp3)) {
				lines(x=c(tmp3[k,"Start"], tmp3[k,"End"]), y=rep(tmp3[k,"SegmentedLog2Ratio"],2), lty=1, lwd=3.5, col="red")
			}
			for (k in 1:(nrow(tmp3)-1)) {
				lines(x=c(tmp3[k,"End"], tmp3[k+1,"Start"]), c(tmp3[k,"SegmentedLog2Ratio"],tmp3[k+1,"SegmentedLog2Ratio"]), lty=1, lwd=1, col="red")
			}
			axis(2, at=c(-2, -1, 0, 1, 2), labels=c(-2, -1, 0, 1, 2), cex.axis=1.25, las=1)
			abline(h=0, col="red")
			cen = max(CytoBand[CytoBand[,"Chromosome"]==j & CytoBand[,"Arm"]=="p","End"])
			abline(v=cen, col="goldenrod3", lty=2, lwd=2.5)
			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")
				}
			}
			mtext(side = 2, text = expression(Log[2]~"Ratio"), line = 3, cex = 1.5)
			box(lwd=1.5)
			screen(figs[2])
			plot(tmp2[,"Position"], tmp2[,"Log2Ratio"], type="n", axes=FALSE, frame=FALSE, xlab="", ylab="", ylim=c(0,2), xlim=c(0,max(CytoBand[CytoBand[,"Chromosome"]==j,"End"])))
			plotIdiogram(chromosome=j, build="hg19", new=FALSE, cex.axis=.6)
			points(cen, 1.1, pch=19, col="goldenrod3", cex=2)
			close.screen(all.screens=TRUE)
			dev.off()
		}
		return(1)
	}
	close(pb)
	if (sum(unlist(res))==length(sampleNames)) {
		cat("\n")
	}
	return(invisible(as.logical(unlist(res))))
}