# TODO: Add comment # # Author: E.Korsching Feb 8, 2012 ############################################################################### plot.genome.inline <- function(x, event.marker=1, bin.num=10, y, pos.col.name="Physical.Position", chr.col.name="Chromosome", xlab="", ylab="", chr.sel=0, pos1.sel=0, pos2.sel=0, marker1=NULL, scale.individually=F, scale.ext.f=0.05, sig.offset=0, loss.th=-0.08, gain.th=0.08, ampli.th=0.45, pl.density=F, plot.type="p", pt.s.size=0.4, colorI=c("green","red","blue"), cex=0.5){ # plot all chr in a line from 1 to n # needs data.frame (x) with rownames and one or more cols (experiments) # and annotation (y) with chromosomal position col and chromosome number col (no characters!) # chr.sel: one or more numbers in vec, pos1.sel/pos2.sel: only one starting and end point / only possible if one chr selected # marker1: put a vertical line (or more) in the graph, NULL: no line # data.frame: col1: physical position, col2: chromosome, col3: description (char, empty "" str are valid) # loss.th, gain.th, ampli.th: threshold values for horizontal lines for loss, gain, amplification # scale.individually: T : every experiment by it's own range # scale.ext.f: extend scale range by a certain factor # sig.offset: to make signals positions which are '0' visible # pl.density: plot density curve on top with signal values >= event.marker value # plot.type: p:points, b:bars -- pt.s.size: plot symbol size [cex] # colorI: Color definition (3: green:loss, red:gain, blue:amplification), other lines per default: black # cex: Character size of the chromosome numbers #ini ny <- ncol(x) exp.names <- names(x) #check if(length(chr.sel)>1 & (pos1.sel!=0 | pos2.sel!=0)){cat("\n error - pos1.sel,pos2.sel info only applicable by one chr ");break} if(sum(is.na(match(chr.sel, unique(y[,chr.col.name]))))>=1){cat("\n error - selected chr not in annotation set ");break} #order anno concerning chr and pos y1 <- y[order(y[,chr.col.name],y[,pos.col.name]),,drop=F] #discard all lines in anno where in chr or pos is a NA y2 <- y1[!is.na(y1[,pos.col.name]) & !is.na(y1[,chr.col.name]),,drop=F] #ini2 if(sum(match(chr.sel,0), na.rm=T)==1){ chr.set <- unique(y2[,chr.col.name])[order(unique(y2[,chr.col.name]))] #always ascending ! }else{ chr.set <- chr.sel[order(chr.sel)] } #check max of anno with max of centromere.human for all chr bb <- 0 for(i in chr.set){ if(max(y2[y2[,chr.col.name]==i,pos.col.name]) > centromere.human[i,"bottom"]){ cat("\n chr ",i," length error in annotation - centromere.human : ",max(y2[y2[,chr.col.name]==i,pos.col.name])," - ",centromere.human[i,"bottom"]) bb <- 1 } } if(bb==1){cat("\n\n Fix length error in centromere.human - backup old centromere.human \n",); break} #order x, resize x and anno if(sum(match(chr.set,0), na.rm=T)==1){ # all chr x1 <- x[row.names(y2),,drop=F] # only order, take complete x }else{ # chr(s) selected if(pos1.sel==0){ # more than one chr bb <- matrix(0,nrow(y2),length(chr.set)) for(i in 1:length(chr.set)){ bb[,i] <- y2[,chr.col.name]==chr.set[i] } bb <- rowSums(bb) if(max(bb)>1){cat("\n error - overlapping selections - check the selection process ");break} y2 <- y2[as.logical(bb),,drop=F] # resize }else{ # one chr y2 <- y2[y2[,chr.col.name]==chr.set & y2[,pos.col.name]>=pos1.sel & y2[,pos.col.name]<=pos2.sel,,drop=F] # resize } x1 <- x[row.names(y2),,drop=F] # resize and order } #add running sum per (next) chr - add in new col to anno data set bb <- 0 bc <- ncol(y2)+1 for(i in chr.set){ y2[y2[,chr.col.name]==i,bc] <- y2[y2[,chr.col.name]==i,pos.col.name] + bb bb <- bb + centromere.human[i,"bottom"] } #ini3 nr <- nrow(x1) if(!scale.individually){ # max/min for all experiments identical if(max(x1)<1){ max.x1 <- 1.1 }else{ max.x1 <- max(x1) } if(min(x1)>-1){ min.x1 <- -1.1 }else{ min.x1 <- min(x1) } } # cat("\n max min ext ",max.x1," ",min.x1," ",min.x1-(abs(min.x1)*scale.ext.f)) #ini4 - plot data f1 <- function(w,sig.offset,col){ #segments if(w[2]>gain.th){tmp <- col[2]} # choose one col from colorI if(w[2]<=gain.th & w[2]>=loss.th){tmp <- "black"} if(w[2]gain.th){tmp <- col[2]} # choose one col from colorI if(w[2]<=gain.th & w[2]>=loss.th){tmp <- "black"} if(w[2]0.45){ segments(0, 0.45, sum(centromere.human[chr.set,"bottom"]), ampli.th, lty=1, lwd=1, col=colorI[3]) # plot amplification line } bb <- 0 for(j in chr.set){ #plot ticks aa <- bb + centromere.human[j,"bottom"] segments(bb+centromere.human[j,"centromere"], min.x1, bb+centromere.human[j,"centromere"], 0, lty=3) # plot chr centromer tick if(j==chr.set[1]){segments(sig.offset, min.x1, sig.offset, 0, lty=2)} # plot start chr axis tick segments(aa, min.x1, aa, 0, lty=2) # plot end chr axis tick # plot chromosome marker1 + description if(!is.null(marker1)){ if(sum(marker1[,2]==j)>0){ marker.tmp <- marker1[marker1[,2]==j, 1] segments(bb+marker.tmp, min.x1/2, bb+marker.tmp, max.x1/2, lty="431313", col="black") marker.tmp2 <- marker1[marker1[,2]==j, 3] text(x=bb+marker.tmp,y=(max.x1/2)+(max.x1/12),marker.tmp2, cex=cex, adj=0.5) } } text(x=aa-((aa-bb)/2),y=min.x1-(abs(min.x1)*scale.ext.f),as.character(j), cex=cex, adj=0.5) # plot chr ID bb <- aa } }else{ # part of a chromosome # plot area plot(x=0,y=0, type="n", xlim=c(pos1.sel,pos2.sel), ylim=c(min.x1-(abs(min.x1)*scale.ext.f),max.x1+(max.x1*scale.ext.f)), xlab=xlab, ylab=ylab, axes=F) text(x=pos2.sel-((pos2.sel-pos1.sel)/2),y=min.x1-(abs(min.x1)*scale.ext.f),as.character(chr.set), cex=cex) # plot chr ID } axis(2) title(exp.names[i], cex=0.7) tmp1 <- cbind(y2[,bc],x1[,i]) if(plot.type=="b"){ apply(tmp1,1,f1,sig.offset=sig.offset,col=colorI) } # plot data [apply-segments,points - factor 17 faster] if(plot.type=="p"){ apply(tmp1,1,f2,sig.offset=sig.offset,col=colorI,cex=pt.s.size) } if(pl.density){ bin.pos.data <- genome.event.density(scale.len=c(0,sum(centromere.human[chr.set,"bottom"])),event=y2[x1[,i]>=event.marker,bc],bin.num=bin.num) #count certain signals plot.genome.event.density(x=bin.pos.data, scaling=T, range.y=c(par()$yaxp[1],par()$yaxp[2]), range.fix=F, new=F) } } } #fix chromosome length error: #centromere.human.old1 <- centromere.human #and set manually to appropriate value #plot.genome.inline(x=d.meth.g.n.1.diff.inv.dcis.1, y=h.methylation.anno.c, pos.col.name="loc", chr.col.name="GeneChromosome", # xlab="chr 1 - 22, X,Y", ylab="signal", chr.sel=0, pos1.sel=0, pos2.sel=0, scale.individually=F, scale.ext.f=0.05, sig.offset=0.05, pl.density=F) #plot.genome.inline(x=d.meth.g.n.1.diff.inv.dcis.1, y=h.methylation.anno.c, pos.col.name="loc", chr.col.name="GeneChromosome", # xlab="chr 1 - 22, X,Y", ylab="signal", chr.sel=c(1,4), pos1.sel=0, pos2.sel=0, scale.individually=F, scale.ext.f=0.05, sig.offset=0.05, pl.density=F) #plot.genome.inline(x=d.meth.g.n.1.diff.inv.dcis.1, y=h.methylation.anno.c, pos.col.name="loc", chr.col.name="GeneChromosome", # xlab="chr 1 - 22, X,Y", ylab="signal", chr.sel=4, pos1.sel=70000000, pos2.sel=140000000, scale.individually=F, scale.ext.f=0.05, sig.offset=0.05, pl.density=F) #plot.genome.inline(x=d.meth.g.n.1.diff.inv.dcis.1, event.marker=1, bin.num=500, y=h.methylation.anno.c, pos.col.name="loc", chr.col.name="GeneChromosome", # xlab="chr 1 - 22, X,Y", ylab="signal", chr.sel=4, pos1.sel=0, pos2.sel=0, scale.individually=F, scale.ext.f=0.05, sig.offset=0.05, pl.density=T) genome.event.density <- function(scale.len,event,bin.num){ # density of events on a ratio scale per interval (bin) # scale.len: lower+upper value, bin.num: number of bins # event: vec of numbers - all or part is element of scale.len #create bin scalar interval <- (scale.len[2]-scale.len[1])/bin.num bin.table <- matrix(0,bin.num,1) for(i in 1:bin.num){ bin.table[i,1] <- scale.len[1] + interval * i } #calculate bins f1 <- function(bin.upper,event,interval,max.point){ bin.lower <- bin.upper - interval if(max.point==bin.upper){ erg <- event[event>=bin.lower & event<=bin.upper] #bin on edge might be slightly larger !OBS! }else{ erg <- event[event>=bin.lower & event