R/plottings.R

In r3fang/SnapATAC: Single Nucleus Analysis Package for ATAC-Seq

#' Plot Barocde Quality Control Distribution
#'
#' @param obj a snap object
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @param ... Arguments passed to hist method.
#' @examples 
#' data(demo.sp);
#' plotBarcode(demo.sp, col="grey", border="grey");
#' @importFrom graphics par hist 
#' @importFrom grDevices pdf dev.off
#' @export
plotBarcode <- function(obj, pdf.file.name, pdf.width, pdf.height, ...) {
  UseMethod("plotBarcode", obj);
}

#' @export
plotBarcode.default <- function(
	obj,
	pdf.file.name=NULL,
	pdf.width=7,
	pdf.height=7,
	...
	){	

	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
		if(nrow(obj@metaData) == 0){stop(paste("obj deos not have meta data", sep=""))}
		# check the input
		barcode = obj@metaData;
	}
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
			
	par(mfrow=c(3,2));
	hist(log(barcode$TN + 1,10), main="log10(Total Fragments)", xlab="log10(Total Fragments + 1)", ...);
	hist(log(barcode$UQ + 1,10), main="log10(UMI)", xlab="log10(UMI + 1)", ...);
	hist((barcode$UM+1)/(barcode$TN+1), main="Mappability Ratio", xlab="Mappability", xlim=c(0, 1), ...);
	hist((barcode$PP+1)/(barcode$UM+1), main="Proper Paired Ratio", xlab="Proper Paired",xlim=c(0, 1), ...);
	hist(1 - (barcode$UQ+1)/(barcode$PP+1), main="Duplicate Rate Ratio", xlab="Duplicate Rate", xlim=c(0, 1), ...);
	hist((barcode$CM+1) / (barcode$UQ+1), main="chrM Rate Ratio", xlab="chrM Rate", xlim=c(0, 1), ...);	
	if(!is.null(pdf.file.name)){
		dev.off()		
	}
	graphics::par(mfrow=c(1,1));
}

#' Plot Bin Coverage Distribution
#'
#' @param obj a snap object
#' @param rm.zeros Remove bins of zero coverage when ploting the coverage distribution [TRUE].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @param ... Arguments passed to hist method.
#' @examples 
#' data(demo.sp);
#' plotBinCoverage(demo.sp, col="grey", border="grey");
#' @importFrom stats sd
#' @importFrom grDevices pdf dev.off
#' @importFrom graphics hist
#' @export
plotBinCoverage <- function(obj, rm.zeros, pdf.file.name, pdf.width, pdf.height, ...) {
  UseMethod("plotBinCoverage", obj);
}

#' @export
plotBinCoverage.default <- function(
	obj, 
	rm.zeros=TRUE,
	pdf.file.name=NULL,
	pdf.width=7, 
	pdf.height=7, 
	...
){	
	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
		data.use = obj@bmat;
		
		if((x=nrow(data.use)) == 0L){
			stop("cell-by-bin matrix is empty, run addBmatToSnap first")
		}
	}
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		grDevices::pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	
	cov = Matrix::colSums(data.use);
	idy = seq(cov);	
	cov = cov[which(cov > 0)];	
	cov = log10(cov + 1);
	cov = (cov - mean(cov)) / sd(cov);
	hist(cov, ...);
	
	if(!is.null(pdf.file.name)){
		grDevices::dev.off()		
	}
	graphics::par(mfrow=c(1,1));
}

