將下三角矩陣轉成完整矩陣--R

下三角矩陣轉成完整矩陣，對角線為0

只有上三角矩陣，比照辦理將upper.tri() 改成lower.tri()

原始數據 fst.csv ，是由 Arlequin或是MEGA等程式產生 Fst之類的距離矩陣，是下三角矩陣。先整理成只保留 column name，row name要先自己刪除，最後就像一般的數據框一樣

例如：

DMd,DMs,GK,HR,MD

0,,,,

0.04223,0,,,

0.00382,0.05202,0,,

0.03624,0.0194,0.02805,0,

0.30585,0.26611,0.28621,0.36128,0

>read.csv("fst.csv", header=T)->fst   #下三角矩陣

> head(fst)  #看一下

      DMd     DMs      GK      HR MD

1 0.00000      NA      NA      NA NA

2 0.04223 0.00000      NA      NA NA

3 0.00382 0.05202 0.00000      NA NA

4 0.03624 0.01940 0.02805 0.00000 NA

5 0.30585 0.26611 0.28621 0.36128  0

※如果原始數據沒有 column 或是row name，可以

colnames(fst) <-c("DMd","DMs","GK","HR","MD")

> makeSymm <- function(fst) {fst[upper.tri(fst)]<-t(fst)[upper.tri(fst)] 

+ return(fst) } 

#函數makeSymm。 將轉置矩陣中 定義為TRUE的數據取出 （t(fst)[upper.tri(fst)]），然後賦值給原始矩陣中定義為TRUE的位置

>upper.tri(fst)  #看一看這些指令在做什麼。還有一個參數 diag=T，可以考慮

      [,1]  [,2]  [,3]  [,4]  [,5]

[1,] FALSE  TRUE  TRUE  TRUE  TRUE

[2,] FALSE FALSE  TRUE  TRUE  TRUE

[3,] FALSE FALSE FALSE  TRUE  TRUE

[4,] FALSE FALSE FALSE FALSE  TRUE

[5,] FALSE FALSE FALSE FALSE FALSE 

> t(fst)[upper.tri(fst)]   

 [1] 0.04223 0.00382 0.05202 0.03624 0.01940 0.02805 0.30585 0.26611 0.28621

[10] 0.36128

> makeSymm(fst) ->fstM  ＃對fst執行函數

> fstM   ＃檢查一下對不對

      DMd     DMs      GK      HR      MD

1 0.00000 0.04223 0.00382 0.03624 0.30585

2 0.04223 0.00000 0.05202 0.01940 0.26611

3 0.00382 0.05202 0.00000 0.02805 0.28621

4 0.03624 0.01940 0.02805 0.00000 0.36128

5 0.30585 0.26611 0.28621 0.36128 0.00000

或是輸出數據與原始數據在Excel相減，下三角 ＝ 0，確定無誤。

>write.csv(fstM,"fst_matrix.csv")

※後續分析，例如 PCoA

>library(ape)

>pcoa(fstM)->pcoaF

>pcoaF

...局部結果

$vectors

          Axis.1        Axis.2        Axis.3

[1,] -0.06028953  1.647973e-02  0.0087032005

[2,] -0.02505095 -2.823545e-02  0.0006659622

[3,] -0.04271856  2.245764e-02 -0.0067301870

[4,] -0.11741318 -1.067173e-02 -0.0024847685

[5,]  0.24547221 -3.019853e-05 -0.0001542072

>biplot(pcoaF)

或是

> plot(pcoaF$vectors[,1],pcoaF$vectors[,2])  #自己挑選要作圖的軸

> plot(pcoaF$vectors[,1],pcoaF$vectors[,2], xlab="PCoA 1", ylab="PCoA 2", pch=16,col="red")  #修一下

> text(pcoaF$vectors[,1],pcoaF$vectors[,2],labels=c("DMd","DMs","GK","HR","MD"),family="Helvetica")  #點太密時，直接把標籤壓在點上，比較不會搞錯，再輸出到Inkscape手工微調。