rm(list=ls())
library(ape)
library(ggtree)
## Registered S3 method overwritten by 'treeio':
## method from
## root.phylo ape
## ggtree v1.16.5 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use ggtree in published research, please cite the most appropriate paper(s):
##
## [36m-[39m Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on phylogeny using ggtree. Molecular Biology and Evolution 2018, accepted. doi: 10.1093/molbev/msy194[36m-[39m Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 2017, 8(1):28-36, doi:10.1111/2041-210X.12628
##
## Attaching package: 'ggtree'
## The following object is masked from 'package:ape':
##
## rotate
tree<-read.tree("http://tonyjhwueng.info/phyrates/ple.nwk")
plot(tree)
vcv(tree)[1,1]
## [1] 25.41057
fulltipname<-tree$tip.label
firstup <- function(x) {
substr(x, 1, 1) <- toupper(substr(x, 1, 1))
x
}
shorttipname<-array(NA,length(tree$tip.label))
for(Index in 1:length(shorttipname)){
shorttipname[Index]<-paste("P.",firstup(strsplit(tree$tip.label,"_")[[Index]][2]),sep = "")
}
shorttipname
## [1] "P.Dorsalis" "P.Ventralis" "P.Angusticlavius"
## [4] "P.Welleri" "P.Punctatus" "P.Wehrlei"
## [7] "P.Websteri" "P.Teyahalee" "P.Cylindraceus"
## [10] "P.Variolatus" "P.Chlorobryonis" "P.Chattahoochee"
## [13] "P.Cheoah" "P.Shermani" "P.Amplus"
## [16] "P.Meridianus" "P.Montanus" "P.Albagula"
## [19] "P.Sequoyah" "P.Ocmulgee" "P.Savannah"
## [22] "P.Grobmani" "P.Kisatchie" "P.Mississippi"
## [25] "P.Kiamichi" "P.Aureolus" "P.Glutinosus"
## [28] "P.Jordani" "P.Metcalfi" "P.Ouachitae"
## [31] "P.Fourchensis" "P.Caddoensis" "P.Kentucki"
## [34] "P.Petraeus" "P.Yonahlossee" "P.Hubrichti"
## [37] "P.Nettingi" "P.Richmondi" "P.Electromorphus"
## [40] "P.Cinereus" "P.Shenandoah" "P.Hoffmani"
## [43] "P.Virginia" "P.Serratus"
tree$tip.label<-shorttipname
ggtree(tree,size=1) + geom_tiplab()+ xlim(NA, 30)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
tree$tip.label<-fulltipname
df<-read.csv("~/Dropbox/FCU/Teaching/Mentoring/2020Fall/MingHanHsu/rcodeMHH/Adams2012-SystBiolData.csv")
head(df)
## X HeadLength BodyWidth Forelimb
## 1 P_albagula 12.320000 8.160000 16.12000
## 2 P_amplus 11.120000 7.580000 16.44000
## 3 P_angusticlavius 7.033333 3.700000 8.10000
## 4 P_aureolus 9.200000 5.314286 12.52857
## 5 P_caddoensis 8.200000 3.600000 9.20000
## 6 P_chattahoochee 10.371429 6.400000 13.67143
?strsplit
spX<-strsplit(as.character(df$X),"_")
spname<-array(NA,length(tree$tip.label))
for(Index in 1:length(tree$tip.label)){
spname[Index]<-paste("Plethodon_", spX[[Index]][2],sep="")
}
spname
## [1] "Plethodon_albagula" "Plethodon_amplus"
## [3] "Plethodon_angusticlavius" "Plethodon_aureolus"
## [5] "Plethodon_caddoensis" "Plethodon_chattahoochee"
## [7] "Plethodon_cheoah" "Plethodon_chlorobryonis"
## [9] "Plethodon_cinereus" "Plethodon_cylindraceus"
## [11] "Plethodon_dorsalis" "Plethodon_electromorphus"
## [13] "Plethodon_fourchensis" "Plethodon_glutinosus"
## [15] "Plethodon_grobmani" "Plethodon_hoffmani"
## [17] "Plethodon_hubrichti" "Plethodon_jordani"
## [19] "Plethodon_kentucki" "Plethodon_kiamichi"
## [21] "Plethodon_kisatchie" "Plethodon_meridianus"
## [23] "Plethodon_metcalfi" "Plethodon_mississippi"
## [25] "Plethodon_montanus" "Plethodon_nettingi"
## [27] "Plethodon_ocmulgee" "Plethodon_ouachitae"
## [29] "Plethodon_petraeus" "Plethodon_punctatus"
## [31] "Plethodon_richmondi" "Plethodon_savannah"
## [33] "Plethodon_sequoyah" "Plethodon_serratus"
## [35] "Plethodon_shenandoah" "Plethodon_shermani"
## [37] "Plethodon_teyahalee" "Plethodon_variolatus"
## [39] "Plethodon_ventralis" "Plethodon_virginia"
## [41] "Plethodon_websteri" "Plethodon_wehrlei"
## [43] "Plethodon_welleri" "Plethodon_yonahlossee"
df$X<-spname
HeadLength<-df$HeadLength
names(HeadLength)<-spname
BodyWidth<-df$BodyWidth
names(BodyWidth)<-spname
#ace(x=HeadLength, phy=tree)
compar.gee(HeadLength ~ BodyWidth, phy = tree)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept) BodyWidth
## 2.865456 1.190700
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept) BodyWidth
## 2.865456 1.190700
## Call: compar.gee(formula = HeadLength ~ BodyWidth, phy = tree)
## Number of observations: 44
## Model:
## Link: identity
## Variance to Mean Relation: gaussian
##
## QIC: 53.55486
##
## Summary of Residuals:
## Min 1Q Median 3Q Max
## -1.96753636 -0.79230097 -0.07935879 0.66661590 2.03557941
##
##
## Coefficients:
## Estimate S.E. t Pr(T > |t|)
## (Intercept) 2.748068 0.46810002 5.870687 5.033548e-05
## BodyWidth 1.214010 0.04096275 29.636915 1.585398e-13
##
## Estimated Scale Parameter: 0.8682678
## "Phylogenetic" df (dfP): 15.2847