sp_density2d.sh

#!/bin/bash

#set -x

usage()
{
cat <<EOF
${txtcyn}

***CREATED BY Chen Tong (chentong_biology@163.com)***

Usage:

$0 options${txtrst}

${bldblu}Function${txtrst}:

This script is used to draw a line or multiple lines using ggplot2.
You can specify whether or not smooth your line or lines.

fileformat for -f (suitable for data extracted from one sample, the
number of columns is unlimited. Column 'Set' is not necessary)
------------------------------------------------------------
Pos	gene1	gene2	gene3	gene4	gene5 ...
time1	1	2000	3000	4000	4000  ...
time2	1	2000	3000	4000	4000  ...
time3	1	2000	3000	4000	4000  ...
time4	1	2000	3000	4000	4000  ...
time5	1	2000	3000	4000	4000  ...
time6	1	2000	3000	4000	4000  ...
-------------------------------------------------------------

${txtbld}OPTIONS${txtrst}:
	-f	Data file (with header line, the first column is the
 		rownames (currently only numerical rownames are supported), 
		tab seperated. First row as column names.)${bldred}[NECESSARY]${txtrst}
	-A	The attribute of x-axis variable.
		[${txtred}Default <numeric>, means X-axis label is number.
		Accept <text> means X-axis label is text.${txtrst}]
	-H	Name for legend variable.
		${bldred}[Default variable, this should only be set when -m is TRUE]${txtrst}
	-J	Name for color variable.
		${bldred}[Default same as -H, this should only be set when -m is TRUE]${txtrst}
	-l	Levels for legend variable
		[${txtred}Default column order, accept a string like
		"'ctcf','h3k27ac','enhancer'"  
		***When -m is used, this default will be ignored too.********* 
	   	${txtrst}]
	-P	Legend position[${txtred}Default none. Accept
		top,bottom,left,right, or c(0.08,0.8).${txtrst}]
	-g	Set Ylim.
		${bldred}[Default null. Accept a string in format like '-2,2'.]${txtrst}
	-L	Levels for x-axis variable, suitable when x-axis is not used
		as a number. 
		[${txtred}Default the order of first column, accept a string like
		"'g','a','j','x','s','c','o','u'"  
	   	This will only be considered when -A is TRUE.
		***When -m is used, this default will be ignored too.********* 
		${txtrst}]
	-o	Smooth your data or not.
		[${txtred}Default FALSE means no smooth. Accept TRUE to smooth
		lines.${txtrst}]
	-O	The smooth method you want to use.
		[${txtred}smoothing method (function) to use,  eg. lm, glm,
		gam, loess,rlm. 
		For datasets with n < 1000 default is 'loess'. 
		For datasets with 1000 or more observations defaults to 'gam'.
		${txtrst}]
	-V	Add vertical lines.${bldred}[Default FALSE, accept a series of
		numbers in following format "c(1,2,3,4,5)" or other
		R code that can generate a vector]${txtrst}
	-D	Add labels to vline.
		${bldred}[Default same as -V
		Accept a series of numbers in following format "c(1,2,3,4,5)" or other R code
		that can generate a vector as labels.
		Or one can give '1' to disallow labels]${txtrst}
	-j	Add horizontal lines.${bldred}[Default FALSE, accept a series of
		numbers in following format "c(1,2,3,4,5)" or other
		R code that can generate a vector]${txtrst}
	-d	Add labels to hline.
		${bldred}[Default same as -j
		Accept a series of numbers in following format "c(1,2,3,4,5)" or other R code
		that can generate a vector as labels.
		Or one can give '1' to disallow labels]${txtrst}
	-I	Manually set the position of xtics.
		${bldred}[Default FALSE,  accept a series of
		numbers in following format "c(1,2,3,4,5)" or other R code
		that can generate a vector to set the position of xtics]${txtrst}
	-b	Manually set the value of xtics when -I is specified.
		${bldred}[Default the content of -I when -I is specified, 
		accept a series of
		numbers in following format "c(1,2,3,4,5)" or other R code
		that can generate a vector to set the position of xtics]${txtrst}
	-X	Display xtics. ${bldred}[Default TRUE]${txtrst}
	-Y	Display ytics. ${bldred}[Default TRUE]${txtrst}
	-R	Rotation angle for x-axis value(anti clockwise)
		${bldred}[Default 0]${txtrst}
	-B	line size.[${txtred}Default 1. Accept a number.${txtrst}]
	-t	Title of picture[${txtred}Default empty title${txtrst}]
	-x	xlab of picture[${txtred}Default empty xlab${txtrst}]
	-y	ylab of picture[${txtred}Default empty ylab${txtrst}]
	-c	Manually set colors for each line.[${txtred}Default FALSE,
		meaning using ggplot2 default.${txtrst}]
	-C	Color for each line.[${txtred}
		When -c is TRUE, one have two options
		1. Supplying a function to generate colors, 
		like "rainbow(11)" or "rainbow(11, alpha=0.6)", 
		    rainbow is an R color palletes, 
		    11 is the number of colors you want to get, 
			0.6 is the alpha value.
		The R palletes include heat.colors, terrain.colors,
		topo.colors, cm.colors.
		2. Supplying a list of colors in given format, 
		ususlly the number of colors should be equal to the number of
		bars like "'red','pink','blue','cyan','green','yellow'" or
		"rgb(255/255,0/255,0/255),rgb(255/255,0/255,255/255),
		 rgb(0/255,0/255,255/255),rgb(0/255,255/255,255/255),
		 rgb(0/255,255/255,0/255),rgb(255/255,255/255,0/255)"
		${txtrst}]
	-s	Scale y axis
		[${txtred}Default null. Accept TRUE. This function is
		depleted. If the supplied number after -S is not 0, this
		parameter is TRUE${txtrst}]
	-F	The formula for facets.[${bldred}Default no facets, 
		+facet_grid(level ~ .) means divide by levels of 'level' vertcally.
		+facet_grid(. ~ level) means divide by levels of 'level' horizontally.
		+facet_grid(lev1 ~ lev2) means divide by lev1 vertically and lev2
		horizontally.
		+facet_wrap(~level, ncol=2) means wrap horizontally with 2
		columns.
		#Pay attention to the single quote for parameters in function for scale.
		Example: +facet_wrap(~Size,ncol=6,scale='free')
		Example: +facet_grid(Size ~ .,scale='free_y')
		${txtrst}]
	-G	If facet is given, you may want to specifize the order of
		variable in your facet, default alphabetically.
		[${txtred}Accept sth like 
		(one level one sentence, separate by';') 
		data\$size <- factor(data\$size, levels=c("l1",
		"l2",...,"l10"), ordered=T) ${txtrst}]
	-v	If scale is TRUE, give the following
		scale_y_log10()[default], coord_trans(y="log10"), or other legal
		command for ggplot2 or simply log2.)${txtrst}]
	-S	A number to add if scale is used.
		[${txtred}Default 0. If a non-zero number is given, -s is
		TRUE.${txtrst}]	
	-p	Other legal R codes for gggplot2 will be given here.
		[${txtres}Begin with '+' ${txtrst}]
	-w	The width of output picture.[${txtred}Default 20${txtrst}]
	-u	The height of output picture.[${txtred}Default 12${txtrst}] 
	-E	The type of output figures.[${txtred}Default pdf, accept
		eps/ps, tex (pictex), png, jpeg, tiff, bmp, svg and wmf)${txtrst}]
	-r	The resolution of output picture.[${txtred}Default 300 ppi${txtrst}]
	-z	Is there a header[${bldred}Default TRUE${txtrst}]
	-e	Execute or not[${bldred}Default TRUE${txtrst}]
	-i	Install depended packages[${bldred}Default FALSE${txtrst}]