#' Visulization
#'
#' @param obj A snap object.
#' @param method Visulization method c("tsne", "umap").
#' @param point.size Point size [1].
#' @param point.shape Point shape type [19].
#' @param point.alpha Point transparancy level [0.8].
#' @param point.color Color of point. Two options, points color by 
#' cell cluster label or sample ID c("cluster", "sample"). 
#' @param text.add Whether to add cluster label text at the centroid of each cluster [TRUE].
#' @param text.size Cluster label text size [5].
#' @param text.color Cluster label text color ["black"].
#' @param text.halo.add Add halo to cluster label text [TRUE].
#' @param text.halo.color Halo color ["white"].
#' @param text.halo.width Halo width [0.5].
#' @param legend.add Add a legend to the plot [FALSE].
#' @param legend.pos Position of the legend c("bottomleft", "bottom", "left", "topleft", "top", "topright", "right", "center").
#' @param legend.text.size Size of the text in the legend [1].
#' @param legend.text.color Color of the text in the legend ["black"].
#' @param down.sample Downsample the original cells to down.sample cells to ovoid large dataset [10,000].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @param ... Arguments passed to plot method.
#'
#' @examples 
#' data(demo.sp);
#' plotViz(
#' 	obj=demo.sp, 
#' 	method="tsne", 
#' 	point.color="sample", 
#' 	text.add=FALSE,
#'	legend.add=TRUE
#' );
#' 
#' @importFrom grDevices pdf dev.off
#' @importFrom methods slot
#' @importFrom scales alpha
#' @importFrom graphics plot text title legend
#' @export
plotViz <- function(obj, 
	method, 
	point.size, 
	point.shape, 
	point.alpha, 
	point.color, 
	text.add, 
	text.size, 
	text.color, 
	text.halo.add, 
	text.halo.color, 
	text.halo.width, 
	legend.add, 
	legend.pos, 
	legend.text.size,
	legend.text.color,
	down.sample, 
	pdf.file.name, 
	pdf.width, 
	pdf.height, 
	...
){
  UseMethod("plotViz", obj);
}

#' @export
plotViz.default <- function(obj, 
		method=c("tsne", "umap"), 
		point.size=1, 
		point.shape=19, 
		point.alpha=0.8, 
		point.color=NULL,
		text.add=TRUE,
		text.size=1, 
		text.color="black",
		text.halo.add=TRUE,
		text.halo.color="white",
		text.halo.width=0.2,
		legend.add=FALSE,
		legend.pos=c("bottomleft", "bottom", "left", "topleft", "top", "topright", "right", "center"),
		legend.text.size=1,
		legend.text.color="black",
		down.sample=10000,
		pdf.file.name=NULL,
		pdf.width=7, 
		pdf.height=7,
		...
){	
	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
	}
		
	if(is.integer(down.sample)){
		stop("down.sample must be an integer");
	}
	
	method = match.arg(method);
	data.use = as.data.frame(slot(obj, method));
	if(method=="tsne"){
		colnames(data.use) = c("TSNE-1", "TSNE-2")
	}else{
		colnames(data.use) = c("UMAP-1", "UMAP-2")
	}
	if((x=nrow(data.use)) == 0L){
		stop("visulization method does not exist, run runViz first!")
	}
	xlims = c(-max(abs(data.use[,1])) * 1.05, max(abs(data.use[,1])) * 1.2);
	ylims = c(-max(abs(data.use[,2])) * 1.05, max(abs(data.use[,2])) * 1.05);
	
	cluster = point.color;
	if(((x=length(cluster)) == 0L) | (is.null(point.color))){
		warning("cluster does not exist, text.add is ignored")
		text.add = FALSE;
	}
	
	if(length(cluster) != 0L){	
		data.use$col = factor(cluster);
	}else{
		data.use$col = factor(1);
	}
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	legend.pos = match.arg(legend.pos);
	down.sample = min(down.sample, ncell);
	idx.ds = sort(sample(seq(ncell), down.sample));
	data.use = data.use[idx.ds,,drop=FALSE]
										
	colPanel = createColorPanel(length(unique(data.use$col)));
	graphics::plot(
				   data.use[,c(1,2)],
		 		   cex=point.size, 
		 		   pch=point.shape, 
		 		   col=scales::alpha(colPanel[factor(data.use$col)], point.alpha),
				   bty="l",
				   font.lab=2,
				   col.axis = 'darkgrey',
				   xlim=xlims,
				   ylim=ylims,
		 		   ...
				   );
  	box(bty="l", lwd=2)
	if(text.add){
		xx = findCentrod(data.use[,c(1,2)], data.use$col);
		textHalo(x=xx[,1], y=xx[,2], labels = xx[,3], col=text.color, bg=text.halo.color, r=text.halo.width, cex=text.size);
  	}
	
	if(legend.add){
		legend.ncol = as.integer(length(levels(data.use$col)) / 25) +1;
		legend(
		  "topright", 
		  legend = levels(factor(data.use$col)), 
		  col = colPanel,
		  pch = point.shape,
		  pt.cex=1, 
		  bty = "n", 
		  cex = legend.text.size, 
		  text.col = legend.text.color,
		  horiz = FALSE,
		  ncol=legend.ncol
		  )		
	}

	if(!is.null(pdf.file.name)){
		dev.off()		
	}
}

