snpPlot <- function(x,range=NULL, yLabels=NULL, zlim=NULL,title=NULL,colRamp=12){ require(colorRamps) if (missing(x)) { stop("Data matrix x missing. Stopped.") } if (!is.matrix(x)) { x <- as.matrix(x) } if (!is.numeric(x)) { stop("Data matrix x must be numeric. Stopped.") } min <- min(x) max <- max(x) if( !is.null(zlim) ){ min <- zlim[1] max <- zlim[2] } if( !is.null(yLabels) ){ yLabels <- c(yLabels) } else { yLabels <- c(1:nrow(x)) } if( !is.null(range) ){ xLabels <- c(range) } else { xLabels <- c(1:ncol(x)) } if( is.null(title) ){ title <- c() } # Red and green range from 0 to 1 while Blue ranges from 1 to 0 ColorRamp <- switch(colRamp, rgb( seq(0,1,length=256),seq(0,1,length=256),seq(1,0,length=256)), rgb( seq(0,1,length=256),seq(0.5,1,length=256),seq(1,0,length=256)), rgb( seq(0,1,length=256),seq(0,0.5,length=256),seq(1,0,length=256)), blue2green(256), green2red(256), blue2yellow(256), cyan2yellow(256), magenta2green(256), matlab.like(256), matlab.like2(256), blue2green2red(256), rgb(t(rep(c(255,255,255),256) - col2rgb(matlab.like2(256))),maxColorValue = 255), rgb(t(rep(c(255,255,255),256) - col2rgb(blue2green2red(256))),maxColorValue = 255), rgb(t(abs(rep(c(0,255,0),256) - col2rgb(matlab.like2(256)))),maxColorValue = 255), rgb(t(abs(rep(c(0,55,0),256) - col2rgb(blue2green2red(256)))),maxColorValue = 255), rgb(t(abs(rep(c(255,0,255),256) - col2rgb(matlab.like2(256)))),maxColorValue = 255), rgb(t(abs(rep(c(255,0,255),256) - col2rgb(blue2green2red(256)))),maxColorValue = 255) ) # ColorRamp <- rgb( seq(0,1,length=256), # Red # seq(0,1,length=256), # Green # seq(1,0,length=256)) # Blue ColorLevels <- seq(min, max, length=length(ColorRamp)) # Reverse Y axis reverse <- nrow(x) : 1 yLabels <- yLabels[reverse] x <- x[reverse,] # Data Map # par(mar = c(3,5,2.5,2)) # par(mar = c(3,5,2.5,1)) # Original: par(mar=c(5.1,4.1,4.1,2.1)) marO <- par("mar") par(mar = c(2.5,6,3,1)) image(1:ncol(x), 1:nrow(x), t(x), col=ColorRamp, xlab="", ylab="", axes=FALSE, zlim=c(min,max)) if( !is.null(title) ){ title(main=title) } xt <- c(1,axTicks(1)) lxt <- length(xt) if (abs(xt[lxt]-ncol(x))>5) { xt <- c(xt,ncol(x)) } axis(BELOW<-1, at=xt, labels=prettyNum(xLabels[xt], big.mark = ","), cex.axis=0.7) axis(LEFT <-2, at=1:length(yLabels), labels=yLabels, las= HORIZONTAL<-1,cex.axis=0.7) par(mar = marO) layout(1) } plotRangeL <- function(Lout,ibb,physPos=NULL,colRamp=12,val=c(0.0,2.0,1.0),chrom=1,count=0,labelsNA=NULL,prange=NULL,labelsNA1 = c("model L") ) { if (missing(Lout)) { stop("Sample SNVs 'Lout' are missing. Stopped.") } if (missing(ibb)) { stop("tagSNVs 'ibb' are missing. Stopped.") } n <- dim(Lout)[2] n1 <- length(ibb) doMis <- FALSE if (n1>4) { n1 <- n1-1 doMis <- TRUE } if ((!is.null(prange))&&(!is.null(physPos))) { if (length(prange)==2) { ma <- length(ibb[[1]]) a1 <- which(physPos>=prange[1]) a2 <- which(physPos<=prange[2]) aa <- intersect(a1,a2) if (length(aa)>1) { ibb <- ibb[[1]][aa] physPos <- physPos[aa] } } } if (is.null(physPos)) { physPos <- ibb*100 } ma <- length(ibb[[1]]) range <- ibb[[1]][1]:ibb[[1]][ma] lr <- ibb[[1]][ma]-ibb[[1]][1]+1 physRanget <- physPos[1]:physPos[ma] lt <- length(physRanget) ltr <- lt%/%(lr-1) ltt <- c(1,(1:(lr-2))*ltr,lt) physRange <- physRanget[ltt] physLength <- round(lt/1000) physPosM <- (physPos[ma] + physPos[1])%/%2 mat <- matrix(val[1],(n+n1+1),lr) mat[(n+1),] <- val[2] for (iB in 1:n1) { ibb_pos <- match(ibb[[iB]],range) mat[(n+1+iB),ibb_pos] <- val[3] } if (doMis) { ibb_pos <- match(ibb[[n1+1]],range) mat[(n+1+3),ibb_pos] <- 0.35*val[3] } ibb_pos <- match(ibb[[1]],range) for (j in 1:n) { lla <- which(Lout[range,j]>=0.1) llb <- setdiff(1:lr,lla) llc <- intersect(ibb_pos,lla) lla <- setdiff(lla,llc) mat[j,lla] <- val[2] mat[j,llb] <- val[1] mat[j,llc] <- val[3] } if (is.null(labelsNA)) { labelsNA <- as.character(1:n) } if (length(labelsNA1)!=n1) { if (n1>1) { ade <- rep("",(n1-1)) labelsNA1 = c("model L",ade) } else { labelsNA1 = c("model L") } } labelsNA <- c(labelsNA,"",labelsNA1) if (count<=0) { snpPlot(mat,range=physRange,yLabels=labelsNA, title=paste("chr: ",chrom," || pos: ",prettyNum(physPosM, big.mark = ",")," || length: ",physLength,"kbp || #SNPs: ",ma," || #Samples: ",n,sep=""),colRamp=colRamp) } else { snpPlot(mat,range=physRange,yLabels=labelsNA, title=paste("count: ",count," || chr: ",chrom," || pos: ",prettyNum(physPosM, big.mark = ",")," || length: ",physLength,"kbp || #SNPs: ",ma," || #Samples: ",n,sep=""),colRamp=colRamp) } } plotAcientGenomes <- function(IBD,filename,ibdC) { samp <- IBD[[ibdC]]$noSamples ibb <- IBD[[ibdC]]$noSnps ibb <- as.integer(sort.int(as.integer(unique(ibb)))) snpPositions <- IBD[[ibdC]]$snpPositions labels_ALL <- IBD[[ibdC]]$labelSamples Lout <- readSamplesSpfabia(X=filename,samples=samp,lowerB=0,upperB=1000.0) if (!exists("annotationAChr1")) { load(file="/seppdata/sepp/linkage/release/annotation/annotationAChr1.Rda") } if (!exists("posA")) { posA <- annotationAChr1[,2] neaA <- annotationAChr1[,24] refA <- annotationAChr1[,4] altA <- annotationAChr1[,5] } mA <- match(snpPositions,posA,nomatch=0) gA <- mA[which(mA>0)] l2 <- length(gA) neaAG <- neaA[gA] altAG <- altA[gA] refAG <- refA[gA] ibbA <- c() for (it in 1:l2) { if (toupper(neaAG[it]) == toupper(altAG[it])) { ibbA <- c(ibbA,ibb[it]) } } if (!exists("neanderHg19Chr1")) { load(file="/seppdata/sepp/linkage/release/annotation/neanderHg19Chr1.Rda") } if (!exists("posN")) { posN <- neanderHg19Chr1$V2 neaN <- neanderHg19Chr1$ad1 refN <- neanderHg19Chr1$V5 altN <- neanderHg19Chr1$V6 } mN <- match(snpPositions,posN,nomatch=0) gN <- mN[which(mN>0)] l2 <- length(gN) neaNG <- neaN[gN] altNG <- altN[gN] refNG <- refN[gN] ibbN <- c() ibbNa <- c() for (it in 1:l2) { AAs <- strsplit(as.character(neaNG[it]),",")[[1]] if (length(AAs)>1) { if (AAs[1] == altNG[it]) { ibbN <- c(ibbN,ibb[it]) } } else { if (AAs[1]=="-") { ibbNa <- c(ibbNa,ibb[it]) } } } if (!exists("denisovaHg19Chr1")) { load(file="/seppdata/sepp/linkage/release/annotation/denisovaHg19Chr1.Rda") } if (!exists("posD")) { posD <- denisovaHg19Chr1$V2 neaD <- denisovaHg19Chr1$ad1 refD <- denisovaHg19Chr1$V5 altD <- denisovaHg19Chr1$V6 } mD <- match(snpPositions,posD,nomatch=0) gD <- mD[which(mD>0)] l2 <- length(gD) neaDG <- neaD[gD] altDG <- altD[gD] refDG <- refD[gD] ibbD <- c() for (it in 1:l2) { AAs <- strsplit(as.character(neaDG[it]),",")[[1]] if (length(AAs)>1) { if (AAs[1] == altDG[it]) { ibbD <- c(ibbD,ibb[it]) } } } ibbL <- list(ibb,ibbA,ibbN,ibbD,ibbNa) plotRangeL(Lout=Lout,ibb=ibbL,physPos=snpPositions,colRamp=12,val=c(0.0,2.0,1.0),chrom=1,count=0,labelsNA=labels_ALL,labelsNA1=c("model L","Ancestor","Neandertal","Denisova")) }