R/plotbafhetero.R

In ndbrown6/CNu: Copy Number Utilities

'plotbafhetero' <- 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/bafh"))) {
		dir.create(paste0(path, "/resu/bafh"))
	}
	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"))
		tCN = subset(tCN, tCN[,"Genotype"]==1)
		col = rainbow_hcl(nrow(tCN), start=30, end=300)
		pdf(file=paste0(path, "/resu/bafh/", sampleNames[i], ".pdf"), height=7*10/7, width=7*20/7)
		par(mar=c(5, 5, 4, 2)+.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)
		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))))
}