#' plotFeatureSingle
#'
#' @param obj A snap object.
#' @param feature.value Feature enrichment value for each cell. Value will be normalized betweeen 0 and 1.
#' @param method Visulization method c("tsne", "umap").
#' @param point.size Point size [1].
#' @param point.shape Point shape type [19].
#' @param down.sample Downsample the original cells to down.sample cells to ovoid large dataset [10,000].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width Width of the graphics region in inches [7].
#' @param pdf.height Height of the graphics region in inches [7].
#' @param quantiles Feature value outside this range will be removed [c(0.01, 0.99)]
#' @param ... Arguments passed to plot method.
#'
#' @examples
#' data(demo.sp);
#' x = rexp(nrow(demo.sp));
#' x = (x-min(x))/(max(x)-min(x));
#' PlotFeatureSingle(
#' 	obj=demo.sp, 
#' 	feature.value=x,
#' 	method="tsne", 
#' 	point.size=1, 
#' 	point.shape=19, 
#' 	down.sample=10000, 
#' 	pdf.file.name=NULL, 
#' 	pdf.width=7, 
#' 	pdf.height=7
#' )
#' 
#' @importFrom grDevices pdf dev.off
#' @importFrom methods slot
#' @importFrom scales alpha
#' @importFrom graphics plot text title legend
#' @importFrom plot3D scatter2D
#' @importFrom viridis viridis
#' @export
plotFeatureSingle <- function(
	obj, 
	feature.value,
	method, 
	point.size, 
	point.shape, 
	down.sample, 
	pdf.file.name, 
	pdf.width, 
	pdf.height, 
	quantiles,
	...
){
  UseMethod("plotFeatureSingle", obj);
}

#' @export
plotFeatureSingle.default <- function(
	obj, 
	feature.value,
	method=c("tsne", "umap"), 
	point.size=1, 
	point.shape=19, 
	down.sample=10000,
	pdf.file.name=NULL,
	pdf.width=7, 
	pdf.height=7,
	quantiles=c(0.01, 0.99),
	...
){	
	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
	}
	
	if(missing(feature.value)){
		stop("feature.value is missing")
	}
	
	if(is.integer(down.sample)){
		stop("down.sample must be an integer");
	}
	
	method = match.arg(method);
	data.use = as.data.frame(slot(obj, method));
	
	if(method=="tsne"){
		colnames(data.use) = c("TSNE-1", "TSNE-2")
		lab.names = c("TSNE-1", "TSNE-2");
	}else{
		colnames(data.use) = c("UMAP-1", "UMAP-2")
		lab.names = c("UMAP-1", "UMAP-2");
	}
	
	if((x=nrow(data.use)) == 0L){
		stop("visulization method does not exist, run runViz first!")
	}
	
	if(nrow(data.use) != length(feature.value)){
		stop("feature.value has different length with number of cells in obj")
	}
	
	quantiles.low = quantile(feature.value, quantiles[1]);
	quantiles.high = quantile(feature.value, quantiles[2]);
	feature.value[feature.value > quantiles.high] = quantiles.high;
	feature.value[feature.value < quantiles.low ] = quantiles.low;
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	down.sample = min(down.sample, ncell);
	idx.ds = sort(sample(seq(ncell), down.sample));
	data.use = data.use[idx.ds,,drop=FALSE];
	feature.value = feature.value[idx.ds];
	xlims = c(-max(abs(data.use[,1])) * 1.05, max(abs(data.use[,1])) * 1.2);
	ylims = c(-max(abs(data.use[,2])) * 1.05, max(abs(data.use[,2])) * 1.05);
	
	plot3D::scatter2D(
				   x=data.use[,1],
				   y=data.use[,2],
				   colvar=feature.value,
		 		   cex=point.size, 
		 		   pch=point.shape, 
				   bty="l",
				   font.lab=2,
				   col.axis = 'darkgrey',
				   xlim=xlims, 
				   ylim=ylims, 
				   xlab=lab.names[1],
				   ylab=lab.names[2],
				   col=viridis(256, option = "D"),
		 		   ...
				   );
  	box(bty="l", lwd=2)
										
	if(!is.null(pdf.file.name)){
		dev.off()		
	}
}

