Some Useful Rscripts

barplot_with_pareto.R

#=============================================================================
#     FileName: plot.R
#         Desc: A short script for hist barplot and pareto cumulated frequency.
#       Author: tanhao
#        Email: tanhao2013@gmail.com
#     HomePage: http://buttonwood.github.io
#      Version: 0.0.1
#   LastChange: 2014-03-20 14:35:48
#      History:
#=============================================================================
#!/usr/local/bin/Rscript --slave

barplot_with_pareto <- function(x, main="", xlab="Indentity(%)", r=100, s=5) {
	op <- par(mar = c(5, 4, 4, 5) + 0.1, las = 2)
	
	bw <- seq(0, r, s)
	r1 <- hist(x,breaks=bw,plot=FALSE)
	print(r1)
	x  <- r1$counts
	y  <- cumsum(x)/sum(x)

	barplot(x, xlab=xlab, ylab="Counts",main=main, col="pink")
	axis(2)
	xn <- paste(seq(0,r-s,s), seq(s,r,s), sep="--")
	names(x) <- xn 

	par(new = T)
	plot(y, type="l", col="blue", axes=FALSE,xlab='', ylab='', main='')
	axis(4)
	par(las=0)
	mtext("Cumulated frequency", side=4, line=3)
	
	axis(1, at=1:length(x), labels=names(x))
	par(op)
}

#args = commandArgs(TRUE)
#print(args)
#print(args[0])
#print(args[1])
#data <- read.table(args[0])
#pdf(args[1])
#data <- read.table("rate.tab.best")
#pdf("best1t1.pdf")
data <- read.table("best.tab")
pdf("best1t1.pdf")
#bw <- seq(0,100,5)
#r1 <- hist(data[,1],breaks=bw,plot=FALSE)
#plot(bw,r1$counts,xlab="Indentity(%)",ylab="Counts(%)",type="l",lty=1,col="blue")
barplot_with_pareto(data[,1],main="Average Indentity Counts", xlab="Average Indentity(%)")
barplot_with_pareto(data[,2],main="Effective Query Indentity Counts", xlab="Effective Query Indentity(%)")
barplot_with_pareto(data[,3],main="Query Coverage Counts", xlab="Query Coverage(%)")
barplot_with_pareto(data[,4],main="Query Aln rate Counts", xlab="Query Aln Rate(%)")
barplot_with_pareto(data[,5],main="Effective Target Indentity Counts", xlab="Effective Target Indentity(%)")
barplot_with_pareto(data[,6],main="Target Coverage Counts",xlab="Target Coverage(%)")
barplot_with_pareto(data[,7],main="Target Aln Rate Counts",xlab="Target Aln Rate(%)")

cnames <- c(
    "Average Indentity(%)",
    "Effective Query Indentity(%)",
    "Query Coverage(%)",
    "Query Aln Rate(%)",
    "Effective Target Indentity(%)",
    "Target Coverage(%)",
    "Target Aln Rate(%)"
)

colnames(data) <- cnames

pairs(data)

#col_table <- rainbow(20)
# rev(rich.colors(20))
library(gplots)
col_table <- rich.colors(20)
iden 	  <- data[,1]
col_iden  <- c(1:length(iden))
l_names   <- paste(seq(0,95,5),seq(5,100,5),sep = "-")

for (i in 1:length(iden)) {
	for (j in 1:length(col_table)) {
		if ( (iden[i] >= 5*(j-1)) && (iden[i] <= 5*j) ){
			col_iden[i] = col_table[j]
		}
	}
}

for (i in 1:length(cnames)) {
	    for (j in i:length(cnames)){
	    		if(i != j){
	    	    	plot(data[,i],data[,j],xlim=c(0,100), ylim=c(0,100),
	    	    		xlab=cnames[i], ylab=cnames[j], col=col_iden)
	    	    	legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))
	    	    }
	    	}
	} 

plot(data[,1], data[,2], xlim=c(0,100), ylim=c(0,100), 
	main="Average Indentity vs Effective Query Indentity", 
	xlab=cnames[1],
    ylab=cnames[2],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))

plot(data[,1], data[,3], xlim=c(0,100), ylim=c(0,100),
    main="Average Indentity vs Query Coverage",
    xlab=cnames[1],
    ylab=cnames[3],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))

plot(data[,1], data[,4], xlim=c(0,100), ylim=c(0,100), 
	main="Average Indentity vs Query Aln Rate", 
	xlab=cnames[1],
    ylab=cnames[4],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,1], data[,5], xlim=c(0,100), ylim=c(0,100), 
	main="Average Indentity vs Effective Target Indentity", 
	xlab=cnames[1],
    ylab=cnames[5],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,1], data[,6], xlim=c(0,100), ylim=c(0,100),
    main="Average Indentity vs Target Coverage",
    xlab=cnames[1],
    ylab=cnames[6],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,1], data[,7], xlim=c(0,100), ylim=c(0,100), 
	main="Average Indentity vs Target Aln Rate", 
	xlab=cnames[1],
    ylab=cnames[7],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,2], data[,5], xlim=c(0,100), ylim=c(0,100), 
	main="Effective Query Indentity vs Effective Target Indentity", 
	xlab=cnames[2],
    ylab=cnames[5],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,3], data[,6], xlim=c(0,100), ylim=c(0,100), 
	main="Query Coverage vs Target Coverage", 
	xlab=cnames[3],
    ylab=cnames[6],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))


plot(data[,4], data[,7], xlim=c(0,100), ylim=c(0,100), 
	main="Query Aln Rate vs Target Aln Rate", 
	xlab=cnames[4],
    ylab=cnames[7],
    col=col_iden
)
legend("topleft",legend=l_names, col=col_table,pch=rep(1,20))



#barplot_with_pareto(data[,1]*data[,2]/10000,
#	main="Indentity Counts Adjust by Query Coverage", xlab="Indentity Adjusted(%)", r=1, s=0.05)

#barplot_with_pareto(data[,1]*data[,3]/10000,
#	main="Indentity Counts Adjust by Query Aln rate", xlab="Indentity Adjusted(%)", r=1, s=0.05)

#barplot_with_pareto(data[,1]*data[,4]/10000,
#	main="Indentity Counts Adjust by Target Coverage", xlab="Indentity Adjusted(%)", r=1, s=0.05)

#barplot_with_pareto(data[,1]*data[,5]/10000,
#	main="Indentity Counts Adjust by Target Aln rate", xlab="Indentity Adjusted(%)", r=1, s=0.05)

dev.off()

#*************  Float like a butterfly! Stay like a buttonwood!  *************