Mercurial > repos > xuebing > sharplab_interval_analysis
comparison genomeview-old2.r @ 20:16ba480adf96
Uploaded
| author | xuebing |
|---|---|
| date | Sat, 31 Mar 2012 08:31:22 -0400 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 19:d325683ec368 | 20:16ba480adf96 |
|---|---|
| 1 | |
| 2 caloffset = function(genome){ | |
| 3 total_len = sum(as.numeric(genome[,2])) | |
| 4 offset = 0 | |
| 5 for (i in 1:nrow(genome)) { | |
| 6 offset = c(offset,offset[i]+genome[i,2]) | |
| 7 } | |
| 8 offset | |
| 9 } | |
| 10 | |
| 11 coverage = function(intervals,genome,offset,resolution) { | |
| 12 | |
| 13 nChr = length(offset) - 1 | |
| 14 total_len = offset[nChr+1] | |
| 15 nbin = as.integer(total_len / resolution) | |
| 16 #cat('nbin=',nbin,'genomelen=',total_len,'\n') | |
| 17 cov = numeric(nbin)#coverage | |
| 18 col = numeric(nbin)#color | |
| 19 for (i in 1:nChr) { | |
| 20 d = x[x[,1]==as.character(genome[i,1]),2:3] | |
| 21 if (nrow(d) > 0){ | |
| 22 #cat('dim(d)=',dim(d),'\n') | |
| 23 d = ceiling((d+offset[i])*nbin/total_len) | |
| 24 for (j in 1:nrow(d)){ | |
| 25 cov[d[j,1]:d[j,2]] = cov[d[j,1]:d[j,2]] + 1 | |
| 26 } | |
| 27 } | |
| 28 col[ceiling(offset[i]*nbin/total_len):ceiling(offset[i]*nbin/total_len)] = i | |
| 29 } | |
| 30 list(nbin=nbin,cov=cov) | |
| 31 } | |
| 32 | |
| 33 # plot coverage | |
| 34 # res = genomeView(x,genome,100000) | |
| 35 plotcov = function(res,genome,offset,title,uselog) { | |
| 36 if (uselog == 'log'){ | |
| 37 res$cov = log10(res$cov+1) | |
| 38 } | |
| 39 ymax = max(res$cov) | |
| 40 par(mar=c(5,5,5,1)) | |
| 41 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)) | |
| 42 xticks = numeric(nrow(genome)) | |
| 43 for (i in 1:nrow(genome)){ | |
| 44 xticks[i] = (offset[i]+offset[i+1])/2*res$nbin/offset[length(offset)] | |
| 45 } | |
| 46 mtext(genome[,1],side=1,at=xticks,adj=1,las=2,col=seq(nrow(genome))) | |
| 47 } | |
| 48 | |
| 49 union_correlation = function(x,y){ | |
| 50 z = x>0 | y>0 | |
| 51 cor(x[z],y[z]) | |
| 52 } |