Differences

This shows you the differences between two versions of the page.

--- population_average_pairwise_difference [2008/01/28 11:40] – created heidi
+++ population_average_pairwise_difference [2008/07/22 13:31] (current) – external edit 127.0.0.1
@@ Line 1: / Line 1: @@
-====== mismatch distribution ======
+====== Population average pairwise difference ======
-  * file: lines_Mismatch.r
+  * file: matrix_PairwiseDifferences.r
-  * document page: 1
+  * document page: 33
 \\
 **code:**
 <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
+</code>