pdf('output.qual.boxplot.pdf')
p0<-rep(c(33,43,45,51,54,58,59,60,61,62,63,64,66,67,69,70,71),times=c(1,2,1,3,1,1,1,1,1,3,1,1,1,2,5,5,10)/1000)
p1<-rep(c(43,45,51,56,57,58,60,61,62,63,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,2,1,1,2,1,3,3,1,8,5,6)/1000)
p2<-rep(c(43,49,51,53,54,56,58,59,60,61,64,65,66,67,69,70,71),times=c(1,1,1,1,2,1,1,1,2,1,1,2,2,2,7,6,8)/1000)
p3<-rep(c(33,39,53,54,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,1,1,1,2,1,1,1,4,4,1,5,5,6)/1000)
p4<-rep(c(33,55,58,59,60,61,62,64,65,66,67,68,69,70,71),times=c(1,1,1,1,5,1,1,4,1,2,5,2,3,3,9)/1000)
p5<-rep(c(33,53,58,60,61,62,64,65,67,68,69,70,71),times=c(1,1,1,2,1,4,3,2,2,3,5,4,11)/1000)
p6<-rep(c(33,40,54,56,58,60,64,66,67,68,69,70,71),times=c(3,2,1,1,1,1,2,2,4,4,4,2,13)/1000)
p7<-rep(c(41,42,43,49,50,51,57,58,62,64,65,66,67,68,69,70,71),times=c(1,1,1,1,2,1,1,1,2,1,1,2,7,1,2,3,12)/1000)
p8<-rep(c(33,39,53,56,58,59,60,62,63,66,67,68,69,70,71),times=c(1,1,1,1,1,1,2,1,2,4,3,2,3,5,12)/1000)
p9<-rep(c(33,40,50,52,53,57,58,60,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,3,1,4,1,1,3,2,5,4,10)/1000)
p10<-rep(c(33,40,53,54,55,58,60,63,64,66,67,68,69,70,71),times=c(2,1,1,2,1,1,2,1,4,1,1,1,7,6,9)/1000)
p11<-rep(c(45,50,51,52,53,56,57,58,60,62,63,64,65,66,67,68,69,70,71),times=c(1,2,1,1,1,1,1,1,3,2,1,1,1,3,1,1,3,5,10)/1000)
p12<-rep(c(33,41,52,53,54,58,59,60,64,65,66,67,69,70,71),times=c(3,1,1,2,1,1,1,3,1,2,3,2,5,6,8)/1000)
p13<-rep(c(33,40,51,53,55,56,58,59,60,63,64,65,66,67,68,69,70,71),times=c(4,1,1,1,1,1,1,1,1,1,3,1,2,2,2,3,5,9)/1000)
p14<-rep(c(33,39,54,56,57,59,60,61,62,63,65,66,67,68,69,70,71),times=c(4,1,3,1,1,1,1,2,1,1,1,2,2,1,7,2,9)/1000)
p15<-rep(c(33,39,40,42,45,50,52,53,57,58,59,60,61,64,66,68,69,70,71),times=c(2,2,1,1,1,1,1,1,1,3,1,1,1,3,1,1,4,5,9)/1000)
p16<-rep(c(33,47,51,52,53,58,59,60,61,64,67,69,70,71),times=c(4,1,1,1,2,3,1,1,2,3,3,7,2,9)/1000)
p17<-rep(c(33,48,50,51,53,54,55,56,58,60,61,63,64,65,66,67,69,70,71),times=c(1,1,1,1,2,1,1,1,2,2,3,2,2,1,2,1,7,2,7)/1000)
p18<-rep(c(33,43,48,51,53,58,59,60,61,63,64,66,67,68,69,70,71),times=c(2,1,1,1,2,2,1,2,2,3,1,1,3,1,7,2,8)/1000)
p19<-rep(c(33,44,47,50,51,52,54,58,59,61,62,64,65,67,69,70,71),times=c(2,1,1,2,1,1,2,1,1,1,1,3,1,1,8,4,9)/1000)
p20<-rep(c(33,46,47,51,54,56,58,59,61,62,63,64,66,67,69,70,71),times=c(1,1,1,1,2,1,1,2,1,2,1,5,5,3,5,2,6)/1000)
p21<-rep(c(33,43,54,55,57,58,62,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,5,2,3,1,3,2,4,5,4,6)/1000)
p22<-rep(c(33,47,51,53,54,57,58,60,62,63,64,65,66,68,69,70,71),times=c(1,1,1,1,1,1,1,1,2,1,5,1,4,3,5,5,6)/1000)
p23<-rep(c(33,42,53,54,55,57,58,62,63,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,2,1,1,5,2,2,3,2,9,3,4)/1000)
p24<-rep(c(33,42,52,54,57,60,61,63,64,65,66,67,69,70,71),times=c(1,1,1,1,1,2,1,1,5,1,6,4,5,4,6)/1000)
p25<-rep(c(33,53,54,57,61,62,63,64,66,67,68,69,70,71),times=c(1,1,1,2,1,1,1,2,4,5,2,9,5,5)/1000)
p26<-rep(c(46,53,54,57,58,60,61,62,64,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,2,5,8,4,1,5,3,5)/1000)
p27<-rep(c(42,43,48,54,56,57,60,61,62,66,67,68,69,70,71),times=c(1,1,1,2,1,1,4,1,2,1,4,1,5,6,9)/1000)
p28<-rep(c(51,55,56,57,60,62,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,2,4,2,2,2,4,10,1,8)/1000)
p29<-rep(c(49,52,56,57,58,60,63,64,65,66,67,69,70,71),times=c(2,1,1,1,2,1,1,3,1,3,6,8,2,8)/1000)
p30<-rep(c(45,47,50,57,61,62,64,66,67,68,69,70,71),times=c(1,1,1,1,1,1,2,6,3,2,8,5,8)/1000)
p31<-rep(c(48,52,53,54,57,58,59,60,61,62,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,1,2,1,4,1,2,3,3,7,3,6)/1000)
p32<-rep(c(43,47,48,54,56,62,64,66,67,68,69,70,71),times=c(1,1,1,1,2,1,2,1,5,3,10,5,7)/1000)
p33<-rep(c(52,55,58,60,61,63,64,68,69,70,71),times=c(1,1,2,1,1,1,5,4,11,5,8)/1000)
p34<-rep(c(42,43,50,56,59,60,63,64,67,68,69,70,71),times=c(1,1,1,1,1,1,1,3,1,4,9,5,11)/1000)
p35<-rep(c(42,53,57,58,60,64,66,68,69,70,71),times=c(1,1,1,2,1,3,2,2,12,7,8)/1000)
p36<-rep(c(48,53,56,61,63,64,66,67,69,70,71),times=c(2,1,1,2,1,1,2,6,7,3,14)/1000)
p37<-rep(c(53,56,60,63,64,66,68,69,70,71),times=c(2,2,1,1,3,2,6,7,8,8)/1000)
p38<-rep(c(41,48,53,57,59,61,62,63,64,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,1,3,3,4,1,4,6,11)/1000)
p39<-rep(c(38,42,51,53,56,58,61,63,64,65,66,67,68,69,70,71),times=c(1,1,1,1,1,1,1,1,2,1,4,2,2,7,3,11)/1000)
p40<-rep(c(53,58,61,62,63,64,66,67,68,69,70,71),times=c(1,1,3,1,2,2,1,2,2,9,4,12)/1000)
p41<-rep(c(48,53,54,57,58,59,60,61,63,64,66,67,68,69,70,71),times=c(1,1,1,1,1,1,2,1,1,1,3,3,1,5,7,9)/1000)
p42<-rep(c(38,49,54,58,59,64,65,66,67,68,69,70,71),times=c(1,1,3,1,1,3,1,2,1,1,7,7,9)/1000)
p43<-rep(c(50,51,62,63,64,65,66,67,69,70,71),times=c(2,1,1,1,3,1,2,4,3,8,12)/1000)
p44<-rep(c(48,54,63,64,65,66,67,68,69,70,71),times=c(1,2,4,1,1,2,2,1,7,8,8)/1000)
p45<-rep(c(50,57,58,59,60,62,64,67,69,70,71),times=c(1,1,1,1,1,1,1,1,10,8,7)/1000)
p46<-rep(c(43,48,54,59,60,64,65,66,67,68,69,70,71),times=c(2,1,1,2,1,2,1,1,4,2,1,6,8)/1000)
p47<-rep(c(49,53,56,61,64,66,67,69,70,71),times=c(1,1,1,1,2,2,2,3,10,7)/1000)
p48<-rep(c(61,64,66,67,68,69,70,71),times=c(2,1,2,2,2,6,5,7)/1000)
p49<-rep(c(33,56,60,64,66,68,69,70,71),times=c(1,1,2,1,1,2,2,5,10)/1000)
p50<-rep(c(33,66,67,68,69,70,71),times=c(1,1,1,2,4,5,7)/1000)
boxplot(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,p26,p27,p28,p29,p30,p31,p32,p33,p34,p35,p36,p37,p38,p39,p40,p41,p42,p43,p44,p45,p46,p47,p48,p49,p50,xlab="Position of Read(5'->3')",ylab="Phred Quality Score",outline=F)
dev.off()


