Intro

Reanalysis of the genome formula data with PERMANOVA. The script follows the method as presented in Boezen et al. 2023, as deposited in Zenodo.

Load the data and necessary libraries

Load the experimental data from Sicard et al. 2013 on FBNSV in V. faba and M. truncatula, underlying Figure 2b. Only the data from aphid transmission are included in the datafile and analysis described here.

# Load libraries
library(magrittr)
library(utils)
library(plyr)
library(dplyr)
library(vegan)

# Load data
setwd("C:/R_Data")
c.data <- read.csv("Dataset_7_RBNSV_Sicard_Fig_2b.csv") %>%
    as.data.frame()
# head(c.data, 20)

Analysis of data with PERMANOVA

We analyze the data, chcecking for an effect of host on the genome formula distance. By modifying the filtering criteria (see below), the script can also perform pairwise comparisons between experiments.

# No filtering of the data needed for this example, unless you want to perform
# pairwise comparisons.
 c.data.f <- c.data
# c.data.f <- c.data %>%
#  filter(Uni_Num %in% c(1, 4))
  
# Generate distances matrix and run PERMANOVA
set.seed(42)
dist.f <- vegdist(c.data.f[,6:13], method = "euclidean") 
per.gf <- adonis2(formula = dist.f ~ Uni_Num, data = c.data.f, permutations = 10^4) 
print(per.gf)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = dist.f ~ Uni_Num, data = c.data.f, permutations = 10^4)
##          Df SumOfSqs      R2      F    Pr(>F)    
## Uni_Num   1   1.0109 0.36577 40.946 9.999e-05 ***
## Residual 71   1.7529 0.63423                     
## Total    72   2.7639 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Check for significant differences in spread between the groups, using a PERMDISP2 procedure
pd.gf <- betadisper(d = dist.f, group = c.data.f$Uni_Num, type = "median", bias.adjust = TRUE, 
           sqrt.dist = FALSE, add = FALSE)
print(pd.gf)
## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = dist.f, group = c.data.f$Uni_Num, type = "median",
## bias.adjust = TRUE, sqrt.dist = FALSE, add = FALSE)
## 
## No. of Positive Eigenvalues: 8
## No. of Negative Eigenvalues: 0
## 
## Average distance to median:
##       1       2       3       4       5 
## 0.08583 0.10001 0.13412 0.07912 0.10113 
## 
## Eigenvalues for PCoA axes:
##     PCoA1     PCoA2     PCoA3     PCoA4     PCoA5     PCoA6     PCoA7     PCoA8 
## 1.911e+00 4.327e-01 2.417e-01 1.183e-01 3.695e-02 1.844e-02 4.510e-03 3.131e-06
ph.pd.gf <- permutest(pd.gf, pairwise = TRUE, permutations = 10^4)
print(ph.pd.gf)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm  Pr(>F)  
## Groups     4 0.037519 0.0093797 2.0823  10000 0.08789 .
## Residuals 68 0.306314 0.0045046                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          1        2        3        4      5
## 1          0.492551 0.061694 0.702330 0.6428
## 2 0.454024          0.287971 0.156784 0.9804
## 3 0.067751 0.296111          0.018698 0.2243
## 4 0.671439 0.163812 0.022303          0.4531
## 5 0.615760 0.975922 0.218734 0.418444
p.adjust(ph.pd.gf$pairwise$permuted, method = "holm")
##       1-2       1-3       1-4       1-5       2-3       2-4       2-5       3-4 
## 1.0000000 0.5552445 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.1869813 
##       3-5       4-5 
## 1.0000000 1.0000000