• 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)