# TODO: Add comment # # Author: E.Korsching 2016 ############################################################################### #### two functions: plot.exon.levels(), genomic.presentation.scale() plot.exon.levels <- function(x, y=NULL, gene.symbol, sig=NULL, logt="", y.name="normalized counts", cex=0.7, lcex=1, legend=T, col.v=c("blue","red")){ # x: data.frame column with exon counts of one gene in condition 1 # y: data.frame column with exon counts of the same gene in condition 2 # if you want to show only one exon series - take same colors and give x+y the same values # grep: gene name needs to be in the beginning of the string and an underscore is added to mark the end of the name # both own row.names with a format: name,underscore,chr,underscore,start,underscore,end : "A1BG-AS1_chr19_58351970_58353044" # gene.symbol: a name (case sensitive) # sig: exon number or vector of exon numbers with differential exons -> different color # logt: "y": log on y axis xlab <- row.names(x) ylab <- row.names(y) x <- x[grepl(pattern=paste("^",gene.symbol,"_",sep=""), x=xlab),,drop=F] y <- y[grepl(pattern=paste("^",gene.symbol,"_",sep=""), x=ylab),,drop=F] xlab <- row.names(x) ylab <- row.names(y) if(xlab[1]!=ylab[1]){ stop("gene/exon names differ") } nrx <- nrow(x) nry <- nrow(y) if(nrx!=nry){ stop("x and y length differ") } ncx <- ncol(x) ncy <- ncol(y) f1 <- function(x){ # not for big jobs - rbind too slow xlen <- length(x) tmp <- c("0","0","0","0") for(i in 1:xlen){ tmp <- rbind(tmp,x[[i]]) } return(tmp[-1,]) } anno <- data.frame(f1(strsplit(xlab, split="_")), stringsAsFactors=F) # build annotation from 'speaking' row names names(anno) <- c("gene","chr","begin","end") # set data types anno <- transform(anno, begin=as.numeric(begin), end=as.numeric(end) ) # sort all ascending o.data <- order(anno[,3]) cat("\n x order changed ? ",o.data,"\n") anno <- anno[o.data,,drop=F] x <- x[o.data,,drop=F] y <- y[o.data,,drop=F] # return(anno) ## start plotting layout( matrix(c(1,2),nrow=2,ncol=1), width=1, height=c(1,1) ) # layout.show(2) ## level plot xlab <- paste("exon ",seq(1,nrx),sep="") # create new x labels nx <- c(1:nrx) # x positions ma.y <- max(x,y) # detect range mi.y <- min(x,y) myoffset <- -ma.y/3 plot(1, 1, log=logt, type="n", xlim=c(0,nrx+1), ylim=c(myoffset,ma.y), xlab="", ylab="", cex=cex, adj=0, axes=F) for(i in 1:ncx){ segments(x0=nx-0.3, y0=x[,i], x1=nx, y1=x[,i], col=col.v[1], lwd=1.2) } for(i in 1:ncy){ segments(x0=nx, y0=y[,i], x1=nx+0.3, y1=y[,i], col=col.v[2], lwd=1.2) } axis(side=1, at=nx, labels=F, pos=-1, cex=cex) text(x=nx, y=myoffset/2, labels=xlab, adj=c(1,1), srt=45, xpd=T, cex=cex) mtext(text=paste(anno[1,2],", exons of the gene: ",gene.symbol,sep=""), side=1, line=1.5, cex=lcex) axis(side=2, at=c(0,round(ma.y/2,0),round(ma.y,0)),las=1, cex=cex) mtext(text=y.name, side=2, line=3, adj=1, cex=lcex) if(legend){ legend(x=0, y=ma.y, legend=c(paste("condition 1 [x]",sep=""), paste("condition 2 [y]",sep="")), density=-1, fill=col.v, yjust=1, border="black", bty="n", cex=cex) } ## transcript plot min.pos <- min(anno[ ,3]) max.pos <- max(anno[ ,4]) # customize scale of genomic scheme anno2 <- genomic.presentation.scale(x=anno, begin=3, end=4) min.npos <- min(anno2[ ,3]) max.npos <- max(anno2[ ,4]) plot(1, 1, type="n", xlim=c(min.npos,max.npos), ylim=c(3,0), xlab="", ylab="", cex=cex, adj=0, axes=F) lines(x=c(min.npos,max.npos),y=c(0.5,0.5)) for(i in 1:nrx){ xb <- c(anno2[i,3], anno2[i,4], anno2[i,4], anno2[i,3]) yb <- c(0, 0, 1, 1) if(!is.null(sig)){ if(i %in% sig){ cole <- "yellow" }else{ cole <- "green" } }else{ cole <- "green" } polygon(xb, yb, density=-1, col=cole, border="black", lwd=1.2) } mtext(text=paste(min.pos," to ",max.pos," [scaled]",sep=""), side=3, line=1.5, cex=lcex) # return() } genomic.presentation.scale <- function(x, begin=3, end=4){ # customize scale of a genomic gene/exon scheme # to fit into a small view # preserve all exonic and separately all intronic relations # but give each group 50% of the space # sum of all exons and sum of all introns = 0.5 each # scaling of all introns and exons each by: exon-length * 0.5 / sum(all exon length) rsp. introns # create again a scaled exon table # x: table of exonic start end positions # begin: column number of start position, end: column number of end position nr <- nrow(x) x.names <- names(x) x[,5] <- x[,4]-x[,3]+1 # length of exons y <- (x[2:nr,3]-x[1:(nr-1),4])-1 # length of introns x.sum <- sum(x[,5]) y.sum <- sum(y) x[,6] <- x[,5]/x.sum # relative exon length y2 <- y/y.sum # relative intron length max.sum <- max(x.sum,y.sum) # normalization to the larger one (more precise) x[,7] <- ceiling(x[,6] * max.sum) # new exon length y3 <- ceiling(y2 * max.sum) # new intron length k <- 1 for(i in 1:nr){ x[i,8] <- k k <- k + x[i,7] x[i,9] <- k if(i!=nr){ k <- k + y3[i] } } x <- x[,-c(3:7)] names(x) <- x.names return(x) } #genomic.presentation.scale(x=ab, begin=3, end=4) # test with first return() uncommented #aa <- c("ab_chr1_47219720_47219722","ab_chr1_167936912_167937008","ab_chr1_167951800_167951861","ab_chr1_167966629_167966721", # "ab_chr1_167974830_167975015","ab_chr1_167987495_167987608","ab_chr1_167991204_167991339","ab_chr1_167993226_167993440", # "ab_chr1_168002482_168002575") #plot.exon.levels(x=data.frame(c(1,0,0,8,33,54,2,0,0),row.names=aa), # y=data.frame(c(1,0,0,8,3,54,2,0,0),row.names=aa), # gene.symbol="ab") #ab <- c("ad_chr1_167931912_167932008","ab_chr1_167936912_167937008","ab_chr1_167951800_167951861","ab_chr1_167966629_167966721", # "ab_chr1_167974830_167975015","ab_chr1_167987495_167987608","ab_chr1_167991204_167991339","ad_chr1_167993226_167993440", # "ab_chr1_168002482_168002575") #plot.exon.levels(x=data.frame(c(1,0,0,8,33,54,2,0,0),c(2,2,0,5,30,60,2,2,0),row.names=ab), # y=data.frame(c(1,0,0,8,3,54,2,0,0),c(2,2,0,5,10,6,2,2,0),row.names=ab), # sig=4, gene.symbol="ab")