ndbrown6/CNu
Copy Number Utilities

Package index
Search the ndbrown6/CNu package
Vignettes
Functions
12
Source code
12
Man pages
12
Home
/
GitHub
/
ndbrown6/CNu
/
R/zzz.R

R/zzz.R
In ndbrown6/CNu: Copy Number Utilities 

.CnuEnv <- new.env()
assign("cex_lE", .85, envir=.CnuEnv)
assign("pch_lE", ".", envir=.CnuEnv)
assign("cex_bE", 2.75, envir=.CnuEnv)
assign("pch_bE", ".", envir=.CnuEnv)

assign("cex_lT", .85, envir=.CnuEnv)
assign("pch_lT", 16, envir=.CnuEnv)
assign("cex_bT", .85, envir=.CnuEnv)
assign("pch_bT", 16, envir=.CnuEnv)


'.AscatDirectoryStructure' <- function(path) {
	dir.create(paste0(path, "/resu/ASCAT/"))
	dir.create(paste0(path, "/resu/ASCAT/Segmented/"))
	dir.create(paste0(path, "/resu/ASCAT/Results/"))
	dir.create(paste0(path, "/resu/ASCAT/ASCAT_Params/"))
}

'.RunAscatAuto' <- function(path) {
	ASCAT_Params = read.csv(file=paste0(path, "/resu/ASCAT/ASCAT_Params/ASCAT_Params.csv"), header=TRUE, sep=",", stringsAsFactors=FALSE)
	sampleNames = ASCAT_Params[,"sample_name"]
	resu = rep(0, length(sampleNames))
	pb = txtProgressBar(min=1, max=length(sampleNames), style=3)
	for (i in 1:length(sampleNames)) {
		setTxtProgressBar(pb, i)
		if (!dir.exists(paste0(path, "/resu/ASCAT/Results/", sampleNames[i]))) {
			dir.create(paste0(path, "/resu/ASCAT/Results/", sampleNames[i]))
		}
		load(paste0(path, "/resu/ASCAT/Segmented/", sampleNames[i], ".RData"))
		tmp2 = try(runASCAT(lrr=tmp$Tumor_LogR,
						baf=tmp$Tumor_BAF,
						lrrsegmented=tmp$Tumor_LogR_segmented,
						bafsegmented=tmp$Tumor_BAF_segmented,
				 		gender=tmp$gender,
				 		SNPpos=tmp$SNPpos,
				 		chromosomes=tmp$chromosomes,
				 		chrnames=tmp$chrnames,
				 		sexchromosomes=tmp$sexchromosomes,
				 		failedqualitycheck=FALSE,
				 		distance = paste0(path, "/resu/ASCAT/Results/", sampleNames[i], "/Heatmap_wfit.pdf"),
				 		copynumberprofile = NULL,
				 		nonroundedprofile = NULL, 
    			 		aberrationreliability = NULL,
    			 		gamma = 1, rho_manual = NA, psi_manual = NA, y_limit = 3, circos = NA))
    	if (!("try-error" %in% is(tmp2))) {
    		save(tmp, tmp2, sCN, tCN, file=paste0(path, "/resu/ASCAT/Results/", sampleNames[i], "/CN_and_Genotype.RData"))
	    	ASCAT_Params[i,"rho"] = tmp2$rho
	    	ASCAT_Params[i,"psi"] = tmp2$psi
	    	ASCAT_Params[i,"purity"] = tmp2$rho
	    	ASCAT_Params[i,"ploidy"] = sum((tmp2$seg_raw[,"nAraw"]+tmp2$seg_raw[,"nBraw"])*(tmp$SNPpos[tmp2$seg_raw[,"end"],"pos"]-tmp$SNPpos[tmp2$seg_raw[,"start"],"pos"])/sum(as.numeric(tmp$SNPpos[tmp2$seg_raw[,"end"],"pos"]-tmp$SNPpos[tmp2$seg_raw[,"start"],"pos"])))
    		write.table(ASCAT_Params, file=paste0(path, "/resu/ASCAT/ASCAT_Params/ASCAT_Params.csv"), col.names=TRUE, row.names=FALSE, sep=",", quote=FALSE)
    		resu[i] = 1
    	}
    }
    close(pb)
    return(invisible(resu))
}

