Statistical analysis

Statistical analysis

CHC data

1. Trim CHC data to the following 19 peaks that represent well-established honey bee CHCs: C_18:1, C_23:1, C₂₃, C₂₄, C_25:1, C₂₅, C_27:1, C₂₇, 11,13 methyl C₂₇, C₂₈, C_29:1, C₂₉, 11,13 methyl C₂₉, z(7) C_31:1, z(9) C_31:1, C₃₁, C_33:2, C_33:1, C₃₃.
2. Calculate the relative proportion of each compound in each sample.
   1. Calculate by dividing the ng of the individual compound by the total ng of all compounds per sample.
   2. Rescale each proportion to a value between 0 and 1: 𝑧𝑖 = (𝑥𝑖 − min(𝑥))/(max(𝑥) − min(𝑥)).
3. Make a data frame, with each column representing a different compound. The first column should be your group designations.

Microbiome data

1. Use ASV table taken from DADA2 analysis. Make sure your groups are correctly labeled. Each column representing a different ASV. The first column should be your group designations.

Both CHC and Microbiome

1. After steps above, you can statistically analyze each dataset.
2. For each dataset, a permutation MANOVA was run using the ADONIS function in the vegan package of R with Bray-Curtis dissimilarity measures. 
   1. Example R code: adonis(C[,2:20] ~ C$group, perm = 10000) #columns 2-20 represent my data, where as column 1 represents my group.

3. Run pairwise comparisons with FDR p-value correction if more than two groups are compared. 

   Example R code for Pairwise PerMANOVA, taken from https://www.researchgate.net/post/How_can_I_do_PerMANOVA_pairwise_contrasts_in_R

   x = C[,2:20] #Include only data in your dataframe, no text.

    

   factors= C$group

   pairwise.adonis <- function(x,factors, sim.method, p.adjust.m)

   {

     library(vegan)

     co = as.matrix(combn(unique(factors),2))

     pairs = c()

     F.Model =c()

     R2 = c()

     p.value = c()

    

     for(elem in 1:ncol(co)){

       ad = adonis(x[factors %in% 

   c(as.character(co[1,elem]),as.character(co[2,elem])),] ~

                        factors[factors %in% 

   c(as.character(co[1,elem]),as.character(co[2,elem]))] , method =sim.method);

       pairs = c(pairs,paste(co[1,elem],'vs',co[2,elem]));

       F.Model =c(F.Model,ad$aov.tab[1,4]);

       R2 = c(R2,ad$aov.tab[1,5]);

       p.value = c(p.value,ad$aov.tab[1,6])

     }

     p.adjusted = p.adjust(p.value,method=p.adjust.m)

     pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)

     return(pairw.res)

   }

   PW.Adonis.fdr=pairwise.adonis(x,factors,sim.method="bray",p.adjust.m = "fdr")

    

4. Visualize data using non-metric multidimensional scaling using Bray-Curtis dissimilarity, and either 2 or three dimensions in order to minimize stress to <0.1. 

   Example R code:

   library(vegan)

   C.2 <- C[,2:20] #Just the data from your dataframe, no text.

   C.NMDS <- metaMDS(C.2, k=3, trymax= 100) #k represents dimensions

   C.NMDS

   stressplot(C.NMDS)

   plot(C.NMDS)

   text(C.NMDS, display= "sites")

5. To test for differences in the number of reads of each ASV (or in amount of each CHC compound) between groups, perform an ANOVA, t-test, non-parametric Kruskal Wallis, or Wilcoxon Rank Sum test, depending on normality and number of groups.

   1. Example R code:

      C.s <- stack(C[,4:22]) #Must stack your data before you can run the test.

      write.csv(C.s, "data.csv") #Write as a .csv and then edit to include your group information

      C.s2 <- read.csv("data.csv", header=T)

      for(i in unique(C.s2$Compound)) {

                  c <- subset(C.s2, C.s2$Compound == i)

                  print(i)

                  print(kruskal.test(Proportion ~ group, data = c))}

Behavioral data

1. Calculate the proportion of bees accepted by guard honey bees for each experimental group. Use these numbers in your plots.
2. Run a Pearson’s chi-square for each experimental group. 
   a.    To do this, you must use raw numbers and make a table with the number accepted and number rejected for each group.
   Example R code: 
   T = read.csv("dataAR.csv", header=T)
   library(MASS)
   tbl = table(T$group, T$AR)
   tbl
   b.    Then run a chi-square using this table
   Example R code:
   chisq.test(tbl)
   c.    Finally, run a pairwise chi-square:
   Example R code:
   FUN = function(i,j){     
        chisq.test(matrix(c(tbl[i,1], tbl[i,2],
                            tbl[j,1], tbl[j,2]),
                   nrow=2,
                   byrow=TRUE))$ p.value 
                  }

Strain comparisons

1. Perform alignments in USEARCH and calculate proportion of SNPs between each comparison.
2. Run Kruskal Wallis between groups, with a Dunn’s test for pairwise comparisons.

   1. Example R code:

      KW <- kruskal.test(T$prop, T$group)

      require(PMCMR)

      print(posthoc.kruskal.dunn.test(x=T$prop, g=T$group, p.adjust="fdr"))

3. Plot using Violin plot.

   1. Example R code:

      ggplot(T, aes(x=group, y = SNPs)) +

        geom_violin() + 

        stat_summary(mapping = aes(x = group, y = SNPs),

        fun.ymin = function(z) { quantile(z,0.25) },

        fun.ymax = function(z) { quantile(z,0.75) },

        fun.y = median) +

        theme_classic()

Category