|
20
|
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 }
|
