AMOVA with unbalanced groups/Adonis permanova

[D-rgz]

Hi I want to do permutational analysis of variance to calculate variation between nematode haplotypes of the same species in samples , i have unbalanced haplotypes in samples and in regions (108 samples, 3 regions, sharing 27 haplotypes). Each sample has a range of between 1 to 9 haplotypes, I have built a pairwise distance matrix (hamming) for all haplotypes. I want to build an R script for performing analysis of molecular variance within samples and between samples in a regions and between regions (I have 3 regions). I guess my challenge is to build the correct dataframe shape ..i tried the nested design making Region and SampleID for each Haplotype so the rows match the dist_matrix but it just give errors, and also i realize that i have to build it using the permute function over samples??? I also want to do an ordination for the different sample to show how they vary between the 3 regions
Your help will be much appreciated.

[gavinsimpson]

You're asking to do at least three different things here. I'm not going to do your analysis for you and I doubt Jari will either. Please focus your question to one thing at a time. Do you want help with PERMANOVA, AMOVA, or designing a permutation test?

[D-rgz]

# Load necessary libraries
install.packages("adegenet")
library(adegenet)
library(adegenet)
library(pegas)
library(tidyr)
library(dplyr)
library("permute")
library(vegan)
library(tidymodels)
library(pegas)
library(ape)

# make sure you have Rtools installed
if (!require('devtools')) install.packages('devtools')
# install sprex development version
devtools::install_github("ericarcher/sprex")
# install strataG latest version
devtools::install_github('ericarcher/strataG', build_vignettes = TRUE)

options(repos = c(
  zkamvar = 'https://zkamvar.r-universe.dev',
  CRAN = 'https://cloud.r-project.org'))

install.packages('strataG')
library(strataG)
# install diveRsity from github
devtools::install_github("kkeenan02/diveRsity")
library(diveRsity)




# Import Arlequin data using adegenet's read.arlequin function
# Replace "path/to/your_file.arp" with the path to your Arlequin file
arp_data <- arlequinRead("C:\\Users\\DR DEO\\Documents\\genetic_div\\Nad4haplotypes.arp")

arp_data2 <- arp2gen("C:\\Users\\DR DEO\\Documents\\genetic_div\\Nad4haplotypes.arp")
a<-as.data.frame(arp_data$data.info$haplotype.definition$list)
View(a)
# Convert the Arlequin data to a genind object
sequences <- a$locus_1

# Convert sequences from string to a list of character vectors
seq_list <- strsplit(sequences, "")

# Convert list of character vectors to a matrix
seq_matrix <- do.call(rbind, seq_list)

# Convert the matrix to a DNAbin object
dnabin_obj <- as.DNAbin(seq_matrix)
# Inspect the dnabin object
str(dnabin_obj)

rownames(dnabin_obj)<-a$id