#' Elbow Plot for Dimentionality Reduction Result
#'
#' @param obj A snap object
#' @param point.size Point size [2].
#' @param point.shape Point shape type [19].
#' @param point.color Point color ["red"].
#' @param point.alpha Point transparancy level [1].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @param ... Arguments passed to plot function. 
#'
#' @examples 
#' data(demo.sp);
#' plotDimReductElbow(demo.sp);
#' 
#' @importFrom grDevices pdf dev.off
#' @importFrom methods slot
#' @importFrom scales alpha
#' @importFrom graphics plot text title
#' @export
plotDimReductElbow <- function(obj, point.size, point.shape, point.color, point.alpha, pdf.file.name, pdf.height, pdf.width, ...){
  UseMethod("plotDimReductElbow", obj);
}

#' @export
plotDimReductElbow.default <- function(
	obj, 
	point.size=1.5,
	point.shape=19,
	point.color="red",
	point.alpha=1,
	pdf.file.name=NULL,
	pdf.height=7,
	pdf.width=7,
	...
){
	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
		if(!isDimReductComplete(obj@smat)){
			stop("obj does not have valid dim.reduct object, run 'runDimReduct' first");			
		}
	}

	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	data.use = data.frame(PC=1:length(obj@smat@sdev), sd=obj@smat@sdev);	
	plot(x=data.use[,1], y=data.use[,2], cex=point.size, pch=point.shape, col=alpha(point.color, point.alpha), xlab="PCs", ylab="Standard Deviation of PCs");
		
	if(!is.null(pdf.file.name)){
		dev.off()
	}
}


#' Pairwise plot for Dimentionality Reduction Result
#'
#' @param obj A snap object
#' @param eigs.dims Eigenvectors to plot [1:30]
#' @param point.size Point size [2].
#' @param point.shape Point shape type [19].
#' @param point.color Point color ["grey"].
#' @param point.alpha Point transparancy level [1].
#' @param down.sample Number of cells to plot. down.sample cells will be randomly downsampled for plot if there are more cells.
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' 
#' @examples 
#' data(demo.sp);
#' plotDimReductPW(demo.sp, eigs.dims=1:10);
#' 
#' @importFrom grDevices pdf dev.off
#' @importFrom methods slot
#' @importFrom scales alpha
#' @importFrom graphics plot text title mtext
#' @export
plotDimReductPW <- function(obj, eigs.dims, point.size, point.color, point.shape, point.alpha, down.sample, pdf.file.name, pdf.height, pdf.width){
  UseMethod("plotDimReductPW", obj);
}

