Mercurial > repos > xuebing > sharplabtool
diff mytools/genomeview.r @ 7:f0dc65e7f6c0
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:59:07 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mytools/genomeview.r Fri Mar 09 19:59:07 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)) + +}