Differences

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

--- haplotype_distance_between_within_two_populations [2007/12/21 13:34] – created heidi
+++ haplotype_distance_between_within_two_populations [2008/07/22 13:31] (current) – external edit 127.0.0.1
@@ Line 70: / Line 70: @@
         axis(2, at = c(1:b), labels=Data[(Row:1),1], cex.axis=0.7)
         box()
         half <- (Row/2) + 0.5
         lines(c(0,half),c(half,half), lwd=2)
         lines(c(half,half),c(0,half), lwd=2)
         mtext(side = 1, at =(Row/4), line = 2, text = "Population 1", cex=0.8)
         mtext(side = 1, at =(3*Row/4), line = 2, text = "Population 2", cex=0.8)
@@ Line 81: / Line 80: @@
         mtext(side = 2, at =(3*Row/4), line = 2, text = "Population 1", cex=0.8)
 </code>
+===== mixed data from several populations =====
+  * file: matrix_HapDistance_between-within_complex.r
+  * document page: 13
+\\
+**code:**
+<code>
+#----read haplotype list----
+Data <- read.table("D:/Heidi/Master/R_Daten/HaplotypeDistance/ListHaplotype_betweenPop.txt", skip=1)
+Row <- nrow(Data)
+Columns <- Row
+#----read data row by row:----
+x <- 0
+n <- 1
+DistanceMatrix <- as.matrix(scan("D:/Heidi/Master/R_Daten/HaplotypeDistance/HapDistanceMatrix_betweenPop.txt",
+                                  what=double(0), skip=x, nlines=1, nmax=n))
+DistanceMatrix <- cbind(DistanceMatrix, matrix(NA, ncol=Columns, nrow=1))
+DistanceMatrix <- DistanceMatrix[,1:Columns]
+n <- n + 1
+x <- x + 1
+while(n<(Row+1)){
+  nextrow <- as.matrix(scan("D:/Heidi/Master/R_Daten/HaplotypeDistance/HapDistanceMatrix_betweenPop.txt",
+                             what=double(0), skip=x, nlines=1, nmax=n))
+  nextrow <- cbind(t(nextrow), matrix(NA, ncol=Columns, nrow=1))
+  nextrow <- nextrow[,1:Columns]
+  DistanceMatrix <- rbind(DistanceMatrix, nextrow)
+  n <- n + 1
+  x <- x + 1
+}
+#----select the data of DistanceMatrix----
+dimnames(DistanceMatrix) <- list(Data[,1], Data[,1])
+#--population1:
+Pop1 <- read.table("D:/Heidi/Master/R_Daten/HaplotypeDistance/HapDistanceMatrix_withinPop1.txt", skip=5)
+Pop1 <- as.character(Pop1[,1])
+DistanceMatrixPop1 <- DistanceMatrix[Pop1, Pop1]
+#--population2:
+Pop2 <- read.table("D:/Heidi/Master/R_Daten/HaplotypeDistance/HapDistanceMatrix_withinPop2.txt", skip=5)
+Pop2 <- as.character(Pop2[,1])
+DistanceMatrixPop2 <- DistanceMatrix[Pop2, Pop2]
+#--between population1/2:
+#whole DistanceMatrix
+wholeDistanceMatrix <- DistanceMatrix
+x <- 1
+while(x <= ncol(wholeDistanceMatrix)){
+  twholeDistanceMatrix <- t(wholeDistanceMatrix)
+  wholeDistanceMatrix[x,(x:ncol(wholeDistanceMatrix))] <- twholeDistanceMatrix[x,(x:ncol(wholeDistanceMatrix))]
+  x <- x + 1
+}
+wholeDistanceMatrixBetween <- wholeDistanceMatrix[Pop2, Pop1]
+#----together:-----
+DistanceMatrixUp <- cbind(DistanceMatrixPop1, matrix(NA, ncol(DistanceMatrixPop1), nrow(DistanceMatrixPop2)))
+DistanceMatrixLo <- cbind(wholeDistanceMatrixBetween, DistanceMatrixPop2)
+DistanceMatrixTogether <- rbind(DistanceMatrixUp, DistanceMatrixLo)
+#----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))
+}
+DistanceMatrixTogether <- rotate270.matrix(DistanceMatrixTogether)
+#----draw matrix----
+library(fields)
+ ColorRamp <- rgb( seq(1,0,length=256),  # Red
+                   seq(1,0,length=256),  # Green
+                   seq(1,1,length=256))  # Blue
+ a <- ncol(DistanceMatrixTogether)
+ b <- nrow(DistanceMatrixTogether)
+ x <- c(1:a)
+ y <- c(1:b)
+ image.plot(x,y,DistanceMatrixTogether, col=ColorRamp, main="haplotype distancematrix between/within populations",
+             xlab="Haplotype", ylab="Haplotype", axes = FALSE)
+        contour(DistanceMatrixTogether, add = TRUE)
+        axis(1, at = c(1:a), labels=c(Pop1,Pop2), cex.axis=0.6)
+        axis(2, at = c(1:b), labels=c(Pop2[NROW(Pop2):1],Pop1[NROW(Pop1):1]), cex.axis=0.6)
+        box()
+        lines(c(0,NROW(Pop1)+0.5),c(NROW(Pop2)+0.5,NROW(Pop2)+0.5), lwd=2)
+        lines(c(NROW(Pop1)+0.5,NROW(Pop1)+0.5),c(0,NROW(Pop2)+0.5), lwd=2)
+        mtext(side = 1, at =(Row/4), line = 2, text = "Population 1", cex=0.8)
+        mtext(side = 1, at =(3*Row/4), line = 2, text = "Population 2", cex=0.8)
+        mtext(side = 2, at =(Row/4), line = 2, text = "Population 2", cex=0.8)
+        mtext(side = 2, at =(3*Row/4), line = 2, text = "Population 1", cex=0.8)
+</code>