#' @export
plotDimReductPW.default <- function(
	obj, 
	eigs.dims=1:50,
	point.size=0.5,
	point.color="grey",
	point.shape=19,
	point.alpha=0.5,
	down.sample=3000,
	pdf.file.name=NULL, 
	pdf.height=7, 
	pdf.width=7
){
	
	if(missing(obj)){
		stop("obj is missing");
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object");
		}
		ncell = nrow(obj);
		if(!isDimReductComplete(obj@smat)){
			stop("dim.reduct is not complete, run 'runDimReduct' first")
		}
	}
		
	down.sample = min(down.sample, ncell);
	idx.ds = sort(sample(seq(ncell), down.sample));
	obj = obj[idx.ds,,drop=FALSE];

	if(max(eigs.dims) > length(obj@smat@sdev)){
		stop(paste("eigs.dims exceeds PCA dimentions ", length(obj@smat@sdev)));
	}
	
	if((x=length(eigs.dims)) > 50L){
		stop("eigs.dims must be within 1:50")
	}
	
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	op <- par(mfrow = c(5,5), oma = c(3,3,1,1) + 0.2, mar = c(0,0,1,1) + 0.2);
	PCA.plot <- split(sort(eigs.dims), ceiling(seq(eigs.dims)/2));
	if((length(x = eigs.dims)  %% 2) == 1){
		PCA.plot = PCA.plot[1:(length(PCA.plot) - 1)]
	}
	
	for(x in PCA.plot){
		data.use = data.frame(obj@smat@dmat[,c(x[1],x[2])]);	
		colnames(data.use) = c("dim1", "dim2");
		plot(x=data.use[,1], 
			 y=data.use[,2],
			 cex=point.size, 
			 col=scales::alpha(point.color, point.alpha),
			 mtext(paste(paste("eigs", x[1]), x[2], sep=" vs "), side=3),
			 yaxt='n', 
			 xaxt="n",
			 xlab="", 
			 ylab=""
		);
	}

	if(!is.null(pdf.file.name)){
		dev.off()		
	}	
	graphics::par(mfrow=c(1,1));
}

#' Plot gene-body accessibility level
#'
#' @param obj A snap object.
#' @param gene.names Name of genes to plot.
#' @param viz.method Visulization method c("tsne", "umap").
#' @param point.size Point size [0.5].
#' @param point.shape Point shape type [19].
#' @param point.color Point color ["blue"].
#' @param background.point If add points as background [TRUE].
#' @param background.point.color Color of background points ["grey"].
#' @param background.point.alpha Transparency level of background points [0.3].
#' @param background.point.size Size of background points [0.5].
#' @param background.point.shape Shape of background points [19].
#' @param low.value Feature value is standarded to 0-1, value less than low.value will be set to low.value.
#' @param high.value Feature value is standarded to 0-1, value greater than high.value will be set to high.value.
#' @param down.sample Number of cells to plot. Cells will be randomly downsampled for this number for plotting.
#' @param seed.use Random seed [10].
#' @param plot.nrow Number of rows in the plot [3].
#' @param plot.ncol Number of columns in the plot [3].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @param ... Arguments passed to plot method.
#' 
#' @examples 
#' data(demo.sp);
#' demo.sp = scaleCountMatrix(
#'	obj=demo.sp,
#'	mat="gmat",
#'	cov=SnapATAC::rowSums(demo.sp, mat="bmat"),
#'	method="RPM"
#'	)
#' gene.names = c(
#'	"Prdm14",   "E330040D14Rik", "Gm17971",       
#'	"Gm17970",	"Defb44-ps",     "Gm7357",
#'	"Gm37265",  "Kctd18",		"Gm37143"
#'	)
#' plotGene(
#'	obj=demo.sp, 
#'	gene.names=gene.names, 
#'	viz.method="tsne"
#'	)
#' @importFrom grDevices pdf dev.off
#' @importFrom methods slot 
#' @importFrom scales alpha
#' @importFrom graphics par points plot 
#' @importFrom grDevices pdf dev.off
#' @export
plotGene <- function(obj, gene.names, viz.method, point.size, point.color, point.shape,  background.point, background.point.color, background.point.alpha, background.point.size, background.point.shape, low.value, high.value, down.sample, seed.use, plot.nrow, plot.ncol, pdf.file.name, pdf.height, pdf.width,...){
  UseMethod("plotGene", obj);
}

