====== haplotype distance between/within two populations ====== * file: matrix_HapDistance_between-within.r * document page: 11 \\ **code:** #----read haplotype list---- Data <- read.table("D:/Heidi/Master/R_Daten/HaplotypeDistance/ListHaplotype_betweenBsp.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_betweenBsp.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_betweenBsp.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 } #----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)) } DistanceMatrix <- rotate270.matrix(DistanceMatrix) #----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(DistanceMatrix) b <- nrow(DistanceMatrix) x <- c(1:a) y <- c(1:b) image.plot(x,y,DistanceMatrix, col=ColorRamp, main="haplotype distance matrixbetween/within populations", xlab="Haplotype", ylab="Haplotype", axes = FALSE) contour(DistanceMatrix, add = TRUE) axis(1, at = c(1:a), labels=Data[,1], cex.axis=0.7) 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) 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) ===== mixed data from several populations ===== * file: matrix_HapDistance_between-within_complex.r * document page: 13 \\ **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)