new file: mixplot.R & mixplotBW.R
[galaxycodebases:main.git] / perl / etc / justonce / sperm / mixplotBW.R
1 DATa <- read.table("bamrsplot.tsv",skip=1)
2 DATb <- read.table("depstat/nrss1m.tsv",skip=3)
3 for(k in c(1:2)) {
4      png(filename = paste0("mixplot.chr",k,".png"),
5           width = 3780, height = 2835, units = "px", pointsize = 96)
6      par(mar=c(2, 4, 2, 0.5))     # c(bottom, left, top, right)
7      layout(rbind(1,2,3), heights=c(7,7,1))
8           data <- subset(DATa, V1==k);
9           ymax <- 400;
10           thelwd = 8
11           plot(data[,2],data[,3]/1000,ylim=c(0,ymax),lwd=thelwd,
12                main=paste0("Chr",k," Reads distribution"),xlab="M",ylab="Reads Count (k)",type="l",col="#999999",xaxs="r",yaxs="r")
13           lines(data[,2],data[,4]/1000,type="l",col="#AAAAAA",lwd=thelwd)
14           lines(data[,2],data[,5]/1000,type="l",col="#CCCCCC",lwd=thelwd)
15           #lines(data[,2],data[,6]/1000,type="l",col="black",lwd=thelwd)
16           lines(data[,2],data[,7]/1000,type="l",col="#222222",lwd=thelwd)
17           lines(data[,2],data[,8]/1000,type="l",col="#444444",lwd=thelwd)
18           lines(data[,2],data[,9]/1000,type="l",col="#555555",lwd=thelwd)
19          
20           lines(data[,2],data[,6]/1000,type="l",col="black",lwd=thelwd)
21           box(lwd=thelwd)
22           axis(1,lwd=thelwd)
23           axis(2,lwd=thelwd)
24
25           data <- subset(DATb, V1==k);
26           X <- data[,24:30]
27           YY <- X[cbind(row(X)[which(!X == 0)], col(X)[which(!X == 0)])]
28           ValueMax <- mean( YY )
29           ymax <- 5*ValueMax
30           MinusRatio <- -25*ceiling(ValueMax/10)
31           plot(data[,2],(data[,10]*MinusRatio),ylim=c(MinusRatio,ymax), main=paste0("Chr",k," Bases distribution"),xlab="",ylab="UnCovered Rate      Base Coverage Depth",type="l",pch=1 ,col="#999999",lwd=thelwd,yaxt='n')
32           lines(data[,2],data[,24],type="l",col="#999999",lwd=thelwd)
33           lines(data[,2],(data[,11]*MinusRatio),type="l",col="#AAAAAA",lwd=thelwd)
34           lines(data[,2],data[,25],type="l",col="#AAAAAA",lwd=thelwd)
35           lines(data[,2],(data[,12]*MinusRatio),type="l",col="#CCCCCC",lwd=thelwd)
36           lines(data[,2],data[,26],type="l",col="#CCCCCC",lwd=thelwd)
37           #lines(data[,2],(data[,13]*MinusRatio),type="l",col="black",lwd=thelwd)
38           #lines(data[,2],data[,27],type="l",col="black",lwd=thelwd)
39           lines(data[,2],(data[,14]*MinusRatio),type="l",col="#222222",lwd=thelwd)
40           lines(data[,2],data[,28],type="l",col="#222222",lwd=thelwd)
41           lines(data[,2],(data[,15]*MinusRatio),type="l",col="#444444",lwd=thelwd)
42           lines(data[,2],data[,29],type="l",col="#444444",lwd=thelwd)
43           lines(data[,2],(data[,16]*MinusRatio),type="l",col="#555555",lwd=thelwd)
44           lines(data[,2],data[,30],type="l",col="#555555",lwd=thelwd)
45
46           lines(data[,2],(data[,13]*MinusRatio),type="l",col="black",lwd=thelwd)
47           lines(data[,2],data[,27],type="l",col="black",lwd=thelwd)
48
49           box(lwd=thelwd)
50           axis(1,lwd=thelwd)
51           #axis(2,lwd=thelwd,labels=NA)
52           axis(2,lwd=thelwd,at=axTicks(2)[axTicks(2)>=0])
53           #TmpTmp <- c(axTicks(2)[axTicks(2)<0], MinusRatio)
54           TmpTmp <- c(.25,.5,.75)*MinusRatio
55           axis(2,lwd=thelwd,at=c(TmpTmp,MinusRatio),labels=c(paste0(100*TmpTmp/MinusRatio,'%'),1),las=1 )
56
57      # setup for no margins on the legend
58      par(mar=c(0, 0, 0, 0))
59      # c(bottom, left, top, right)
60      plot.new()
61      legend('center','groups',c("Sperm23","Sperm24","Sperm28","Donor","SpermS01","SpermS02","SpermS03"), lty=seq(1,1,length=7),lwd=24,col=c("#999999","#AAAAAA","#CCCCCC","black","#222222","#444444","#555555"),bty="n",horiz=TRUE)
62      dev.off()
63 }