I am trying to compare two distance matrices with R. I would like to
create a XY plot of these matrices and do some linear regression on
it. But, I am a bit new to R, so i have a few questions (I searched in
the documentation with no success).
The first problem is loading a distance matrix into R. This matrix is
the output of a the Phylip program Protdist and lookes like this:
I tried with the scan() function to load the files, but with no
success. How should i load in these files? ....
you can separately load each matrix with two text files.
require(ade4)
mat1 <- read.table("mat1.txt")
nam1 <- mat1[,1]
mat1 <- mat1[,-1]
row.names(mat1) <- names(mat1) <- nam1
mat2 <- read.table("mat2.txt")
nam2 <- mat2[,1]
mat2 <- mat2[,-1]
row.names(mat2) <- names(mat2) <- nam2
dist1 <- mat2dist(mat1)
dist2 <- mat2dist(mat2)
# with mat1:
# n_crassa 0.000000 0.690737 0.895257 0.882576 2.365386
# c_neufor 0.690737 0.000000 0.956910 0.979988 2.103041
# a_thaliana 0.895257 0.956910 0.000000 1.003668 2.724847
# pompep 0.882576 0.979988 1.003668 0.000000 2.065202
# s_cerevis 2.365386 2.103041 2.724847 2.065202 0.000000
# and mat2:
# n_crassa 0.000000 0.739560 0.933986 0.861644 2.207467
# c_neufor 0.739560 0.000000 0.988779 0.925168 1.941141
# a_thaliana 0.933986 0.988779 0.000000 1.007803 2.415320
# pompep 0.861644 0.925168 1.007803 0.000000 2.394490
# s_cerevis 2.207467 1.941141 2.415320 2.394490 0.000000
I think that it?s the more simple solution, NOT the optimal solution. Maybe,
there is an interface to read Phylip file in Bioconductor project(?).
To transform matrix to dist, you can use the function mat2dist of the library
ade4 (see example)
To compare 2 distances matrices, there are many possibility ! For example, if
the distances matrices are Euclidian, you can used directly principal
coordinate analyses (dudi.pco), and co-inertia analysis or procuste
analysis,etc ...
# examples:
# Data preparation
require(ade4)
mat1 <- read.table("mat1.txt")
nam1 <- mat1[,1]
mat1 <- mat1[,-1]
row.names(mat1) <- names(mat1) <- nam1
mat2 <- read.table("mat2.txt")
nam2 <- mat2[,1]
mat2 <- mat2[,-1]
row.names(mat2) <- names(mat2) <- nam2
? mat2dist
dist1 <- mat2dist(mat1)
dist2 <- mat2dist(mat2)
dist1
dist2
is.euclid(dist1)
is.euclid(dist2)
# cool the distances matrices are euclidian :)
# to compare the 2 matrices
# example 1 : mantel test
? mantel.randtest
mt1 <- mantel.randtest(dist1,dist2,nrepet=10)
plot(mt1)
# Example 2: coefficient of vectorial correlation (Escoufier, 1973)
? RVdist.randtest
RV1 <- RVdist.randtest(dist1,dist2,nrepet=10)
plot(RV1)
# Example 3: coinertia analysis
?dudi.pco
# ?cmdscale
?coinertia
?randtest.coinertia
pco1 <- dudi.pco(dist1,nf=3,scannf=F)
pco2 <- dudi.pco(dist2,nf=3,scannf=F)
co1 <- coinertia(pco1,pco2,nf=3,scannf=F)
plot(co1)
testco1 <-randtest.coinertia(co1,nrepet=10)
plot(testco1)
# Example 4: procuste analysis
? procuste
? procuste.randtest
pco1 <- dudi.pco(dist1,nf=3,scannf=F)
pco2 <- dudi.pco(dist2,nf=3,scannf=F)
proc1 <- procuste(pco1$tab,pco2$tab)
plot(proc1)
testproc1 <-procuste.randtest(pco1$tab,pco2$tab,nrepet=10)
plot(testproc1)
....
the choice depend to your outcome.
hopes this help
P.BADY