diff genomeview.r @ 11:b7f1d9f8f3bc

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/genomeview.r	Sat Mar 10 07:59:27 2012 -0500
@@ -0,0 +1,70 @@
+
+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)
+    cov = numeric(nbin)#coverage
+    col = numeric(nbin)#color
+    for (i in 1:nChr) {
+        d = x[x[,1]==as.character(genome[i,1]),2:3]
+        d = ceiling(((d[,1]+d[,2])/2+offset[i])*nbin/total_len)
+        t = table(d)
+		pos = as.numeric(row.names(t)) 
+        cov[pos] = cov[pos] + as.numeric(t)
+        col[pos] = i
+    }
+    list(nbin=nbin,cov=cov,col=col)
+}
+
+# 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)
+    #print(ymax)
+	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])
+}
+
+
+heatmap2 = function(x,scale,sym,margins,labRow,labCol){
+    h = heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow,labCol=labCol)
+    x = x[h$rowInd,h$colInd]
+    tx = numeric(0)
+    ty = numeric(0)
+    txt = character(0)
+    for (i in 1:nrow(x)){
+        for (j in 1:ncol(x)){
+            tx <- c(tx,i)
+            ty <- c(ty,ncol(x)-j+1)
+            txt <- c(txt,format(x[i,j],digits=2,nsmall=2))
+        }    
+    }
+    #heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow[h$rowInd],labCol=labCol[h$colInd],add.expr=text(1:4,1:4,1:4))
+	cat(dim(tx))
+	text(tx,ty,txt)
+	heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow[h$rowInd],labCol=labCol[h$colInd],add.expr=text(tx,ty,txt))
+
+}