data <- c(
  "NAD4.1008 =  Muleba: ASV_001, ASV_002, ASV_004",
  "NAD4.1064 =  Muleba: ASV_001, ASV_002, ASV_004, ASV_006",
  "NAD4.1103 =  Muleba: ASV_001, ASV_002",
  "NAD4.1125 =  Muleba: ASV_001, ASV_003",
  "NAD4.1128 =  Muleba: ASV_001, ASV_002",
  "NAD4.1131 =  Muleba: ASV_002, ASV_006, ASV_007",
  "NAD4.1135 =  Muleba: ASV_001, ASV_002, ASV_006, ASV_016",
  "NAD4.1144 =  Muleba: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_025",
  "NAD4.1157 =  Muleba: ASV_001, ASV_002, ASV_006, ASV_025",
  "NAD4.1230 =  Muleba: ASV_001, ASV_002, ASV_004, ASV_006",
  "NAD4.1237 =  Muleba: ASV_001, ASV_002, ASV_004, ASV_005, ASV_027",
  "NAD4.1263 =  Muleba: ASV_002, ASV_004",
  "NAD4.1304 =  Muleba: ASV_001, ASV_002, ASV_011",
  "NAD4.1305 =  Muleba: ASV_001, ASV_006",
  "NAD4.1314 =  Muleba: ASV_001, ASV_002, ASV_004",
  "NAD4.1381 =  Muleba: ASV_001, ASV_002, ASV_024",
  "NAD4.1387 =  Muleba: ASV_001, ASV_002, ASV_004, ASV_027",
  "NAD4.1420 =  Muleba: ASV_002",
  "NAD4.1434 =  Muleba: ASV_001, ASV_002, ASV_004, ASV_024",
  "NAD4.1512 =  Muleba: ASV_001, ASV_002, ASV_004 ",
  "NAD4.819  =  Muleba: ASV_001, ASV_002, ASV_003, ASV_004, ASV_006",
  "NAD4.833  =  Muleba: ASV_001, ASV_002, ASV_003, ASV_004, ASV_006, ASV_025",
  "NAD4.834  =  Muleba: ASV_001, ASV_002, ASV_004, ASV_005",
  "NAD4.933  =  Muleba: ASV_001, ASV_002, ASV_004, ASV_017",
  "NAD4.950  =  Muleba: ASV_001, ASV_003, ASV_017",
  "NAD4.977  =  Muleba: ASV_001, ASV_003, ASV_007, ASV_011",
  "NAD4.988  =  Muleba: ASV_002, ASV_004, ASV_007",
  "NAD4.996  =  Muleba: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007",
  "NAD4.44   =  Buchosa: ASV_001, ASV_002",
  "NAD4.148  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.151  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_011, ASV_024, ASV_025, ASV_027",
  "NAD4.172  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007",
  "NAD4.195  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_007",
  "NAD4.316  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004, ASV_015, ASV_025" ,
  "NAD4.326  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_006, ASV_007, ASV_012",
  "NAD4.335  =  Buchosa: ASV_001, ASV_002, ASV_005, ASV_006, ASV_007, ASV_014, ASV_017, ASV_025",
  "NAD4.338  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.341  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_006, ASV_007, ASV_026",
  "NAD4.342  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_006, ASV_007, ASV_026",
  "NAD4.345  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007",
  "NAD4.363  =  Buchosa: ASV_001, ASV_002, ASV_004, ASV_015",
  "NAD4.373  =  Buchosa: ASV_001, ASV_002, ASV_013",
  "NAD4.385  =  Buchosa: ASV_001, ASV_002, ASV_003",
  "NAD4.389  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004, ASV_015, ASV_027",
  "NAD4.391  =  Buchosa: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007",
  "NAD4.1517 =  Kaliua: ASV_001, ASV_002",
  "NAD4.1520 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.1534 =  Kaliua: ASV_001, ASV_002, ASV_007",
  "NAD4.1538 =  Kaliua: ASV_001, ASV_002, ASV_023",
  "NAD4.1546 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.1571 =  Kaliua: ASV_001",
  "NAD4.1572 =  Kaliua: ASV_001",
  "NAD4.1576 =  Kaliua: ASV_001, ASV_002, ASV_004",
  "NAD4.1577 =  Kaliua: ASV_003",
  "NAD4.1586 =  Kaliua: ASV_002, ASV_004",
  "NAD4.1589 =  Kaliua: ASV_001, ASV_002, ASV_005",
  "NAD4.1651 =  Kaliua: ASV_001, ASV_004, ASV_008",
  "NAD4.1702 =  Kaliua: ASV_001, ASV_002, ASV_018",
  "NAD4.1703 =  Kaliua: ASV_001, ASV_002, ASV_004",
  "NAD4.1705 =  Kaliua: ASV_002, ASV_003, ASV_008, ASV_022",
  "NAD4.1754 =  Kaliua: ASV_001, ASV_002, ASV_012, ASV_026",
  "NAD4.1757 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_007",
  "NAD4.1828 =  Kaliua: ASV_001, ASV_002, ASV_007, ASV_010, ASV_016, ASV_018",
  "NAD4.1848 =  Kaliua: ASV_001, ASV_002, ASV_003",
  "NAD4.1873 =  Kaliua: ASV_001, ASV_002, ASV_004",
  "NAD4.1875 =  Kaliua: ASV_001, ASV_002, ASV_004",
  "NAD4.1883 =  Kaliua: ASV_001",
  "NAD4.1884 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_010",
  "NAD4.1889 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_007",
  "NAD4.1892 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_006, ASV_018",
  "NAD4.1893 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.1906 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_027",
  "NAD4.1933 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_017, ASV_024",
  "NAD4.1944 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_008, ASV_010, ASV_022",
  "NAD4.1947 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_006, ASV_007",
  "NAD4.1948 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.1974 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_010, ASV_017, ASV_025",
  "NAD4.2013 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_011",
  "NAD4.2017 =  Kaliua: ASV_001, ASV_017",
  "NAD4.2023 =  Kaliua: ASV_002, ASV_003, ASV_020",
  "NAD4.2030 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004",
  "NAD4.2031 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_010, ASV_020, ASV_025",
  "NAD4.2040 =  Kaliua: ASV_002, ASV_005, ASV_010",
  "NAD4.2041 =  Kaliua: ASV_001, ASV_002",
  "NAD4.2044 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_010, ASV_017, ASV_027",
  "NAD4.2045 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_019",
  "NAD4.2046 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_007, ASV_009, ASV_025",
  "NAD4.2047 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_027",
  "NAD4.2048 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_006, ASV_019",
  "NAD4.2195 =  Kaliua: ASV_001, ASV_002, ASV_008",
  "NAD4.2207 =  Kaliua: ASV_001, ASV_008, ASV_023",
  "NAD4.2212 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_012",
  "NAD4.2224 =  Kaliua: ASV_001, ASV_002",
  "NAD4.2226 =  Kaliua: ASV_001, ASV_002, ASV_009",
  "NAD4.2227 =  Kaliua: ASV_001, ASV_002, ASV_003",
  "NAD4.2228 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_007, ASV_012, ASV_015, ASV_017",
  "NAD4.2236 =  Kaliua: ASV_001, ASV_002, ASV_007",
  "NAD4.2239 =  Kaliua: ASV_002, ASV_013",
  "NAD4.2245 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_005, ASV_006, ASV_022, ASV_023, ASV_026",
  "NAD4.2365 =  Kaliua: ASV_001, ASV_002",
  "NAD4.2373 =  Kaliua: ASV_001, ASV_002",
  "NAD4.2377 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_005, ASV_007",
  "NAD4.2422 =  Kaliua: ASV_001, ASV_002, ASV_021",
  "NAD4.2437 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_007, ASV_009",
  "NAD4.2439 =  Kaliua: ASV_001, ASV_002, ASV_004",
  "NAD4.2440 =  Kaliua: ASV_001, ASV_002, ASV_004, ASV_005, ASV_010",
  "NAD4.2450 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_004, ASV_006, ASV_007, ASV_008, ASV_009, ASV_010",
  "NAD4.2451 =  Kaliua: ASV_001, ASV_002, ASV_003, ASV_010"
)

