Duck Phyloaccc
<<<<<<< HEAD2026-07-27
=======2026-01-12
>>>>>>> 50c7f3404db70be3474bff72ea1e16567c81bd1eNote
Following includes codes used to do shape analysis of duck bones and then run phyloacc-c
1 Introduction
We combined surface scans and CT/μCT scans to create a database of skeletal elements for three hindlimb bones (femur, tibiotarsus, and tarsometatarsus) along with the four forelimb bone (humerus, radius, ulna, and carpometacarpus) for 65 species of ducks. Our sampling for the species and the bone scans for those species are presented in the following table.
| Images | Scientific_Name | Common_Name | Voucher_Info | Genomes_NCBI | Femur | Tibiotarsus | Tarsometatarsus | Humerus | Radius | Ulna | Carpometacarpus |
|---|---|---|---|---|---|---|---|---|---|---|---|
|
Aix galericulata | Mandarin Duck | USNM 622444 | GCA_024635365.1 | X | X | X | X | X | X | X |
|
Aix sponsa | Wood Duck | MCZ:Orn:347504 | GCA_032270905.1 | X | X | X | X | X | X | X |
|
Amazonetta brasiliensis | Brazilian Teal | USNM 631047 | X | X | X | X | X | X | X | |
|
Anas acuta | Northern Pintail | MCZ:Orn:341895 | GCA_963932015.1 | X | X | X | ||||
|
Anas aucklandica | Auckland Islands Teal | UMZC:Vertebrates:12/Ana/42/a/8-10 | X | X | X | X | X | X | X | |
|
Anas crecca | Green-winged Teal | MCZ:Orn:347036 | X | X | X | X | X | |||
|
Anas erythrorhyncha | Red-billed Duck | NHMUK:Zoo:S/1952.1.159 | X | X | X | X | ||||
|
Anas laysanensis | Laysan Duck | MCZ:Orn:347544 | X | X | X | |||||
|
Anas platyrhynchos | Mallard | MCZ:Orn:347156 | GCA900411745.1 | X | X | X | X | X | X | X |
|
Anas rubripes | American Black Duck | MCZ:Orn:337845 | X | X | X | X | X | X | ||
|
Anas zonorhyncha | Eastern Spot-billed Duck | USNM 633481 | GCA_002224875.1 | X | X | X | X | X | X | X |
|
Anhima cornuta | Horned Screamer | MCZ:Orn:346993 | GCA_013399475.1 | X | X | X | X | X | X | |
|
Anser albifrons | Greater White-fronted Goose | NHMUK:Zoo:1895.2.6.8 | X | X | X | X | X | X | X | |
|
Anser caerulescens | Snow Goose | MCZ:Orn:341885 | X | X | X | X | X | X | ||
|
Anser cygnoides | Swan Goose | USNM 19637 | GCA_026259575.1 | X | X | X | X | X | ||
|
Anser erythropus | Lesser White-fronted Goose | MCZ:Orn:340330 | X | X | X | |||||
|
Anser fabalis | Taiga Bean-Goose | MCZ:Orn:340262 | GCA_002592135.1 | X | X | X | X | X | ||
|
Anser indicus | Bar-headed Goose | USNM 557494 | GCA_025583725.1 | X | X | X | X | X | X | |
|
Anseranas semipalmata | Magpie Goose | NHMUK:Zoo:1852.7.22.1 | GCA_013399115.1 | X | X | X | X | X | ||
|
Asarcornis scutulata | White-winged Duck | USNM 612610 | GCA_013398475.1 | X | X | X | X | X | X | X |
|
Aythya americana | Redhead | MCZ:Orn:346918 | X | X | X | X | X | X | ||
|
Aythya fuligula | Tufted Duck | USNM 641850 | GCA_009819795.1 | X | X | X | X | X | X | X |
|
Branta bernicla | Brant | MCZ:Orn:336993 | X | X | X | X | X | X | X | |
|
Branta canadensis | Canada Goose | MCZ:Orn:346738 | GCA_006130075.1 | X | X | X | X | X | X | X |
|
Branta leucopsis | Barnacle Goose | MCZ:Orn:346732 | X | X | X | X | X | X | X | |
|
Bucephala clangula | Common Goldeneye | MCZ:Orn:347155 | X | X | X | X | X | X | X | |
|
Cairina moschata | Muscovy Duck | MCZ:Orn:348425 | GCA_018104995.1 | X | X | X | X | X | X | |
|
Cereopsis novaehollandiae | Cape Barren Goose | MCZ:Orn:347095 | X | X | X | X | X | X | X | |
|
Chauna chavaria | Northern Screamer | MCZ:Orn:340307 | X | X | X | X | ||||
|
Chenonetta jubata | Maned Duck | MCZ:Orn:339522 | X | X | X | |||||
|
Chloephaga hybrida | Kelp Goose | MCZ:Orn:340225 | X | X | X | X | X | X | ||
|
Chloephaga picta | Upland Goose | MCZ:Orn:342847 | X | X | X | X | ||||
|
Chloephaga poliocephala | Ashy-headed Goose | NHMUK:Zoo:S/1963.10.1 | X | X | X | X | X | X | ||
|
Chloephaga rubidiceps | Ruddy-headed Goose | MCZ:Orn:343209 | X | X | X | X | X | X | ||
|
Clangula hyemalis | Long-tailed Duck | MCZ:Orn:346497 | GCA_029619115.1 | X | X | X | X | X | ||
|
Cygnus columbianus | Tundra Swan | MCZ:Orn:343048 | X | X | X | X | ||||
|
Cygnus olor | Mute Swan | NHMUK:Zoo:1854.6.25.2 | GCA_009769625.2 | X | X | X | X | X | X | X |
|
Dendrocygna autumnalis | Black-bellied Whistling-Duck | MCZ:Orn:335880 | X | X | X | X | X | X | X | |
|
Dendrocygna bicolor | Fulvous Whistling-Duck | MCZ:Orn:347071 | X | X | X | X | X | X | X | |
|
Heteronetta atricapilla | Black-headed Duck | NHMUK:Zoo:S/2002.47.6 | GCA_011075105.1 | X | X | X | X | X | X | |
|
Histrionicus histrionicus | Harlequin Duck | NHMUK:Zoo:1843.6.23.23 | X | X | X | X | X | X | ||
|
Lophodytes cucullatus | Hooded Merganser | MCZ:Orn:347140 | X | X | X | X | X | X | ||
|
Lophonetta specularioides | Crested Duck | USNM 490992 | GCA_032165235.1 | X | X | X | X | X | X | X |
|
Malacorhynchus membranaceus | Pink-eared Duck | USNM 553596 | X | X | X | X | X | X | X | |
|
Mareca strepera | Gadwall | MCZ:Orn:341872 | GCA_035583915.1 | X | X | X | X | X | ||
|
Mareca sibilatrix | Chiloe Wigeon | NHMUK:Zoo:S/1985.72.2 | GCA_034783215.1 | X | X | X | X | X | X | X |
|
Marmaronetta angustirostris | Marbled Teal | USNM 500299 | X | X | X | X | X | X | X | |
|
Melanitta fusca | Velvet Scoter | MCZ:Orn:348703 | GCA_029620185.1 | X | X | X | X | X | ||
|
Melanitta perspicillata | Surf Scoter | MCZ:Orn:347151 | X | X | X | X | X | |||
|
Merganetta armata | Torrent Duck | MCZ:Orn:345094 | X | X | X | X | X | X | X | |
|
Mergellus albellus | Smew | MCZ:Orn:348520 | X | X | X | X | X | X | X | |
|
Mergus merganser | Common Merganser | MCZ:Orn:340318 | X | X | ||||||
|
Mergus serrator | Red-breasted Merganser | MCZ:Orn:343636 | X | X | X | X | X | X | X | |
|
Nettapus auritus | African Pygmy-Goose | MCZ:Orn:341374 | GCA_011076525.1 | X | X | X | ||||
|
Oxyura jamaicensis | Ruddy Duck | NHMUK:Zoo:S/1985.70.4 | GCA_011077185.1 | X | X | X | X | X | X | |
|
Oxyura ferruginea | Andean Duck | MCZ:Orn:341437 | X | X | X | X | X | |||
|
Sarkidiornis melanotos | Knob-billed Duck | USNM 488140 | X | X | X | X | X | X | X | |
|
Somateria mollissima | Common Eider | MCZ:Orn:336958 | GCA_951411735.1 | X | X | X | X | X | X | |
|
Spatula clypeata | Northern Shoveler | MCZ:Orn:347105 | X | X | X | X | X | X | ||
|
Spatula discors | Blue-winged Teal | NHMUK:Zoo:S/2001.41.5 | GCA_032165275.1 | X | X | X | X | X | X | X |
|
Stictonetta naevosa | Freckled Duck | USNM 555594 | GCA_011074415.1 | X | X | X | X | X | X | X |
|
Tachyeres brachypterus | Falkland Steamer-Duck | MCZ:Orn:342206 | GCA_032357605.1 | X | X | X | X | X | X | X |
|
Tachyeres patachonicus | Flying Steamer-Duck | USNM 491013 | X | X | X | X | X | X | X | |
|
Tachyeres pteneres | Flightless Steamer-Duck | USNM 490930 | X | X | X | X | X | X | X | |
|
Tadorna tadorna | Common Shelduck | MCZ:Orn:347538 | X | X | X | X | X | X | X |
2 Getting the samples
- Surface Scanner: For the surface scanning we used EinScan SE Desktop scanner. After calibration, we put individual bones on the pedestal and ran the algorithm in the program EXScan S v3.1.2.0. Meshes were saved as .obj files.
- μCT scans: For the μCT scanning we used the μCT scanner at the Shared Use Facility at the Harvard Museum of Comparative Zoology. Slices were saved as .tif files.
- Morphosource scans: For several species, we used material available in the Morphosource public database. We either downloaded the segmented bone mesh or downloaded the raw CT files and used Slicer to segment out the necessary bones.
After all the scans were collected, all scans were roughly aligned at the same coordinate space using the Align function in Meshmixer. For CT scans, we removed any internal structure by using the Select function to select the outer surface (setting the Crease Angle Thres option to >70 that prevents the selection from moving to the internal surface through holes) and using the Invert function to switch to internal surface, which were then deleted. We used Make Solid function to close any holes on the outside mesh, with a comparable Solid Accuracy and Mesh Density to avoid over smoothing of scans. Finally we exported this scan as a .obj file.
3 Landmarking
3.1 Pseudolandmarking
For landmarking, we used the pseudolandmarking (Rolfe, Davis, et al. 2021) approach implemented through PseudoLMGenerator function of SlicerMorph (Rolfe, Pieper, et al. 2021) in Slicer. We used the Mallard Duck (Anas platyrhynchos) for the reference specimen to generate the primary landmark for all bones except the tibiotarsus (where the mallard specimen lacks the fibula) where we used the closely related Spot-billed Duck (Anas zonorhyncha). Spacing tolerance for each bone was adjusted to apply ~400 pseudolandmarks per bone.
| Bone | Spatial Tolerance | Number of Pseudolandmarks | Point Density Adjustment |
|---|---|---|---|
| Femur | 2.7 | 459 | 1.7 |
| Tibiotarsus | 2.0 | 480 | 1.9 |
| Tarsometatarsus | 2.8 | 431 | 1.7 |
| Humerus | 2.4 | 458 | 1.9 |
| Radius | 1.8 | 432 | 2.6 |
| Ulna | 2.2 | 452 | 2.0 |
| Carpometacarpus | 2.6 | 432 | 2.6 |
3.2 ALPACA
To transfer pseudolanmarks from one species to another, we used the single-template ALPACA framework (Porto et al. 2021). We adjusted the point density adjustment to create meshes of ~4000 points as recommended for ALPACA (values in table above). We used a BCPD acceleration for processing the samples. We used a batch processing to generate landmarks for bones of all the species.
3.3 Phylogeny
We used the phylogeny from () for this study. For a few taxa missing in the tree that we had scans for, we appended a branch next to the known closest relative species.
duck_data <- xlsx::read.xlsx("Ducks_MCZ.xlsx", sheetIndex = 1, header = T)
duck_data <- duck_data[-66,]
#get body size
duckbm <- log10(duck_data$Body.mass.g)
names(duckbm) <- duck_data$Scientific_Name
#get flight data
flight <- as.factor(duck_data$Flight)
names(flight) <- duck_data$Scientific_Name
flight## Aix_galericulata Aix_sponsa
## 1 1
## Amazonetta_brasiliensis Anas_acuta
## 1 1
## Anas_aucklandica Anas_crecca
## 0 1
## Anas_erythrorhyncha Anas_laysanensis
## 1 1
## Anas_platyrhynchos Anas_rubripes
## 1 1
## Anas_zonorhyncha Anhima_cornuta
## 1 1
## Anser_albifrons Anser_caerulescens
## 1 1
## Anser_cygnoides Anser_erythropus
## 1 1
## Anser_fabalis Anser_indicus
## 1 1
## Anseranas_semipalmata Asarcornis_scutulata
## 1 1
## Aythya_americana Aythya_fuligula
## 1 1
## Branta_bernicla Branta_canadensis
## 1 1
## Branta_leucopsis Bucephala_clangula
## 1 1
## Cairina_moschata Cereopsis_novaehollandiae
## 1 1
## Chauna_chavaria Chenonetta_jubata
## 1 1
## Chloephaga_hybrida Chloephaga_picta
## 1 1
## Chloephaga_poliocephala Chloephaga_rubidiceps
## 1 1
## Clangula_hyemalis Cygnus_columbianus
## 1 1
## Cygnus_olor Dendrocygna_autumnalis
## 1 1
## Dendrocygna_bicolor Heteronetta_atricapilla
## 1 1
## Histrionicus_histrionicus Lophodytes_cucullatus
## 1 1
## Lophonetta_specularioides Malacorhynchus_membranaceus
## 1 1
## Mareca_strepera Mareca_sibilatrix
## 1 1
## Marmaronetta_angustirostris Melanitta_fusca
## 1 1
## Melanitta_perspicillata Merganetta_armata
## 1 1
## Mergellus_albellus Mergus_merganser
## 1 1
## Mergus_serrator Nettapus_auritus
## 1 1
## Oxyura_jamaicensis Oxyura_ferruginea
## 1 1
## Sarkidiornis_melanotos Somateria_mollissima
## 1 1
## Spatula_clypeata Spatula_discors
## 1 1
## Stictonetta_naevosa Tachyeres_brachypterus
## 1 0
## Tachyeres_patachonicus Tachyeres_pteneres
## 1 0
## Tadorna_tadorna
## 1
## Levels: 0 1#get diving data
diving <- as.factor(duck_data$Diving)
names(diving) <- duck_data$Scientific_Name
diving## Aix_galericulata Aix_sponsa
## Dabbling Dabbling
## Amazonetta_brasiliensis Anas_acuta
## Dabbling Occasional_diving
## Anas_aucklandica Anas_crecca
## Occasional_diving Dabbling
## Anas_erythrorhyncha Anas_laysanensis
## Dabbling Dabbling
## Anas_platyrhynchos Anas_rubripes
## Occasional_diving Diving
## Anas_zonorhyncha Anhima_cornuta
## Dabbling Walking
## Anser_albifrons Anser_caerulescens
## Walking Walking
## Anser_cygnoides Anser_erythropus
## Walking Surface_swimmer
## Anser_fabalis Anser_indicus
## Surface_swimmer Surface_swimmer
## Anseranas_semipalmata Asarcornis_scutulata
## Surface_swimmer Occasional_diving
## Aythya_americana Aythya_fuligula
## Diving Diving
## Branta_bernicla Branta_canadensis
## Dabbling Occasional_diving
## Branta_leucopsis Bucephala_clangula
## Surface_swimmer Diving
## Cairina_moschata Cereopsis_novaehollandiae
## Surface_swimmer Walking
## Chauna_chavaria Chenonetta_jubata
## Walking Surface_swimmer
## Chloephaga_hybrida Chloephaga_picta
## Surface_swimmer Walking
## Chloephaga_poliocephala Chloephaga_rubidiceps
## Walking Walking
## Clangula_hyemalis Cygnus_columbianus
## Diving Dabbling
## Cygnus_olor Dendrocygna_autumnalis
## Dabbling Dabbling
## Dendrocygna_bicolor Heteronetta_atricapilla
## Occasional_diving Diving
## Histrionicus_histrionicus Lophodytes_cucullatus
## Diving Diving
## Lophonetta_specularioides Malacorhynchus_membranaceus
## Dabbling Dabbling
## Mareca_strepera Mareca_sibilatrix
## Dabbling Dabbling
## Marmaronetta_angustirostris Melanitta_fusca
## Dabbling Diving
## Melanitta_perspicillata Merganetta_armata
## Diving Diving
## Mergellus_albellus Mergus_merganser
## Diving Diving
## Mergus_serrator Nettapus_auritus
## Diving Dabbling
## Oxyura_jamaicensis Oxyura_ferruginea
## Diving Diving
## Sarkidiornis_melanotos Somateria_mollissima
## Dabbling Diving
## Spatula_clypeata Spatula_discors
## Occasional_diving Dabbling
## Stictonetta_naevosa Tachyeres_brachypterus
## Dabbling Occasional_diving
## Tachyeres_patachonicus Tachyeres_pteneres
## Diving Occasional_diving
## Tadorna_tadorna
## Dabbling
## Levels: Dabbling Diving Occasional_diving Surface_swimmer Walkingduck_tree <- read.tree("duck_tree.tre")
# Root at the Anhimidae
duck_tree <- root(duck_tree, outgroup = c("Anhima_cornuta", "Chauna_chavaria", "Chauna_torquata"))
# Add missing species
duck_tree <- bind.tip(duck_tree, "Anas_aucklandica", where = which(duck_tree$tip.label=="Anas_chlorotis"), edge.length = 0.0143477498/4, position = 0.0143477498/4)
duck_tree <- bind.tip(duck_tree, "Anas_zonorhyncha", where = which(duck_tree$tip.label=="Anas_platyrhynchos"), edge.length = 0.0022507147/4, position = 0.0022507147/4)
duck_tree <- bind.tip(duck_tree, "Anser_erythropus", where = which(duck_tree$tip.label=="Anser_albifrons"), edge.length = 0.0017958047/4, position = 0.0017958047/4)
duck_tree <- bind.tip(duck_tree, "Asarcornis_scutulata", where = which(duck_tree$tip.label=="Cairina_moschata"), edge.length = 0.0067085060/4, position = 0.0067085060/4)
duck_tree$tip.label[which(duck_tree$tip.label == "Malacorhyncus_membranaceus")] <- "Malacorhynchus_membranaceus"
duck_tree <- bind.tip(duck_tree, "Mergus_merganser", where = which(duck_tree$tip.label=="Mergus_serrator"), edge.length = 0.0027206995/4, position = 0.0027206995/4)
duck_tree <- bind.tip(duck_tree, "Oxyura_ferruginea", where = which(duck_tree$tip.label=="Oxyura_jamaicensis"), edge.length = 0.0024525321/4, position = 0.0024525321/4)
duck_tree <- bind.tip(duck_tree, "Spatula_discors", where = which(duck_tree$tip.label=="Spatula_cyanoptera"), edge.length = 0.0027173372/4, position = 0.0027173372/4)
duck_tree <- bind.tip(duck_tree, "Anas_crecca", where = which(duck_tree$tip.label=="Anas_flavirostris"), edge.length = 0.0032631134/4, position = 0.0032631134/4)
duck_tree_orig <- duck_tree3.4 GPA and PCA
To conduct a Generalized Procrustes Analysis we used the gpagen function of the geommorph package. We first summarized the information for all the specimens of a single bone using the GPA function in SlicerMorph to allow for easier upload into R work environment. For pca we used gm.prcomp with phylogeny from above.
3.4.1 Femur
# Load summarized data for all landmark files
logfile <- parser("gpa/Femur/analysis.log", forceLPS = T)
# Run GPA
femur.gpa <- gpagen(logfile$LM, print.progress = F)
femurshape <- femur.gpa$coords
dimnames(femurshape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(femurshape)[[3]])
femursize <- log10(femur.gpa$Csize)
names(femursize) <- dimnames(femurshape)[[3]]
# Subset phylogenetic tree to keep only taxa for which scans exist
duck_tree_femur <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Femur %in% "X"])
dimnames(femur.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(femur.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
femur_pca <- gm.prcomp(femur.gpa$coords, phy = duck_tree_femur)
# pc1 (25.87%), pc2 (19.36%), pc3 (7.10%), pc4 (5.34%), pc5 (5.01%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(femur_pca$x[,1]), PC1=femur_pca$x[,1], PC2=femur_pca$x[,2], PC3=femur_pca$x[,3])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Femur")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
femurshape <- femurshape[,,duck_tree_femur$tip.label]
femursize <- femursize[duck_tree_femur$tip.label]
femurbm <- duckbm[duck_tree_femur$tip.label]
femflight <- flight[duck_tree_femur$tip.label]
femdiving <- diving[duck_tree_femur$tip.label]
femallom1 <- procD.pgls(femurshape ~ femurbm, phy = duck_tree_femur, SS.type = "II", iter = 999)
summary(femallom1)
femallom2 <- procD.pgls(femurshape ~ femursize, phy = duck_tree_femur,SS.type = "II", iter = 999)
summary(femallom2)
femallom3 <- procD.pgls(femurshape ~ femurbm + femursize, phy = duck_tree_femur, SS.type = "II", iter = 999)
summary(femallom3)
femur.mod <- procD.pgls(femurshape ~ femurbm + femursize + femdiving, phy = duck_tree_femur, SS.type = "II", iter = 999)
summary(femur.mod)
femur.mod <- procD.pgls(femurshape ~ femurbm + femursize + femdiving + femflight, phy = duck_tree_femur, SS.type = "II", iter = 999)
summary(femur.mod)
#convergence
### find convergence test for 3D data
femdivelist <- names(femdiving[femdiving == "Diving"])
femdivconv <- convrat(duck_tree_femur, femur_pca$x[,1:5], femdivelist)
femdivconvsig <- convratsig(duck_tree_femur, femur_pca$x[,1:5], femdivelist, 10)
femdivconvsig
## convergence in flightless ducks
femflightlesslist <- names(femflight[femflight == "0"])
femflightlessconv <- convrat(femurtree, femur_pca$x[,1:5], femflightlesslist)
femflightlessconv
femflightconvsig <- convratsig(femurtree, femur_pca$x[,1:5], femflightlesslist, 10)
femflightconvsig
## Convergence test with PC axes
femconv <- conv.map(femdivelist, scores = femur_pca$x, pcs = femur_pca$rotation, mshape = femurmean, focal = "Somateria_mollissima" )
femconvtest <- convnum(duck_tree_femur, femur_pca$x[,1:5], femdivelist)We also generated warped meshes to visualize bone structure relevant to the different PC axes.
ref <- mshape(femur.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_femur.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_femur_new.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "femur_ref")
mesh2ply(g2, "femur_pc1_min")
mesh2ply(g3, "femur_pc1_max")
mesh2ply(g4, "femur_pc2_min")
mesh2ply(g5, "femur_pc2_max")3.4.2 Tibiotarsus
logfile <- parser("gpa/Tibiotarsus/analysis.log", forceLPS = T)
tibiotarsus.gpa <- gpagen(logfile$LM, print.progress = F)
tibioshape <- tibiotarsus.gpa$coords
dimnames(tibioshape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(tibioshape)[[3]])
tibiosize <- log10(tibiotarsus.gpa$Csize)
names(tibiosize) <- dimnames(tibioshape)[[3]]
duck_tree_tibio <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Tibiotarsus %in% "X"])
dimnames(tibiotarsus.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(tibiotarsus.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
tibio_pca <- gm.prcomp(tibiotarsus.gpa$coords, phy = duck_tree_tibio)
# pc1 (29.38%), pc2 (11.31%), pc3 (10.16%), pc4 (7.18%), pc5 (4.51%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(tibio_pca$x[,1]), PC1=tibio_pca$x[,1], PC2=tibio_pca$x[,2], PC3=tibio_pca$x[,3], PC4=tibio_pca$x[,4], PC5=tibio_pca$x[,5])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Tibiotarsus")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC3)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
tibioshape <- tibioshape[,,duck_tree_tibio$tip.label]
tibiosize <- tibiosize[duck_tree_tibio$tip.label]
tibiobm <- duckbm[duck_tree_tibio$tip.label]
tibioflight <- flight[duck_tree_tibio$tip.label]
tibiodiving <- diving[duck_tree_tibio$tip.label]
tibioallom1 <- procD.pgls(tibioshape ~ tibiobm, phy = duck_tree_tibio, SS.type = "II", iter = 999)
summary(tibioallom1)
tibioallom2 <- procD.pgls(tibioshape ~ tibiosize, phy = duck_tree_tibio,SS.type = "II", iter = 999)
summary(tibioallom2)
tibioallom3 <- procD.pgls(tibioshape ~ tibiobm + tibiosize, phy = duck_tree_tibio, SS.type = "II", iter = 999)
summary(tibioallom3)
tibio.mod <- procD.pgls(tibioshape ~ tibiobm + tibiosize + tibiodiving, phy = duck_tree_tibio, SS.type = "II", iter = 999)
summary(tibio.mod)
tibio.mod <- procD.pgls(tibioshape ~ tibiobm + tibiosize + tibiodiving + tibioflight, phy = duck_tree_tibio, SS.type = "II", iter = 999)
summary(tibio.mod)
#convergence
### find convergence test for 3D data
tibiodivelist <- names(tibiodiving[tibiodiving == "Diving"])
tibiodivconv <- convrat(duck_tree_tibio, tibio_pca$x[,1:5], tibiodivelist)
tibiodivconvsig <- convratsig(duck_tree_tibio, tibio_pca$x[,1:5], tibiodivelist, 10)
tibiodivconvsig
## convergence in flightless ducks
tibioflightlesslist <- names(tibioflight[tibioflight == "0"])
tibioflightlessconv <- convrat(tibiotree, tibio_pca$x[,1:5], tibioflightlesslist)
tibioflightlessconv
tibioflightconvsig <- convratsig(tibiotree, tibio_pca$x[,1:5], tibioflightlesslist, 10)
tibioflightconvsig
## Convergence test with PC axes
tibioconv <- conv.map(tibiodivelist, scores = tibio_pca$x, pcs = tibio_pca$rotation, mshape = tibiomean, focal = "Somateria_mollissima" )
tibioconvtest <- convnum(duck_tree_tibio, tibio_pca$x[,1:5], tibiodivelist)ref <- mshape(tibiotarsus.gpa$coords)
ref_obj <- read.ply("Anas_zonorhyncha_USNM633481_tibiotarsus.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_zonorhyncha_USNM633481_tibiotarsus.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp3$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp3$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "tibiotarsus_ref")
mesh2ply(g2, "tibiotarsus_pc1_min")
mesh2ply(g3, "tibiotarsus_pc1_max")
mesh2ply(g4, "tibiotarsus_pc3_min")
mesh2ply(g5, "tibiotarsus_pc3_max")3.4.3 Tarsometatarsus
logfile <- parser("gpa/Tarsometatarsus/analysis.log", forceLPS = T)
tarsometatarsus.gpa <- gpagen(logfile$LM, print.progress = F)
tarsoshape <- tarsometatarsus.gpa$coords
dimnames(tarsoshape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(tarsoshape)[[3]])
tarsosize <- log10(tarsometatarsus.gpa$Csize)
names(tarsosize) <- dimnames(tarsoshape)[[3]]
duck_tree_tarso <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Tarsometatarsus %in% "X"])
dimnames(tarsometatarsus.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(tarsometatarsus.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
tarso_pca <- gm.prcomp(tarsometatarsus.gpa$coords, phy = duck_tree_tarso)
# pc1 (40.49%), pc2 (10.91%), pc3 (7.74%), pc4 (5.44%), pc5 (4.45%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(tarso_pca$x[,1]), PC1=tarso_pca$x[,1], PC2=tarso_pca$x[,2], PC3=tarso_pca$x[,3], PC4=tarso_pca$x[,4], PC5=tarso_pca$x[,5])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Tarsometatarsus")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
tarsoshape <- tarsoshape[,,duck_tree_tarso$tip.label]
tarsosize <- tarsosize[duck_tree_tarso$tip.label]
tarsobm <- duckbm[duck_tree_tarso$tip.label]
tarsoflight <- flight[duck_tree_tarso$tip.label]
tarsodiving <- diving[duck_tree_tarso$tip.label]
tarsoallom1 <- procD.pgls(tarsoshape ~ tarsobm, phy = duck_tree_tarso, SS.type = "II", iter = 999)
summary(tarsoallom1)
tarsoallom2 <- procD.pgls(tarsoshape ~ tarsosize, phy = duck_tree_tarso,SS.type = "II", iter = 999)
summary(tarsoallom2)
tarsoallom3 <- procD.pgls(tarsoshape ~ tarsobm + tarsosize, phy = duck_tree_tarso, SS.type = "II", iter = 999)
summary(tarsoallom3)
tarso.mod <- procD.pgls(tarsoshape ~ tarsobm + tarsosize + tarsodiving, phy = duck_tree_tarso, SS.type = "II", iter = 999)
summary(tarso.mod)
tarso.mod <- procD.pgls(tarsoshape ~ tarsobm + tarsosize + tarsodiving + tarsoflight, phy = duck_tree_tarso, SS.type = "II", iter = 999)
summary(tarso.mod)
#convergence
### find convergence test for 3D data
tarsodivelist <- names(tarsodiving[tarsodiving == "Diving"])
tarsodivconv <- convrat(duck_tree_tarso, tarso_pca$x[,1:5], tarsodivelist)
tarsodivconvsig <- convratsig(duck_tree_tarso, tarso_pca$x[,1:5], tarsodivelist, 10)
tarsodivconvsig
## convergence in flightless ducks
tarsoflightlesslist <- names(tarsoflight[tarsoflight == "0"])
tarsoflightlessconv <- convrat(tarsotree, tarso_pca$x[,1:5], tarsoflightlesslist)
tarsoflightlessconv
tarsoflightconvsig <- convratsig(tarsotree, tarso_pca$x[,1:5], tarsoflightlesslist, 10)
tarsoflightconvsig
## Convergence test with PC axes
tarsoconv <- conv.map(tarsodivelist, scores = tarso_pca$x, pcs = tarso_pca$rotation, mshape = tarsomean, focal = "Somateria_mollissima" )
tarsoconvtest <- convnum(duck_tree_tarso, tarso_pca$x[,1:5], tarsodivelist)ref <- mshape(tarsometatarsus.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_tarsometatarsus.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_tarsometatarsus.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "tarsometatarsus_ref")
mesh2ply(g2, "tarsometatarsus_pc1_min")
mesh2ply(g3, "tarsometatarsus_pc1_max")
mesh2ply(g4, "tarsometatarsus_pc2_min")
mesh2ply(g5, "tarsometatarsus_pc2_max")3.4.4 Humerus
logfile <- parser("gpa/Humerus/analysis.log", forceLPS = T)
humerus.gpa <- gpagen(logfile$LM, print.progress = F)
humershape <- humerus.gpa$coords
dimnames(humershape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(humershape)[[3]])
humersize <- log10(humerus.gpa$Csize)
names(humersize) <- dimnames(humershape)[[3]]
duck_tree_humer <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Humerus %in% "X"])
duck_tree_humer <- root(duck_tree_humer, outgroup = c("Anhima_cornuta", "Chauna_chavaria"))
dimnames(humerus.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(humerus.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
humer_pca <- gm.prcomp(humerus.gpa$coords, phy = duck_tree_humer)
# pc1 (24.13%), pc2 (12.16%), pc3 (9.49%), pc4 (8.66%), pc5 (6.59%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(humer_pca$x[,1]), PC1=humer_pca$x[,1], PC2=humer_pca$x[,2], PC3=humer_pca$x[,3])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Humerus")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
humershape <- humershape[,,duck_tree_humer$tip.label]
humersize <- humersize[duck_tree_humer$tip.label]
humerbm <- duckbm[duck_tree_humer$tip.label]
humerflight <- flight[duck_tree_humer$tip.label]
humerdiving <- diving[duck_tree_humer$tip.label]
humerallom1 <- procD.pgls(humershape ~ humerbm, phy = duck_tree_humer, SS.type = "II", iter = 999)
summary(humerallom1)
humerallom2 <- procD.pgls(humershape ~ humersize, phy = duck_tree_humer,SS.type = "II", iter = 999)
summary(humerallom2)
humerallom3 <- procD.pgls(humershape ~ humerbm + humersize, phy = duck_tree_humer, SS.type = "II", iter = 999)
summary(humerallom3)
humer.mod <- procD.pgls(humershape ~ humerbm + humersize + humerdiving, phy = duck_tree_humer, SS.type = "II", iter = 999)
summary(humer.mod)
humer.mod <- procD.pgls(humershape ~ humerbm + humersize + humerdiving + humerflight, phy = duck_tree_humer, SS.type = "II", iter = 999)
summary(humer.mod)
#convergence
### find convergence test for 3D data
humerdivelist <- names(humerdiving[humerdiving == "Diving"])
humerdivconv <- convrat(duck_tree_humer, humer_pca$x[,1:5], humerdivelist)
humerdivconvsig <- convratsig(duck_tree_humer, humer_pca$x[,1:5], humerdivelist, 10)
humerdivconvsig
## convergence in flightless ducks
humerflightlesslist <- names(humerflight[humerflight == "0"])
humerflightlessconv <- convrat(humertree, humer_pca$x[,1:5], humerflightlesslist)
humerflightlessconv
humerflightconvsig <- convratsig(humertree, humer_pca$x[,1:5], humerflightlesslist, 10)
humerflightconvsig
## Convergence test with PC axes
humerconv <- conv.map(humerdivelist, scores = humer_pca$x, pcs = humer_pca$rotation, mshape = humermean, focal = "Somateria_mollissima" )
humerconvtest <- convnum(duck_tree_humer, humer_pca$x[,1:5], humerdivelist)ref <- mshape(humerus.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_humerus.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_humerus_new.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "humerus_ref")
mesh2ply(g2, "humerus_pc1_min")
mesh2ply(g3, "humerus_pc1_max")
mesh2ply(g4, "humerus_pc2_min")
mesh2ply(g5, "humerus_pc2_max")3.4.5 Radius
logfile <- parser("gpa/Radius/analysis.log", forceLPS = T)
radius.gpa <- gpagen(logfile$LM, print.progress = F)
radiishape <- radius.gpa$coords
dimnames(radiishape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(radiishape)[[3]])
radiisize <- log10(radius.gpa$Csize)
names(radiisize) <- dimnames(radiishape)[[3]]
duck_tree_radii <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Radius %in% "X"])
dimnames(radius.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(radius.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
radii_pca <- gm.prcomp(radius.gpa$coords, phy = duck_tree_radii)
# pc1 (28.75%), pc2 (18.06%), pc3 (10.83%), pc4 (9.17%), pc5 (5.86%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(radii_pca$x[,1]), PC1=radii_pca$x[,1], PC2=radii_pca$x[,2], PC3=radii_pca$x[,3], PC4=radii_pca$x[,4], PC5=radii_pca$x[,5])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Radius")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
radiishape <- radiishape[,,duck_tree_radii$tip.label]
radiisize <- radiisize[duck_tree_radii$tip.label]
radiibm <- duckbm[duck_tree_radii$tip.label]
radiiflight <- flight[duck_tree_radii$tip.label]
radiidiving <- diving[duck_tree_radii$tip.label]
radiiallom1 <- procD.pgls(radiishape ~ radiibm, phy = duck_tree_radii, SS.type = "II", iter = 999)
summary(radiiallom1)
radiiallom2 <- procD.pgls(radiishape ~ radiisize, phy = duck_tree_radii,SS.type = "II", iter = 999)
summary(radiiallom2)
radiiallom3 <- procD.pgls(radiishape ~ radiibm + radiisize, phy = duck_tree_radii, SS.type = "II", iter = 999)
summary(radiiallom3)
radii.mod <- procD.pgls(radiishape ~ radiibm + radiisize + radiidiving, phy = duck_tree_radii, SS.type = "II", iter = 999)
summary(radii.mod)
radii.mod <- procD.pgls(radiishape ~ radiibm + radiisize + radiidiving + radiiflight, phy = duck_tree_radii, SS.type = "II", iter = 999)
summary(radii.mod)
#convergence
### find convergence test for 3D data
radiidivelist <- names(radiidiving[radiidiving == "Diving"])
radiidivconv <- convrat(duck_tree_radii, radii_pca$x[,1:5], radiidivelist)
radiidivconvsig <- convratsig(duck_tree_radii, radii_pca$x[,1:5], radiidivelist, 10)
radiidivconvsig
## convergence in flightless ducks
radiiflightlesslist <- names(radiiflight[radiiflight == "0"])
radiiflightlessconv <- convrat(radiitree, radii_pca$x[,1:5], radiiflightlesslist)
radiiflightlessconv
radiiflightconvsig <- convratsig(radiitree, radii_pca$x[,1:5], radiiflightlesslist, 10)
radiiflightconvsig
## Convergence test with PC axes
radiiconv <- conv.map(radiidivelist, scores = radii_pca$x, pcs = radii_pca$rotation, mshape = radiimean, focal = "Somateria_mollissima" )
radiiconvtest <- convnum(duck_tree_radii, radii_pca$x[,1:5], radiidivelist)ref <- mshape(radius.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_radius.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_radius.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "radius_ref")
mesh2ply(g2, "radius_pc1_min")
mesh2ply(g3, "radius_pc1_max")
mesh2ply(g4, "radius_pc2_min")
mesh2ply(g5, "radius_pc2_max")3.4.6 Ulna
logfile <- parser("gpa/Ulna/analysis.log", forceLPS = T)
ulna.gpa <- gpagen(logfile$LM, print.progress = F)
ulnashape <- ulna.gpa$coords
dimnames(ulnashape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(ulnashape)[[3]])
ulnasize <- log10(ulna.gpa$Csize)
names(ulnasize) <- dimnames(ulnashape)[[3]]
duck_tree_ulna <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Ulna %in% "X"])
dimnames(ulna.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(ulna.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
ulna_pca <- gm.prcomp(ulna.gpa$coords, phy = duck_tree_ulna)
# pc1 (28.35%), pc2 (19.04%), pc3 (14.70%), pc4 (7.03%), pc5 (5.61%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(ulna_pca$x[,1]), PC1=ulna_pca$x[,1], PC2=ulna_pca$x[,2], PC3=ulna_pca$x[,3], PC4=ulna_pca$x[,4], PC5=ulna_pca$x[,5])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Ulna")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
ulnashape <- ulnashape[,,duck_tree_ulna$tip.label]
ulnasize <- ulnasize[duck_tree_ulna$tip.label]
ulnabm <- duckbm[duck_tree_ulna$tip.label]
ulnaflight <- flight[duck_tree_ulna$tip.label]
ulnadiving <- diving[duck_tree_ulna$tip.label]
ulnaallom1 <- procD.pgls(ulnashape ~ ulnabm, phy = duck_tree_ulna, SS.type = "II", iter = 999)
summary(ulnaallom1)
ulnaallom2 <- procD.pgls(ulnashape ~ ulnasize, phy = duck_tree_ulna,SS.type = "II", iter = 999)
summary(ulnaallom2)
ulnaallom3 <- procD.pgls(ulnashape ~ ulnabm + ulnasize, phy = duck_tree_ulna, SS.type = "II", iter = 999)
summary(ulnaallom3)
ulna.mod <- procD.pgls(ulnashape ~ ulnabm + ulnasize + ulnadiving, phy = duck_tree_ulna, SS.type = "II", iter = 999)
summary(ulna.mod)
ulna.mod <- procD.pgls(ulnashape ~ ulnabm + ulnasize + ulnadiving + ulnaflight, phy = duck_tree_ulna, SS.type = "II", iter = 999)
summary(ulna.mod)
#convergence
### find convergence test for 3D data
ulnadivelist <- names(ulnadiving[ulnadiving == "Diving"])
ulnadivconv <- convrat(duck_tree_ulna, ulna_pca$x[,1:5], ulnadivelist)
ulnadivconvsig <- convratsig(duck_tree_ulna, ulna_pca$x[,1:5], ulnadivelist, 10)
ulnadivconvsig
## convergence in flightless ducks
ulnaflightlesslist <- names(ulnaflight[ulnaflight == "0"])
ulnaflightlessconv <- convrat(ulnatree, ulna_pca$x[,1:5], ulnaflightlesslist)
ulnaflightlessconv
ulnaflightconvsig <- convratsig(ulnatree, ulna_pca$x[,1:5], ulnaflightlesslist, 10)
ulnaflightconvsig
## Convergence test with PC axes
ulnaconv <- conv.map(ulnadivelist, scores = ulna_pca$x, pcs = ulna_pca$rotation, mshape = ulnamean, focal = "Somateria_mollissima" )
ulnaconvtest <- convnum(duck_tree_ulna, ulna_pca$x[,1:5], ulnadivelist)ref <- mshape(ulna.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_ulna.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_ulna.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "ulna_ref")
mesh2ply(g2, "ulna_pc1_min")
mesh2ply(g3, "ulna_pc1_max")
mesh2ply(g4, "ulna_pc2_min")
mesh2ply(g5, "ulna_pc2_max")3.4.7 Carpometacarpus
logfile <- parser("gpa/Carpometacarpus/analysis.log", forceLPS = T)
carpometacarpus.gpa <- gpagen(logfile$LM, print.progress = F)
carposhape <- carpometacarpus.gpa$coords
dimnames(carposhape)[[3]] <- sub("^([^_]+_[^_]+).*", "\\1", dimnames(carposhape)[[3]])
carposize <- log10(carpometacarpus.gpa$Csize)
names(carposize) <- dimnames(carposhape)[[3]]
duck_tree_carpo <- keep.tip(duck_tree_orig, duck_data$Scientific_Name[duck_data$Carpometacarpus %in% "X"])
dimnames(carpometacarpus.gpa$coords)[[3]] <- unlist(lapply(str_split(dimnames(carpometacarpus.gpa$coords)[[3]], "_"), FUN = function(x){paste(x[1], x[2], sep="_")}))
carpo_pca <- gm.prcomp(carpometacarpus.gpa$coords, phy = duck_tree_carpo)
# pc1 (16.79%), pc2 (15.00%), pc3 (14.15%), pc4 (8.73%), pc5 (8.11%)
asp_ratio <- 1.618
duck_pca <- data.frame(species=names(carpo_pca$x[,1]), PC1=carpo_pca$x[,1], PC2=carpo_pca$x[,2], PC3=carpo_pca$x[,3], PC4=carpo_pca$x[,4], PC5=carpo_pca$x[,5])
duck_pca <- merge(duck_pca, duck_data[,c("Scientific_Name", "Images", "Carpometacarpus")], by.x = "species", by.y = "Scientific_Name", all.x = T, all.y = F)
g1 <- ggplot(duck_pca, aes(PC1, PC2)) +
geom_image(aes(image=Images), size=0.075, by = "width", asp = asp_ratio)+
theme_classic()+
theme(aspect.ratio = 1/asp_ratio)
g1
Then we check for convergence
carposhape <- carposhape[,,duck_tree_carpo$tip.label]
carposize <- carposize[duck_tree_carpo$tip.label]
carpobm <- duckbm[duck_tree_carpo$tip.label]
carpoflight <- flight[duck_tree_carpo$tip.label]
carpodiving <- diving[duck_tree_carpo$tip.label]
carpoallom1 <- procD.pgls(carposhape ~ carpobm, phy = duck_tree_carpo, SS.type = "II", iter = 999)
summary(carpoallom1)
carpoallom2 <- procD.pgls(carposhape ~ carposize, phy = duck_tree_carpo,SS.type = "II", iter = 999)
summary(carpoallom2)
carpoallom3 <- procD.pgls(carposhape ~ carpobm + carposize, phy = duck_tree_carpo, SS.type = "II", iter = 999)
summary(carpoallom3)
carpo.mod <- procD.pgls(carposhape ~ carpobm + carposize + carpodiving, phy = duck_tree_carpo, SS.type = "II", iter = 999)
summary(carpo.mod)
carpo.mod <- procD.pgls(carposhape ~ carpobm + carposize + carpodiving + carpoflight, phy = duck_tree_carpo, SS.type = "II", iter = 999)
summary(carpo.mod)
#convergence
### find convergence test for 3D data
carpodivelist <- names(carpodiving[carpodiving == "Diving"])
carpodivconv <- convrat(duck_tree_carpo, carpo_pca$x[,1:5], carpodivelist)
carpodivconvsig <- convratsig(duck_tree_carpo, carpo_pca$x[,1:5], carpodivelist, 10)
carpodivconvsig
## convergence in flightless ducks
carpoflightlesslist <- names(carpoflight[carpoflight == "0"])
carpoflightlessconv <- convrat(carpotree, carpo_pca$x[,1:5], carpoflightlesslist)
carpoflightlessconv
carpoflightconvsig <- convratsig(carpotree, carpo_pca$x[,1:5], carpoflightlesslist, 10)
carpoflightconvsig
## Convergence test with PC axes
carpoconv <- conv.map(carpodivelist, scores = carpo_pca$x, pcs = carpo_pca$rotation, mshape = carpomean, focal = "Somateria_mollissima" )
carpoconvtest <- convnum(duck_tree_carpo, carpo_pca$x[,1:5], carpodivelist)ref <- mshape(carpometacarpus.gpa$coords)
ref_obj <- read.ply("Anas_platyrhynchos_MCZ347156_carpometacarpus.ply", ShowSpecimen = F)
ref_lnd <- readland.fcsv("Anas_platyrhynchos_MCZ347156_carpometacarpus.fcsv")
warped_ref <- warpRefMesh(ref_obj, ref_lnd, ref, color = "ivory", centered = FALSE)
g2 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$min, mesh=warped_ref, method = "surface")
g3 <- plotRefToTarget(ref, pca$shapes$shapes.comp1$max, mesh=warped_ref, method = "surface")
g4 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$min, mesh=warped_ref, method = "surface")
g5 <- plotRefToTarget(ref, pca$shapes$shapes.comp2$max, mesh=warped_ref, method = "surface")
mesh2ply(warped_ref, "carpometacarpus_ref")
mesh2ply(g2, "carpometacarpus_pc1_min")
mesh2ply(g3, "carpometacarpus_pc1_max")
mesh2ply(g4, "carpometacarpus_pc2_min")
mesh2ply(g5, "carpometacarpus_pc2_max")4 PhyloAcc
We did both a binary version of PhyloAcc with diving vs non-diving as traits (not for flightlessness as we only had 3 species, 2 with genomes), and a continuous version using PhyloAcc-C.
4.1 Genomic Data
4.1.1 Cactus
We aligned 27 duck genomes (GenBank numbers in table) and also included the chicken genomes using progressive cactus. The tree used for the alignment is shown below.
(((((((((Aix_galericulata:0.0051473029,Aix_sponsa:0.0053483288):0.0017143373,(Cairina_moschata:0.0016771265,Asarcornis_scutulata:0.0016771265):0.0050313795):0.0030667213,((((Lophonetta_specularioides:0.0100951666,Tachyeres_brachypterus:0.0055556181):0.0028194957,((Anas_acuta:0.0042863335,(Anas_platyrhynchos:0.00124933445,Anas_zonorhyncha:0.00124933445):0.00374800335):0.001201751,(Mareca_sibilatrix:0.0037871541,Mareca_strepera:0.0037867301):0.0014935713):0.0011642544):0.0004165566,Spatula_discors:0.007856085):0.0023604916,Aythya_fuligula:0.009531404):0.0012462578):0.0009209619,(Melanitta_fusca:0.0064818964,(Clangula_hyemalis:0.0070627445,Somateria_mollissima:0.0073593856):0.0003851771):0.0014142496):0.0173506804,((((Anser_fabalis:0.0020284789,Anser_cygnoides:0.001826199):0.0008868183,Anser_indicus:0.0025060375):0.0036566335,Branta_canadensis:0.006571712):0.0029378261,Cygnus_olor:0.0091439003):0.0085572595):0.000654519,(((Heteronetta_atricapilla:0.0213298248,Stictonetta_naevosa:0.0132137813):0.0004721905,Oxyura_jamaicensis:0.0155460859):0.0065766967,Nettapus_auritus:0.0249469325):0.0034281317):0.0328791421,Anseranas_semipalmata:0.0308188616):0.01,Anhima_cornuta:0.0500106949):0.01,Gallus_gallus:0.0308188616);To run cactus:
#!/bin/bash
#SBATCH -c 1
#SBATCH -t 99-00:00
#SBATCH -p holy-smokes
#SBATCH --nodelist=holy7c0908
#SBATCH --mem=1000
#SBATCH -o ./logs/cac_%j.out
#SBATCH -e ./logs/cac_%j.err
cd /n/holyscratch01/informatics/sshakya/ducks/
source activate PhyloAcc
singularity exec --cleanenv cactus_v2.7.1-gpu.sif cactus-prepare duck_aln.txt --outDir snk_output --jobStore snk_temp --gpu
snakemake -p -s ~/GenomeAnnotation/genome_alignment2/cactus-gpu-snakemake/cactus_gpu.smk --configfile duck_aln.yaml --profile ~/GenomeAnnotation/genome_alignment2/cactus-gpu-snakemake/profiles/slurm_profile/4.1.2 Neutral Rates
To get neutral rates, we used the chicken annotations. See conserved_rates for steps to get neutral rates.
4.1.3 Conserved elements alignment
We used the coordinates from conserved elements produced in Shakya et al. for getting alignments of conserved elements. See elements_fasta for steps to get conserved elements.
4.2 Morphological input data
4.2.1 Femur
duck_pca_genome <- duck_pca[which(duck_pca$species %in% duck_data$Scientific_Name[!is.na(duck_data$Genomes)]),]
duck_tree_genome <- keep.tip(duck_tree, duck_pca_genome$species)
duck_pca_genome$species <- factor(duck_pca_genome$species, levels = duck_tree_genome$tip.label)
duck_pca_genome <- duck_pca_genome[order(duck_pca_genome$species),]
p <- ggtree(duck_tree_genome) +
xlim(NA,0.2)+
geom_tiplab(size=3)+
geom_facet(panel = "Femur_PC1", geom = geom_col, data = duck_pca_genome,
mapping = aes(x = PC1), orientation='y', width=0.6)+
geom_tiplab(aes(image = c(duck_pca_genome$Images, rep("NA", 24))), geom="image",
offset=0.1, align = T)
facet_widths(p, c(0.75, 0.25))4.2.2 Tarsometatarsus
duck_pca_genome <- duck_pca[which(duck_pca$species %in% duck_data$Scientific_Name[!is.na(duck_data$Genomes)]),]
duck_tree_genome <- keep.tip(duck_tree, duck_pca_genome$species)
duck_pca_genome$species <- factor(duck_pca_genome$species, levels = duck_tree_genome$tip.label)
duck_pca_genome <- duck_pca_genome[order(duck_pca_genome$species),]
p <- ggtree(duck_tree_genome) +
xlim(NA,0.2)+
geom_tiplab(size=3)+
geom_facet(panel = "tarsometatarsus_PC1", geom = geom_col, data = duck_pca_genome,
mapping = aes(x = PC1), orientation='y', width=0.6)+
geom_tiplab(aes(image = c(duck_pca_genome$Images, rep("NA", 22))), geom="image",
offset=0.1, align = T)
facet_widths(p, c(0.75, 0.25))4.2.3 Humerus
duck_pca_genome <- duck_pca[which(duck_pca$species %in% duck_data$Scientific_Name[!is.na(duck_data$Genomes)]),]
duck_tree_genome <- keep.tip(duck_tree, duck_pca_genome$species)
duck_pca_genome$species <- factor(duck_pca_genome$species, levels = duck_tree_genome$tip.label)
duck_pca_genome <- duck_pca_genome[order(duck_pca_genome$species),]
p <- ggtree(duck_tree_genome) +
xlim(NA,0.2)+
geom_tiplab(size=3)+
geom_facet(panel = "humerus_PC1", geom = geom_col, data = duck_pca_genome,
mapping = aes(x = PC1), orientation='y', width=0.6)+
geom_tiplab(aes(image = c(duck_pca_genome$Images, rep("NA", 26))), geom="image",
offset=0.1, align = T)
facet_widths(p, c(0.75, 0.25))4.2.4 Carpometacarpus
duck_pca_genome <- duck_pca[which(duck_pca$species %in% duck_data$Scientific_Name[!is.na(duck_data$Genomes)]),]
duck_tree_genome <- keep.tip(duck_tree, duck_pca_genome$species)
duck_pca_genome$species <- factor(duck_pca_genome$species, levels = duck_tree_genome$tip.label)
duck_pca_genome <- duck_pca_genome[order(duck_pca_genome$species),]
p <- ggtree(duck_tree_genome) +
xlim(NA,0.2)+
geom_tiplab(size=3)+
geom_facet(panel = "carpometacarpus_PC1", geom = geom_col, data = duck_pca_genome,
mapping = aes(x = PC1), orientation='y', width=0.6)+
geom_tiplab(aes(image = c(duck_pca_genome$Images, rep("NA", 19))), geom="image",
offset=0.1, align = T)
facet_widths(p, c(0.75, 0.25))4.3 Binary
We ran PhyloAcc with the species tree model as such. We ran each separately for the three different cases-macrochromosomes, microchromosomes, and Z chromosome.
First we modified the input fasta files to account for missing taxa:
# 28 species in desired order
species_list <- c(
"Aix_galericulata", "Aix_sponsa", "Cairina_moschata",
"Asarcornis_scutulata", "Lophonetta_specularioides", "Tachyeres_brachypterus",
"Anas_acuta", "Anas_platyrhynchos", "Anas_zonorhyncha",
"Mareca_sibilatrix", "Mareca_strepera", "Spatula_discors",
"Aythya_fuligula", "Melanitta_fusca", "Clangula_hyemalis",
"Somateria_mollissima", "Anser_fabalis", "Anser_cygnoides",
"Anser_indicus", "Branta_canadensis", "Cygnus_olor",
"Heteronetta_atricapilla", "Stictonetta_naevosa", "Oxyura_jamaicensis",
"Nettapus_auritus", "Anseranas_semipalmata", "Anhima_cornuta",
"Gallus_gallus"
)
# Function to process one fasta file
add_missing_species <- function(infile, outfile = infile) {
aln <- readDNAStringSet(infile)
# Alignment length (assumes all sequences in file same length)
aln_widths <- width(aln)
if (length(unique(aln_widths)) != 1L) {
stop("Sequences in ", infile, " do not all have the same length.")
}
aln_len <- unique(aln_widths)
present <- names(aln)
missing <- setdiff(species_list, present)
if (length(missing) > 0) {
# Create N-only sequences for missing species
n_string <- paste(rep("N", aln_len), collapse = "")
n_seqs <- DNAStringSet(rep(n_string, length(missing)))
names(n_seqs) <- missing
# Combine and reorder according to species_list
aln_all <- c(aln, n_seqs)
aln_all <- aln_all[species_list[species_list %in% names(aln_all)]]
} else {
aln_all <- aln
}
writeXStringSet(aln_all, filepath = outfile, format = "fasta")
}
## Process all fasta files in a directory
## Set your directory path here:
getwd()
fasta_dir <- "./ducks/macros"
#fasta_dir <- "./ducks/zchrom"
#fasta_dir <- "./ducks/micros"
fasta_files <- list.files(fasta_dir, pattern = "\\.(fa|fasta)$", full.names = TRUE)
fasta_files[[1]]
f<-1
for (f in 1:length(fasta_files)) {
message("Processing: ", f)
add_missing_species(fasta_files[[f]])
}PhyloAcc set up
# set up batch runs of phyloacc
FAS=./ducks/macros
MOD=ducks_4d_macro_non-conserved.mod
OUT=./ducks/phyloacc_macro_out
phyloacc -d $FAS -m $MOD -t "Aythya_fuligula;Clangula_hyemalis;Heteronetta_atricapilla;Melanitta_fusca;Oxyura_jamaicensis;Somateria_mollissima" \
-g "Gallus_gallus" -r st -burnin 1000 -mcmc 10000 -n 1 -batch 10 -j 20 -part "edwards" -o $OUT --overwrite Running PhyloAcc
CFGLIST=list_paths_duck_macros_cfgs.txt #generated by phyloacc command
OFFSET=0
LINENUM=$((SLURM_ARRAY_TASK_ID + OFFSET))
CFG=$(awk -v n="$LINENUM" 'FNR == n { print $1; exit }' "$CFGLIST")
LOG=$(awk -v n="$LINENUM" 'FNR == n { print $2; exit }' "$CFGLIST")
OUT=./ducks/phyloacc_macro_out/phyloacc-job-files/phyloacc-output
PhyloAcc-ST $CFG &> $OUT/$LOGPostprocessing
TAR=./ducks/phyloacc_macro_out
OUT=./ducks/post_macro_results
phyloacc_post.py -i $TAR -o $OUT --overwriteThe same approach is repeated for the microchromosome and the Z chromosome using the respective mod file.
4.3.1 Interpreting results
macro_results <- fread("post_macro_results/results/elem_lik.txt.gz")
micro_results <- fread("post_micros_results/results/elem_lik.txt.gz")
Z_results <- fread("post_Z_results/results/elem_lik.txt.gz")
all_results <- rbind(macro_results, micro_results, Z_results)
all_results$color <- "gray50"
all_results$color[which(all_results$logbf1 > 10 & all_results$logbf2 > 1)] <- "darkorange"
table(all_results$color)##
## darkorange gray50
## 951 812467# get all elements that are accelerated and count the number of clades it is accelerated in
clade_map <- c(
"Aythya_fuligula" = "Aythya",
"Clangula_hyemalis" = "Mergini",
"Heteronetta_atricapilla" = "Heteronetta",
"Melanitta_fusca" = "Mergini",
"Oxyura_jamaicensis" = "Oxyura",
"Somateria_mollissima" = "Mergini"
)
accel_elements <- subset(all_results, color %in% c("darkorange"))
# count the number of clades
accel_elements <- accel_elements %>%
mutate(
clades = map(accel.lineages.m1, function(x) {
taxa <- str_split(x, "-", simplify = FALSE)[[1]] |>
str_split(",", simplify = FALSE) |>
unlist()
sort(unique(unname(clade_map[taxa])))
}),
n_clades = map_int(clades, length)
)
table(accel_elements$n_clades)##
## 0 1 2 3 4
## 2 580 227 115 27
accel_elements <- subset(accel_elements, n_clades >= 2)
phy2 <- ggplot(accel_elements, aes(x = clades)) +
geom_bar() +
scale_x_upset(order_by = "freq") +
ggtitle("B)") +
labs(
x = "Accelerated clades",
y = "Number of accelerated elements"
) +
theme_classic()
phy2
all_data <- fread("conserved_elements_all.txt") #from Shakya et al. 2025
final_data_filtered <- subset(all_data, id %in% all_results$original.id)
accel_elements <- subset(accel_elements, n_clades >= 3)
final_data_accelerated <- subset(final_data_filtered, id %in% accel_elements$original.id)
fwrite(final_data_filtered[,1:4], "loci_all.txt", sep = "\t", col.names = F, row.names = F, quote = F)
fwrite(final_data_accelerated[,1:4], "loci_accl.txt", sep = "\t", col.names = F, row.names = F, quote = F)# after getting closest genes
closest_all <- fread("loci_all_genes.txt", header = F)
closest_all$gene <- str_extract(closest_all$V13, "gene=[a-zA-Z0-9_]+")
closest_all$gene <- gsub("^[A-Za-z_]+=", "", closest_all$gene)
closest_accl <- subset(closest_all, V4 %in% accel_elements$original.id)
#Organism 9031 is Gallus gallus
#Dataset GO:0008150 is GO terms in biological processes group
accl_enrich <- rba_panther_enrich(genes = unique(closest_accl$gene),
organism = 9031,
annot_dataset = "GO:0008150")
accl_enrich_table <- accl_enrich$resultGO enrichment by sampling
nums <- 142 #number of accelerated CNEEs
enrich_out <- data.frame(term.id=character())
gene_out <- data.frame(term.id=character())
unique_genes <- c()
for (i in seq(1,1000)){
rand_ces <- sample(all_data$id, nums)
closest <- subset(closest_all, V4 %in% rand_ces)
unique_genes <- c(unique_genes, length(unique(closest$gene)))
enrich <- rba_panther_enrich(genes = unique(closest$gene),
organism = 9031,
annot_dataset = "GO:0008150")
try({
print(paste("Run", i, sep=" "))
colnames(enrich$result)[1] <- paste(colnames(enrich$result)[1], i, sep = "_")
colnames(enrich$result)[2] <- paste(colnames(enrich$result)[2], i, sep = "_")
enrich_out <- merge(enrich_out, enrich$result[,c(paste("fold_enrichment",i, sep = "_"), "term.id")], by="term.id", all=T)
gene_out <- merge(gene_out, enrich$result[,c(paste("number_in_list",i, sep = "_"), "term.id")], by="term.id", all=T)
}, silent = T)
}
fwrite(enrich_out, paste("enrich_out_", nums, ".txt", sep=""), sep = "\t", col.names = T, row.names = F, quote = F)
fwrite(gene_out, paste("genes_out_", nums, ".txt", sep=""), sep = "\t", col.names = T, row.names = F, quote = F)gene_out <- fread("genes_out_142.txt", header=TRUE)
get_pgo <- function(GO, value, distr){
if (GO %in% distr$term.id){
vals <- as.numeric(distr[distr$term.id == GO, -1])
vals <- vals[!is.na(vals)]
return((sum(vals >= value)+1)/(length(vals)+1))
} else {
return(0)
}
}
accl_enrich_table$pgo <- mapply(
get_pgo,
accl_enrich_table$term.id,
accl_enrich_table$number_in_list,
MoreArgs=list(distr=gene_out)
)
gene_out$median <- MatrixGenerics::rowMedians(as.matrix(gene_out[,2:ncol(gene_out)]), useNames = F)
accl_enrich_table <- merge(
accl_enrich_table,
gene_out[,c("term.id","median")],
by="term.id",
all.x=TRUE
)
tot <- median(unique_genes) # 125, 322
accl_enrich_table$median[accl_enrich_table$median == 0] <- 1e-5
accl_enrich_table$enrich <-
(accl_enrich_table$number_in_list / 119) /
(accl_enrich_table$median / 125)
all_sigs <- subset(
accl_enrich_table,
pgo < 0.05 &
number_in_reference >= 5 &
median >= 1
)
write.table(all_sigs, "GO_analysis_142.txt", row.names = F, col.names = T, sep = "\t", quote = F)4.4 Continuous
We ran PhyloAcc-C model for PC1 of each of the four selected bones. To run PhyloAcc-C, we used the following code:
library(coda)
library(parallel)
suppressMessages(library(sm))
library(argparse)
library(stringi)
library(ape)
ELL_do_mcmc = function(row_i) {
mclapply(1:args$runs, function(i) {
pac3::run_model(tree=tree, alignment=(all_alignments[[row_i]]),
trait_tbl=trait_tbl,
transition_matrix=Q, stationary_distro=PI, n_iter=args$iter,
inner_iter=500,
Phi_prior=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
prior_r2=c(4,0.1),
prior_r3=c(10,0.4))
}, mc.cores=1) # N.B. set to 1 = lapply, so do parallelism at element level
}
ELL_make_mcmclist = function(mcmc_outs, n_burnin=n_burnin, n_thin=1) {
cs <- lapply(mcmc_outs, function(mcmc_out) {
the_trace = with(mcmc_out$mcmc_trace,
data.frame(lv =v_trace[,1],
lb2=v_trace[,2],
lb3=v_trace[,3]))
as.mcmc(the_trace, start=(n_burnin+1), thin=n_thin)
})
as.mcmc.list(cs)
}
do_ith_element = function(i) {
print(paste("Working on element ",i," of ",length(all_alignments)," of batch ",ceiling(start_index/args$batch_size) , sep=""))
the_mcmc_outs = ELL_do_mcmc(i)
the_trace = do.call(rbind, lapply(1:args$runs, function(i) the_mcmc_outs[[i]]$mcmc_trace$v_trace[-(1:args$burnin),])) # needs the rbind otherwise don't get correct quantiles
the_mcmc_list = ELL_make_mcmclist(the_mcmc_outs)
rhat = gelman.diag(the_mcmc_list)$mpsrf
the_density = sm.density(the_trace[,2:3], model="Normal", eval.grid=F, eval.points=matrix(data=c(0,0), byrow=T, nrow=1, ncol=2), display="none")
the_density2= sm.density(the_trace[,3]-the_trace[,2], model="Normal", eval.grid=F, eval.points=0, display="none")
qlv = as.numeric(quantile(the_trace[,1], prob=c(.1,.5,.9)))
qlb2 = as.numeric(quantile(the_trace[,2], prob=c(.1,.5,.9)))
qlb3 = as.numeric(quantile(the_trace[,3], prob=c(.1,.5,.9)))
qrat = as.numeric(quantile(the_trace[,3]- the_trace[,2], prob=c(.1,.5,.9)))
list( lv.10=qlv [1], lv.50=qlv [2], lv.90=qlv [3],
lb2.10=qlb2[1], lb2.50=qlb2[2], lb2.90=qlb2[3],
lb3.10=qlb3[1], lb3.50=qlb3[2], lb3.90=qlb3[3],
rat.10=qrat[1], rat.50=qrat[2], rat.90=qrat[3],
rhat=as.numeric(rhat),
dens0=as.numeric(the_density$estimate),
densr=as.numeric(the_density2$estimate))
}
do_these_elements = function(indices) {
mclapply(indices, function(i) {
the_res = do_ith_element(i)
the_res$NAME = names(all_alignments[i])
the_res
}, mc.cores=args$cores)
}
# Create Argument Parser
parser <- ArgumentParser(description = "Run PhyloAcc-c")
# Add arguments
parser$add_argument("--mod_file", type="character", help = "Path to mod file with tree and neutral rates")
parser$add_argument("--output", type="character", help = "Output file name")
parser$add_argument("--output_folder", default="./", type="character", help = "Output folder name (default = ./)")
parser$add_argument("--fasta", type="character", help = "Fasta folder")
parser$add_argument("--fasta_list", type="character", help = "Optional list to only include from fasta")
parser$add_argument("--traits", type="character", help = "Tab-delimited file with species and trait")
parser$add_argument("--batch_size", type="integer", default=100, help = "How many elements to run and write at one time")
parser$add_argument("--runs", type="integer", default=3, help = "Number of runs per element (default=3)")
parser$add_argument("--cores", type="integer", default=1, help = "How many cores to run in parallel")
parser$add_argument("--iter", type="integer", default=10000, help = "Number of iterations per element (default=10000)")
parser$add_argument("--burnin", type="integer", default=1000, help = "Number of iterations of burnin to discard (default=1000)")
parser$add_argument("--max_missing", type="double", default=0.5, help = "Fraction of taxa need to be present in alignment (default=0.5)")
# Parse arguments
args <- parser$parse_args()
print(args)
# Parse mod file
lines <- readLines(args$mod_file)
for (line_index in seq_along(lines)) {
# Split the line by ":"
parts <- stri_split_fixed(lines[line_index], ":", n=2)[[1]]
# Extract variable name and IDs
variable <- trimws(parts[1]) # Trim whitespace
if (variable == "BACKGROUND") {
PI <- as.numeric(unlist(strsplit(trimws(parts[2]), "\\s+")))
}
if (variable == "RATE_MAT"){
mat_vals <- strsplit(lines[(line_index+1):(line_index+4)], "\\s+")
Q <- do.call(rbind, mat_vals)
Q <- Q[,-1]
Q <- matrix(as.numeric(Q), nrow=4, byrow=T)
line_index <- line_index + 4
}
if (variable == "TREE"){
tree = read.tree(text = trimws(parts[2]))
tree$node.label <- paste("Anc", 0:Nnode(tree), sep="")
}
}
# Get trait data
trait_tbl <- read.table(args$traits)
colnames(trait_tbl) <- c("NAME", "Y")
tree <- keep.tip(tree, trait_tbl$NAME)
trait_tbl$NAME <- factor(trait_tbl$NAME, levels = tree$tip.label)
trait_tbl <- trait_tbl[order(trait_tbl$NAME),]
row.names(trait_tbl) <- 1:nrow(trait_tbl)
alignment_loc <- args$fasta
fasta_files <- list.files(alignment_loc, full.names = T, include.dirs = F)
fasta_basenames <- tools::file_path_sans_ext(basename(fasta_files))
if (!is.null(args$fasta_list)){
keep_fasta <- read.table(args$fasta_list)
fasta_files <- fasta_files[which(fasta_basenames %in% keep_fasta$V1)]
fasta_basenames <- fasta_basenames[which(fasta_basenames %in% keep_fasta$V1)]
}
#some basic housekeeping
start_run <- T
if (!file.exists(args$output_folder)) {
# If it doesn't exist, create the output folder
dir.create(args$output_folder)
}
# Loop in batches of batch size
for (start_index in seq(1, length(fasta_files), by = args$batch_size)){
all_alignments <- list()
for (each_file in start_index:(start_index + args$batch_size - 1)){
element <- fasta_basenames[each_file]
alignment <- read.dna(fasta_files[each_file], format = "fasta", as.matrix = T, as.character = T)
if ((nrow(alignment) > ((length(tree$tip.label))*args$max_missing))) {
alignment <- alignment[rownames(alignment) %in% trait_tbl$NAME,]
alignment[alignment=="*"] <- "-"
alignment[alignment=="n"] <- "-"
all_same <- all(apply(alignment, 2, function(x) length(unique(x))) == 1)
missing_taxa <- setdiff(tree$tip.label, rownames(alignment))
if (length(missing_taxa > 0)){
for (miss in missing_taxa){
missing_row <- matrix(rep("-", ncol(alignment)), nrow=1)
rownames(missing_row) <- miss
alignment <- rbind(alignment, missing_row)
}
}
if (!all_same){
all_alignments <- append(all_alignments, setNames(list(alignment), element))
}
}
}
if(length(all_alignments) > 0){
system.time(yyy <- do_these_elements(1:length(all_alignments)))
tbl_to_write = do.call(rbind, yyy)
out_fn = paste(args$output_folder, "/", args$output, ".tsv", sep="")
if (file.exists(out_fn) & start_run){
write.table(tbl_to_write, out_fn, sep="\t", row.names=F, quote=F, col.names=T)
start_run <- F
}
write.table(tbl_to_write, out_fn, sep="\t", row.names=F, quote=F, col.names=F, append=T)
}
}This script can then be run as such:
scaffold="NC_052572.1"
Rscript run_phyloaccc.r --mod_file ducks_4d_Z_non-conserved.mod --output ${scaffold} --output_folder phylo_out_carpometacarpus --trait duck_carpometacarpus_pc1.txt --fasta fasta/fasta_Z/${scaffold}/ --fasta_list bed_ranges/${scaffold}.txt --batch_size 200 --cores 40 --iter 10000 --burnin 10004.4.1 Parsing outputs
4.4.1.1 Femur
4.4.1.2 Humerus
4.4.1.3 Tarsometatarsus
4.4.1.4 Carpometacarpus
4.4.2 Individual runs
For highly significant loci, we can run phyloacc independently to generate more plots showing parameter mixing and trees with rates drawn on it. An Rdata file can also generated if you want to load that run. The code you can use is here:
library(coda)
library(parallel)
suppressMessages(library(sm))
library(argparse)
library(stringi)
library(ape)
library(pac3)
# Create Argument Parser
parser <- ArgumentParser(description = "Run PhyloAcc-c")
# Add arguments
parser$add_argument("--mod_file", type="character", help = "Path to mod file with tree and neutral rates")
parser$add_argument("--fasta", type="character", help = "Fasta file")
parser$add_argument("--traits", type="character", help = "Tab-delimited file with species and trait")
parser$add_argument("--iter", type="integer", default=10000, help = "Number of iterations per element (default=10000)")
parser$add_argument("--burnin", type="integer", default=1000, help = "Number of iterations of burnin to discard (default=1000)")
parser$add_argument("--max_missing", type="double", default=0.5, help = "Fraction of taxa need to be present in alignment (default=0.5)")
parser$add_argument("--outfolder", type="character", default="./", help="output folder path")
# Parse arguments
args <- parser$parse_args()
# Parse mod file
lines <- readLines(args$mod_file)
trait_file <- args$traits
alignment_loc <- args$fasta
for (line_index in seq_along(lines)) {
# Split the line by ":"
parts <- stri_split_fixed(lines[line_index], ":", n=2)[[1]]
# Extract variable name and IDs
variable <- trimws(parts[1]) # Trim whitespace
if (variable == "BACKGROUND") {
PI <- as.numeric(unlist(strsplit(trimws(parts[2]), "\\s+")))
}
if (variable == "RATE_MAT"){
mat_vals <- strsplit(lines[(line_index+1):(line_index+4)], "\\s+")
Q <- do.call(rbind, mat_vals)
Q <- Q[,-1]
Q <- matrix(as.numeric(Q), nrow=4, byrow=T)
line_index <- line_index + 4
}
if (variable == "TREE"){
tree = read.tree(text = trimws(parts[2]))
tree$node.label <- paste("Anc", 0:Nnode(tree), sep="")
}
}
# Get trait data
trait_tbl <- read.table(trait_file)
colnames(trait_tbl) <- c("NAME", "Y")
tree <- keep.tip(tree, trait_tbl$NAME)
trait_tbl$NAME <- factor(trait_tbl$NAME, levels = tree$tip.label)
trait_tbl <- trait_tbl[order(trait_tbl$NAME),]
row.names(trait_tbl) <- 1:nrow(trait_tbl)
element <- tools::file_path_sans_ext(basename(alignment_loc))
alignment <- read.dna(alignment_loc, format = "fasta", as.matrix = T, as.character = T)
if ((nrow(alignment) > ((length(tree$tip.label))*args$max_missing))) {
alignment <- alignment[rownames(alignment) %in% trait_tbl$NAME,]
alignment[alignment=="*"] <- "-"
alignment[alignment=="n"] <- "-"
all_same <- all(apply(alignment, 2, function(x) length(unique(x))) == 1)
missing_taxa <- setdiff(tree$tip.label, rownames(alignment))
if (length(missing_taxa > 0)){
for (miss in missing_taxa){
missing_row <- matrix(rep("-", ncol(alignment)), nrow=1)
rownames(missing_row) <- miss
alignment <- rbind(alignment, missing_row)
}
}
}
save_model <- run_model(tree=tree, alignment=alignment,
trait_tbl=trait_tbl,
transition_matrix=Q, stationary_distro=PI, n_iter=args$iter,
inner_iter=500,
Phi_prior=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
#Phi_prior=c(2.0,1.0,1.0,50.0,1.0,1.0),
#prior_lv=c(0.0,3.0),
#prior_lb2=c(0,0.1),
#prior_lb3=c(0,0.1),
prior_r2=c(4,0.1),
#prior_r2=c(8,0.04),
prior_r3=c(10,0.4))
r_trace = save_model$mcmc_trace$r_trace
v_trace = save_model$mcmc_trace$v_trace
burnin = args$burnin
elem = strsplit(basename(args$fasta), "\\.")[[1]][1]
pdf(paste(args$outfolder, "/", elem, ".pdf", sep=""), paper="letter")
plot.new()
text(.1,.1, paste("sig^2", c("10%", "50%", "90%"), round(quantile(v_trace[-burnin,1], prob=c(.1,.5,.9)), 2), collapse = ";", sep = " "), pos=4, cex=0.75)
text(.1,.2, paste("B3:B2", c("10%", "50%", "90%"), round(quantile(v_trace[-burnin,3] - v_trace[-burnin,2], prob=c(.1,.5,.9)), 2), collapse = ";", sep=" "), pos=4, cex=0.75)
text(.1,.3, paste("r2", c("10%", "50%", "90%"), round(quantile(r_trace[-burnin,2], probs=c(.1,.5,.9)), 2), collapse = ";", sep=" "), pos=4, cex=0.75)
text(.1,.4, paste("r3", c("10%", "50%", "90%"), round(quantile(r_trace[-burnin,3], probs=c(.1,.5,.9)), 2), collapse = ";", sep=" "), pos=4, cex=0.75)
text(.1,.5, elem)
plot_z(save_model, drop_n = burnin)
plot_mixing(save_model, ith=100)
invisible(capture.output(dev.off()))
#save.image(file = paste(args$outfolder, "/", elem, ".rdata", sep=""))4.4.3 Generate skeletal set
The following python code was used to generte the skeletal element GO set
import os
import csv
import gzip
HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
DATA = os.path.join(ROOT, "data")
RESULTS = os.path.join(ROOT, "results")
# ---- seed terms (high-level; DAG closure pulls in all descendants) ----------
# Two profiles are computed per gene:
# is_skeletal INCLUSIVE (recall-oriented; skeleton-wide phenotype)
# limb_appendicular STRICT (limb/appendicular skeleton only; precision)
# The inclusive flag is the default because skeleton-wide phenotype annotations
# are what catch genuine limb/skeletal genes (e.g. FGF18) whose mouse phenotype
# is recorded generically rather than as "limb bone". The strict flag lets the
# enrichment test be re-run as a sensitivity analysis (drops craniofacial-only
# genes such as BCAS3).
# Developmental GO processes (functional evidence; bone/cartilage/limb focused).
GO_SEEDS = {
"GO:0001501": "skeletal system development",
"GO:0060173": "limb development",
"GO:0048736": "appendage development",
"GO:0051216": "cartilage development",
"GO:0060348": "bone development",
"GO:0001503": "ossification",
"GO:0001649": "osteoblast differentiation",
"GO:0002062": "chondrocyte differentiation",
"GO:0061448": "connective tissue development",
"GO:0030509": "BMP signaling pathway",
}
# Inclusive (skeleton-wide + whole-limb) phenotype roots. MP nodes are the
# high-level system categories HMD actually uses.
MP_SEEDS = {
"MP:0005390": "skeleton phenotype",
"MP:0005371": "limbs/digits/tail phenotype",
}
HP_SEEDS = {
"HP:0000924": "Abnormality of the skeletal system",
"HP:0040064": "Abnormality of limbs",
}
# Strict limb/appendicular-skeleton subset (for the limb_appendicular sub-flag).
GO_LIMB = {
"GO:0060173": "limb development",
"GO:0048736": "appendage development",
}
# MGI is system-level (see header): the limb system category is the finest that
# matches HMD, and it correctly captures limb regulators like FGF18.
MP_LIMB = {"MP:0005371": "limbs/digits/tail phenotype"}
# HPO is term-level, so the limb tier can use bone-specific HPO nodes (keeping
# craniofacial-only genes such as BCAS3 out of the strict tier).
HP_LIMB = {
"HP:0011844": "Abnormal appendicular skeleton morphology",
"HP:0002813": "Abnormal limb bone morphology",
}
def load_obo_dag(path):
"""Return (alt2id, parents) where parents[id] = set of is_a/part_of targets."""
alt2id, parents = {}, {}
cur, obsolete = None, False
op = gzip.open if path.endswith(".gz") else open
with op(path, "rt") as fh:
for line in fh:
line = line.rstrip("\n")
if line == "[Term]":
cur, obsolete = None, False
elif line.startswith("id: "):
cur = line[4:].strip()
parents.setdefault(cur, set())
elif line.startswith("alt_id: ") and cur:
alt2id[line[8:].strip()] = cur
elif line.startswith("is_obsolete: true"):
obsolete = True
if cur:
parents.pop(cur, None)
cur = None
elif cur and not obsolete and line.startswith("is_a: "):
parents[cur].add(line[6:].split("!")[0].strip())
elif cur and not obsolete and line.startswith("relationship: part_of "):
parents[cur].add(line[len("relationship: part_of "):].split("!")[0].strip())
return alt2id, parents
def closure(seeds, parents):
"""All terms that are the seeds or descend (is_a/part_of) from them."""
children = {}
for c, ps in parents.items():
for p in ps:
children.setdefault(p, set()).add(c)
out, stack = set(seeds), list(seeds)
while stack:
t = stack.pop()
for c in children.get(t, ()):
if c not in out:
out.add(c)
stack.append(c)
return out
def norm_terms(terms, alt2id):
return {alt2id.get(t, t) for t in terms}
def load_chicken_go():
"""chicken GeneID -> set(GO)."""
p = os.path.join(DATA, "GCF_016699485.2_gene_ontology.gaf.gz")
d = {}
with gzip.open(p, "rt") as fh:
for line in fh:
if line.startswith("!"):
continue
f = line.rstrip("\n").split("\t")
if len(f) < 5:
continue
d.setdefault(f[1], set()).add(f[4])
return d
def load_goa_human():
"""human SYMBOL(upper) -> set(GO)."""
p = os.path.join(DATA, "goa_human.gaf.gz")
d = {}
with gzip.open(p, "rt") as fh:
for line in fh:
if line.startswith("!"):
continue
f = line.rstrip("\n").split("\t")
if len(f) < 5:
continue
d.setdefault(f[2].upper(), set()).add(f[4])
return d
def load_orthologs():
"""chicken GeneID -> human GeneID (NCBI gene_orthologs, hub=human 9606)."""
p = os.path.join(DATA, "gene_orthologs.gz")
c2h = {}
with gzip.open(p, "rt") as fh:
next(fh, None)
for line in fh:
f = line.rstrip("\n").split("\t")
if len(f) < 5:
continue
tax, gid, _, otax, ogid = f[0], f[1], f[2], f[3], f[4]
if tax == "9606" and otax == "9031":
c2h[ogid] = gid
elif otax == "9606" and tax == "9031":
c2h[gid] = ogid
return c2h
def load_gene_info(path, gz=True):
"""GeneID -> symbol ; symbol(upper) -> GeneID."""
id2sym, sym2id = {}, {}
op = gzip.open if gz else open
with op(path, "rt") as fh:
header = fh.readline().lstrip("#").rstrip("\n").split("\t")
gi, si = header.index("GeneID"), header.index("Symbol")
for line in fh:
f = line.rstrip("\n").split("\t")
id2sym[f[gi]] = f[si]
sym2id[f[si].upper()] = f[gi]
return id2sym, sym2id
def load_mgi_skeletal(mp_closure):
"""human symbols / GeneIDs whose mouse ortholog has a skeletal/limb MP."""
p = os.path.join(DATA, "HMD_HumanPhenotype.rpt")
syms, ids = set(), set()
with open(p) as fh:
for line in fh:
f = line.rstrip("\n").split("\t")
if len(f) < 5:
continue
hsym, hgid, mps = f[0], f[1], f[4]
terms = {t.strip() for t in mps.split(",") if t.strip()}
if terms & mp_closure:
syms.add(hsym.upper())
if hgid.isdigit():
ids.add(hgid)
return syms, ids
def load_hpo_skeletal(hp_closure):
"""human GeneIDs / symbols with a skeletal/limb HPO annotation."""
p = os.path.join(DATA, "genes_to_phenotype.txt")
ids, syms = set(), set()
with open(p) as fh:
header = fh.readline().rstrip("\n").split("\t")
gi = header.index("ncbi_gene_id"); si = header.index("gene_symbol"); hi = header.index("hpo_id")
for line in fh:
f = line.rstrip("\n").split("\t")
if len(f) <= hi:
continue
if f[hi] in hp_closure:
ids.add(f[gi]); syms.add(f[si].upper())
return ids, syms
def gene_sources(gid, symU, hsym_eff, hgid_eff, chick_go, go_alt, go_c,
human_go, mgi_sym, mgi_id, hpo_id, hpo_sym):
"""Which evidence sources flag this gene under a given closure profile."""
s = []
if norm_terms(chick_go.get(gid, set()), go_alt) & go_c:
s.append("chickenGO")
if norm_terms(human_go.get(hsym_eff, set()), go_alt) & go_c:
s.append("humanGO")
if hsym_eff in mgi_sym or (hgid_eff and hgid_eff in mgi_id):
s.append("MGI")
if (hgid_eff and hgid_eff in hpo_id) or hsym_eff in hpo_sym:
s.append("HPO")
return s
def main():
print("Building ontology DAG closures ...")
go_alt, go_par = load_obo_dag(os.path.join(DATA, "go-basic.obo"))
mp_alt, mp_par = load_obo_dag(os.path.join(DATA, "mp.obo"))
hp_alt, hp_par = load_obo_dag(os.path.join(DATA, "hp.obo"))
go_c = closure(GO_SEEDS, go_par); mp_c = closure(MP_SEEDS, mp_par); hp_c = closure(HP_SEEDS, hp_par)
go_l = closure(GO_LIMB, go_par); mp_l = closure(MP_LIMB, mp_par); hp_l = closure(HP_LIMB, hp_par)
print(f" inclusive: GO {len(go_c)} MP {len(mp_c)} HP {len(hp_c)} | "
f"limb: GO {len(go_l)} MP {len(mp_l)} HP {len(hp_l)}")
print("Loading annotations ...")
chick_go = load_chicken_go()
human_go = load_goa_human()
c2h = load_orthologs()
h_id2sym, _ = load_gene_info(os.path.join(DATA, "Homo_sapiens.gene_info.gz"))
chick_id2sym, _ = load_gene_info(os.path.join(DATA, "Gallus_gallus.gene_info.gz"))
mgi_sym, mgi_id = load_mgi_skeletal(mp_c) # inclusive
hpo_id, hpo_sym = load_hpo_skeletal(hp_c)
mgi_sym_l, mgi_id_l = load_mgi_skeletal(mp_l) # limb-restricted
hpo_id_l, hpo_sym_l = load_hpo_skeletal(hp_l)
print(f" chicken genes={len(chick_id2sym)} | chicken->human orthologs={len(c2h)}")
out = os.path.join(RESULTS, "skeletal_genes.tsv")
src_counts = {k: 0 for k in ("chickenGO", "humanGO", "MGI", "HPO")}
n_skel = n_limb = 0
with open(out, "w", newline="") as fh:
w = csv.writer(fh, delimiter="\t")
w.writerow(["chicken_gene_id", "chicken_symbol", "human_ortholog_symbol",
"is_skeletal", "sources", "limb_appendicular"])
for gid, sym in sorted(chick_id2sym.items(), key=lambda x: x[1]):
symU = sym.upper()
hgid = c2h.get(gid)
hsym = h_id2sym.get(hgid, "").upper() if hgid else ""
hsym_eff = hsym or symU # symbol fallback when no NCBI ortholog
hgid_eff = hgid or ""
sources = gene_sources(gid, symU, hsym_eff, hgid_eff, chick_go, go_alt,
go_c, human_go, mgi_sym, mgi_id, hpo_id, hpo_sym)
limb = gene_sources(gid, symU, hsym_eff, hgid_eff, chick_go, go_alt,
go_l, human_go, mgi_sym_l, mgi_id_l, hpo_id_l, hpo_sym_l)
is_skel = bool(sources)
if is_skel:
n_skel += 1
for s in sources:
src_counts[s] += 1
if limb:
n_limb += 1
w.writerow([gid, sym, h_id2sym.get(hgid, "") if hgid else "",
int(is_skel), ";".join(sources), int(bool(limb))])
print(f"\nWrote {out}")
n = len(chick_id2sym)
print(f" is_skeletal (inclusive): {n_skel}/{n} ({100*n_skel/n:.1f}%)")
print(f" limb_appendicular (strict): {n_limb}/{n} ({100*n_limb/n:.1f}%)")
print(" inclusive by source (non-exclusive):", src_counts)
print("\n verification:")
tsv = {r["chicken_symbol"]: r for r in csv.DictReader(open(out), delimiter="\t")}
for g in ["IGF2", "SALL4", "SHOX", "RUNX2", "SOX5", "FGF18", "GLI3", "TBX5",
"GDF5", "BCAS3", "INS", "CILP2", "GAPDH"]:
r = tsv.get(g)
if r:
print(f" {g:8s} skeletal={r['is_skeletal']} limb={r['limb_appendicular']} "
f"via={r['sources'] or '-'}")
else:
print(f" {g:8s} (not in chicken gene_info)")
if __name__ == "__main__":
main()
4.4 Parsing outputs
4.4.1 Femur
4.4.2 Humerus
4.4.3 Tarsometatarsus
4.4.4 Carpometacarpus
References
Porto, Arthur, Sara Rolfe, and A. Murat Maga. 2021. “ALPACA: A
Fast and Accurate Computer Vision Approach for Automated Landmarking of
Three-Dimensional Biological Structures.” Methods in Ecology
and Evolution 12 (11): 2129–44. https://doi.org/https://doi.org/10.1111/2041-210X.13689.
Rolfe, Sara, Christopher Davis, and A. Murat Maga. 2021.
“Comparing Semi-Landmarking Approaches for Analyzing
Three-Dimensional Cranial Morphology.” American Journal of
Physical Anthropology 175 (1): 227–37. https://doi.org/https://doi.org/10.1002/ajpa.24214.
Rolfe, Sara, Steve Pieper, Arthur Porto, et al. 2021.
“SlicerMorph: An Open and Extensible Platform to Retrieve,
Visualize and Analyse 3D Morphology.” Methods in Ecology and
Evolution 12 (10): 1816–25. https://doi.org/https://doi.org/10.1111/2041-210X.13669.