EOF
}

file=
title=''
melted='FALSE'
xlab=''
ylab=''
xvariable='id'
variable='variable'
color_variable='variable'
level=""
x_level=""
x_type='numeric'
scaleY='FALSE'
y_add=0
scaleY_x='scale_y_log10()'
header='TRUE'
execute='TRUE'
ist='FALSE'
uwid=20
vhig=8
res=300
ext='pdf'
par=''
legend_pos='none'
smooth='FALSE'
smooth_method='auto'
line_size=1
xtics='TRUE'
xtics_angle=0
ytics='TRUE'
color='FALSE'
color_v=''
custom_vline='NULL'
custom_vanno='NULL'
custom_hline='NULL'
custom_hanno='NULL'
facet=''
facet_o=''
xtics_pos=0
xtics_value=0
Ylim='NULL'


while getopts "hf:m:a:A:b:I:t:x:l:j:J:d:F:g:G:H:P:L:y:V:D:c:C:B:X:Y:R:w:u:r:o:O:s:S:p:z:v:e:E:i:" OPTION
do
	case $OPTION in
		h)
			usage
			exit 1
			;;
		f)
			file=$OPTARG
			;;
		m)
			melted=$OPTARG
			;;
		a)
			xvariable=$OPTARG
			;;
		A)
			x_type=$OPTARG
			;;
		g)
			Ylim=$OPTARG
			;;
		H)
			variable=$OPTARG
			;;
		V)
			custom_vline=$OPTARG
			;;
		D)
			custom_vanno=$OPTARG
			;;
		j)
			custom_hline=$OPTARG
			;;
		J)
			color_variable=$OPTARG
			;;
		d)
			custom_hanno=$OPTARG
			;;
		I)
			xtics_pos=$OPTARG
			;;
		b)
			xtics_value=$OPTARG
			;;
		t)
			title=$OPTARG
			;;
		x)
			xlab=$OPTARG
			;;
		l)
			level=$OPTARG
			;;
		F)
			facet=$OPTARG
			;;
		G)
			facet_o=$OPTARG
			;;
		P)
			legend_pos=$OPTARG
			;;
		B)
			line_size=$OPTARG
			;;
		c)
			color=$OPTARG
			;;
		C)
			color_v=$OPTARG
			;;
		X)
			xtics=$OPTARG
			;;
		R)
			xtics_angle=$OPTARG
			;;
		Y)
			ytics=$OPTARG
			;;
		L)
			x_level=$OPTARG
			;;
		p)
			par=$OPTARG
			;;
		y)
			ylab=$OPTARG
			;;
		w)
			uwid=$OPTARG
			;;
		u)
			vhig=$OPTARG
			;;
		r)
			res=$OPTARG
			;;
		E)
			ext=$OPTARG
			;;
		o)
			smooth=$OPTARG
			;;
		O)
			smooth_method=$OPTARG
			;;
		s)
			scaleY=$OPTARG
			;;
		S)
			y_add=$OPTARG
			;;
		v)
			scaleY_x=$OPTARG
			;;
		z)
			header=$OPTARG
			;;
		e)
			execute=$OPTARG
			;;
		i)
			ist=$OPTARG
			;;
		?)
			usage
			exit 1
			;;
	esac