# Step 2: Parsing the data
library(stringr)
parsed_data <- str_match(data, "(.*?) = (.*?): (.*)")

# Step 3: Create a data frame
result <- data.frame(
  SampleID = parsed_data[, 2],
  Region = parsed_data[, 3],
  ASVs = parsed_data[, 4],
  stringsAsFactors = FALSE
)



df_long <- result %>%
  separate_rows(ASVs, sep = ", ") %>%
  rename(ASV = ASVs) %>%
  arrange(Region, SampleID, ASV)

transform_asv_column <- function(df, asv_column='ASV') {
  df[[asv_column]] <- paste0("Hap_", as.numeric(sub("ASV_", "", df[[asv_column]])))
  return(df)
}

# Applying the function to the dataframe
df_long_transformed <- transform_asv_column(df_long)

# Displaying the transformed dataframe
print(df_long_transformed)

pop_structure<-df_long_transformed

pop_structure$Region <- as.factor(pop_structure$Region)
pop_structure$SampleID <- as.factor(pop_structure$SampleID)


dist_matrix <- dist.dna(dnabin_obj, model = "N", pairwise.deletion = TRUE)

#### Distance matrix has 27 rows corresponding to haplotypes



pop_nested<-pop_structure %>% group_by(ASV) %>% nest() %>%  arrange(as.numeric(sub("Hap_", "", ASV)))

