#' @title plotGene
#' @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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.