'.RunAscatNoAuto' <- function(path) {
	ASCAT_Params = read.csv(file=paste0(path, "/resu/ASCAT/ASCAT_Params/ASCAT_Params.csv"), header=TRUE, sep=",", stringsAsFactors=FALSE)
	sampleNames = ASCAT_Params[,"sample_name"]
	resu = rep(0, length(sampleNames))
	pb = txtProgressBar(min=1, max=length(sampleNames), style=3)
	for (i in 1:length(sampleNames)) {
		setTxtProgressBar(pb, i)
		if (!dir.exists(paste0(path, "/resu/ASCAT/Results/", sampleNames[i]))) {
			dir.create(paste0(path, "/resu/ASCAT/Results/", sampleNames[i]))
		}
		load(paste0(path, "/resu/ASCAT/Segmented/", sampleNames[i], ".RData"))
		tmp2 = try(runASCAT(lrr=tmp$Tumor_LogR,
						baf=tmp$Tumor_BAF,
						lrrsegmented=tmp$Tumor_LogR_segmented,
						bafsegmented=tmp$Tumor_BAF_segmented,
				 		gender=tmp$gender,
				 		SNPpos=tmp$SNPpos,
				 		chromosomes=tmp$chromosomes,
				 		chrnames=tmp$chrnames,
				 		sexchromosomes=tmp$sexchromosomes,
				 		failedqualitycheck=FALSE,
				 		distance = paste0(path, "/resu/ASCAT/Results/", sampleNames[i], "/Heatmap_wfit.pdf"),
				 		copynumberprofile = NULL,
				 		nonroundedprofile = NULL, 
    			 		aberrationreliability = NULL,
    			 		gamma = 1, rho_manual = ASCAT_Params[i,"rho"], psi_manual = ASCAT_Params[i,"psi"], y_limit = 3, circos = NA))
    	if (!("try-error" %in% is(tmp2))) {
 		   	save(tmp, tmp2, sCN, tCN, file=paste0(path, "/resu/ASCAT/Results/", sampleNames[i], "/CN_and_Genotype.RData"))
		    ASCAT_Params[i,"purity"] = tmp2$rho
 	    	ASCAT_Params[i,"ploidy"] = sum((tmp2$seg_raw[,"nAraw"]+tmp2$seg_raw[,"nBraw"])*(tmp$SNPpos[tmp2$seg_raw[,"end"],"pos"]-tmp$SNPpos[tmp2$seg_raw[,"start"],"pos"])/sum(as.numeric(tmp$SNPpos[tmp2$seg_raw[,"end"],"pos"]-tmp$SNPpos[tmp2$seg_raw[,"start"],"pos"])))
	    	write.table(ASCAT_Params, file=paste0(path, "/resu/ASCAT/ASCAT_Params/ASCAT_Params.csv"), col.names=TRUE, row.names=FALSE, sep=",", quote=FALSE)
	    	resu[i] = 1
	    }
    }
    close(pb)
    return(invisible(resu))
}

'.GapDirectoryStructure' <- function(path) {
	dir.create(paste0(path, "/resu/GAP/"))
	dir.create(paste0(path, "/resu/GAP/Raw/"))
	dir.create(paste0(path, "/resu/GAP/Segmented/"))
	dir.create(paste0(path, "/resu/GAP/Results/"))
	dir.create(paste0(path, "/resu/GAP/GAP_Params/"))
}