pdf('output.qual.heatmap.pdf')
qual=c(1,0,0,0,0,0,0,0,0,0,2,0,1,0,0,0,0,0,3,0,0,1,0,0,0,1,1,1,1,3,1,1,0,1,2,0,5,5,10,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,1,1,1,0,1,2,1,1,2,1,3,3,1,8,5,6,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,1,2,0,1,0,1,1,2,1,0,0,1,2,2,2,0,7,6,8,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,1,1,1,1,2,1,1,1,4,4,1,5,5,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,5,1,1,0,4,1,2,5,2,3,3,9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,2,1,4,0,3,2,0,2,3,5,4,11,3,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,2,0,2,4,4,4,2,13,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,1,2,1,0,0,0,0,0,1,1,0,0,0,2,0,1,1,2,7,1,2,3,12,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,2,0,1,2,0,0,4,3,2,3,5,12,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,1,3,0,1,0,0,0,4,1,1,3,2,5,4,10,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,2,1,0,0,1,0,2,0,0,1,4,0,1,1,1,7,6,9,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,2,1,1,1,0,0,1,1,1,0,3,0,2,1,1,1,3,1,1,3,5,10,3,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,2,1,0,0,0,1,1,3,0,0,0,1,2,3,2,0,5,6,8,4,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,1,0,1,1,1,0,0,1,3,1,2,2,2,3,5,9,4,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,1,1,0,1,1,2,1,1,0,1,2,2,1,7,2,9,2,0,0,0,0,0,2,1,0,1,0,0,1,0,0,0,0,1,0,1,1,0,0,0,1,3,1,1,1,0,0,3,0,1,0,1,4,5,9,4,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,2,0,0,0,0,3,1,1,2,0,0,3,0,0,3,0,7,2,9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,2,1,1,1,0,2,0,2,3,0,2,2,1,2,1,0,7,2,7,2,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,2,0,0,0,0,2,1,2,2,0,3,1,0,1,3,1,7,2,8,2,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,2,1,1,0,2,0,0,0,1,1,0,1,1,0,3,1,0,1,0,8,4,9,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,2,0,1,0,1,2,0,1,2,1,5,0,5,3,0,5,2,6,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,1,0,1,5,0,0,0,2,0,3,1,3,2,4,5,4,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,1,1,0,1,0,2,1,5,1,4,0,3,5,5,6,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,2,0,0,0,1,1,5,2,2,3,2,9,3,4,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,2,1,0,1,5,1,6,4,0,5,4,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,2,0,0,0,1,1,1,2,0,4,5,2,9,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,2,0,5,0,8,4,1,5,3,5,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,0,2,0,1,1,0,0,4,1,2,0,0,0,1,4,1,5,6,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,1,0,2,0,4,2,2,2,4,10,1,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,1,0,0,0,1,1,2,0,1,0,0,1,3,1,3,6,0,8,2,8,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,1,1,0,2,0,6,3,2,8,5,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,2,1,0,4,1,2,3,3,7,3,6,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,1,0,2,0,0,0,0,0,1,0,2,0,1,5,3,10,5,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,2,0,1,1,0,1,5,0,0,0,4,11,5,8,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,1,0,0,1,3,0,0,1,4,9,5,11,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,2,0,1,0,0,0,3,0,2,0,2,12,7,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,1,0,0,1,0,0,0,0,2,0,1,1,0,2,6,0,7,3,14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,2,0,0,0,1,0,0,1,3,0,2,0,6,7,8,8,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,1,1,3,0,3,4,1,4,6,11,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,2,1,4,2,2,7,3,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,3,1,2,2,0,1,2,2,9,4,12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0,0,1,1,1,2,1,0,1,1,0,3,3,1,5,7,9,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,3,0,0,0,1,1,0,0,0,0,3,1,2,1,1,7,7,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0,0,0,1,1,3,1,2,4,0,3,8,12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,2,0,0,0,0,0,0,0,0,4,1,1,2,2,1,7,8,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,1,1,0,1,0,1,0,0,1,0,10,8,7,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,2,1,0,0,0,2,1,1,4,2,1,6,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,1,0,0,2,0,2,2,0,3,10,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,1,0,2,2,2,6,5,7,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,2,0,0,0,1,0,1,0,2,2,5,10,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,2,4,5,7)
mat=matrix(qual,ncol=51,byrow=F)
Lab.palette <- colorRampPalette(c("blue", "orange", "red3","red2","red1","red"), space = "rgb",interpolate=c('spline'))
heatmap(mat,Rowv=NA,Colv=NA,xlab="Position of Read",ylab="Phred Quality Score",labRow=seq(from=33,to=71),col = Lab.palette(256),scale="none" )
dev.off()