Mercurial > repos > xuebing > sharplabtool
view tools/mytools/genomeview-old2.r @ 1:cdcb0ce84a1b
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:45:15 -0500 |
parents | 9071e359b9a3 |
children |
line wrap: on
line source
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]) }