====== Population average pairwise difference ====== * file: matrix_PairwiseDifferences.r * document page: 33 \\ **code:** #----read data------------------------------------------------------------------ Data <- read.table("D:/Heidi/Master/R_Daten/SummaryStatistics/pairwise_differences.txt", skip=10, row.names=1) DataMatrix <- as.matrix(Data) #----UnderMatrix---- n <- 2 x <- 1 UnderMatrix <- matrix(NA, ncol=ncol(DataMatrix), nrow=1) while(n<=nrow(DataMatrix)){ nextrow <- DataMatrix[n,1:x] nextrow <- cbind(t(nextrow), matrix(NA, ncol=ncol(DataMatrix)-x, nrow=1)) UnderMatrix <- rbind(UnderMatrix, nextrow) n <- n+1 x <- x+1 } #----UpperMatrix---- n <- 2 x <- 1 UpperMatrix <- matrix(NA, ncol=x, nrow=1) UpperMatrix <- cbind(UpperMatrix, t(DataMatrix[x,n:ncol(DataMatrix)])) n <- n+1 x <- x+1 while(n<=nrow(DataMatrix)){ nextrow <- matrix(NA, ncol=x, nrow=1) nextrow <- cbind(nextrow, t(DataMatrix[x,n:ncol(DataMatrix)])) UpperMatrix <- rbind(UpperMatrix, nextrow) n <- n+1 x <- x+1 } UpperMatrix <- rbind(UpperMatrix, matrix(NA, ncol=ncol(DataMatrix), nrow=1)) #----DiagonalMatrix---- n <- 1 DiagonalMatrix <- DataMatrix[n,n] DiagonalMatrix <- cbind(DiagonalMatrix, matrix(NA, ncol=ncol(DataMatrix)-n, nrow=1)) n <- n+1 x <- 1 while(n<=nrow(DataMatrix)){ nextrow <- matrix(NA, ncol=x, nrow=1) nextrow <- cbind(nextrow, t(DataMatrix[n,n])) nextrow <- cbind(nextrow, matrix(NA, ncol=ncol(DataMatrix)-n, nrow=1)) DiagonalMatrix <- rbind(DiagonalMatrix, nextrow) n <- n+1 x <- x+1 } #----plot----------------------------------------------------------------------- a <- ncol(DataMatrix) b <- nrow(DataMatrix) x <- c(1:a) y <- c(1:b) ColorRamp <- rgb( seq(1,0,length=256), # Red seq(1,0,length=256), # Green seq(1,1,length=256)) # Blue ColorRamp2 <- rgb( seq(1,0,length=256), # Red seq(1,0.6,length=256), # Green seq(1,0,length=256)) # Blue ColorRamp3 <- rgb( seq(1,1,length=256), # Red seq(1,0.5,length=256), # Green seq(1,0,length=256)) # Blue #----Mirror matrix (left-right)---- mirror.matrix <- function(x) { xx <- as.data.frame(x); xx <- rev(xx); xx <- as.matrix(xx); xx; } #----Rotate matrix 270 clockworks---- rotate270.matrix <- function(x) { mirror.matrix(t(x)) } UnderMatrix <- rotate270.matrix(UnderMatrix) UpperMatrix <- rotate270.matrix(UpperMatrix) DiagonalMatrix <- rotate270.matrix(DiagonalMatrix) #-----draw matrix--------------------------------------------------------------- #png("D:/Heidi/Master/R_Graphiken/matrix_PairwiseDifferences.png", width = 680, height = 660) #----devide plot region in 4 parts---- def.par <- par(no.readonly = TRUE) # save default, for resetting... layout(rbind(c(1,2), c(1,3), c(1,4)), heights=rbind(c(2,1), c(2,1), c(2,1)), respect=rbind(c(0,1), c(0,0), c(0,0))) #-----draw matrixe plot---- op <- par(mar=c(7.1, 5.1, 9.1, 2.1)) image(x,y, UnderMatrix, col=ColorRamp, xlab="", ylab="", axes = FALSE) image(x,y, UpperMatrix, col=ColorRamp2, xlab="", ylab="", axes = FALSE, add=TRUE) image(x,y, DiagonalMatrix, col=ColorRamp3, xlab="", ylab="", axes = FALSE, add=TRUE) axis(1, at = c(1:a)) axis(2, at = c(1:b), labels=c(b:1)) box() mtext(text="Population", side=1, line=2.5) mtext(text="Population", side=2, line=2.5) mtext(text="Population average pairwise difference", line=4.5, cex=1.2, font=2) par(op) #----draw legends---- #--upper legend-- op2 <- par(mar=c(0.6, 0, 3.6, 7.6)) UpperLegend <- seq(from=0, to=max(UpperMatrix, na.rm=TRUE), length=100) UpperLegend <- as.matrix(UpperLegend) image(y=UpperLegend, t(UpperLegend), col=ColorRamp2, axes=FALSE) axis(side=4, las=2) mtext(text="Average number of pairwise differences between populations (PiXY)", side=4, line=4.5, cex=0.75) box() par(op2) #--diagonal legend-- op3 <- par(mar=c(2.1, 0, 2.1, 7.6)) DiagonalLegend <- seq(from=0, to=max(DiagonalMatrix, na.rm=TRUE), length=100) DiagonalLegend <- as.matrix(DiagonalLegend) image(y=DiagonalLegend, t(DiagonalLegend), col=ColorRamp3, axes=FALSE) axis(side=4, las=2) mtext(text="Average number of pairwise differences within population (PiX)", side=4, line=4.5, cex=0.75) box() par(op3) #--under legend-- op4 <- par(mar=c(3.6, 0, 0.6, 7.6)) UnderLegend <- seq(from=0, to=max(UnderMatrix, na.rm=TRUE), length=100) UnderLegend <- as.matrix(UnderLegend) image(y=UnderLegend, t(UnderLegend), col=ColorRamp, axes=FALSE) axis(side=4, las=2) mtext(text="Corrected average pairwise difference(PiXY-(PiX+PiY)/2)", side=4, line=4.5, cex=0.75) box() par(op4) #dev.off() par(def.par) #reset to default