#' @export
plotGene.default <- function(
	obj, 
	gene.names,
	viz.method=c("tsne", "umap"),
	point.size=0.5,
	point.color="red",
	point.shape=19,
	background.point=TRUE,
	background.point.color="grey",
	background.point.alpha=0.3,
	background.point.size=0.5,
	background.point.shape=19,
	low.value=0.0,
	high.value=1.0,
	down.sample=5000,
	seed.use=10,
	plot.nrow=3,
	plot.ncol=3,
	pdf.file.name=NULL, 
	pdf.height=7, 
	pdf.width=7,
	...
){
	
	if(missing(obj)){
		stop("obj is missing")
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object")
		}
		ncell = nrow(obj);
		data.use = obj@gmat;
		if((x=nrow(data.use)) == 0L){
			stop("gmat is empty, run addGamtToSnap frist")
		}
		if((x = nrow(data.use))!= ncell){
			stop("gmat has different number of rows from cell number, add gmat again")
		}
	}
	
	if(missing(gene.names)){
		stop("gene.names is missing")
	}else{
		if((x=length(gene.names)) > 9L){
			warning("plotGene works best with less than 9 genes in on plot");
		}
		if(any(!(gene.names %in% colnames(data.use)))){
			stop(paste(gene.names[which(!(gene.names %in% colnames(data.use)))], "does not exist in cell-by-gene matrix"))
		};	
	}

	viz.method = match.arg(viz.method);
	viz.use = methods::slot(obj, viz.method);
	
	if((x = nrow(viz.use)) == 0L){
		stop("visulization matrix is empty, run runViz first")
	}
	
	
	if(down.sample < ncell){
		set.seed(seed.use);
		idx = sort(sample(seq(ncell), down.sample));
		data.use = data.use[idx,];
		viz.use = viz.use[idx,];
	}
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	ndim = length(gene.names);
	
	op <- par(mfrow = c(plot.nrow,plot.ncol), oma = c(3,3,1,1) + 0.2, mar = c(0,0,1,1) + 0.2);
	
	for(i in seq(ndim)){
		y = data.use[,which(colnames(data.use) == gene.names[i])];
		y = pmin(1, y/quantile(y[which(y > 0)], 0.99));
		y[y < low.value] = 0;
		y[y > high.value] = 1;
		if(background.point){
			plot(viz.use, 
	   			 main=gene.names[i],
				 col=scales::alpha(background.point.color, background.point.alpha), 
				 cex=background.point.size,
				 pch=background.point.shape,
		   		 yaxt='n', 
		   		 xaxt="n",
		   		 xlab="", 
		   		 ylab="",
				 ...
				 );	
		 	points(viz.use, 
		 		   col=alpha(point.color, y), 
		 		   cex=point.size,
		 		   pch=point.shape
		 		 );	
		}else{
			plot(viz.use, 
	   			 main=gene.names[i],
				 col=alpha(point.color, y), 
				 cex=point.size,
				 pch=point.shape,
		   		 yaxt='n', 
		   		 xaxt="n",
		   		 xlab="", 
		   		 ylab="",
				 ...
				 );	
		}
	}
	
	if(!is.null(pdf.file.name)){
		dev.off()		
	}
	par(mfrow=c(1,1));
}

#' Feature Enrichment Boxplot
#'
#' @param obj A snap object.
#' @param feature Feature enrichment value for each cell.
#' @param outline If 'outline' is not true, the outliers are not drawn (as points whereas S+ uses lines).
#' @param ylab Name of ylab.
#' @param main Main title.
#' @param add.point If 'add.point' is true, scatter points will be added to the top of the boxplot.
#' @param point.size Point size [0.5].
#' @param point.shape Point shape type [19].
#' @param point.alpha Point transparancy level [0.9].
#' @param pdf.file.name pdf file name to save the plot [NULL].
#' @param pdf.width the width of the graphics region in inches [7].
#' @param pdf.height the height of the graphics region in inches [7].
#' @importFrom grDevices pdf dev.off
#' @importFrom graphics stripchart boxplot
#' @importFrom methods slot 
#' @importFrom scales alpha
#' @importFrom methods is
#' 
#' @examples
#' data(demo.sp);
#' boxPlotFeature(demo.sp, rnorm(nrow(demo.sp)));
#' 
#' @export
boxPlotFeature <- function(obj, feature, outline, ylab, main, add.point,  point.size, point.shape, point.alpha, pdf.file.name, pdf.height, pdf.width){
  UseMethod("boxPlotFeature", obj);
}

