comparison genomeview-old2.r @ 20:16ba480adf96

Uploaded

20:16ba480adf96

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])
}