# Date : 06.07.2020 # Author: Eberhard Korsching ############################################################################### # Title: scale up exons (and scale down introns) # for illustration purpose # # proportionaly means: # if the gene reading frame sequence is 100% and all the exons are taking 20% / introns 80% # then we would tune by a factor # that all exons (proportionaly) have e.g. 60% and introns 40% # # note 1 : the resulting graph is no longer a true image of the genomic structure/scale # # note 2 : the graph legend or the graph itself needs to state this fact for the reader # # note 3 : raw data from EBI Ensemble web site e.g. # https://www.ensembl.org/Homo_sapiens/Transcript/Summary?db=core;g=ENSG00000106086;r=7:30028685-30084663;t=ENST00000440706 # left panel 'summary', select one transcript, left panel 'exon', wait until loaded, download full table (on the right) # delete last column (CSV format) # No. Exon / Intron Start End Start Phase End Phase Length # 5' upstream sequence # 1 ENSE00001884360 29268268 29268349 - 1 82 # Intron 1-2 29268350 29272215 3866 gene.exon.illustration <- function(file, outFile, ze=1, lwd.s=0.2, lwd.m=0.2, genelabel="gene", othertext="", posMarker=NULL, pmtype="", lineMarker=NULL, col=c("#CC6600","#FF9900","#0099FF","#33CC00","black"), cexD=0.6, cexA=1){ # schematic graph of a gene structure # input: file: path+file name of EBI ensemble gene export data set (see presentation) # ze: emphasize factor for exons # lwd.s: line width scheme, lwd.m: line width marker line # genelabel: gene name or else as title # othertext: discription beside from/to positions # posMarker: above exons: one or more positions as vertical lines in the image, pmtype: marker type: "": none, "a": arrow, "s": line # lineMarker: below exons: one or more lines: vector with 2,4,6, ..., 12 entries (max) # col: colors: [1] exon, [2] intron, [3] marker, [4] line, [5] title line, exon numbers, footer line # cexD: font size for subtitle, cexA: font size for title and exon numbers # output: (path and) file name # constants ylim <- c(0,1.5) # c(0.15,0.85) y1 <- 0.4 y2 <- 0.6 y3 <- 0.46 # = y5 ! y4 <- 0.54 n1 <- c(0.75,1.05) # exon numbers n2 <- c(0.77,0.85) # marker lines y start, end n3 <- c(0.3,0.025) # range lines y start, step is.even <- function(x){x%%2 == 0} # ini require(gtools) is.odd <- function(x){ x %% 2 != 0 } # import # 1 header line x <- read.table(file=file, header=T, sep=",", dec=".", quote="", stringsAsFactors=F) # dim xr <- nrow(x) # delete unwanted lines x <- x[-c(1,xr),,drop=F] # strip "," in character strings of numbers and convert to numeric type x[,3] <- as.numeric(gsub(pattern=",", replacement="", x=x[,3], fixed=T)) x[,4] <- as.numeric(gsub(pattern=",", replacement="", x=x[,4], fixed=T)) x[,7] <- as.numeric(gsub(pattern=",", replacement="", x=x[,7], fixed=T)) # get/customize core data xr <- nrow(x) # update if(x[1,3]>x[xr,4]){ # flip matrix vertically if on reverse strand flip <- T x.new <- x k <- 0 for(i in 1:xr){ x.new[i,] <- x[(xr-k),] k <- k+1 } xnam <- names(x.new) x <- x.new[,c(1,2,4,3,5,6,7)] # flip additionally column 3 and 4 start end names(x) <- xnam # correct start end names }else{ flip <- F } FromTo.seq <- c(x[1,3],x[xr,4]) # start end cat("\n from-to: ",FromTo.seq) len.std <- x[ ,7] # get length vector len.seq <- sum(len.std) # length seq cat("\n length seq: ",len.seq) # --Marker (1) if(!is.null(posMarker)){ posMarker.relative <- posMarker-FromTo.seq[1]+1 # make absolute position a relative position for gene len.sections.std <- create.sections(len.std) # sub function to create section marker vector } if(!is.null(lineMarker)){ nl <- length(lineMarker) if(nl>12 | is.odd(nl)){ stop("lineMarker: wrong number of items") } lineMarker.relative <- lineMarker-FromTo.seq[1]+1 # make absolute position a relative position len.sections.std <- create.sections(len.std) # sub function to create section marker vector } # proportions (in old length - return new length values) comp.new <- proportion(x=len.std, ze=ze) len.new <- comp.new[[1]] i.ze <- comp.new[[2]] # --Marker (2) if(!is.null(posMarker)){ posMarker.rel.new <- translate.Marker(pMr=posMarker.relative, lss=len.sections.std, lst=len.std, lne=len.new, ze=ze, i.ze=i.ze) } if(!is.null(lineMarker)){ lineMarker.rel.new <- translate.Marker(pMr=lineMarker.relative, lss=len.sections.std, lst=len.std, lne=len.new, ze=ze, i.ze=i.ze) } # create coordinates d1.row <- (xr+1)/2 d1 <- matrix(0, d1.row, 5) # y1,y2: exon box height, y3=y5, y4: intron triangle (all constants) dimnames(d1)[[2]] <- c("x1","x2","x3","x4","x5") # x1: exon start, x2: exon end, x3: intron start, x4: intron mid point, x5: intron end (all variable) tmp <- 0 count <- 1 for(i in seq(1,xr,2)){ if(i!=xr){ d1[count,c(1,2)] <- c(tmp, (tmp+len.new[i])) tmp <- tmp + len.new[i] d1[count,c(3,5)] <- c(tmp, (tmp+len.new[i+1])) tmp <- tmp + len.new[i+1] count <- count + 1 }else{ d1[count,c(1,2)] <- c(tmp, (tmp+len.new[i])) tmp <- 0 count <- 0 } } # dependend value of d1 for(i in 1:d1.row){ if(i!=d1.row){ d1[i,4] <- d1[i,3] + (d1[i,5]-d1[i,3])/2 } } # adjust gene direction if(flip){ tx <- paste("from",FromTo.seq[1],"to",FromTo.seq[2],othertext," <- reverse strand ",sep=" ") }else{ tx <- paste("from",FromTo.seq[1],"to",FromTo.seq[2],othertext," forward strand -> ",sep=" ") } # graph png(filename=outFile, width=1200, height=600, units="px", bg="white", res=600) # pointsize=12, p.mar <- par()$mar par(mar=c(2,0.1,2,0.1)) plot(0,0, type="n", xlab="", ylab="", xlim=c(-(sum(len.new)/50),sum(len.new)), ylim=ylim, axes=F, col=col[5]) # footer color # draw boxes and connectors - because of the order of drawings: two loops for(i in 1:(d1.row-1)){ segments( x0=d1[i,c(3,4)], y0=c(y3,y4), x1=d1[i,c(4,5)], y1=c(y4,y3), lwd=lwd.s, col=col[2] ) # intron lines } for(i in 1:d1.row){ polygon( x=d1[i,c(1,2,2,1)], y=c(y1,y1,y2,y2), border=NA, col=col[1] ) # exons boxes } # annotate the graph en1 <- d1[,1]+(d1[,2]-d1[,1])/2 en1.l <- length(en1) en1.even <- is.even(en1.l) if(en1.even){ en1.f <- en1.l/2 v1 <- rep(n1,en1.f) }else{ en1.f <- ceiling(en1.l/2) v1 <- rep(n1,en1.f) v1 <- v1[-(en1.l+1)] } if(flip){ text(x=en1, y=v1, col=col[6], labels=as.character(d1.row:1), cex=cexA, srt=30) # exon numbers }else{ text(x=en1, y=v1, col=col[6], labels=as.character(1:d1.row), cex=cexA, srt=30) # exon numbers } mtext(text=genelabel, side=3, col=col[5], cex=cexA) # gene name if(ze>1){ # exon number & note mtext(text=paste(d1.row," exons (Note: exons enlarged)",sep=""), side=1, line=0.5, col=col[5], cex=cexD) }else{ mtext(text=paste(d1.row," exons",sep=""), side=1, line=0.5, col=col[5], cex=cexD) } mtext(text=tx, side=1, line=0, col=col[5], cex=cexD) # any marker to set? note: scaling broadens the size of one base in the case of exons if(!is.null(posMarker)){ if(ze==1){ ze <- 0 } # no adjustment if(pmtype=="a"){ # -ze/2: a scaling adjustment, visible in the tests example, nearly not visible >100 bases pml <- length(posMarker.rel.new) points(x=posMarker.rel.new-ze/2, y=rep(n2[2],pml), pch=25, col=col[5], bg=col[3], lwd=lwd.m, cex=cexD) }else if(pmtype=="s"){ segments(x0=posMarker.rel.new-ze/2, y0=n2[2],y1=n2[1], col=col[3], lwd=lwd.m) } } if(!is.null(lineMarker)){ # if(ze==1){ ze <- 0 } j <- 1 k <- n3[1] for(i in 1:(nl/2)){ # segments(x0=lineMarker.rel.new[j]-ze/2, x1=lineMarker.rel.new[(j+1)]-ze/2, y0=k, col=col[4], lwd=lwd) segments(x0=lineMarker.rel.new[j], x1=lineMarker.rel.new[(j+1)], y0=k, col=col[4], lwd=lwd.m) k <- k-n3[2] j <- j+2 } } # cleanup par(mar=p.mar) dev.off() # if other functions want to add something: return the positions return(d1) } # tests # #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1, lwd.s=1, genelabel="gene", othertext="", posMarker=NULL, pmtype="", lineMarker=NULL, # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1.2, lwd.s=1, genelabel="gene", othertext="", posMarker=NULL, pmtype="", lineMarker=NULL, # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) # #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1, lwd.s=1, genelabel="gene", othertext="", posMarker=NULL, pmtype="", lineMarker=NULL, # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1, lwd.s=1, genelabel="gene", othertext="", posMarker=c(1), pmtype="s", lineMarker=c(2,3), # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) # #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1, lwd.s=1, genelabel="gene", othertext="", posMarker=c(1), pmtype="s", lineMarker=c(2,3, 5,8, 2,9), # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) #a <- gene.exon.illustration(file="/home/korschi/on1/data/eclipseR/0functions/chromosome/gene.exon.illustration.testdata.csv", outFile="test1.png", # ze=1.2, lwd.s=1, genelabel="gene", othertext="", posMarker=c(1), pmtype="s", lineMarker=c(2,3, 5,8, 2,9), # col=c("blue","red","darkgreen","blue","orange3","black"), cexD=0.6, cexA=1) ## examples #a <- gene.exon.illustration(file="projectwork/CDH1.csv", outFile="test1.png", ze=1, genelabel="CDH1", othertext="+ strand") #a <- gene.exon.illustration(file="projectwork/CDH1.csv", outFile="test2.png", ze=8, genelabel="CDH1", othertext="+ strand") # #a <- gene.exon.illustration(file="projectwork/COL8A1.csv", outFile="test3.png", ze=1, genelabel="COL8A1", othertext="+ strand") #a <- gene.exon.illustration(file="projectwork/COL8A1.csv", outFile="test4.png", ze=8, genelabel="COL8A1", othertext="+ strand") # #save.image() # sub functions create.sections <- function(len.std){ # create a category vector len <- length(len.std) cat.std <- NULL for(i in 1:len){ cat.std <- c(cat.std,rep(i,times=len.std[i])) } return(cat.std) } translate.Marker <- function(pMr, lss, lst, lne, ze, i.ze){ # calulate new marker position for the scaled graph # pMr: vector relative marker positions # lss: vector sections standard # lst: vector length elements standard # lne: vector length elements new len.pMr <- length(pMr) # maybe more than one position pMrrf <- vector("numeric",len.pMr) # final marker positions # cat("\n pMr :",pMr) sec.M <- lss[pMr] # pMr associated section pointer cat("\n std positions in sections :",sec.M) pMrr <- vector("numeric",len.pMr) for(i in 1:len.pMr){ # relative per section if(sec.M[i]>1){ #cat("\n pMr ",pMr[i]," sum ",sum( lst[1:(sec.M[i]-1)] )) pMrr[i] <- pMr[i] - sum( lst[1:(sec.M[i]-1)] ) }else{ pMrr[i] <- pMr[i] } } # cat("\n pMrr ",pMrr) # cat("\n lst ",lst) # cat("\n lne ",lne) for(i in 1:len.pMr){ # new relative section position if(odd(sec.M[i])){ rpn <- pMrr[i] * ze # the arrow head has now the width of ze }else{ rpn <- pMrr[i] * i.ze # the arrow head has now the width of i.ze } if(sec.M[i]>1){ pMrrf[i] <- rpn + sum( lne[1:(sec.M[i]-1)] ) # new relative position in gene range }else{ pMrrf[i] <- rpn } } # cat("\n pMrrf ",pMrrf) return(pMrrf) } proportion <- function(x, ze){ # x: length column , ze: ze factor # calculate new proportions based on ze and original length of exon / introns # length of the parameter/argument x len <- length(x) # vector structure at this point: exon, intron, exon, intron, ... , exon sum.exon <- sum(x[seq(1,len,2)]) sum.intron <- sum(x[seq(2,len,2)]) sum.ei <- sum.exon + sum.intron p.exon <- sum.exon / sum.ei p.intron <- sum.intron / sum.ei # cat("\n input : exon :",sum.exon," intron :",sum.intron," sum : ",sum.ei) # for each exon / intron value n.ei <- vector("numeric", len) # new exon length depending on ze n.ei[seq(1,len,2)] <- x[seq(1,len,2)] * ze # new sum of exons sum.exon.new <- sum(n.ei[seq(1,len,2)]) if(sum.exon.new>=(0.98*sum.ei)){ cat("\n"); stop(paste("stop - zoom factor too high, ze: ",ze,sep="")) } # new sum introns sum.intron.new <- sum.ei - sum.exon.new # new intron scaling factor i.ze <- sum.intron.new / sum.intron # new intron lenth n.ei[seq(2,len,2)] <- x[seq(2,len,2)] * i.ze sum.intron.new.control <- sum(n.ei[seq(2,len,2)]) # cat("\n output : exon :",sum.exon.new," intron :",sum.intron.new," sum : ",sum.exon.new+sum.intron.new) # cat("\n intron control :",sum.intron.new.control,"\n") cat("\n i.ze :",i.ze) return(list(n.ei,i.ze)) }