annotate genomeview.r @ 12:2f4ea569f048

Uploaded
author xuebing
date Sat, 10 Mar 2012 08:10:44 -0500
parents b7f1d9f8f3bc
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
2 caloffset = function(genome){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
3 total_len = sum(as.numeric(genome[,2]))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
4 offset = 0
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
5 for (i in 1:nrow(genome)) {
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
6 offset = c(offset,offset[i]+genome[i,2])
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
7 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
8 offset
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
9 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
10
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
11 coverage = function(intervals,genome,offset,resolution) {
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
12
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
13 nChr = length(offset) - 1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
14 total_len = offset[nChr+1]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
15 nbin = as.integer(total_len / resolution)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
16 cov = numeric(nbin)#coverage
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
17 col = numeric(nbin)#color
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
18 for (i in 1:nChr) {
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
19 d = x[x[,1]==as.character(genome[i,1]),2:3]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
20 d = ceiling(((d[,1]+d[,2])/2+offset[i])*nbin/total_len)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
21 t = table(d)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
22 pos = as.numeric(row.names(t))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
23 cov[pos] = cov[pos] + as.numeric(t)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
24 col[pos] = i
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
25 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
26 list(nbin=nbin,cov=cov,col=col)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
27 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
28
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
29 # plot coverage
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
30 # res = genomeView(x,genome,100000)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
31 plotcov = function(res,genome,offset,title,uselog) {
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
32 if (uselog == 'log'){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
33 res$cov = log10(res$cov+1)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
34 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
35 ymax = max(res$cov)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
36 #print(ymax)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
37 par(mar=c(5,5,5,1))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
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))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
39 xticks = numeric(nrow(genome))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
40 for (i in 1:nrow(genome)){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
41 xticks[i] = (offset[i]+offset[i+1])/2*res$nbin/offset[length(offset)]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
42 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
43 mtext(genome[,1],side=1,at=xticks,adj=1,las=2,col=seq(nrow(genome)))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
44 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
45
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
46 union_correlation = function(x,y){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
47 z = x>0 | y>0
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
48 cor(x[z],y[z])
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
49 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
50
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
51
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
52 heatmap2 = function(x,scale,sym,margins,labRow,labCol){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
53 h = heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow,labCol=labCol)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
54 x = x[h$rowInd,h$colInd]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
55 tx = numeric(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
56 ty = numeric(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
57 txt = character(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
58 for (i in 1:nrow(x)){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
59 for (j in 1:ncol(x)){
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
60 tx <- c(tx,i)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
61 ty <- c(ty,ncol(x)-j+1)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
62 txt <- c(txt,format(x[i,j],digits=2,nsmall=2))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
63 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
64 }
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
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))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
66 cat(dim(tx))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
67 text(tx,ty,txt)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
68 heatmap(x,scale=scale,sym=sym,margins=margins,labRow=labRow[h$rowInd],labCol=labCol[h$colInd],add.expr=text(tx,ty,txt))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
69
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
70 }