done

if [ -z $file ]; then
	usage
	exit 1
fi

if test ${y_add} -ne 0; then
	scaleY="TRUE"
fi

if test "${color_variable}" == "variable"; then
	color_variable=${variable}
fi

mid='.density2d'

if test "${smooth}" == 'TRUE'; then
	mid=${mid}'.smooth'
fi

cat <<END >${file}${mid}.r

if ($ist){
	install.packages("ggplot2", repo="http://cran.us.r-project.org")
	install.packages("reshape2", repo="http://cran.us.r-project.org")
	install.packages("grid", repo="http://cran.us.r-project.org")
}
library(ggplot2)
library(reshape2)
library(grid)

if(! $melted){

	data <- read.table(file="${file}", sep="\t", header=$header,
	row.names=1, check.names=F, quote="")
	data_rownames <- rownames(data)
	data_colnames <- colnames(data)

	data_mean <- rowMeans(data)
	data_mean <- data.frame(${xvariable}=data_rownames, data_mean=data_mean)
	data_mean\$${xvariable} <- as.numeric(levels(data_mean\$${xvariable}))[data_mean\$${xvariable}]

	data2 <- cbind(${xvariable}=rownames(data), data)
	data_m <- melt(data2, id="${xvariable}")
	data_m\$${xvariable} <- as.numeric(levels(data_m\$${xvariable}))[data_m\$${xvariable}]

	#data\$${xvariable} <- data_rownames
	#data_m <- melt(data, id.vars=c("${xvariable}"))
} else {
	data_m <- read.table(file="$file", sep="\t",
	header=$header, check.names=F, quote="")
}

