Mercurial > repos > xuebing > sharplabtool
comparison genomeview.r @ 11:b7f1d9f8f3bc
Uploaded
author | xuebing |
---|---|
date | Sat, 10 Mar 2012 07:59:27 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
10:1558594a3c2e | 11:b7f1d9f8f3bc |
---|---|
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 cov = numeric(nbin)#coverage | |
17 col = numeric(nbin)#color | |
18 for (i in 1:nChr) { | |
19 d = x[x[,1]==as.character(genome[i,1]),2:3] | |
20 d = ceiling(((d[,1]+d[,2])/2+offset[i])*nbin/total_len) | |
21 t = table(d) | |
22 pos = as.numeric(row.names(t)) | |
23 cov[pos] = cov[pos] + as.numeric(t) | |
24 col[pos] = i | |
25 } | |
26 list(nbin=nbin,cov=cov,col=col) | |
27 } | |
28 | |
29 # plot coverage | |
30 # res = genomeView(x,genome,100000) | |
31 plotcov = function(res,genome,offset,title,uselog) { | |
32 if (uselog == 'log'){ | |
33 res$cov = log10(res$cov+1) | |
34 } | |
35 ymax = max(res$cov) | |
36 #print(ymax) | |
37 par(mar=c(5,5,5,1)) | |
38 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)) | |
39 xticks = numeric(nrow(genome)) | |
40 for (i in 1:nrow(genome)){ | |
41 xticks[i] = (offset[i]+offset[i+1])/2*res$nbin/offset[length(offset)] | |
42 } | |
43 mtext(genome[,1],side=1,at=xticks,adj=1,las=2,col=seq(nrow(genome))) | |
44 } | |
45 | |
46 union_correlation = function(x,y){ | |
47 z = x>0 | y>0 | |
48 cor(x[z],y[z]) | |
49 } | |
50 | |
51 | |
52 heatmap2 = function(x,scale,sym,margins,labRow,labCol){ | |
53 h = heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow,labCol=labCol) | |
54 x = x[h$rowInd,h$colInd] | |
55 tx = numeric(0) | |
56 ty = numeric(0) | |
57 txt = character(0) | |
58 for (i in 1:nrow(x)){ | |
59 for (j in 1:ncol(x)){ | |
60 tx <- c(tx,i) | |
61 ty <- c(ty,ncol(x)-j+1) | |
62 txt <- c(txt,format(x[i,j],digits=2,nsmall=2)) | |
63 } | |
64 } | |
65 #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)) | |
66 cat(dim(tx)) | |
67 text(tx,ty,txt) | |
68 heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow[h$rowInd],labCol=labCol[h$colInd],add.expr=text(tx,ty,txt)) | |
69 | |
70 } |