Symbiodiniaceae ITS2 analyses
Ryan Eckert – ryan.j.eckert@gmail.com
version: April 3, 2020
About this document
All analyses preformed with R version 3.6.1.
This is the code accompanies the publication Eckert RJ, Reaume AM, Sturm AB, Studivan MS and Voss JD (2020) Depth influences Symbiodiniaceae associations among Montastraea cavernosa corals on the Belize Barrier Reef. Front. Microbiol. 11:518. doi: 10.3389/fmicb.2020.00518. Here you will find all the code to repeat the statistical analyses performed for this manuscript. All of the accompanying data can be found on my GitHub.
Background
Here we sampled Montastraea cavernosa colonies from 4 different depth zones across 4 sites on the Belize Barrier Reef. We extracted holobiont DNA using a modified CTAB protocol and amplified the ITS2 region of algal symbiont (Family Symbiodiniaceae) DNA for sequencing on Illumina MiSeq (2 x 300 bp). We used the SymPortal analysis frame work to identify Symbiodiniaceae ITS2 type profiles within the coral host to examine variation with depth.
If you download my entire accompanying github directory you should be able to re-run these analyses by following along with the code chunks in R Studio. If you download the code separtely or you are using this pipeline on your own data, you may need to change the working directory to where the associated files are housed (i.e. setwd("~/path/to/directory/with/data")).
The data used for this analysis are calculated with SymPortal . The raw Symbiodiniaceae ITS2 sequences obtained from Montastraea cavernosa samples can be found in the NCBI SRA under project number PRJNA579363. Hopefully you are able to follow along with this file and find it useful to use with your own data!
Basic setup of R environment
Loading required packages
For the following analyses we will require the use of a number of different R packages. Most of these can be sourced from CRAN, but a couple need to be downloaded from GitHub or BioConducter. We can use the following code to quickly load in the packages and install any packages not previously installed in the R console.
if (!require("pacman")) install.packages("pacman")
pacman::p_load("dplyr", "edgeR", "ggplot2", "MCMC.OTU", "pairwiseAdonis", "rcartocolor", "RColorBrewer", "Redmonder", "reshape2", "vegan")
pacman::p_load_gh("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
if (!require("edgeR")){BiocManager::install("edgeR", update = FALSE)
library(edgeR)}
options("scipen" = 10)ITS2 sequences
Loading ITS2 sequence data into R
First, we need to load in the data from SymPortal analysis and clean it up in R. We can order sites from north to south so later data can be plotted shallow to deep; north to south. To do this we can set “Depth_zone” and “Sample_site” as factors and define the order of the factors.
its2Seq = read.delim("62_20190310_DBV_2019-03-11_01-11-25.167036.seqs.absolute.clean.txt", header = TRUE, check.names = FALSE)
head(its2Seq)
its2MetaData = read.delim("CBC_MCAV_sampling_metadata.txt", header = TRUE, check.names = FALSE)
head(its2MetaData)
its2Seq = cbind(its2Seq[1], its2MetaData[,2:3], its2Seq[,c(2:length(its2Seq))])
colnames(its2Seq)[3] = "Depth_zone"
head(its2Seq)
its2Seq$Depth_zone = factor(its2Seq$Depth_zone, levels = c("10", "16", "25", "35"))
levels(its2Seq$Depth_zone)
its2Seq$Sample_site = factor(its2Seq$Sample_site, levels(its2Seq$Sample_site)[c(4, 2, 3, 1)])
its2Seq = its2Seq[order(its2Seq$Sample_site, its2Seq$Depth_zone), ]
levels(its2Seq$Sample_site)
head(its2Seq)## sample_name Sample_site Depth_zone noName Clade A noName Clade B
## 107 CBC_107 TR 10 0 0
## 108 CBC_108 TR 10 0 0
## 109 CBC_109 TR 10 0 0
## 110 CBC_110 TR 10 0 0
## 111 CBC_111 TR 10 0 0
## 112 CBC_112 TR 10 0 0
## noName Clade C noName Clade D noName Clade E noName Clade F noName Clade G
## 107 204 0 0 0 0
## 108 15750 0 0 0 0
## 109 1101 0 0 0 0
## 110 13491 0 0 0 0
## 111 11425 0 0 0 0
## 112 1429 0 0 0 0
## noName Clade H noName Clade I A13 A4a A4 50815_A 71474_A 45275_A 71448_A
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71475_A 29923_A 71501_A B18b 71441_B B21 40228_B 44807_B B1 71508_B 71509_B
## 107 0 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 6 0 0
## 111 0 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0 0
## 71490_B 71510_B 71511_B 71514_B 71515_B 71516_B B19 40229_B 71512_B 71517_B
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0
## 71518_B 71513_B 71519_B 71520_B 71521_B 71522_B 71523_B 71367_B 71499_B
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71524_B 71525_B 71368_B 71403_B 71526_B 71527_B 71528_B 71402_B 71491_B
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 29623_B 71420_B 71369_B 71500_B 71407_B 71408_B 71547_B 71506_B C3 C3de
## 107 0 0 0 0 0 0 0 0 1729 906
## 108 0 0 0 0 0 0 0 0 92525 30348
## 109 0 0 0 0 0 0 0 0 7587 1269
## 110 0 0 0 0 0 0 0 0 59532 42690
## 111 0 0 0 0 0 0 0 0 52304 39964
## 112 0 0 0 0 0 0 0 0 7106 4657
## C3bb C21ae C3an C3s C3dk C3cd C3dm 71372_C 71371_C 71374_C C3dn 71376_C
## 107 549 207 467 114 103 0 0 37 26 42 0 30
## 108 14205 14597 8941 6955 5609 0 0 3131 1788 1487 0 0
## 109 713 428 601 286 224 0 0 64 67 66 0 0
## 110 20025 18009 10091 8928 7498 4643 0 3452 2498 2165 2303 0
## 111 19281 16220 9980 8003 7245 0 0 2571 1769 1616 0 1540
## 112 1846 1609 1078 611 629 0 0 265 248 209 0 132
## 12353_C C21 71394_C 71377_C C3b 71442_C 71459_C 71460_C 71458_C 37807_C
## 107 0 114 0 0 0 0 0 0 0 0
## 108 3864 0 0 0 1925 0 1189 1260 939 0
## 109 411 99 0 0 222 0 59 56 69 0
## 110 0 0 0 0 0 0 958 1159 1289 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 99 96 107
## 71397_C 71431_C 71444_C 71392_C 2301_C 71375_C 71386_C 71387_C 71440_C
## 107 0 34 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 1471 0 0 0 0 0 0 0 0
## 112 163 0 0 0 0 0 0 0 0
## 23653_C 71395_C 71445_C 71413_C 71391_C 71406_C 71443_C 11189_C C1 71379_C
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 880 0 0
## 109 56 0 0 0 0 79 0 99 0 0
## 110 0 0 0 0 0 0 0 0 0 0
## 111 0 0 0 1320 0 0 0 0 0 0
## 112 0 0 0 110 0 0 0 0 0 0
## 71393_C 71469_C 71411_C 71381_C 71382_C 49051_C 71380_C 25182_C 3031_C
## 107 0 17 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 1138 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71477_C 71543_C 20889_C 71373_C 71423_C 71422_C 71428_C C3i 71404_C 71446_C
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 1212 0 0 0 0 0 0 0
## 109 0 0 75 0 0 0 0 0 0 0
## 110 0 0 1128 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0
## 71426_C 71388_C 71432_C 23920_C 71436_C 71409_C 71430_C 71384_C 71437_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 7689_C 71414_C 71502_C 71486_C 71396_C 71412_C 71468_C 25087_C 20939_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71435_C 71424_C 71439_C 71390_C 71539_C 25491_C 71505_C C3n 71496_C 71389_C
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0
## 71434_C 71465_C 71415_C 71450_C 71370_C 5436_C 71504_C 71489_C 71548_C
## 107 0 0 0 0 18 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71503_C 44082_C 71529_C 71530_C 71451_C 71433_C 71410_C 3461_C 71452_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71487_C 71454_C 71449_C 71544_C 71537_C 71400_C 71495_C 71457_C 71461_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71455_C 71533_C C1c 71462_C 71470_C 45080_C 71484_C 71485_C 71538_C 22650_C
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0
## 71545_C 71463_C 71478_C 71534_C 71481_C 71479_C 71438_C 71456_C 71492_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71542_C 44826_C 71464_C 71416_C 71401_C 71540_C 71480_C 71476_C 71493_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71497_C 71494_C 71532_C 71398_C 71531_C 71541_C 23413_C 71453_C C3.7
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71466_C 71482_C 71488_C 71405_C 71467_C 71549_C 71425_C 71498_C 71383_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## C3ad 71546_C 71472_C 19660_C 1813_C 71483_C 892_C 71421_C 9930_C 71419_C
## 107 0 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0 0
## 71447_C 71417_C 22245_C 71473_C 10133_C 71535_C 71429_C 71378_C 71385_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 842 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71418_C 71471_F 71536_G
## 107 0 0 0
## 108 0 0 0
## 109 0 0 0
## 110 0 0 0
## 111 0 0 0
## 112 0 0 0Purging outlier sequences and normalizing reads
Here we remove low abundance (< 0.01%) sequences and normalize sequence counts with weighted trimmed mean of M-values (TMM; Robinson and Oshlack 2010). This helps to account for disparity in sequencing depth across libraries.
First we purge sequences and transpose the data to work with edgeR
goods = purgeOutliers(its2Seq, count.columns = 4:length(its2Seq), otu.cut = 0.0001, sampleZcut = -5)## [1] "samples with counts below z-score -5 :"
## character(0)
## [1] "zscores:"
## named numeric(0)
## [1] "OTUs passing frequency cutoff 0.0001 : 81"its2SeqTransposed = t(goods[, 4:length(goods[1, ])])
its2SeqList = DGEList(counts = its2SeqTransposed)
head(its2SeqList$samples)Now we can use TMM normalization in edgeR
its2SeqNorm = calcNormFactors(its2SeqList, method = "TMM")
head(its2SeqNorm$samples)
its2TMM = t(cpm(its2SeqNorm, normalized.lib.sizes = TRUE))
its2SeqNorm = cbind(its2Seq[,c(1:3)], its2TMM)
colnames(its2SeqNorm)[4] = "Clade C"head(its2SeqNorm)## sample_name Sample_site Depth_zone Clade C C3 C3de C3bb
## 107 CBC_107 TR 10 44129.25 374017.0 195985.8 118759.60
## 108 CBC_108 TR 10 63533.52 373234.2 122420.0 57301.19
## 109 CBC_109 TR 10 146557.87 1009931.5 168920.9 94909.87
## 110 CBC_110 TR 10 55295.48 244003.4 174973.2 82076.34
## 111 CBC_111 TR 10 52913.45 242239.4 185088.2 89297.52
## 112 CBC_112 TR 10 57065.48 283770.0 185972.0 73717.90
## C21ae C3an C3s C3dk C3cd C3dm 71372_C 71371_C
## 107 44778.21 101021.37 24660.46 22280.94 0.00 0 8003.834 5624.316
## 108 58882.46 36066.87 28055.60 22626.00 0.00 0 12630.061 7212.567
## 109 56972.54 80001.16 38070.44 29817.41 0.00 0 8519.259 8918.599
## 110 73813.38 41359.92 36593.14 30732.01 19030.24 0 14148.691 10238.537
## 111 75120.88 46221.11 37064.89 33554.30 0.00 0 11907.262 8192.900
## 112 64253.57 43048.70 24399.59 25118.40 0.00 0 10582.472 9903.596
## 71374_C C3dn 71376_C 12353_C C21 71394_C 71377_C C3b
## 107 9085.434 0.000 6489.595 0.00 24660.46 0 0 0.000
## 108 5998.371 0.000 0.000 15586.89 0.00 0 0 7765.208
## 109 8785.486 0.000 0.000 54709.62 13178.23 0 0 29551.179
## 110 8873.672 9439.292 0.000 0.00 0.00 0 0 0.000
## 111 7484.300 0.000 7132.316 0.00 0.00 0 0 0.000
## 112 8346.176 0.000 5271.269 0.00 0.00 0 0 0.000
## 71442_C 71459_C 71460_C 71458_C 37807_C 71397_C 71431_C 71444_C
## 107 0 0.000 0.000 0.000 0.000 0.000 7354.875 0
## 108 0 4796.277 5082.682 3787.808 0.000 0.000 0.000 0
## 109 0 7853.692 7454.351 9184.826 0.000 0.000 0.000 0
## 110 0 3926.549 4750.386 5283.216 0.000 0.000 0.000 0
## 111 0 0.000 0.000 0.000 0.000 6812.751 0.000 0
## 112 0 0.000 3953.452 3833.650 4272.923 6509.218 0.000 0
## 71392_C 2301_C 71375_C 71386_C 71387_C 71440_C 23653_C 71395_C 71445_C
## 107 0 0 0 0 0 0 0.000 0 0
## 108 0 0 0 0 0 0 0.000 0 0
## 109 0 0 0 0 0 0 7454.351 0 0
## 110 0 0 0 0 0 0 0.000 0 0
## 111 0 0 0 0 0 0 0.000 0 0
## 112 0 0 0 0 0 0 0.000 0 0
## 71413_C 71391_C 71406_C 71443_C 11189_C C1 71379_C 71393_C 71469_C
## 107 0.000 0 0.00 0 0.000 0 0 0 3677.437
## 108 0.000 0 0.00 0 3549.809 0 0 0 0.000
## 109 0.000 0 10515.96 0 13178.229 0 0 0 0.000
## 110 0.000 0 0.00 0 0.000 0 0 0 0.000
## 111 6113.414 0 0.00 0 0.000 0 0 0 5270.503
## 112 4392.724 0 0.00 0 0.000 0 0 0 0.000
## 71411_C 71381_C 71382_C 49051_C 71380_C 25182_C 3031_C 71477_C 71543_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 20889_C 71373_C 71423_C 71422_C 71428_C C3i 71404_C 71446_C 71426_C
## 107 0.000 0 0 0 0 0 0 0 0
## 108 4889.056 0 0 0 0 0 0 0 0
## 109 9983.506 0 0 0 0 0 0 0 0
## 110 4623.327 0 0 0 0 0 0 0 0
## 111 0.000 0 0 0 0 0 0 0 0
## 112 0.000 0 0 0 0 0 0 0 0
## 71388_C 71432_C 23920_C 71436_C 71409_C 71430_C 71384_C 71437_C 7689_C
## 107 0 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0 0
## 71502_C 71486_C 71412_C 71435_C 71439_C 71505_C 71434_C 71548_C
## 107 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0Calculation of ITS2 sequence relative abundances
Now we can calculate the relative abundance of each ITS2 sequences per sample. This will allow us to view the assemblages in a faceted barplot.
colOrder = order(colSums(its2SeqNorm[4:length(its2SeqNorm[1,])]), decreasing = FALSE) + 3
its2SeqPerc = cbind(its2SeqNorm[,c(1:3)], its2SeqNorm[,c(colOrder)])
its2SeqPerc$sum = apply(its2SeqPerc[, c(4:length(its2SeqPerc[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
its2SeqPerc = cbind(its2SeqPerc[, c(1:3)], (its2SeqPerc[, c(4:(ncol(its2SeqPerc)-1))] / its2SeqPerc$sum))
Now a quick sanity check. If this worked the sum of each row should = 100% (i.e. “1”).
apply(its2SeqPerc[, c(4:(ncol(its2SeqPerc)))], 1, function(x) {
sum(x, na.rm = T)
})## 107 108 109 110 111 112 113 114 116 117 118 119 120 121 156 93 94 95 96 97
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 98 99 100 101 102 103 104 105 106 137 122 123 124 125 126 127 128 129 130 131
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 132 133 134 135 136 226 227 228 229 230 231 232 233 234 235 236 237 238 239 29
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 30 31 32 33 34 35 36 37 38 39 40 41 42 43 17 18 19 20 21 22
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 23 24 25 26 27 28 44 45 46 1 2 3 4 5 6 7 8 9 10 11
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 12 13 14 15 16 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 47 48 49 50
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 51 52 53 54 55 56 57 58 59 60 92 61 62 63 64 65 66 67 68 69
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 70 71 72 73 74 75 198 199 200 201 202 203 204 205 206 207 208 209 210 240
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 241 115 152 153 154 155 157 158 159 160 161 162 163 164 165 166 148 149 150 151
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 176 177 178 179 180 181 182 183 187 188 189 145 146 147 167 168 169 170 171 172
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 173 174 175 184 185 186 138 139 140 141 142 143 144 190 191 192 193 194 195 196
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 197
## 1Everything adds up to 1, this is good! The code works.
I added an additional column to sort better for the stacked barplot. This was just a work around to get the facet_grid() function to play nice with our data. I added a coulumn “barPlotOrder” and for each population I filled in a series 1:n foreach sample at each Site:Depth combo, so now there’s no large blank expanses on the plot.
sampleCounts = plyr::count(its2SeqPerc, c('Sample_site','Depth_zone'))
meltedList = melt(lapply(sampleCounts$freq,function(x){c(1:x)}))
its2SeqPerc$barPlotOrder = meltedList$value
its2SeqPerc=its2SeqPerc[c(1,ncol(its2SeqPerc),2:(ncol(its2SeqPerc)-1))]
head(its2SeqPerc)
gssSeq = otuStack(its2SeqPerc, count.columns = c(5:length(its2SeqPerc[1, ])),
condition.columns = c(1:4))[1:19521,] # remove summ rows
levels(gssSeq$otu)
levels(gssSeq$Depth_zone)
levels(gssSeq$Depth_zone) = c("10 m", "16 m", "25 m", "35 m")
levels(gssSeq$Sample_site)
levels(gssSeq$Sample_site) = c("Tobacco Reef", "Raph's Wall", "South Reef", "Glover's Reef")
levels(gssSeq$Depth_zone)
levels(gssSeq$Sample_site)Consruct ITS2 sequence barplot
colorCount = length(c(5:length(its2SeqPerc[1,])))
getPalette = colorRampPalette(redmonder.pal(8, "qPBI"), bias = 1.7)
its2SeqPlotA = ggplot(gssSeq, aes(x = barPlotOrder, y = count, fill = factor(otu))) +
geom_bar(position = "stack", stat = "identity", color = "black",
size = 0.25) +
ylab("Proportion") +
scale_fill_manual(values=rev(getPalette(colorCount)))+
labs(fill = expression(paste(italic("ITS2"), " sequence"))) +
guides(fill = guide_legend(ncol = 9, reverse = TRUE)) +
facet_grid(Depth_zone ~ Sample_site, scales = "free_x") + #faceting plots by Depth and Site
theme_bw()
its2SeqPlot = its2SeqPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 12),
legend.position = "bottom",
legend.justification = "left",
legend.direction = "horizontal",
legend.title = element_text(color = "black", size = 12, hjust = 0.5, angle = 90),
legend.text = element_text(color = "black", size = 10),
legend.key = element_blank(),
legend.key.size = unit(0.4,"line"),
legend.background = element_blank(),
panel.border = element_blank(),
panel.background = element_rect(fill = "white"),
plot.background = element_blank(),
strip.text.x = element_text(size = 12),
strip.text.y = element_text(size = 12),
strip.background = element_rect(fill = "white", size = 0.9)
)
its2SeqPlotSave the barplot
ggsave("../figures/Fig3.eps", plot = its2SeqPlot, width = 8.5, height = 8.5, unit = "in", dpi = 600)
ITS2 type profiles
We can now look at the ITS2 type profiles predicted by SymPortal.
Prepare ITS2 type profile data
Similar to what we did with the sequence data
its2Profs = read.delim("62_20190310_DBV_2019-03-11_01-11-25.167036.profiles.absolute.clean.txt", header = TRUE, check.names = FALSE)
head(its2Profs)
its2Profs = cbind(its2Profs[1], its2MetaData[,2:3], its2Profs[,c(2:length(its2Profs))])
colnames(its2Profs)[3] = "Depth_zone"
head(its2Profs)
its2Profs$Depth_zone = factor(its2Profs$Depth_zone, levels = c("10", "16", "25", "35"))
levels(its2Profs$Depth_zone)
its2Profs$Sample_site = factor(its2Profs$Sample_site, levels(its2Profs$Sample_site)[c(4, 2, 3, 1)])
its2Profs = its2Profs[order(its2Profs$Sample_site, its2Profs$Depth_zone), ]
head(its2Profs)
sampleCounts = plyr::count(its2Profs, c('Sample_site','Depth_zone'))
meltedList = reshape2::melt(lapply(sampleCounts$freq,function(x){c(1:x)}))
its2Profs$barPlotOrder = meltedList$value
its2Profs=its2Profs[c(1,ncol(its2Profs),2:(ncol(its2Profs)-1))]head(its2Profs)## Sample barPlotOrder Sample_site Depth_zone A4/A4a B18b
## 107 CBC_107 1 TR 10 0 0
## 108 CBC_108 2 TR 10 0 0
## 109 CBC_109 3 TR 10 0 0
## 110 CBC_110 4 TR 10 0 0
## 111 CBC_111 5 TR 10 0 0
## 112 CBC_112 6 TR 10 0 0
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk
## 107 4075 0
## 108 173180 0
## 109 11108 0
## 110 166773 0
## 111 152997 0
## 112 17536 0
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an C21 C3-C3b C3cd/C3-C3dn C3cd
## 107 0 0 0 0 0
## 108 0 0 0 0 0
## 109 0 0 0 0 0
## 110 0 0 0 0 0
## 111 0 0 0 0 0
## 112 0 0 0 0 0
## C3-C3de-C3an-C3bb-C21ae-C3s C3dm C1-C21 C3
## 107 0 0 0 0
## 108 0 0 0 0
## 109 0 0 0 0
## 110 0 0 0 0
## 111 0 0 0 0
## 112 0 0 0 0Normalization of ITS2 type profile reads
Similar to what we did with the ITS2 sequences, but we won’t purge any low abundance reads, since we only have 13 total ITS2 type profiles.
its2ProfsTransposed = t(its2Profs[, 5:length(its2Profs[1, ])])
its2ProfsList = DGEList(counts = its2ProfsTransposed)
head(its2ProfsList$samples)
its2ProfsNorm = calcNormFactors(its2ProfsList, method = "TMM")
head(its2ProfsNorm$samples)
its2TMM = t(cpm(its2ProfsNorm, normalized.lib.sizes = TRUE))
its2ProfsNorm = cbind(its2Profs[,c(1:4)], its2TMM)head(its2ProfsNorm)## Sample barPlotOrder Sample_site Depth_zone A4/A4a B18b
## 107 CBC_107 1 TR 10 0 0
## 108 CBC_108 2 TR 10 0 0
## 109 CBC_109 3 TR 10 0 0
## 110 CBC_110 4 TR 10 0 0
## 111 CBC_111 5 TR 10 0 0
## 112 CBC_112 6 TR 10 0 0
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk
## 107 994029.3 0
## 108 994029.3 0
## 109 994029.3 0
## 110 994029.3 0
## 111 994029.3 0
## 112 994029.3 0
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an C21 C3-C3b C3cd/C3-C3dn C3cd
## 107 0 0 0 0 0
## 108 0 0 0 0 0
## 109 0 0 0 0 0
## 110 0 0 0 0 0
## 111 0 0 0 0 0
## 112 0 0 0 0 0
## C3-C3de-C3an-C3bb-C21ae-C3s C3dm C1-C21 C3
## 107 0 0 0 0
## 108 0 0 0 0
## 109 0 0 0 0
## 110 0 0 0 0
## 111 0 0 0 0
## 112 0 0 0 0Preparing ITS2 type profiles for plotting
colOrder2 = order(colSums(its2ProfsNorm[5:length(its2ProfsNorm[1,])]), decreasing = TRUE) + 4
its2ProfsPerc = cbind(its2ProfsNorm[,c(1:4)],its2ProfsNorm[,c(colOrder2)])
its2ProfsPerc$sum = apply(its2ProfsPerc[, c(5:length(its2ProfsPerc[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
its2ProfsPerc = cbind(its2ProfsPerc[, c(1:4)], (its2ProfsPerc[, c(5:(ncol(its2ProfsPerc)-1))]
/ its2ProfsPerc$sum))
head(its2ProfsPerc)## Sample barPlotOrder Sample_site Depth_zone
## 107 CBC_107 1 TR 10
## 108 CBC_108 2 TR 10
## 109 CBC_109 3 TR 10
## 110 CBC_110 4 TR 10
## 111 CBC_111 5 TR 10
## 112 CBC_112 6 TR 10
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk
## 107 1 0
## 108 1 0
## 109 1 0
## 110 1 0
## 111 1 0
## 112 1 0
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an C3dm C3-C3de-C3an-C3bb-C21ae-C3s
## 107 0 0 0
## 108 0 0 0
## 109 0 0 0
## 110 0 0 0
## 111 0 0 0
## 112 0 0 0
## C3cd/C3-C3dn C3-C3b C3cd C3 C21 C1-C21 A4/A4a B18b
## 107 0 0 0 0 0 0 0 0
## 108 0 0 0 0 0 0 0 0
## 109 0 0 0 0 0 0 0 0
## 110 0 0 0 0 0 0 0 0
## 111 0 0 0 0 0 0 0 0
## 112 0 0 0 0 0 0 0 0apply(its2ProfsPerc[, c(5:(ncol(its2ProfsPerc)))], 1, function(x) {
sum(x, na.rm = T)
})## 107 108 109 110 111 112 113 114 116 117 118 119 120 121 156 93 94 95 96 97
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 98 99 100 101 102 103 104 105 106 137 122 123 124 125 126 127 128 129 130 131
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 132 133 134 135 136 226 227 228 229 230 231 232 233 234 235 236 237 238 239 29
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 30 31 32 33 34 35 36 37 38 39 40 41 42 43 17 18 19 20 21 22
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 23 24 25 26 27 28 44 45 46 1 2 3 4 5 6 7 8 9 10 11
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 12 13 14 15 16 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 47 48 49 50
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 51 52 53 54 55 56 57 58 59 60 92 61 62 63 64 65 66 67 68 69
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 70 71 72 73 74 75 198 199 200 201 202 203 204 205 206 207 208 209 210 240
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 241 115 152 153 154 155 157 158 159 160 161 162 163 164 165 166 148 149 150 151
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 176 177 178 179 180 181 182 183 187 188 189 145 146 147 167 168 169 170 171 172
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 173 174 175 184 185 186 138 139 140 141 142 143 144 190 191 192 193 194 195 196
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 197
## 1Again, everything looks good!
gssProf = otuStack(its2ProfsPerc, count.columns = c(5:length(its2ProfsPerc[1, ])),
condition.columns = c(1:4))[1:3133,] # remove summ rows
levels(gssProf$otu)
levels(gssProf$Depth_zone)
levels(gssProf$Depth_zone) = c("10 m", "16 m", "25 m", "35 m")
levels(gssProf$Sample_site)
levels(gssProf$Sample_site) = c("Tobacco Reef", "Raph's Wall", "South Reef", "Glover's Reef")
levels(gssProf$Depth_zone)
levels(gssProf$Sample_site)Consruct ITS2 type profile barplot
colorCount2 = length(c(4:length(its2ProfsPerc[1,]))) +1
getPalette2 = colorRampPalette(redmonder.pal(8, "qPBI"))
its2ProfsPlotA = ggplot(gssProf, aes(x = barPlotOrder, y = count, fill = factor(otu))) +
geom_bar(position = "stack", stat = "identity", color = "black", size = 0.25) +
ylab("Proportion") +
scale_fill_manual(values = c(getPalette2(colorCount2)[2:12],
carto_pal(n = 7, "Tropic")[c(6,7)]))+
labs(fill = expression(paste(italic("ITS2"), " type profile"))) +
guides(fill = guide_legend(ncol = 2, reverse = FALSE)) +
facet_grid(Depth_zone ~ Sample_site, scales = "free_x") + #faceting plots by Depth and Site
theme_bw()
its2ProfsPlot = its2ProfsPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 12),
legend.position = "bottom",
legend.title = element_text(color = "black", size = 12, hjust = 0.5, angle = 90),
legend.text = element_text(color = "black", size = 10),
legend.key = element_blank(),
legend.key.size = unit(0.75,"line"),
legend.background = element_blank(),
panel.border = element_blank(),
panel.background = element_rect(fill = "white"),
plot.background = element_blank(),
strip.text.x = element_text(size = 12),
strip.text.y = element_text(size = 12),
strip.background = element_rect(fill = "white", size = 0.9)
)
its2ProfsPlotSave ITS2 type profile barplot
ggsave("../figures/Fig4.eps", plot = its2ProfsPlot, width = 8.5, height = 8.5, unit = "in", dpi = 600)
Statistical analyses on ITS2 type profiles
Cheking dispersion with PERMDISP
Using betadisper() in vegan to look at multivariate homogeneity of dispersion (PERMDISP) between sites and depths. This is using Bray-Curtis dissimilarity.
set.seed(694) #setting seed allows repetition of randomized processes
its2dispS = betadisper(vegdist(its2ProfsNorm[, c(5:ncol(its2ProfsNorm))]), its2ProfsNorm$Sample_site)
anova(its2dispS)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.283 0.094301 0.5837 0.6263
## Residuals 237 38.291 0.161566No significant effect of Site on betadiversity.
set.seed(694)
its2dispD = betadisper(vegdist(its2ProfsNorm[, c(6:ncol(its2ProfsNorm))]), its2ProfsNorm$Depth_zone)
anova(its2dispD)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 4.926 1.64192 11.565 0.0000004215 ***
## Residuals 237 33.648 0.14198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depth does significantly affect beta diversity.
Permutation test for pairwise comparisons
Follow up with a permutation test to see where differences occur.
set.seed(694)
its2PermTest = permutest(its2dispD, permutations = 9999, pairwise = T, model = "full", )
its2PermTest##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 4.926 1.64192 11.565 9999 0.0001 ***
## Residuals 237 33.648 0.14198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## 10 16 25 35
## 10 0.887100000 0.000400000 0.0001
## 16 0.890151321 0.000100000 0.0001
## 25 0.000118583 0.000078099 0.7267
## 35 0.000045284 0.000029815 0.733416414its2PermTest$statistic## Overall (F) 10-16 (t) 10-25 (t) 10-35 (t) 16-25 (t) 16-35 (t)
## 11.5648579 -0.1384088 3.9796639 4.2362895 4.0923225 4.3451550
## 25-35 (t)
## 0.3413910Running PERMANOVA in R
Now let’s see how different communities are from each other with PERMANOVA. We will utilize the adonis() function in vegan. We will use Bray-Curtis similarity for our distance matrix and run a total 0f 9,99999 permutations, and test the effects of Site, Depth, and the interaction between Site and Depth. Dispersion is heteroschedastic, but PERMANOVA is robust to deviations in homgeneity of variance ( Anderson and Walsh, 2013)
set.seed(694)
its2Adonis = adonis(its2ProfsNorm[, c(5:ncol(its2ProfsNorm))] ~ Depth_zone * Sample_site,
data = its2ProfsNorm, permutations = 9999, method = "bray")
its2Adonis##
## Call:
## adonis(formula = its2ProfsNorm[, c(5:ncol(its2ProfsNorm))] ~ Depth_zone * Sample_site, data = its2ProfsNorm, permutations = 9999, method = "bray")
##
## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Depth_zone 3 3.331 1.11039 6.8801 0.07925 0.0001 ***
## Sample_site 3 0.705 0.23497 1.4559 0.01677 0.1711
## Depth_zone:Sample_site 9 1.687 0.18740 1.1611 0.04012 0.2675
## Residuals 225 36.313 0.16139 0.86386
## Total 240 42.036 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PERMANOVA reveals that Depth has a signifcant effect on Symbiodiniaceae associations in our M. cavernosa samples.
Pairwise PERMANOVA for multiple comparisons
Since we found that Depth was a significant factor in our PERMANOVA we can now use pairwise PERMANOVA to reveal where differences occur across depth. This utilizes the package pairwiseAdonis, where we will again use Bray-Curtis similarity and 9,99999 permutations. We also have added false discovery rate (FDR) corrections since we are perfoming multiple comparisons.
set.seed(694)
its2PWAdonis = pairwise.adonis(its2ProfsNorm[, c(5:ncol(its2ProfsNorm))],
factors = its2ProfsNorm$Depth_zone,
sim.method = "bray", p.adjust.m = "BH", perm = 9999)
its2PWAdonis## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 10 vs 16 1 0.264839474 0.9837600 0.008199110 0.3709 0.44508
## 2 10 vs 25 1 1.491477086 9.0843560 0.070375344 0.0005 0.00075 **
## 3 10 vs 35 1 1.597044155 10.0260431 0.078312528 0.0002 0.00060 **
## 4 16 vs 25 1 1.573437103 9.4051037 0.073245560 0.0004 0.00075 **
## 5 16 vs 35 1 1.730759227 10.6563156 0.083476603 0.0001 0.00060 **
## 6 25 vs 35 1 0.006758822 0.1196029 0.001012558 1.0000 1.00000We see that again see differences between our deeper (25 + 35 m) and shallower (10 + 16 m) samples.
PERMANOVA without 35 m samples
First we need to remove the deep samples from the dataframe. We will use our dataframe of profiles that haven’t been normalized yet. This way we can calculate the normalization based on only the samples we are keeping in the analysis.
its2Profs2 = subset(its2Profs, !Depth_zone=="35")
its2Profs2[] = lapply(its2Profs2, function(x) if(is.factor(x)) factor(x) else x)
summary(its2Profs2)head(its2Profs2)## Sample barPlotOrder Sample_site Depth_zone A4/A4a B18b
## 107 CBC_107 1 TR 10 0 0
## 108 CBC_108 2 TR 10 0 0
## 109 CBC_109 3 TR 10 0 0
## 110 CBC_110 4 TR 10 0 0
## 111 CBC_111 5 TR 10 0 0
## 112 CBC_112 6 TR 10 0 0
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk
## 107 4075 0
## 108 173180 0
## 109 11108 0
## 110 166773 0
## 111 152997 0
## 112 17536 0
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an C21 C3-C3b C3cd/C3-C3dn C3cd
## 107 0 0 0 0 0
## 108 0 0 0 0 0
## 109 0 0 0 0 0
## 110 0 0 0 0 0
## 111 0 0 0 0 0
## 112 0 0 0 0 0
## C3-C3de-C3an-C3bb-C21ae-C3s C3dm C1-C21 C3
## 107 0 0 0 0
## 108 0 0 0 0
## 109 0 0 0 0
## 110 0 0 0 0
## 111 0 0 0 0
## 112 0 0 0 0
Normalize samples again, as above.
its2ProfsTransposed2 = t(its2Profs2[, 5:length(its2Profs2[1, ])])
its2ProfsList2 = DGEList(counts = its2ProfsTransposed2)
head(its2ProfsList2$samples)
its2ProfsNorm2 = calcNormFactors(its2ProfsList2, method = "TMM")
head(its2ProfsNorm2$samples)
its2TMM2 = t(cpm(its2ProfsNorm2, normalized.lib.sizes = TRUE))
its2ProfsNorm2 = cbind(its2Profs2[,c(1:4)], its2TMM2)head(its2ProfsNorm2)## Sample barPlotOrder Sample_site Depth_zone A4/A4a B18b
## 107 CBC_107 1 TR 10 0 0
## 108 CBC_108 2 TR 10 0 0
## 109 CBC_109 3 TR 10 0 0
## 110 CBC_110 4 TR 10 0 0
## 111 CBC_111 5 TR 10 0 0
## 112 CBC_112 6 TR 10 0 0
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk
## 107 992101.4 0
## 108 992101.4 0
## 109 992101.4 0
## 110 992101.4 0
## 111 992101.4 0
## 112 992101.4 0
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an C21 C3-C3b C3cd/C3-C3dn C3cd
## 107 0 0 0 0 0
## 108 0 0 0 0 0
## 109 0 0 0 0 0
## 110 0 0 0 0 0
## 111 0 0 0 0 0
## 112 0 0 0 0 0
## C3-C3de-C3an-C3bb-C21ae-C3s C3dm C1-C21 C3
## 107 0 0 0 0
## 108 0 0 0 0
## 109 0 0 0 0
## 110 0 0 0 0
## 111 0 0 0 0
## 112 0 0 0 0
Run PERMANOVA on normalized subset of data
set.seed(694)
its2AdonisNo35 = adonis(its2ProfsNorm2[, c(5:ncol(its2ProfsNorm2))] ~ Depth_zone * Sample_site,
data = its2ProfsNorm2, permutations = 9999, method = "bray")
its2AdonisNo35##
## Call:
## adonis(formula = its2ProfsNorm2[, c(5:ncol(its2ProfsNorm2))] ~ Depth_zone * Sample_site, data = its2ProfsNorm2, permutations = 9999, method = "bray")
##
## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Depth_zone 2 2.222 1.11097 5.6348 0.05840 0.0005 ***
## Sample_site 3 0.873 0.29109 1.4764 0.02295 0.1576
## Depth_zone:Sample_site 6 1.432 0.23869 1.2106 0.03764 0.2537
## Residuals 170 33.518 0.19716 0.88100
## Total 181 38.045 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depth is still has a significant effect on community structure
Let’ see where differences lie across depth using pairwise PERMANOVA again.
set.seed(694)
its2PWAdonisNo35 = pairwise.adonis(its2ProfsNorm2[, c(5:ncol(its2ProfsNorm2))],
factors = its2ProfsNorm2$Depth_zone,
sim.method = "bray", p.adjust.m = "BH", perm = 9999)
its2PWAdonisNo35## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 10 vs 16 1 0.2648395 0.983760 0.00819911 0.3709 0.37090
## 2 10 vs 25 1 1.4914771 9.084356 0.07037534 0.0005 0.00075 **
## 3 16 vs 25 1 1.5734371 9.405104 0.07324556 0.0002 0.00060 **The differences are still between 25 m and the shallower sites (10 + 16 m). This gives us confidence that temporal differences are not a cofactor with depth.
SIMPER test between deep and shallow populations
Similarity percentage test (SIMPER) will show us the ITS2 type profiles that contribute the most to the dissimilarity between depth zones
depths = its2ProfsNorm$Depth_zone
levels(depths) = c("Shallow","Shallow","Deep","Deep")
depths
its2ProfsNorm$depth = depths
its2ProfsNorm = cbind(its2ProfsNorm[c(1:4)],its2ProfsNorm[ncol(its2ProfsNorm)],its2ProfsNorm[c(5:(ncol(its2ProfsNorm)-1))])
head(its2ProfsNorm)its2SimperD = simper(sqrt(its2ProfsNorm[, c(6:ncol(its2ProfsNorm))]), its2ProfsNorm$depth)
summary(its2SimperD)##
## Contrast: Shallow_Deep
##
## average sd ratio ava
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk 0.1739207 0.233859 0.74370 659.1803
## C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk 0.0735606 0.176809 0.41604 115.3565
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an 0.0716481 0.174440 0.41073 131.6755
## C3-C3de-C3an-C3bb-C21ae-C3s 0.0123967 0.077750 0.15944 24.7193
## C3cd/C3-C3dn 0.0101081 0.063408 0.15941 23.0803
## C3-C3b 0.0087577 0.059117 0.14814 19.4279
## C21 0.0060702 0.033354 0.18199 14.2738
## C3dm 0.0056642 0.044231 0.12806 20.3637
## C3cd 0.0053058 0.035963 0.14753 15.2092
## C3 0.0047578 0.037056 0.12839 11.5608
## C1-C21 0.0019356 0.015391 0.12576 4.3170
## A4/A4a 0.0002471 0.002707 0.09128 0.6077
## B18b 0.0001825 0.001999 0.09128 0.3718
## avb cumsum
## C3-C3de-C3bb-C21ae-C3an-C3s-C3dk 938.85 0.4643
## C3/C3dm-C3de-C3bb-C21ae-C3cd-C3an-C3s-C3dk 41.54 0.6607
## C3/C3cd-C3de-C3dn-C3bb-C21ae-C3an 16.62 0.8520
## C3-C3de-C3an-C3bb-C21ae-C3s 0.00 0.8851
## C3cd/C3-C3dn 0.00 0.9121
## C3-C3b 0.00 0.9355
## C21 0.00 0.9517
## C3dm 0.00 0.9668
## C3cd 0.00 0.9810
## C3 0.00 0.9937
## C1-C21 0.00 0.9989
## A4/A4a 0.00 0.9995
## B18b 0.00 1.0000
## Permutation: free
## Number of permutations: 0