if (${y_add} != 0){
	data_m\$value <- data_m\$value + ${y_add}
	if ("${scaleY_x}" == "log2") {
		data_m\$value <- log2(data_m\$value)
	}
}

if ("${level}" != ""){
	level_i <- c(${level})
	data_m\$${variable} <- factor(data_m\$${variable}, levels=level_i)
} else if(! $melted){
	data_m\$${variable} <- factor(data_m\$${variable}, levels=data_colnames,
	ordered=T)
}

#if (${x_type}){
#	if ("${x_level}" != ""){
#		x_level <- c(${x_level})
#		data_m\$${xvariable} <- factor(data_m\$${xvariable},levels=x_level)
#	} else if(! $melted){
#		data_m\$${xvariable} <- factor(data_m\$${xvariable},levels=data_rownames,ordered=TRUE)
#	}
#} 

if ("${x_level}" != ""){
	x_level <- c(${x_level})
	data_m\$${xvariable} <- factor(data_m\$${xvariable},levels=x_level)
} else if(! $melted){
	data_m\$${xvariable} <- factor(data_m\$${xvariable},levels=data_rownames,ordered=TRUE)
}

if ("${x_type}" == "numeric"){
	data_m\$${xvariable} <- 
		as.numeric(levels(data_m\$${xvariable}))[data_m\$${xvariable}]
}

#if(! ${color}){
#	data_m\$${variable} <- 
#		as.numeric(levels(data_m\$${variable}))[data_m\$${variable}]
#}

${facet_o}

p = ggplot(data_m, aes($xvariable, value)) + 
	stat_density_2d(aes(alpha=..level.., group=1)) 
	#geom_density_2d(aes(color="grey")) 

p = p + xlab("$xlab") + ylab("$ylab") + theme_bw() + labs(title="${title}") +
	theme(legend.title=element_blank(),
   	panel.grid.major = element_blank(), panel.grid.minor = element_blank())


#p <- p + expand_limits(y = 0)
#p <- p + scale_y_continuous(expand=c(0, 0))

p <- p + theme(legend.key=element_blank()) 

#legend.background = element_rect(colour='white'))

#legend.background = element_rect(fill = "white"), legend.box=NULL, 
#legend.margin=unit(0,"cm"))

if (${smooth}){
	if ("${line_size}" != ""){
		#p <- p + stat_smooth(method="${smooth_method}", se=FALSE,
		p = p + stat_smooth(data=data_mean, mapping=aes(x=${xvariable}, 
			y=data_mean, colour="blue", group=1), method="${smooth_method}", se=F, 
			size=${line_size})
	}else{
		#p = p + geom_smooth(data=data_mean, mapping=aes(x=${xvariable}, 
		p = p + stat_smooth(data=data_mean, mapping=aes(x=${xvariable}, 
			y=data_mean, colour="blue", group=1), method="${smooth_method}", se=F) 
	}	
}else{
	if ("${line_size}" != ""){
		p <- p + geom_line(data=data_mean, mapping=aes(x=${xvariable}, y=data_mean), 
			colour="blue", size=${line_size}, alpha=0.6) 
	}else{
		p <- p + geom_line(data=data_mean, mapping=aes(x=${xvariable}, y=data_mean), 
			colour="blue", alpha=0.6) 
	}
}

if($scaleY & ("$scaleY_x" != "log2")){
	p <- p + $scaleY_x
}

