====== haplotype distance between/within populations and groups ====== * file: matrix_multiplePlots.r * document page: 25 \\ **code:** #----read data------------------------------------------------------------------ #----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) #----whole DistanceMatrix both Populations together:---- Pop_together <- c(Pop1, Pop2) wholeDistanceMatrixTogether <- wholeDistanceMatrix[Pop_together, Pop_together[c(length(Pop_together)):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)) } DistanceMatrixTogether <- rotate270.matrix(DistanceMatrixTogether) #----draw matrix---------------------------------------------------------------- windows() #new graphic window 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) png("D:/Heidi/Master/R_Graphiken/matrix_HapDistance_between-within Pop and Groups.png", width = 680, height = 660) op <- par(mfrow=c(2,2), oma=c(0, 0 ,3, 0)) #----matrix one---- op1 <- par(mfg=c(1,1),mar=c(0.2, 5.2, 4.2, 0.2)) image(x,y,DistanceMatrixTogether, col=ColorRamp, ylab="Group 1", font.lab=2, axes = FALSE) contour(DistanceMatrixTogether, add = TRUE) axis(2, at = c(1:b), labels=c(Pop2[NROW(Pop2):1],Pop1[NROW(Pop1):1]), cex.axis=0.6, las=2) box(bty="l") 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 = 2, at =(3.2*Row/4), line = 2.7, text = "Population 1", cex=0.7) mtext(side = 2, at =(Row/4), line = 2.7, text = "Population 2", cex=0.7) par(op1) #----matrix two---- op2 <- par(mfg=c(2,1), mar=c(4.2, 5.2, 0.2, 0.2)) image(x,y, wholeDistanceMatrixTogether, col=ColorRamp, xlab="Group 1", ylab="Group 2", font.lab=2, axes = FALSE) contour(DistanceMatrixTogether, add = TRUE) axis(1, at = c(1:a), labels=c(Pop_together), cex.axis=0.6, las=2) axis(2, at = c(1:b), labels=c(Pop2[NROW(Pop2):1],Pop1[NROW(Pop1):1]), cex.axis=0.6, las=2) box() lines(c(0,NROW(Pop1)+ NROW(Pop2)+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(Pop1)+ NROW(Pop2)+0.5), lwd=2) mtext(side = 1, at =(Row/4), line = 2.5, text = "Population 1", cex=0.7) mtext(side = 1, at =(3.2*Row/4), line = 2.5, text = "Population 2", cex=0.7) mtext(side = 2, at =(3.2*Row/4), line = 2.7, text = "Population 3", cex=0.7) mtext(side = 2, at =(Row/4), line = 2.7, text = "Population 4", cex=0.7) par(op2) #----matrix three---- op3 <- par(mfg=c(2,2), mar=c(4.2, 0.2, 0.2, 6.2)) image(x,y,DistanceMatrixTogether, col=ColorRamp, xlab="Group 2", font.lab=2, axes = FALSE) contour(DistanceMatrixTogether, add = TRUE) axis(1, at = c(1:a), labels=c(Pop1,Pop2), cex.axis=0.6, las=2) box(bty="l") 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.5, text = "Population 3", cex=0.7) mtext(side = 1, at =(3.2*Row/4), line = 2.5, text = "Population 4", cex=0.7) par(op3) #----legend---- maximum <- max(DistanceMatrixTogether, wholeDistanceMatrixTogether, DistanceMatrixTogether , na.rm = TRUE) Legend <- seq(from=0, to=maximum, length=100) Legend <- as.matrix(Legend) op4 <- par(mfg=c(1,2), mar=c(0, 15.6, 3.2, 6.2) ) image(y=Legend, t(Legend), col=ColorRamp, axes=FALSE) axis(side=4, las=2) box() title("haplotype distance matrix between/within populations and groups", line=0, outer=TRUE) par(op4) #----reset par parameter---- par(op) dev.off()