#' @export 
boxPlotFeature.default <- function(
	obj, 
	feature, 
	outline=FALSE, 
	ylab=NULL, 
	main=NULL, 
	add.point=TRUE, 
	point.size=0.2, 
	point.shape=19, 
	point.alpha=0.5, 
	pdf.file.name=NULL, 
	pdf.height=7, 
	pdf.width=7
){
	if(missing(obj)){
		stop("obj is missing")
	}else{
		if(!is(obj, "snap")){
			stop("obj is not a snap object")
		}
		ncell = nrow(obj);
		data.use = obj@cluster;
		if((x=length(data.use)) == 0L){
			stop("cluster is empty, run runCluster frist")
		}
		if((x = length(data.use))!= ncell){
			stop("gmat has different number of rows from cell number, add gmat again")
		}
	}
	
	if(length(data.use) != length(feature)){
		stop("feature has different length with cell number");
	}
	
	
	if(!is.null(pdf.file.name)){
		if(file.exists(pdf.file.name)){
			warning("pdf.file already exists");
			file.remove(pdf.file.name);
		}else{
			if(!file.create(pdf.file.name)){
				stop("cannot create pdf.file, not a directory")				
			}
			file.remove(pdf.file.name);
		}	
		grDevices::pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
	}
	
	colPanel = createColorPanel(length(unique(obj@cluster)));	
	
	data.use = data.frame(y=feature, cluster=data.use);
	graphics::boxplot(y ~ cluster, data=data.use, las=3, outline=outline, xlab=NULL, mian=main, ylab=ylab, col = colPanel[data.use$cluster]);
	if(add.point){
		graphics::stripchart(y ~ cluster,  data=data.use, vertical = TRUE, method = "jitter", add = TRUE, pch = point.shape, cex=point.size, col = alpha(colPanel[data.use$cluster], point.alpha))		
	}
	
	if(!is.null(pdf.file.name)){
		grDevices::dev.off()		
	}
	graphics::par(mfrow=c(1,1));
}


######' Feature Enrichment scatterplot
######'
######' @param obj A snap object.
######' @param method Visulziation method c("tsne", "umap").
######' @param feature Feature enrichment value for each cell.
######' @param point.size Point size [0.5].
######' @param point.shape Point shape type [19].
######' @param point.color Point color ["blue"].
######' @param pdf.file.name pdf file name to save the plot [NULL].
######' @param pdf.width the width of the graphics region in inches [7].
######' @param pdf.height the height of the graphics region in inches [7].
######' @param ... Arguments passed to plot function.
######' @importFrom grDevices pdf dev.off
######' @importFrom graphics stripchart boxplot
######' @importFrom methods slot 
######' @importFrom cowplot plot_grid
######' @importFrom scales alpha
######' @export
#####scatterPlotFeature <- function(obj, method, feature, point.size, point.shape, pdf.file.name, pdf.height, pdf.width, ...){
#####  UseMethod("scatterPlotFeature", obj);
#####}
#####
#####scatterPlotFeature.default <- function(
#####	obj,
#####	method=c("tsne", "umap"),
#####	point.size=1,
#####	point.shape=19,
#####	pdf.file.name=NULL,
#####	pdf.height=7,
#####	pdf.width=7,
#####	...
#####){
#####	
#####}
#####