'.RunGapAuto' <- function(path) {
	GAP_Params = read.csv(file=paste0(path, "/resu/GAP/GAP_Params/GAP_Params.csv"), header=TRUE, sep=",", row.names=1, stringsAsFactors=FALSE)
	sampleNames = rownames(GAP_Params)
	resu = rep(0, length(sampleNames))
	pb = txtProgressBar(min=1, max=length(sampleNames), style=3)
	for (i in 1:length(sampleNames)) {
		setTxtProgressBar(pb, i)
		if (!dir.exists(paste0(path, "/resu/GAP/Results/", sampleNames[i]))) {
			dir.create(paste0(path, "/resu/GAP/Results/", sampleNames[i]))
		}
		load(paste0(path, "/resu/GAP/Segmented/", sampleNames[i], ".RData"))
		tmp = CBS_breakpoints_and_GAPconstruction(datL=Tumor_LogR_segmented, datA=Tumor_BAF_segmented, SNP_annot=SNP_annot, germHomoThr=0.99)
		notNAs = sort(intersect(which(!is.na(tmp[,GAP::.GapEnv$c_lrr]+tmp[,GAP::.GapEnv$c_baf])),which(tmp[,GAP::.GapEnv$c_chr]<23)))
		resuWfit = try(CopyNumber_and_Genotype(tmp[notNAs,GAP::.GapEnv$c_lrr],tmp[notNAs,GAP::.GapEnv$c_baf],tmp[notNAs,GAP::.GapEnv$c_len],showGAPTemplate=FALSE))
		if (("try-error" %in% is(resuWfit))) {
			resuWfit = try(CopyNumber_and_Genotype(tmp[notNAs,GAP::.GapEnv$c_lrr]+rnorm(length(notNAs), 0, .01),tmp[notNAs,GAP::.GapEnv$c_baf]+abs(rnorm(length(notNAs), 0, .01)),tmp[notNAs,GAP::.GapEnv$c_len],showGAPTemplate=FALSE))
		}
		if (!("try-error" %in% is(resuWfit))) {
			GAP_Params[sampleNames[i],"p_BAF"] = resuWfit$coeff[1]
			GAP_Params[sampleNames[i],"q_LRR"] = resuWfit$coeff[2]
			GAP_Params[sampleNames[i],"X2copy_LRR"] = resuWfit$Centr_CL[5,3]
			
			pdf(file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Preliminary_pattern_wfit.pdf"))
			GAPplot_preliminary(tmp, resuWfit$Centr_CL, resuWfit$coeff[1], resuWfit$coeff[2], resuWfit$NETM, sampleName=sampleNames[i])
			dev.off()
			tmp2 = cbind(tmp, matrix(0,dim(tmp)[1],16))
			colnames(tmp2)[11:26] = GAP::.GapEnv$ColNamesTMP[11:26]
			tmp2[,GAP::.GapEnv$c_nC] = getCNwithPrecision(tmp2[,GAP::.GapEnv$c_lrr], resuWfit$Centr_CL,Precision=1/8)
			tmp2 = getGTwithRank(tmp2,resuWfit$Centr_CL,resuWfit$NETM)
			tmp2[,GAP::.GapEnv$c_CN] = tmp2[,GAP::.GapEnv$c_CN1]
			tmp2[,GAP::.GapEnv$c_BA] = tmp2[,GAP::.GapEnv$c_BA1]
			tmp3 = getConfidence(tmp2)
			tmp3 = shrinkReprTMP(tmp3)
			tmp3 = updateMeans_FinalPolish(tmp3, Tumor_BAF_segmented, Tumor_LogR_segmented, germHomozyg.mBAF.thr=0.99)
			PosSNPstart = SNP_annot[tmp3[,2],4]
			PosSNPend = SNP_annot[tmp3[,3],4]
			tmp4 = cbind(tmp3[,1:3], PosSNPstart, PosSNPend, tmp3[,4:ncol(tmp3)])
			tmp4[,"Chr"] = floor(tmp4[,"Chr"])
			tmp4 = tmp4[,c("Chr", "PosSNPstart", "PosSNPend", "LRR", "MajorBAF", "Copy", "Mallele"),drop=FALSE]

			if (!dir.exists(paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))) {
				dir.create(path=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))
			}
			load(paste0(path, "/resu/GAP/Raw/", sampleNames[i], ".RData"))
			wD = getwd()
			setwd(paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))
			summaryPlots_perChrom(tmp3, sampleName=sampleNames[i], SNP_annot, Tumor_BAF, Tumor_LogR, resuWfit$Centr_CL)
			setwd(wD)
			pdf(file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Pattern_Final.pdf"))
			GAPplot(tmp3, sampleName=sampleNames[i], SNP_annot, resuWfit$Centr_CL, resuWfit$NETM, resuWfit$coeff[1], resuWfit$coeff[2])
			dev.off()
			jpeg(filename=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Summary.jpeg"), width=1920, height=2400, pointsize=24)
			info = summaryPlots(tmp3, sampleName=sampleNames[i], SNP_annot=SNP_annot, datA=Tumor_BAF, datL=Tumor_LogR, resuWfit$Centr_CL, resuWfit$NETM, p_BAF=as.numeric(GAP_Params[sampleNames[i],"p_BAF"]), q_LRR=as.numeric(GAP_Params[sampleNames[i],"q_LRR"]))
			dev.off()
			## info = summaryStats(tmp3, SNP_annot=SNP_annot)
			GAP_Params[sampleNames[i], 4:7] = as.numeric(info)[5:8]
			write.csv(GAP_Params, file=paste0(path, "/resu/GAP/GAP_Params/GAP_Params.csv"))
			save(tmp, tmp2, tmp3, tmp4, notNAs, resuWfit, file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/CopyNumber_and_Genotype.RData"))
			resu[i] = 1
		}
	}
	close(pb)
	return(invisible(resu))
}

'.RunGapNoAuto' <- function(path) {
	GAP_Params = read.csv(file=paste0(path, "/resu/GAP/GAP_Params/GAP_Params.csv"), header=TRUE, sep=",", row.names=1, stringsAsFactors=FALSE)
	sampleNames = rownames(GAP_Params)
	resu = rep(0, length(sampleNames))
	pb = txtProgressBar(min=1, max=length(sampleNames), style=3)
	for (i in 1:length(sampleNames)) {
		setTxtProgressBar(pb, i)
		if (!dir.exists(paste0(path, "/resu/GAP/Results/", sampleNames[i]))) {
			dir.create(paste0(path, "/resu/GAP/Results/", sampleNames[i]))
		}
		load(paste0(path, "/resu/GAP/Segmented/", sampleNames[i], ".RData", sep=""))
		load(paste0(path, "/resu/GAP/Results/", sampleNames[i], "/CopyNumber_and_Genotype.RData"))
		resuWOfit = try(CopyNumber_and_GenotypeWOfitting(tmp[notNAs,GAP::.GapEnv$c_lrr], tmp[notNAs,GAP::.GapEnv$c_baf], tmp[notNAs,GAP::.GapEnv$c_len], showGAPTemplate=FALSE, as.numeric(GAP_Params[sampleNames[i],"p_BAF"]), as.numeric(GAP_Params[sampleNames[i],"q_LRR"]), c(2,as.numeric(GAP_Params[sampleNames[i],"X2copy_LRR"]))))
		if (("try-error" %in% is(resuWOfit))) {
			resuWOfit = try(CopyNumber_and_GenotypeWOfitting(tmp[notNAs,GAP::.GapEnv$c_lrr]+rnorm(length(notNAs), 0, .01),tmp[notNAs,GAP::.GapEnv$c_baf]+abs(rnorm(length(notNAs), 0, .01)),tmp[notNAs,GAP::.GapEnv$c_len],showGAPTemplate=FALSE,as.numeric(GAP_Params[sampleNames[i],"p_BAF"]), as.numeric(GAP_Params[sampleNames[i],"q_LRR"]), c(2,as.numeric(GAP_Params[sampleNames[i],"X2copy_LRR"]))))
		}
		if (!("try-error" %in% is(resuWOfit))) {
			pdf(file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Preliminary_pattern_wofit.pdf"))
			GAPplot_preliminary(tmp, resuWOfit$Centr_CL, resuWOfit$coeff[1], resuWOfit$coeff[2], resuWOfit$NETM, sampleName=sampleNames[i])
			dev.off()
			tmp2 = cbind(tmp, matrix(0,dim(tmp)[1],16))
			colnames(tmp2)[11:26] = GAP::.GapEnv$ColNamesTMP[11:26]
			tmp2[,GAP::.GapEnv$c_nC] = getCNwithPrecision(tmp2[,GAP::.GapEnv$c_lrr], resuWOfit$Centr_CL,Precision=1/8)
			tmp2 = getGTwithRank(tmp2,resuWOfit$Centr_CL,resuWOfit$NETM)
			tmp2[,GAP::.GapEnv$c_CN] = tmp2[,GAP::.GapEnv$c_CN1]
			tmp2[,GAP::.GapEnv$c_BA] = tmp2[,GAP::.GapEnv$c_BA1]
			tmp3 = getConfidence(tmp2)
			tmp3 = shrinkReprTMP(tmp3)
			tmp3 = updateMeans_FinalPolish(tmp3, Tumor_BAF_segmented, Tumor_LogR_segmented, germHomozyg.mBAF.thr=0.99)
			PosSNPstart = SNP_annot[tmp3[,2],4]
			PosSNPend = SNP_annot[tmp3[,3],4]
			tmp4 = cbind(tmp3[,1:3], PosSNPstart, PosSNPend, tmp3[,4:ncol(tmp3)])
			tmp4[,"Chr"] = floor(tmp4[,"Chr"])
			tmp4 = tmp4[,c("Chr", "PosSNPstart", "PosSNPend", "LRR", "MajorBAF", "Copy", "Mallele"),drop=FALSE]

			if (!dir.exists(paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))) {
				dir.create(path=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))
			}
			load(paste(path, "/resu/GAP/Raw/", sampleNames[i], ".RData", sep=""))
			wD = getwd()
			setwd(paste0(path, "/resu/GAP/Results/", sampleNames[i], "/Summary_per_Chrom/"))
			summaryPlots_perChrom(tmp3, sampleName=sampleNames[i], SNP_annot, Tumor_BAF, Tumor_LogR, resuWOfit$Centr_CL)
			setwd(wD)
			pdf(file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Pattern_Final.pdf"))
			GAPplot(tmp3, sampleName=sampleNames[i], SNP_annot, resuWOfit$Centr_CL, resuWOfit$NETM, resuWOfit$coeff[1], resuWOfit$coeff[2])
			dev.off()
			jpeg(filename=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/GAP_Summary.jpeg"), width=1920, height=2400, pointsize=24)
			info = summaryPlots(tmp3, sampleName=sampleNames[i], SNP_annot=SNP_annot, datA=Tumor_BAF, datL=Tumor_LogR, resuWOfit$Centr_CL, resuWOfit$NETM, p_BAF=as.numeric(GAP_Params[sampleNames[i],"p_BAF"]), q_LRR=as.numeric(GAP_Params[sampleNames[i],"q_LRR"]))
			dev.off()
			## info = summaryStats(tmp3, SNP_annot=SNP_annot)
			GAP_Params[sampleNames[i], 4:7] = as.numeric(info)[5:8]
			write.csv(GAP_Params, file=paste0(path, "/resu/GAP/GAP_Params/GAP_Params.csv"))
			save(tmp, tmp2, tmp3, tmp4, notNAs, resuWOfit, file=paste0(path, "/resu/GAP/Results/", sampleNames[i], "/CopyNumber_and_Genotype.RData"))
			resu[i] = 1
		}
	}
	close(pb)
	return(invisible(resu))
}

'.AbsoluteCopyNumber' <- function(rho, psi, gamma=.55, x) {
	return(invisible(((((2^(x/gamma))*(rho*psi+(1-rho)*2)) - ((1-rho)*2))/rho)))
}
ndbrown6/CNu documentation built on May 27, 2019, 1:09 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close