pop_nested2<-pop_structure %>% group_by(SampleID) %>% nest() %>% arrange(SampleID)



# Perform AMOVA
permanova_result <- adonis2(dist_matrix ~ Region + SampleID, data = pop_structure ,permutations = 999)

##Error message
This is the error i get: Error in qr.fitted(qrhs, G) : 
  'qr' and 'y' must have the same number of rows

[gavinsimpson]

Well, pop-structure has 440 rows:

# A tibble: 440 × 3
   SampleID    Region     ASV
   <chr>       <chr>      <chr>
 1 "NAD4.148 " " Buchosa" Hap_1
 2 "NAD4.148 " " Buchosa" Hap_2
 3 "NAD4.148 " " Buchosa" Hap_3
 4 "NAD4.148 " " Buchosa" Hap_4
 5 "NAD4.151 " " Buchosa" Hap_1
 6 "NAD4.151 " " Buchosa" Hap_2
 7 "NAD4.151 " " Buchosa" Hap_4
 8 "NAD4.151 " " Buchosa" Hap_11
 9 "NAD4.151 " " Buchosa" Hap_24
10 "NAD4.151 " " Buchosa" Hap_25
# ℹ 430 more rows
# ℹ Use `print(n = ...)` to see more rows

while you claim dist_matrix has 27 rows (and surely 27 columns if we expand the distance matrix to a full matrix instead of a "dist" class object that it should be?). 440 and 27 are not the same.

[D-rgz]

Yes that is true, is there a workaround you suggest to get adosnis2() call to work, or to match rows in LHS and RHS

[gavinsimpson]

The work around is to get a data frame with 1 row per Haplotype if that's what you say the distance matrix represents. In which case you would have a model where the distances between haplotypes are decomposed into a component due to Region and SampleID.

[jarioksa]

That's really messy. I wouldn't touch anything with that many library() calls and then asking about things in vegan. How did you come up with that kind of workflow? Copy-and-paste? A reproducible example would be much more focussed.

[D-rgz]

Hi Jari and Gavin. Sorry that the other library calls put you off, They shouldn't have been there. I wrote the workflow by myself and I am working with real data. I want to calculate the genetic variation in nematode mtDNA haplotypes I have my distance matrix (pairwise for all the haplotypes) which is a pairwise distance generated by the hamming distance. On my RHS I have I have three levels. 1. Within sample variability ;Samples where I have the haplotypes, 2. Between samples variability; Samples in regions and 3 genetic variation between Regions. The observations are unequal due to uneven distribution of haplotypes within samples and samples within regions. This gives me a challenge in setting up and calling the adonis2() call because of the factors df structure. I have seen of the nested df can match the rows but results in an error message about the column/factor. So the question is am I on the right path??

[gavinsimpson]

No; you need to provide a data frame with as many rows as there are rows / cols in the dissimilarity matrix. I barely know what a haplotype is let alone the details of the workflow you are using to get to the point where you are trying to call adonis2(). But you can't provide nested DFs (from dplyr) or anything like that. You need to create a data frame with the correct numbers of rows (27) one per observation in your dissimilarity object.

[D-rgz]

Thank you so much Gavin I was able to run the adonis2() call after matching the the LHS and RHS rows. It turns out I just had work the LHS dist_matrix(with redundant rows for haplotypes) and of course the RHS df. Now I just have to figure out restricted permutations and I am done. Thank you so much.