R/plotlog2andbafs.R

In ndbrown6/CNu: Copy Number Utilities

'plotlog2andbafs' <- function(path=NULL, platform="exome", 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/log2baf"))) {
		dir.create(paste0(path, "/resu/log2baf"))
	}
	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)
			}
		}
		load(paste0(path, "/resu/mad/", sampleNames[i], ".RData"))
		for (j in 1:nrow(sCN)) {
			index = which(tCN[,"Chromosome"]==sCN[j,"Chromosome"] & tCN[,"Position"]>=sCN[j,"Start"] & tCN[,"Position"]<=sCN[j,"End"] & tCN[,"Genotype"]==0)
			n = length(index)
			sCN[j,"N"] = sCN[j,"N"] - n
		}
		sCN = subset(sCN, sCN[,"N"]>0)
		tCN = subset(tCN, tCN[,"Genotype"]==1)
		sCN[,"N"] = cumsum(sCN[,"N"])
		col = rep("grey80", nrow(tCN))
		col[(tCN[,"Chromosome"] %% 2)==0] = "grey70"
		pdf(file=paste0(path, "/resu/log2baf/", sampleNames[i], ".pdf"), height=7*13/7, width=7*20/7)
		par(mar=c(5, 5, 4, 2)+.1, mfrow=c(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)
	
		plot(tCN[,"BAF"], type="p", pch=pch, cex=cex, col=col, axes=FALSE, frame=TRUE, xlab="", ylab="", main="", ylim=c(0,1))
		points(1-tCN[,"BAF"], type="p", pch=pch, cex=cex, col=col)
		for (j in 2:nrow(sCN)) {
			lines(x=c(sCN[j-1,"N"], sCN[j,"N"]), y=rep(sCN[j,"SegmentedBAF"],2), lty=1, lwd=4, col="red")
			lines(x=c(sCN[j-1,"N"], sCN[j,"N"]), y=1-rep(sCN[j,"SegmentedBAF"],2), lty=1, lwd=4, col="red")
		}
		lines(x=c(1, sCN[1,"N"]), y=rep(sCN[1,"SegmentedBAF"],2), lty=1, lwd=4, col="red")
		lines(x=c(1, sCN[1,"N"]), y=1-rep(sCN[1,"SegmentedBAF"],2), lty=1, lwd=4, 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 = "BAF", line = 3, cex = 1.5)
		abline(v=1, col="goldenrod3")
		abline(h=0.5, 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)
		dev.off()
		return(1)
	}
	close(pb)
	if (sum(unlist(res))==length(sampleNames)) {
		cat("\n")
	}
	return(invisible(as.logical(unlist(res))))
}