R/PlotGene.R

In scottzijiezhang/MeRIPtools: Analyze MeRIP-seq data and QTL calling

#' @title old coverage plot function.
#' @param IP_BAM The bam files for IP samples
#' @param INPUT_BAM The bam files for INPUT samples
#' @param size.IP The size factor for IP libraries
#' @param size.INPUT The size factor for INPUT libraries
#' @param geneName The name (as defined in gtf file) of the gene you want to plot
#' @param geneModel The gene model generated by gtfToGeneModel() function
#' @param libraryType "opposite" for mRNA stranded library, "same" for samll RNA library
#' @param GTF gtf annotation as GRanges object. Can be obtained by GTF <- rtracklayer::import("xxx.gtf",format = "gtf")
#' @export
## the main function to plot m6A-seq on one group of data
plotGene <- function(IP_BAM, INPUT_BAM, size.IP, size.INPUT, geneName, geneModel, libraryType = "opposite", center = mean ,GTF,ZoomIn=NULL){
  IP.cov <- getAveCoverage(geneModel= geneModel,bamFiles = IP_BAM,geneName = geneName,size.factor = size.IP, libraryType = libraryType, center = center, ZoomIn = ZoomIn)
  INPUT.cov <- getAveCoverage(geneModel= geneModel,bamFiles = INPUT_BAM,geneName = geneName,size.factor = size.INPUT, libraryType = libraryType, center = center,ZoomIn = ZoomIn)
  cov.data <- data.frame(IP=IP.cov,Input=INPUT.cov,genome_location=as.numeric(names(IP.cov) ) )
  yscale <- max(IP.cov,INPUT.cov)
  p1 <- "ggplot(data = cov.data,aes(genome_location))+geom_line(aes(y=Input,colour =\"Input\"))+geom_line(aes(y=IP,colour=\"IP\"))+labs(y=\"normalized coverage\")+scale_x_continuous(breaks = round(seq(min(cov.data$genome_location), max(cov.data$genome_location), by = ((max(cov.data$genome_location)-min(cov.data$genome_location))/10) ),1))+
    theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank(), axis.line = element_line(colour = \"black\"))"

  p2 <- .getGeneModelAnno(geneModel,geneName,GTF,ZoomIn)
  p <- paste(p1,p2,sep = "+")
  eval(parse( text = p ))
}

## the main function to plot m6A-seq on two group of data
#' @title plotGenePair
#' @description plot tow groups of samples in the same figure
#' @param Ctl_IP_BAM The bam files for Control IP samples
#' @param Ctl_INPUT_BAM The bam files for Control INPUT samples
#' @param Treat_IP_BAM The bam files for treated IP samples
#' @param Treat_INPUT_BAM The bam files for treated INPUT samples
#' @param Ctl_size.IP The size factor for IP libraries
#' @param Ctl_size.INPUT The size factor for INPUT libraries
#' @param Treat_size.IP The size factor for IP libraries
#' @param Treat_size.INPUT The size factor for INPUT libraries
#' @param geneName The name (as defined in gtf file) of the gene you want to plot
#' @param geneModel The gene model generated by gtfToGeneModel() function
#' @param libraryType "opposite" for mRNA stranded library, "same" for samll RNA library
#' @export
plotGenePair <- function(Ctl_IP_BAM,Ctl_INPUT_BAM,Treat_IP_BAM,Treat_INPUT_BAM,Ctl_size.IP,Ctl_size.INPUT,Treat_size.IP,Treat_size.INPUT,geneName,geneModel, libraryType = "ooposite",center = mean, GTF ,ZoomIn=NULL){
  Ctl_IP.cov <- getAveCoverage(geneModel= geneModel,bamFiles = Ctl_IP_BAM,geneName = geneName,size.factor = Ctl_size.IP, libraryType = libraryType,center = center, ZoomIn = ZoomIn)
  Ctl_INPUT.cov <- getAveCoverage(geneModel= geneModel,bamFiles = Ctl_INPUT_BAM,geneName = geneName,size.factor = Ctl_size.INPUT,libraryType = libraryType, center = center , ZoomIn = ZoomIn)
  Treat_IP.cov <- getAveCoverage(geneModel= geneModel,bamFiles = Treat_IP_BAM,geneName = geneName,size.factor = Treat_size.IP, libraryType = libraryType, center = center,ZoomIn = ZoomIn)
  Treat_INPUT.cov <- getAveCoverage(geneModel= geneModel,bamFiles = Treat_INPUT_BAM,geneName = geneName,size.factor = Treat_size.INPUT, libraryType = libraryType, center = center,ZoomIn = ZoomIn)
  cov.data <- data.frame(Ctl_IP=Ctl_IP.cov, Ctl_Input = Ctl_INPUT.cov,
                         Treat_IP=Treat_IP.cov, Treat_Input = Treat_INPUT.cov,
                         genome_location=as.numeric(names(Ctl_IP.cov) ) )
  yscale <- max(Ctl_IP.cov,Ctl_INPUT.cov,Treat_IP.cov,Treat_INPUT.cov)
  p1 <- "ggplot(data = cov.data,aes(genome_location))+geom_line(aes(y=Ctl_Input,colour =\"Ctl Input\"))+geom_line(aes(y=Treat_IP,colour=\"Treat IP\"))+geom_line(aes(y=Treat_Input,colour =\"Treat Input\"))+geom_line(aes(y=Ctl_IP,colour=\"Ctl IP\"))+labs(y=\"normalized coverage\")+scale_x_continuous(breaks = round(seq(min(cov.data$genome_location), max(cov.data$genome_location), by = ((max(cov.data$genome_location)-min(cov.data$genome_location))/10) ),1))+
    theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank(), axis.line = element_line(colour = \"black\"))"
  p2 <- .getGeneModelAnno(geneModel,geneName,GTF,ZoomIn)
  p <- paste(p1,p2,sep = "+")
  eval(parse( text = p ))

}


#' @title plotGeneMonster
#' @param readsOut The data list from countReads and other analysis.
#' @param geneName The gene symbol to be ploted.
#' @param GTF The GRanges object containing gtf annotation.
#' @param ZoomIn c(start,end) The coordinate to zoom in at the gene to be ploted.
#' @export
plotGeneMonster <- function(readsOut, geneName, libraryType = "opposite", center = "mean", GTF, ZoomIn = NULL){
  if("X" %in% names(readsOut) ){
    X <- readsOut$X
    plotGenePair(Ctl_IP_BAM = readsOut$bamPath.ip[X == unique(X)[1]],
                 Ctl_INPUT_BAM = readsOut$bamPath.input[X == unique(X)[1]],
                 Treat_IP_BAM = readsOut$bamPath.ip[X == unique(X)[2]],
                 Treat_INPUT_BAM = readsOut$bamPath.input[X == unique(X)[2]],
                 Ctl_size.IP = readsOut$sizeFactor$ip[X == unique(X)[1]],
                 Ctl_size.INPUT = readsOut$sizeFactor$input[X == unique(X)[1]],
                 Treat_size.IP = readsOut$sizeFactor$ip[X == unique(X)[2]],
                 Treat_size.INPUT = readsOut$sizeFactor$input[X == unique(X)[2]],
                 geneName = geneName,
                 geneModel = readsOut$geneModel,
                 libraryType = libraryType,center = center,GTF = GTF,ZoomIn = ZoomIn )
  }else{
    plotGene(IP_BAM = readsOut$bamPath.ip,
             INPUT_BAM =  readsOut$bamPath.input,
             size.IP = readsOut$sizeFactor$ip,
             size.INPUT = readsOut$sizeFactor$input,
             geneName = geneName,
             geneModel = readsOut$geneModel,
             libraryType = libraryType,center = center,GTF = GTF,ZoomIn = ZoomIn)
  }
}