-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathplotPMD.v2.R
82 lines (64 loc) · 3.59 KB
/
plotPMD.v2.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
###### Companion script for plotting the results of PMDtools (Skoglund, Northoff, Shunkov, Dervianko, Paabo, Krause, Jakobsson, 2014, PNAS)
###### Contact: [email protected]
######
###### Usage example:
###### samtools view mybam.bam | python pmdtools.py --deamination > PMD_temp.txt
###### R CMD BATCH plotPMD.R
###### cp PMD_plot.pdf mybam_PMD_plot.pdf
pdf("PMD_plot.frag.pdf",width=12,height=4,useDingbats=FALSE)
par(mfrow=c(1,3))
data<-read.table("PMD_temp.txt",header= TRUE)
read_position <-data$z
plot(read_position,read_position,xlab="Distance from 5' end of sequence read",ylab="Mismatch frequency",cex.axis=1,cex.lab=1.5,cex=3,type="n",ylim=c(0,0.50),xlim=c(0,max(read_position)))
lines(read_position,data$CG5,col="black",lwd=2)
lines(read_position,data$CA5,col="black",lwd=2)
lines(read_position,data$GT5,col="black",lwd=2)
lines(read_position,data$GC5,col="black",lwd=2)
lines(read_position,data$AG5,col="black",lwd=2)
lines(read_position,data$AT5,col="black",lwd=2)
lines(read_position,data$AC5,col="black",lwd=2)
lines(read_position,data$TC5,col="black",lwd=2)
lines(read_position,data$TG5,col="black",lwd=2)
lines(read_position,data$TA5,col="black",lwd=2)
lines(read_position,data$CG_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$CA_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$GT_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$GC_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$AG_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$AT_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$AC_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$TC_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$TG_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$TA_CpG_5,col="black",lwd=2,lty=2)
lines(read_position,data$GA_CpG_5,col="blue",lwd=2,lty=2)
lines(read_position,data$CT_CpG_5,col="red",lwd=2,lty=2)
lines(read_position,data$GA5,col="blue",lwd=2)
lines(read_position,data$CT5,col="red",lwd=2)
plot(read_position,read_position,xlab="Distance from 3' end of sequence read",ylab="Mismatch frequency",cex.axis=1,cex.lab=1.5,cex=3,type="n",ylim=c(0,0.50),xlim=c(max(read_position),0))
lines(read_position,data$CG3,col="black",lwd=2)
lines(read_position,data$CA3,col="black",lwd=2)
lines(read_position,data$GT3,col="black",lwd=2)
lines(read_position,data$GC3,col="black",lwd=2)
lines(read_position,data$AG3,col="black",lwd=2)
lines(read_position,data$AT3,col="black",lwd=2)
lines(read_position,data$AC3,col="black",lwd=2)
lines(read_position,data$TC3,col="black",lwd=2)
lines(read_position,data$TG3,col="black",lwd=2)
lines(read_position,data$TA3,col="black",lwd=2)
lines(read_position,data$CG_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$CA_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$GT_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$GC_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$AG_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$AT_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$AC_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$TC_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$TG_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$TA_CpG_3,col="black",lwd=2,lty=2)
lines(read_position,data$GA_CpG_3,col="blue",lwd=2,lty=2)
lines(read_position,data$CT_CpG_3,col="red",lwd=2,lty=2)
lines(read_position,data$GA3,col="blue",lwd=2)
lines(read_position,data$CT3,col="red",lwd=2)
plot.new()
legend("center",c("C->T","G->A","CpG->TpG","CpG->CpA","All others"),lty=c(1,1,2,2,1),lwd=c(2,2,2,2,2),col=c("red","blue","red","blue","black"),bty="b",cex=2,bg="white")
dev.off()