if(${color}){
	p <- p + scale_color_manual(values=c(${color_v}))
} 
#else {
	#p <- p + scale_colour_brewer()
	#p <- p + scale_colour_discrete(limits=levels(data_m\$${variable}),
	#breaks=unique(data_m\$${variable}))
#}

if ("$xtics" == "FALSE"){
	p <- p + theme(axis.text.x=element_blank(), axis.ticks.x = element_blank())
}else{
	if (${xtics_angle} != 0){
		if (${xtics_angle} == 90){
			p <- p + theme(axis.text.x=
			  element_text(angle=${xtics_angle},hjust=1, vjust=0.5))
		}else if (${xtics_angle} == 45){
			p <- p + theme(axis.text.x=
			  element_text(angle=${xtics_angle},hjust=1, vjust=1))
		} else {
			p <- p + theme(axis.text.x=
			  element_text(angle=${xtics_angle},hjust=0.5, vjust=0.5))
		}
	}
}
if ("$ytics" == "FALSE"){
	p <- p + theme(axis.text.y=element_blank())
}


top='top'
botttom='bottom'
left='left'
right='right'
none='none'
legend_pos_par <- ${legend_pos}

p <- p + theme(legend.position=legend_pos_par)

xtics_pos <- ${xtics_pos}
xtics_value <- ${xtics_value}

if(length(xtics_pos) > 1){
	if(length(xtics_value) <= 1){
		xtics_value <- xtics_pos
	}
	p <- p + scale_x_continuous(breaks=xtics_pos, labels=xtics_value)
}

#custom_vline_coord <- ${vline}
#if(length(custom_vline_coord) > 1){
#	p <- p + geom_vline(xintercept=custom_vline_coord, 
#	linetype="dotted" )
#}

custom_vline_coord <- ${custom_vline}
custom_vline_anno <- ${custom_vanno}
if ("${custom_vanno}" == "NULL" ){
	custom_vline_anno <- custom_vline_coord
}

if(is.vector(custom_vline_coord)){
	p <- p + geom_vline(xintercept=custom_vline_coord,
	linetype="dotted", size=0.5)
	if(is.vector(custom_vline_anno)){
		ymax_range <- ggplot_build(p)\$panel\$ranges[[1]]\$y.range
		ymax_v <- ymax_range[2]
		if ("${facet}" != ""){
			ymax_v = 0
		}
		p <- p + annotate("text", x=custom_vline_coord, y=ymax_v,
		label=custom_vline_anno, hjust=0)
	}
}

custom_hline_coord <- ${custom_hline}
custom_hline_anno <- ${custom_hanno}
if ("${custom_hanno}" == "NULL" ){
	custom_hline_anno <- custom_hline_coord
}

if(is.vector(custom_hline_coord)){
	p <- p + geom_hline(yintercept=custom_hline_coord,
	linetype="dotted", size=0.5)
	if(is.vector(custom_hline_anno)){
		#xmax_range <- ggplot_build(p)\$panel\$ranges[[1]]\$x.range
		#xmax_v <- xmax_range[2]
		#if ("${facet}" != ""){
		#	xmax_v = 0
		#}
		p <- p + annotate("text", y=custom_hline_coord, x=0,
		label=custom_hline_anno, vjust=0, hjust=0)
	}
}



p <- p ${facet}

if ("${Ylim}" != "NULL") {
	p <- p + coord_cartesian(ylim=c(${Ylim}))
}

p <- p${par}


#png(filename="${file}${mid}.png", width=$uwid, height=$vhig,
#res=$res)

ggsave(p, filename="${file}${mid}.${ext}", dpi=$res, width=$uwid,
height=$vhig, units=c("cm"))
#postscript(file="${file}${mid}.eps", onefile=FALSE, horizontal=FALSE, 
#paper="special", width=10, height = 12, pointsize=10)
#dev.off()
END

if [ "$execute" == "TRUE" ]; then
	Rscript ${file}${mid}.r
if [ "$?" == "0" ]; then /bin/rm -f ${file}${mid}.r; fi
fi