############################################################
###' Project Gene-body Accessibility to T-SNE plot
###'
###' @param obj A snap object.
###' @param gene.names Name of genes to plot.
###' @param viz.method Project method c("tsne", "umap").
###' @param point.size Point size [0.5].
###' @param point.shape Point shape type [19].
###' @param point.color Point color ["blue"].
###' @param point.alpha Point transparancy level [0.9].
###' @param down.sample Number of cells to plot. Cells will be randomly downsampled for this number for plotting.
###' @param pdf.file.name pdf file name to save the plot [NULL].
###' @param pdf.width the width of the graphics region in inches [7].
###' @param pdf.height the height of the graphics region in inches [7].
###' @param ... Arguments passed to ggplot.
###' @importFrom grDevices pdf dev.off
###' @import ggplot2 
###' @importFrom methods slot 
###' @importFrom cowplot plot_grid
###' @export
##plotGene <- function(obj, gene.names, viz.method, point.size, point.color, point.shape, point.alpha, down.sample, pdf.file.name, pdf.height, pdf.width,...){
##  UseMethod("plotGene", obj);
##}
##
###' @export
##plotGene.default <- function(
##	obj, 
##	gene.names,
##	viz.method=c("tsne", "umap"),
##	point.size=0.5,
##	point.color="blue",
##	point.shape=19,
##	point.alpha=0.9,
##	low.value=0.0,
##	high.value=1.0,
##	down.sample=10000,
##	pdf.file.name=NULL, 
##	pdf.height=7, 
##	pdf.width=7,
##	...
##){
##	ncell = nrow(obj);
##	if(missing(obj)){
##		stop("obj is missing");
##	}else{
##		if(!is.snap(obj)){
##			stop("obj is not a snap object")
##		}
##		if((x=nrow(obj@gmat)) == 0L){
##			stop("gmat is empty, run addGamt frist")
##		}
##		data.use = obj@gmat;
##		if((x = nrow(data.use))!= ncell){
##			stop("gmat has different number of rows from cell number")
##		}
##	}
##
##	viz.use = methods::slot(obj, viz.method);
##	if((x = nrow(viz.use))==0L){
##		stop("visulization method does not exist, run runViz first")
##	}
##	
##	if(missing(gene.names)){
##		stop("gene.names is missing")
##	}else{
##		if((x=length(gene.names)) > 9L){
##			warning("plotGene can only plot up to 9 genes in on plot");
##			gene.names = gene.names[1:9];
##		}
##		if(any(!(gene.names %in% colnames(data.use)))){
##			stop(paste(gene.names[which(!(gene.names %in% colnames(data.use)))], "does not exist in cell-by-gene matrix"))
##		};	
##	}
##	
##	if(!is.null(pdf.file.name)){
##		if(file.exists(pdf.file.name)){
##			warning("pdf.file already exists");
##			file.remove(pdf.file.name);
##		}else{
##			if(!file.create(pdf.file.name)){
##				stop("cannot create pdf.file, not a directory")				
##			}
##			file.remove(pdf.file.name);
##		}	
##		pdf(pdf.file.name,width=pdf.width,height=pdf.height); 
##	}
##
##	plt.list = lapply(as.list(seq(length(gene.names))), function(i){
##		idy = which(colnames(data.use)==gene.names[i]);
##		y = data.use[,idy];
##		y[y==0] = NA;
##		y = pmin(1, y / quantile(y, 0.98, na.rm=TRUE));
##		y[y > high.value] = high.value;
##		y[y < low.value] = low.value;
##		data.use = data.frame(viz.use[,1], viz.use[,2], y);
##		colnames(data.use) = c("tsne1", "tsne2", "value");
##		data.use.pos = data.use[which(y > 0),];
##		sp2 <- ggplot() + 
##			geom_point(data=data.use, aes(x=data.use$tsne1, y=data.use$tsne2), color="#D3D3D3", size=point.size) +
##			geom_point(data=data.use.pos, aes(x=data.use.pos$tsne1, y=data.use.pos$tsne2, color= data.use.pos$value), size=point.size, shape=point.shape, alpha=point.alpha) +
##			scale_colour_gradient2() + 
##  		  	theme(
##				panel.grid.major = element_blank(), 
##				panel.grid.minor = element_blank(),
##				panel.background = element_blank(), 
##				legend.position='none', 
##				axis.line = element_line(colour = "black", linetype = "solid")
##				) +			
##			ggtitle(label=gene.names[i]) +
##		    xlab("Dim-1") + ylab("Dim-2")
##		return(sp2)
##	})
##	
##    if (length(x = plt.list) > 9) {
##      nCol <- 4
##    } else {
##      nCol <- min(length(x = plt.list), 3)
##    }
##	
##	plots.combined <- plot_grid(plotlist = plt.list, ncol = nCol);
##	if(!is.null(pdf.file.name)){
##		dev.off()		
##	}else{
##		print(plots.combined)
##	}	
##}