perl/etc/justonce/sperm/mixplot.r

   1 #!/usr/bin/env littler
   2
   3 COLORS <- c("red","orange","#FF1EEA","black","#02CC7B","blue","#3399CC")
   4 SampleIDs <- c("Sperm23","Sperm24","Sperm28","Donor","SpermS01","SpermS02","SpermS03")
   5 COLORS <- c("gray","black","red","orange","#FF1EEA","#02CC7B","blue","#3399CC")
   6 SampleIDs <- c("MDA-Blood","MALBAC-Blood","MDA-Sperm23","MDA-Sperm24","MDA-Sperm28","MALBAC-SpermS01","MALBAC-SpermS02","MALBAC-SpermS03")
   7 #COLORS <- c("red","blue","green","black","orange","purple","gray")
   8 theLen <- length(SampleIDs)
   9
  10 ChrCount <- 22
  11 #ChrCount <- 1
  12
  13 DATa <- read.table("bamrsplot8.tsv.gz",skip=1)
  14 DATb <- read.table("rss1m.tsv.gz",skip=3)
  15 toplot <- function(value=3) {
  16         for(k in 1:ChrCount) {
  17 #               png(filename = paste0("mixplot.chr",k,".png"),
  18 #                 width = 3780, height = 2835, units = "px", pointsize = 96)
  19                 tiff(filename = paste0("mixplot.chr",k,".tiff"), compression="lzw",
  20                         width = 2049, height = 1800, units = "px", pointsize = 52)
  21                 par(mar=c(2, 4, 2, 0.5))     # c(bottom, left, top, right)
  22                 layout(rbind(1,2,3), heights=c(7,7,1))
  23                 data <- subset(DATa, V1==k);
  24                 ymax <- 800;
  25                 thelwd = value
  26                 plot(data[,2],data[,3]/1000,ylim=c(0,ymax),lwd=thelwd,
  27                      main=paste0("Chr",k," Reads distribution"),xlab="M",ylab="Reads Count (k)",type="l",col=COLORS[1],xaxs="r",yaxs="r")
  28
  29                 for (cc in 2:length(SampleIDs) ) {
  30                         lines(data[,2],data[,2+cc]/1000,type="l",col=COLORS[cc],lwd=thelwd)
  31                 }
  32
  33                 box(lwd=thelwd)
  34                 axis(1,lwd=thelwd)
  35                 axis(2,lwd=thelwd)
  36
  37                 chrk <- paste0('chr',k)
  38                 data <- subset(DATb, V1==chrk);
  39                 X <- data[,(3+3*theLen):(2+4*theLen)]
  40                 YY <- X[cbind(row(X)[which(!X == 0)], col(X)[which(!X == 0)])]
  41                 ValueMax <- mean( YY )
  42                 ymax <- 3*ValueMax
  43                 MinusRatio <- -25*ceiling(ValueMax/10)
  44                 plot(data[,2],(data[,3+theLen]*MinusRatio),ylim=c(MinusRatio,ymax), main=paste0("Chr",k," Bases distribution"),xlab="",ylab="UnCovered Rate          Base Coverage Depth",type="l",pch=1 ,col=COLORS[1],lwd=thelwd,yaxt='n')
  45                 lines(data[,2],data[,3+3*theLen],type="l",col=COLORS[1],lwd=thelwd)
  46                 for (cc in 2:length(SampleIDs) ) {
  47                         lines(data[,2],(data[,2+theLen+cc]*MinusRatio),type="l",col=COLORS[cc],lwd=thelwd)
  48                         lines(data[,2],data[,2+3*theLen+cc],type="l",col=COLORS[cc],lwd=thelwd)
  49                 }
  50                 box(lwd=thelwd)
  51                 axis(1,lwd=thelwd)
  52                 #axis(2,lwd=thelwd,labels=NA)
  53                 axis(2,lwd=thelwd,at=axTicks(2)[axTicks(2)>=0])
  54                 #TmpTmp <- c(axTicks(2)[axTicks(2)<0], MinusRatio)
  55                 TmpTmp <- c(.25,.5,.75)*MinusRatio
  56                 axis(2,lwd=thelwd,at=c(TmpTmp,MinusRatio),labels=c(paste0(100*TmpTmp/MinusRatio,'%'),1),las=1 )
  57
  58                 # setup for no margins on the legend
  59                 par(mar=c(0, 0, 0, 0))
  60                 # c(bottom, left, top, right)
  61                 plot.new()
  62                 legend('center','groups',SampleIDs, lty=seq(1,1,length=theLen),lwd=16,col=COLORS,bty="n",horiz=TRUE,cex=0.56)
  63                 dev.off()
  64         }
  65 }
  66
  67 toplot(3)