hic_matrixplot updates

[robinweide]

• [x] title-option
• [x] option to add names of samples (@teunbrand : this is now in there, right?)
• [ ] calculate the derivatives (add option to limit the distance-range to compute it on)
• [ ] pch and rect mode in bed-plotting

plotPCH <- function(dat, chrom, start, end, col = "#999999", pch = 6, revPCH = NULL, revCOL = NULL, cex = 1){
  subDat = dat[dat[,1] == chrom & dat[,2] >= start & dat[,3] <= end, ]
  subDat$center = apply(subDat[, 2:3], 1, mean)
  subDat$PCH = pch
  if(!is.null(revPCH)){
    subDat[subDat[,6] == '-', "PCH"] = revPCH
  }
  
  subDat$COL = col
  if(!is.null(revCOL)){
    subDat[subDat[,6] == '-', "COL"] = revCOL
  }
  
  
  plot( x = subDat$center, y = rep(0.25, nrow(subDat)), xlim = c(start, end), 
        ylim = c(0,0.5), pch = subDat$PCH, col = subDat$COL, 
        bg = subDat$COL, axes = T, ylab =NA, xlab = NA, cex = cex)
}

• [ ] annotate matrix (decisive diagonals)

annotate.hic.matrixplot = function(exp1, exp2 = NULL, chrom = CHROM, arrowArea = NULL, crossing0 = NULL, percentiles = c(25,75), SKIP = NULL){
  
  SWITCHPOINT = c()
  
################################################################################
# quantile-lines
################################################################################
  if(is.null(exp2)){
    # only percentages
    if(!is.null(crossing0)){
      message('crossing0 has been set to ',crossing0,', but no exp2 was given.')
    }
    MAT = select.subset(exp1, chrom = chrom, start = -Inf, end = Inf)$z
    
    TOTAL = sum(MAT)
    
    MMAT = melt(MAT)
    MMAT$D = abs(MMAT$Var1 - MMAT$Var2)
    MMAT = tapply(MMAT$value, INDEX = MMAT$D, sum)
    MMAT = cumsum(MMAT)
    
    if(percentiles[1] > 1){
      percentiles <- percentiles/100
    }
    QS = percentiles*TOTAL
      #quantile(MMAT, percentiles)
    
    for(q in QS){
      
      
      point = approx(x = MMAT, y = as.numeric(names(MMAT)), xout = q)$y
      abline(v = point)
      SWITCHPOINT = c(SWITCHPOINT, point)

    }

    
  }
  
################################################################################
# diff-lines
################################################################################
  if(!is.null(exp2)){
    # only percentages
    if(is.null(crossing0)){
      message('crossing0 has not been set to TRUE, but exp2 was given.')
    }
    
    DIFF <- select.subset(exp2, chrom = chrom, start = -Inf, end = Inf)$z -
      select.subset(exp1, chrom = chrom, start = -Inf, end = Inf)$z
    DIFF = melt(DIFF)
    DIFF$D = abs(DIFF$Var1 - DIFF$Var2)
    DIFF <- DIFF[DIFF$D < 50, ]
    DIFF = tapply(DIFF$value, INDEX = DIFF$D, mean)
    
    SIGNCHANGE = abs(diff(sign(DIFF)))
    SIGNCHANGE = as.numeric(names(SIGNCHANGE[SIGNCHANGE != 0]))
    
    
    for(i in SIGNCHANGE){
      if(i < 5){
        i = i + (5-i)
      }
      point = round(approx(y = as.numeric(names(DIFF[(i - 5):(i+5)])), x = DIFF[(i -5):(i+5)], xout = 0)$y)
      SWITCHPOINT = c(SWITCHPOINT, point)
    }
    
  }

################################################################################
# wrangling data
################################################################################  
  SWITCHPOINT = exp1$RES *  SWITCHPOINT
  if(!is.null(SKIP)){
    SWITCHPOINT <- SWITCHPOINT[-SKIP]
  }
  STARTS = NULL
  if(is.null(arrowArea )){
    centromere = exp1$CENTROMERES[exp1$CENTROMERES$V1 == chrom,]
    if(is.null(centromere)){
      message('No arrowArea given and no centromere-slot found in exp1.\n',
              'Randomly placing arrows...')
      
      STARTS = quantile(exp1$ABS[exp1$ABS$V1 == chrom,3],
                        seq(0, 1, length.out = length(SWITCHPOINT)+2))
      STARTS = STARTS[-1]
      STARTS = STARTS[-length(STARTS)]
      STARTS = unname(STARTS)
      
    } else {
      if(length(SWITCHPOINT) == 1){
        STARTS = seq(centromere$V2, centromere$V3, length.out = 3 )
        STARTS = STARTS[-1]
        STARTS = STARTS[-length(STARTS)]
      } else if(length(SWITCHPOINT) > 1){
        STARTS = seq(centromere$V2, centromere$V3, length.out = 2+length(SWITCHPOINT) )
        STARTS = STARTS[-1]
        STARTS = STARTS[-length(STARTS)]
      }
    }
  } else {
    arrowArea = arrowArea[arrowArea[,1] == chrom,]
    if(nrow(arrowArea ) == 0){
      message('arrowArea does not contain this chromosome.\n',
              'Randomly placing arrows...')
      
      STARTS = quantile(exp1$ABS[exp1$ABS$V1 == chrom,3],
                    seq(0, 1, length.out = length(SWITCHPOINT)+2))
      STARTS = STARTS[-1]
      STARTS = STARTS[-length(STARTS)]
      STARTS = unname(STARTS)
      
    } else{
      if(length(SWITCHPOINT) == 1){
        STARTS = seq(arrowArea[1,2], arrowArea[1,3], length.out = 3 )
        STARTS = STARTS[-1]
        STARTS = STARTS[-length(STARTS)]
      } else if(length(SWITCHPOINT) > 1){
        STARTS = seq(arrowArea[1,2], arrowArea[1,3], length.out = 2+length(SWITCHPOINT) )
        STARTS = STARTS[-1]
        STARTS = STARTS[-length(STARTS)]
      }
    }
  }

################################################################################
# plotting data
################################################################################    
  
 
  library(shape)
  abline(a = 0, b = 1)
  for(i in SWITCHPOINT){
    if(length(ARROWSHIFT) == 1 & length(SWITCHPOINT > 1)){
      ARROWSHIFT = rep(ARROWSHIFT, length(SWITCHPOINT))
    }
    #abline(a = i, b = 1, lty = which(i == SWITCHPOINT))
    abline(a = -i, b = 1, lty = which(i == SWITCHPOINT)+1)
    
    P1 = (sqrt((i**2) + (i**2)))
    
    locs = c(P1,P1, P1+(i/2), P1-(i/2))+ ( STARTS[match(i ,SWITCHPOINT)] - P1)
    
    Arrows(x0 =locs[1], y0 = locs[2], x1 = locs[3], y1 = locs[4],lty = match(i ,SWITCHPOINT) +1,
           code = 2, arr.type	= 'triangle', arr.adj = 1, arr.length = .2,arr.width = .12)
    Arrows(x0 =locs[1], y0 = locs[2], x1 = locs[3], y1 = locs[4],lty = 1, lcol = NULL,
           code = 2, arr.type	= 'triangle', arr.adj = 1, arr.length = .2,arr.width = .12)
  }
  legend('topright', col = 'black', lty = seq(1:length(SWITCHPOINT))+1, legend =  paste0(round(SWITCHPOINT/1e6,2), "Mb"))
}

[teunbrand]

The names of samples are there for the insulation_/compartment_matrixplot functions, but it'd be trivial to do it too for the hic_matrixplot. I'll include in one of the next PRs

[robinweide]

it'd be trivial to do it too for the hic_matrixplot

Not really, since you can have two contacts-objects. Optional names would make it handy to quickly see which is which.

[teunbrand]

With trivial I mean it might be 5 lines of code, 7 at max. I'm not saying it isn't handy :D