library(vegan) library(haplotypes) library(ape) # read in FASTA file as.cytb <- read.fas("BB_SN_CYTB_alignment.fasta") # assign populations (or read in table) samnames <- as.cytb@seqnames samnames Amova_factors <- read.csv("Amova_factors.csv") ordered_pops <- Amova_factors$Abbrev[match(samnames, Amova_factors$ID)] ordered_morphs <- Amova_factors$Morph[match(samnames, Amova_factors$ID)] ordered_regions <- Amova_factors$Region[match(samnames, Amova_factors$ID)] ordered_fac <- paste(ordered_pops,ordered_morphs,sep="_") table(ordered_fac) # check sample sizes # remove the small small samples and ignore UNIDENTIFIED keepers <- ordered_fac=="COS_BB" | ordered_fac=="COS_SN" | ordered_fac=="OBB_SN" | ordered_fac=="PTH_SN" |ordered_fac=="PTH_BB" | ordered_fac=="RF_BB" | ordered_fac=="RF_SN" cytb.k <- as.dna(as.cytb[keepers,]) fact.k <- ordered_fac[keepers] pops.k <- ordered_pops[keepers] morp.k <- ordered_morphs[keepers] # PhiST estimates amova.1 <- pairPhiST(x=cytb.k, population=fact.k,indels="m",nperm=10000,negatives=FALSE, subset=NULL, showprogbar=TRUE) amova.1 # distance based redundancy analysis dist.k <- haplotypes::distance(cytb.k, indels="missing") cytb.capN <- capscale(dist.k~pops.k+morp.k) # model selection drop1(cytb.capN, test="perm", permutations = how(nperm=99999)) # marginal effects anova(cytb.capN,by="margin", permutations = how(nperm=99999)) anova(capscale(dist.k~morp.k), permutations = how(nperm=99999)) anova(capscale(dist.k~pops.k), permutations = how(nperm=99999))