cis_matrix_body_multi.R

#!/bin/env Rscript

#'
#'  fixed global variables
#'

# load the libraries
source(file.path(library_path, "square_rotations.R"))
source(file.path(library_path, "bed_readers.R"))
source(file.path(library_path, "plot_helpers.R"))
source(file.path(library_path, "normalize.R"))

source("https://bioconductor.org/biocLite.R")
options(bitmapType='cairo')

#biocLite("GenomicRanges")

if(organism == "mm10"){

  if ( !(require("TxDb.Mmusculus.UCSC.mm10.knownGene")) ) {
    biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
  }
  library(TxDb.Mmusculus.UCSC.mm10.knownGene)
  txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene

} else if(organism == "hg19"){

  if ( !(require("TxDb.Hsapiens.UCSC.hg19.knownGene")) ) {
    biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
  }
  library(TxDb.Hsapiens.UCSC.hg19.knownGene)
  txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

} else if(organism == "mm9"){

  if ( !(require("TxDb.Mmusculus.UCSC.mm9.knownGene")) ) {
    biocLite("TxDb.Mmusculus.UCSC.mm9.knownGene")
  }
  library(TxDb.Mmusculus.UCSC.mm9.knownGene)
  txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene

} else {
  stop("Organsim not found")
}

require(plyr, quietly = TRUE)
require(RColorBrewer, quietly = TRUE)
require(gridExtra, quietly = TRUE)
require(grid, quietly = TRUE)
##require(ggplot2 ,quietly = TRUE)
require(ggpmisc, quietly = TRUE)
if ( !(require("ggbio")) ) {
  biocLite("ggbio")
}

require(devtools)
devtools::install_version("ggplot2", version = "2.2.1", repos = "http://cran.us.r-project.org")


ggplot2::theme_set(ggplot2::theme_bw())

# avoid genomic locations to be displayed in scientific notations
options(scipen=10)

#'
#' Functions for the analysis
#'
#'
#'

rotate_data <- function(d) {
  derived <- d
  m <- square_rotation(d$x, d$y)
  m <- scale_square_rotation(m[,1], m[,2])
  derived$x <- m[,1]
  derived$y <- m[,2]
  derived
}

#'
#' Main analysis loop
#'

#' Load the datasets
#'
message(paste("Reading data from", infile))
data <- read_paired_bed(infile)

#' Define the region if not previously determined
#'
if(is.null(region)){
  if(length(unique(c(data[,1], data[,4]))) != 1){
    stop(paste("Cannot infer region:", length(unique(data[,1])), "chromosomes found"))
  }
  region <- list(
    "chr" = data[1,1],
    "start" = min(c(data[,2], data[,5])),
    "end" = max(c(data[,3], data[,6])))
}


#' Select the data
#'
message(paste("Selecting data"))
tf.region <- select_region(data, region)
# tf.self <- is_self_ligation(data)
# tf.same_chr <- ndata$chr_a == ndata$chr_b
# tf.distant <- distance(ndata) > 1000

# combine the datasets
tf <- tf.region # & !tf.self # & tf.same_chr & tf.distant

# assign the dataset to plot
d <- data[tf,]


#' Prepare the plot
#'
message(paste("Transforming data"))
# explicitly define the plot area with a rectangle
bg <- background_data(region)
rotated_bg <- rotate_data(bg)
rotated_bg$y[rotated_bg$y < 0] <- 0
rotated_bg$panel <- "interactions"

#' Get the polygon shapes for the matrix plot
#'

square_data <- bedpe_to_polygon(d)

rotated_data <- rotate_data(square_data)
rotated_data$y[rotated_data$y < 0] <- 0
rotated_data$panel <- "interactions"

#' Plot the data
#'

if(colour_range[2] == 0){
  colour_range[2] = mean(data$signal)*2
}

# set the colour code if not previously defined
#if(is.null(colour_range)){
#  colour_range <- range(d$signal, na.rm=TRUE)
#}

wh <- GenomicRanges::GRanges(region[["chr"]], IRanges::IRanges(as.integer(region[["start"]]),as.integer(region[["end"]])))
tset <- ggbio::autoplot(txdb , which = wh, genome.axis = FALSE, names.expr = "gene_id")
p2 <- tset@ggplot


# get the colour palette
palette <- c(
  rgb(0,0,255,maxColorValue=255),
  rgb(152,20,18,maxColorValue=255),
  rgb(240,48,0,maxColorValue=255),
  rgb(255,188,39,maxColorValue=255),
  rgb(124,252,0,maxColorValue=255))

# triangle plot
message(paste("Constructing plot"))
p <- ggplot2::ggplot()
p <- p + geom_polygon(aes(x=x, y=y, group=group), fill="black", data=rotated_bg)
if(add_borders_to_elements){
  p <- p + geom_polygon(aes(x=x, y=y, group=group, fill=signal, colour=signal), data=rotated_data, size=0.2)
  p <- p + scale_colour_gradientn(limits=colour_range, colours=palette, na.value = "green", guide = "colourbar")
} else {
  p <- p + geom_polygon(aes(x=x, y=y, group=group, fill=signal), data=rotated_data)
}
p <- p + xlab(region[["chr"]]) + ylab("Interaction Matrix") + ggtitle(outfile)
p <- p + scale_y_continuous(labels=NULL)
p <- p + scale_fill_gradientn(limits=colour_range, colours=palette, na.value = "green")
p <- p + theme(
  panel.grid = element_blank(),
  panel.background = element_rect(fill = "white"),
  legend.position = "top",
  axis.ticks.y = element_blank(),
  axis.ticks.length=unit(0.005, "cm"))
p <- p + coord_fixed(xlim = c(region[["start"]], region[["end"]]))
p

TAD_data <- data.table::fread(tadfile, header = TRUE, sep = " ")
p3 <- ggplot2::ggplot(data = TAD_data, aes(x=coordinates,y=TADscore)) + geom_line() +
  labs(x = "", y = "TAD Boundary Score") +
  scale_y_continuous(breaks = c(0,1,2,3)) +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank(), axis.ticks.length=unit(-0.25, "cm"), axis.ticks.margin=unit(0.75, "cm"), text = element_text(size=9))


x_min = region[["start"]]
x_max = region[["end"]]
p1 <- p + scale_x_continuous(limits=c(x_min, x_max))
p2 <- p2 + scale_x_continuous(limits=c(x_min, x_max)) + labs(x = "Genomic coordinate (bp)", y = "Genes") + theme(axis.ticks.length=unit(0.05, "cm"), axis.ticks.margin=unit(0.2, "cm"))
#p3 <- p3 + scale_x_continuous(limits=c(x_min, x_max))

g1 <- ggplotGrob(p1)
g2 <- ggplotGrob(p2)
g3 <- ggplotGrob(p3)
g <- rbind(g1, g3 ,g2, size="first")
wm <- unit.pmax(g1$widths, g2$widths, g3$widths)
g$widths <- wm + unit(10, "cm")

png(outfile, width =  8, height = 8,units = "in",res = 700)
#grid.draw(g)
grid.arrange(
    arrangeGrob(g1,g3,g2,nrow=3,heights=c(.8,.2,.2))
    )
dev.off()