-
Notifications
You must be signed in to change notification settings - Fork 1
/
chrView-mean.R
124 lines (99 loc) · 3.8 KB
/
chrView-mean.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#chrView-mean plot
chrvwlist=read.table("chrView_mean_list.txt",sep='\t', stringsAsFactors=FALSE)
n_rain=nrow(chrvwlist)
cl <- rainbow(n_rain)
viewfile1=read.table(chrvwlist[1,2],sep='\t',header=T)
ym_cg=max(viewfile1$meanCG)*1.5*100
ym_chg=max(viewfile1$meanCHG)*1.5*100
ym_chh=max(viewfile1$meanCHH)*1.5*100
ys_cg=min(viewfile1$meanCG)*1.5*100
ys_chg=min(viewfile1$meanCHG)*1.5*100
ys_chh=min(viewfile1$meanCHH)*1.5*100
#CG
outfile=paste("./chrView_mean_CG.png",sep="")
bitmap(outfile,type="png16m",height=2,width=8,res=400)
plot(NA,ylab=expression(paste(Delta,' CG methylation levels (%)')),xlab="Chromosome",col='white',pch=19,type='p',cex.main = 1,frame = FALSE,ylim=c(ys_cg,ym_cg),xlim=c(0,nrow(viewfile1)),xaxt="n",cex=1.5)
#axis(2,at=seq(ys_cg,ym_cg))
for (i in 1:nrow(chrvwlist)){
name=chrvwlist[i,1]
file=chrvwlist[i,2]
viewfile=read.table(file,sep='\t',header=T)
chr=unique(viewfile$chromosome)
for (j in chr){
sub=viewfile[viewfile$chromosome==j,]
lines(sub$Position,sub$meanCG*100,lwd=0.5,col= cl[i])
lines(sub$Position,vector("numeric", length(sub$Position)),type='h',lwd=0.5,col='black') #draw zero line
#last row
abline(v=sub[nrow(sub),1]+0.5,col="grey",lwd=0.5)
#label
labsite=sub[as.integer(nrow(sub)/2),1]
lab=sub$chromosome[1]
sta_p=sub[1,1]
end_p=sub[(nrow(sub)),1]
axis(1,at=labsite ,label=lab,tick=F)
axis(1,at=sta_p ,label=FALSE)
axis(1,at=end_p ,label=FALSE)
xx<-c(sta_p,end_p)
axis(1,at=xx,labels=F)
}
}
legend("topright", legend = t(chrvwlist[1]), col = t(cl), lwd = 1,cex = 0.5)
dev.off()
#CHG
outfile=paste("./chrView_mean_CHG.png",sep="")
bitmap(outfile,type="png16m",height=2,width=8,res=400)
plot(NA,ylab=expression(paste(Delta,' CHG methylation levels (%)')),xlab="Chromosome",col='white',pch=19,type='p',cex.main = 1,frame = FALSE,ylim=c(ys_chg,ym_chg),xlim=c(0,nrow(viewfile1)),xaxt="n",cex=1.5)
for (i in 1:nrow(chrvwlist)){
name=chrvwlist[i,1]
file=chrvwlist[i,2]
viewfile=read.table(file,sep='\t',header=T)
chr=unique(viewfile$chromosome)
for (j in chr){
sub=viewfile[viewfile$chromosome==j,]
lines(sub$Position,sub$meanCHG*100,lwd=0.5,col= cl[i])
lines(sub$Position,vector("numeric", length(sub$Position)),type='h',lwd=0.5,col='black') #draw zero line
#last row
abline(v=sub[nrow(sub),1]+0.5,col="grey",lwd=0.5)
#label
labsite=sub[as.integer(nrow(sub)/2),1]
lab=sub$chromosome[1]
sta_p=sub[1,1]
end_p=sub[(nrow(sub)),1]
axis(1,at=labsite ,label=lab,tick=F)
axis(1,at=sta_p ,label=FALSE)
axis(1,at=end_p ,label=FALSE)
xx<-c(sta_p,end_p)
axis(1,at=xx,labels=F)
}
}
legend("topright", legend = t(chrvwlist[1]), col = t(cl), lwd = 1,cex = 0.5)
dev.off()
#CHH
outfile=paste("./chrView_mean_CHH.png",sep="")
bitmap(outfile,type="png16m",height=2,width=8,res=400)
plot(NA,ylab=expression(paste(Delta,' CHH methylation levels (%)')),xlab="Chromosome",col='white',pch=19,type='p',cex.main = 1,frame = FALSE,ylim=c(ys_chh,ym_chh),xlim=c(0,nrow(viewfile1)),xaxt="n",cex=1.5)
for (i in 1:nrow(chrvwlist)){
name=chrvwlist[i,1]
file=chrvwlist[i,2]
viewfile=read.table(file,sep='\t',header=T)
chr=unique(viewfile$chromosome)
for (j in chr){
sub=viewfile[viewfile$chromosome==j,]
lines(sub$Position,sub$meanCHH*100,lwd=0.5,col= cl[i])
lines(sub$Position,vector("numeric", length(sub$Position)),type='h',lwd=0.5,col='black') #draw zero line
#last row
abline(v=sub[nrow(sub),1]+0.5,col="grey",lwd=0.5)
#label
labsite=sub[as.integer(nrow(sub)/2),1]
lab=sub$chromosome[1]
sta_p=sub[1,1]
end_p=sub[(nrow(sub)),1]
axis(1,at=labsite ,label=lab,tick=F)
axis(1,at=sta_p ,label=FALSE)
axis(1,at=end_p ,label=FALSE)
xx<-c(sta_p,end_p)
axis(1,at=xx,labels=F)
}
}
legend("topright", legend = t(chrvwlist[1]), col = t(cl), lwd = 1,cex = 0.5)
dev.off()