--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/genomeview-old2.r	Sat Mar 31 08:31:22 2012 -0400
@@ -0,0 +1,52 @@
+
+caloffset = function(genome){
+    total_len = sum(as.numeric(genome[,2]))
+    offset = 0
+    for (i in 1:nrow(genome)) {
+        offset = c(offset,offset[i]+genome[i,2])        
+    }
+    offset    
+}
+
+coverage = function(intervals,genome,offset,resolution) {
+
+    nChr = length(offset) - 1
+    total_len = offset[nChr+1]
+    nbin = as.integer(total_len / resolution)
+    #cat('nbin=',nbin,'genomelen=',total_len,'\n')
+    cov = numeric(nbin)#coverage
+    col = numeric(nbin)#color
+    for (i in 1:nChr) {
+        d = x[x[,1]==as.character(genome[i,1]),2:3]
+		if (nrow(d) > 0){
+			#cat('dim(d)=',dim(d),'\n')
+			d = ceiling((d+offset[i])*nbin/total_len)
+			for (j in 1:nrow(d)){
+				cov[d[j,1]:d[j,2]] = cov[d[j,1]:d[j,2]] + 1
+			}
+		}
+        col[ceiling(offset[i]*nbin/total_len):ceiling(offset[i]*nbin/total_len)] = i
+    }
+    list(nbin=nbin,cov=cov)
+}
+
+# plot coverage
+# res = genomeView(x,genome,100000)
+plotcov = function(res,genome,offset,title,uselog) {
+	if (uselog == 'log'){
+		res$cov = log10(res$cov+1)
+	}
+    ymax = max(res$cov)
+	par(mar=c(5,5,5,1))
+    plot(seq(length(res$cov)),res$cov,type='h',cex=0.1,cex.axis=2,cex.lab=2,cex.main=3,col=res$col,xaxt='n',main=title,xlab='chromosome',ylab='coverage',frame.plot=F,ylim=c(0,ymax))
+    xticks = numeric(nrow(genome))
+    for (i in 1:nrow(genome)){
+       xticks[i] = (offset[i]+offset[i+1])/2*res$nbin/offset[length(offset)]
+    }
+    mtext(genome[,1],side=1,at=xticks,adj=1,las=2,col=seq(nrow(genome)))
+}
+
+union_correlation = function(x,y){
+	z = x>0 | y>0
+	cor(x[z],y[z])
+}