功能性狀的譜系保守性分析是架起進化歷程與生態過程的重要橋梁，是很多高階分析的基礎。譜系保守性是很多基於譜系的預測模型的基本前提，比如生態位模型、達爾文歸化假說、預適應假說。去除譜系保守性的影響，是跨物種進行性狀相關性分析的前提（側重於適應性的協同進化分析，而非系統發育分析）。
本次分享於2021年4月1日18:30在N2-146進行。
示例數據
Species H_N0 H_N2 H_N5 H_N10 H_N20 H_N50 Apocynum venetum 16.42 30.52 37.38 39.98 32.75 34.63 Cynanchum chinense 21.5 45.35 72.96 85.39 82.53 82 Limonium sinense 6.41 10.63 14.34 17.08 17.13 15.86 Chloris virgata 37.06 53.33 72.92 81.41 83.06 78.25 Setaria viridis 31.18 56.84 68.36 73.67 69.77 58.54 Suaeda glauca 24.54 38.52 41.47 45.49 46.29 43.86 Suaeda salsa 18.93 31.84 39.89 41.46 40.96 32.93 Phragmites australis 34.22 44.26 48.48 54.84 54.23 54.49 Xanthium strumarium 17.44 33.25 43 46.7 43.34 30.93 Xanthium sibiricum 16.34 23.03 33.4 43.43 44.03 37.32 Oenothera biennis 6.32 13.14 16.93 16.84 18.26 14.6 Chenopodium ficifolium 11.03 29.57 48.94 39.71 38.36 32.28 Salsola tragus 21.92 37.82 37.02 37.53 34.97 33.21 Polygonum sibiricum 13.92 22.59 27.43 23.57 21.05 22.59 Leymus mollis 30.89 39.74 44.83 38.86 38.86 42.49 Leonurus japonicus 1.98 6.45 10.12 11.72 13.31 8.73 Lactuca indica 10.22 21.82 29.01 28.99 26.33 23.14
示例代碼
#2021年4月1日 劉樂樂
#組會分享:功能性狀的譜系信號與譜系獨立比較的分析方法
#數據集來自祁魯玉：黃河三角洲濕地植物不同物種對氮沉降響應實驗
#V.PhyloMaker包的安裝需要借助devtools::install_github,如下
library("devtools") #提供安裝V.PhyloMaker包和plantlist物種名錄所需要的工具
install_github("jinyizju/V.PhyloMaker")#安裝V.PhyloMaker
install_github("helixcn/plantlist")#安裝plantlist
####載入需要的軟件包####
library("picante")
library("V.PhyloMaker")
library("plantlist") #物種名錄，匹配分類地位（界門綱目科屬種）
#載入數據集
trait <- read.table("data.txt", header = TRUE, sep = "\t", row.names = 1) #第一行為性狀名稱，第一列為物種名；務必保證沒有多余空格
####建立譜系樹####
#搜索物種的屬和科
species_list <- rownames(trait)
sp_lis <- subset(TPL(species_list),select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
colnames(sp_lis) <- c("species", "genus", "family")
#建樹
sp_tree <- phylo.maker(sp.lis=sp_lis) #需要運行一小段時間（幾分鍾），需要耐心等待
sp_tr <- sp_tree$scenario.3
plot(sp_tr)
write.tree(sp_tr, "sp_tree.newick") #寫入文件存儲
####譜系保守性分析####
row.names(trait) <- gsub(" ", "_", rownames(trait)) #物種名與譜系樹的label保持一致
#單個性狀 phylosignal\Kcalc
trait1 <- trait[["H_N0"]]
names(trait1) <- rownames(trait)
phylosignal(trait1,sp_tr)
Kcalc(trait1, sp_tr)
#多個性狀multiPhylosignal\Kcalc
multiPhylosignal(trait, sp_tr)
apply(trait,2,Kcalc, sp_tr)
#校對排序后，順序一致，無需進行此步驟
data_matched <- match.phylo.data(phy = sp_tr, data = trait)
multiPhylosignal(data_matched$data, data_matched$phy) #結果同上
#可視化
trait_full <- trait
trait_full$sp <- row.names(trait_full)
color.plot.phylo(sp_tr,df = trait_full, trait = "H_N0", taxa.names = "sp", num.breaks = 3)
#trait_full$repsonse <- (trait_full$H_N0 + trait_full$H_N50)/(trait_full$H_N0)
#multiPhylosignal(trait_full[,-7], sp_tr)
#color.plot.phylo(sp_tr,df = trait_full, trait = "repsonse", taxa.names = "sp", num.breaks = 3)
####譜系獨立比較###
#譜系獨立的相關性分析
trait1_pic <- pic(trait1,sp_tr) #對單個性狀計算pic值
trait_pic <- apply(trait, 2, pic, sp_tr)#對所有性狀計算pic值
cor.test(trait[[1]],trait[[5]])
cor.test(trait_pic[,1],trait_pic[,5])
par(mfrow=c(1,2))
corrplot::corrplot(cor(trait), type = "upper")
corrplot::corrplot(cor(trait_pic), type = "upper")
R語言及R studio使用小技巧
- 查看數據結構 str、class、head
- 查看幫助文件?、help、F1
- 運行例子 example
- 查看原始文獻 citation
- 查看源代碼 直接輸入函數名
- 需要學習的軟件包 ade4、ape、picante、vegan
- 快捷鍵 ctrl+Enter
PPT及代碼、數據文件可在百度網盤下載
鏈接：https://pan.baidu.com/s/1RUfmZAMkF1b9vuM3EwA7DQ
提取碼：phyl
復制這段內容后打開百度網盤手機App，操作更方便哦