segment.plot <-function( file="/data/www/POOLING/ARABIDOPSIS/FOUNDER/combined.segments.txt", pdffile="segments.pdf", genome=FALSE) { segs=read.table(file,header=TRUE) if (!is.null(pdffile)) pdf(pdffile) par(mfrow=c(5,1)) par(mar=c(3,4,1,2)) if ( genome ) { by (segs, segs$magic, seg.plot.genome ) } else { by (segs, segs$magic, seg.plot ) } if (!is.null(pdffile)) dev.off() } seg.plot <- function( seg ) { print(seg) by( seg, seg$chr, seg.plot.chr ) } seg.plot.chr <- function ( seg ) { print(seg) labels = sub( ".GBS.MERGE", "", seg$acc) y=-log10(seg$error.bp+1.0e-10) plot( c(0), c(0), xlim=c(0,max(seg$to.bp/1.0e6+1)), ylim=c(0,8), xlab="Mb", ylab="log error.bp", t="n") abline( v=seg$from.bp/1.0e6, col="grey") segments( seg$from.bp/1.0e6, y, seg$to.bp/1.0e6, y ) diploid = grep( "^\\*", labels, perl=TRUE, value=FALSE) print("diploid") print(diploid) if ( length(diploid)>0) { xleft = seg$from.bp/1.0e6 xright = seg$to.bp/1.0e6 ybot = rep(0,length(seg$from.bp)) ytop = y print(xleft[diploid]) rect(xleft[diploid], ybot[diploid], xright[diploid], ytop[diploid], col=rep("orange",len=length(diploid)), border=rep(NA,length(diploid))) } yy=c(1,2,3) if ( length(yy) > nrow(seg) ) length(yy) = nrow(seg) text( (seg$from.bp/1.0e6 +seg$to.bp/1.0e6)/2, yy, labels, cex=0.5) print((seg$from.bp/1.0e6 +seg$to.bp/1.0e6)/2) text( 4, 5, paste(seg$magic[1], "chr", seg$chr[1]), col="red") } seg.plot.genome <- function ( seg ) { len = aggregate ( seg$to.bp, by=list(chr=seg$chr), max, simplify = TRUE ) len = c( 0, len$x ) len = cumsum(len) off = len[seg$chr] seg$to.bp = seg$to.bp + off seg$from.bp = seg$from.bp + off labels = sub( ".GBS.MERGE", "", seg$acc) y=-log10(seg$error.bp+1.0e-10) plot( c(0), c(0), xlim=c(0,len[5]/1e6), ylim=c(0,8), xlab="Mb", ylab="log error.bp", t="n") abline( v=seg$from.bp/1.0e6, col="grey") abline( v=len/1.0e6, col="red") segments( seg$from.bp/1.0e6, y, seg$to.bp/1.0e6, y ) diploid = grep( "^\\*", labels, perl=TRUE, value=FALSE) print("diploid") print(diploid) if ( sum(diploid)>0) { xleft = seg$from.bp/1.0e6 xright = seg$to.bp/1.0e6 ybot = 0 ytop = y print(xleft) rect(xleft[diploid], ybot[diploid], xright[diploid], ytop[diploid], col="orange", border=NA) } yy=c(1,2,3) if ( length(yy) > nrow(seg) ) length(yy) = nrow(seg) text( (seg$from.bp/1.0e6 +seg$to.bp/1.0e6)/2, yy, labels, cex=0.5) text( 4, 5, seg$magic[1], col="red") }