Hello everyone !!, Does anyone know how I could analyze the adonis result between months (3 months, July August and September) (PAIRWAISE). The way I am writing this function it gives me the general result of all the months between them. I would like to know the result for months, I mean, July vs. August … July vs. September … and September vs. August.
I try an example I found in internet, but i have an error.
Thank you very much in advance.
Family <-tibble::tribble(
~SampleID, ~K00001, ~K00002, ~K00003, ~K00004, ~K00005, ~K00006,
"1A", 13.0336555, 3.653006048, 13.04832335, 8.669081695, 6.041009978, 0,
"1B", 12.13227078, 6.680420722, 12.30629898, 7.718240989, 8.972651202, 0,
"1C", 12.78444157, 7.824621007, 13.78559758, 9.426997918, 7.568067839, 0,
"1D", 12.93280152, 5.388255506, 13.17357399, 8.668526973, 8.58688218, 0,
"1E", 13.08432509, 5.939290059, 13.05932398, 11.18908222, 8.446807541, 0,
"1F", 12.47032632, 6.970425447, 12.66560006, 9.160178569, 8.586800547, 4.819935626,
"1G", 11.33655907, 5.842294244, 12.09470987, 8.685505314, 9.120690025, 0,
"1H", 12.30842204, 7.56865924, 12.38609878, 7.53779754, 9.444540341, 0,
"2A", 12.25561947, 7.076298974, 12.61309663, 9.081991185, 9.887316299, 2.378671111,
"2B", 11.85658338, 7.653630409, 12.52657494, 8.830352407, 10.08033857, 1.70473684,
"2C", 11.95815429, 5.90722311, 12.36858022, 8.191704706, 9.404747474, 0,
"2D", 12.87295808, 6.181453869, 12.95347313, 7.986927795, 9.340857256, 3.274422907,
"2E", 11.94014957, 5.397316451, 12.7671218, 10.25986242, 8.457825932, 1.110387432,
"2F", 11.79951691, 5.195683999, 12.30144343, 9.435310077, 8.36402902, 0,
"2G", 12.00707972, 6.803597144, 12.71764902, 7.463750858, 9.396296995, 1.10248341,
"2H", 12.5255573, 6.932400736, 12.29855001, 9.675249175, 8.920238803, 3.028271518,
"3A", 12.45882927, 7.611950315, 13.19517448, 8.54723483, 8.871506326, 0,
"3B", 12.00341574, 6.636813413, 13.57137712, 9.476711646, 12.80345854, 0,
"3C", 12.36333134, 4.779227864, 12.56892439, 9.050359862, 5.437461793, 0,
"3D", 12.49152829, 3.881814891, 12.77408748, 8.722490226, 6.851908854, 0,
"3E", 12.39599087, 4.663595488, 12.60508686, 10.60246315, 8.609217558, 0,
"3F", 12.14556933, 6.660872819, 12.48013856, 8.743426142, 7.889679545, 0,
"3G", 12.28639492, 7.117319696, 12.50821694, 8.792408539, 8.878452135, 0,
"3H", 12.05141242, 0, 12.12624923, 4.706533368, 6.254181176, 0,
"4A", 12.39411802, 7.363030556, 12.92729675, 9.415335942, 10.01629272, 1.899641342,
"4B", 11.72855753, 6.221484527, 12.06265014, 8.147175269, 9.295721831, 0,
"4C", 11.85159423, 4.814870357, 12.37250077, 8.016191473, 9.101592483, 0,
"4D", 12.36932793, 6.746061065, 12.70038007, 8.94971237, 9.059720765, 3.368514377,
"4E", 12.28833533, 6.314508898, 12.56820378, 9.216232833, 9.167274079, 0,
"4F", 12.48315568, 0, 13.13202175, 6.642548082, 6.642548082, 0,
"4G", 11.73677502, 6.99428693, 13.05790619, 9.034167857, 9.078354963, 0,
"4H", 12.65268079, 6.785221536, 12.75793551, 8.428070266, 8.812131529, 2.969183234,
"5A", 12.63665764, 5.935521954, 13.28677307, 11.18397148, 8.991509376, 1.951137284,
"5B", 12.9991414, 10.30021555, 13.11517142, 9.148430273, 7.556732093, 0,
"5C", 12.57027162, 5.473590602, 12.37850125, 6.177036375, 7.070435329, 0,
"5D", 11.79382425, 5.540972902, 11.99192094, 8.817703989, 8.415103078, 0,
"5E", 12.03346216, 5.732678664, 13.02960547, 8.746581504, 7.193830755, 0,
"5F", 12.50034004, 6.182918754, 12.56139466, 8.003014868, 8.027762921, 2.443230677,
"5G", 12.30434541, 5.794658392, 12.62777046, 10.44189419, 8.607324415, 0,
"5H", 12.36028824, 6.392542036, 13.41421241, 10.58265869, 8.901848002, 0,
"6A", 12.07518393, 6.435691102, 12.65088245, 8.756632411, 9.519701447, 0,
"6B", 12.43066661, 6.662331413, 12.54069806, 9.347405718, 8.221225728, 0,
"6C", 12.15101277, 6.498225504, 12.70122853, 8.616561285, 9.255914589, 0,
"6D", 12.39475709, 6.66365064, 12.48544659, 8.689841666, 8.698921901, 0,
"6E", 12.3668262, 5.670789363, 12.75953042, 9.422226064, 9.143546483, 0,
"6F", 12.10031439, 6.49757301, 12.44068945, 9.032618798, 8.866041842, 4.654607502,
"6G", 12.29657346, 4.89164898, 13.4304044, 10.92116395, 10.04963159, 1.805013449,
"6H", 11.75859504, 6.335989638, 12.26422847, 9.430388899, 9.200432438, 1.347210803,
"7A", 13.10225898, 7.373853726, 13.19425085, 7.870772426, 8.572795475, 0,
"7B", 12.32636318, 10.03043578, 13.74347983, 6.780744441, 7.467614988, 0,
"7C", 11.6664146, 5.00279211, 12.9842409, 5.464068787, 8.73306237, 0,
"7D", 12.07228615, 6.845364574, 12.66794829, 8.421938557, 10.52932006, 0,
"7E", 12.50937335, 6.985645067, 12.61842509, 8.519299757, 8.894788117, 2.672016242,
"7F", 12.79339894, 6.26851676, 13.02000185, 8.262813761, 8.536659017, 0,
"7G", 12.3867073, 6.339619888, 12.53625712, 10.82896132, 9.397085458, 0,
"7H", 12.44311511, 6.810681638, 13.33144933, 10.35281178, 8.868325835, 0,
"8A", 12.08042873, 6.031668994, 12.58745259, 8.305011098, 9.424533784, 1.782478691,
"8B", 12.27036791, 6.990144415, 12.4638995, 8.528526113, 8.394702511, 1.886720728,
"8C", 11.85862564, 5.833141626, 12.18109761, 7.055536614, 8.566309467, 0,
"8D", 12.61676172, 7.196268668, 12.81073221, 9.081522398, 9.946205911, 0,
"8E", 12.18550064, 7.746967458, 12.51012043, 9.208914404, 9.286629993, 0,
"8F", 12.36429396, 6.508823649, 12.57208181, 8.252325758, 8.098781497, 3.199707024,
"8G", 12.11750417, 5.971224156, 12.49570696, 9.820972793, 9.493208103, 0,
"8H", 12.84640058, 6.247528728, 12.92866438, 8.254912934, 8.651908969, 0,
"9A", 12.68413781, 9.529080155, 13.45030264, 9.315379622, 7.0375159, 0,
"9B", 12.78379921, 9.059947548, 13.04820921, 8.447154365, 6.666105425, 0,
"9C", 12.27320298, 4.988801212, 12.87314207, 8.945563907, 7.998911726, 0,
"9D", 12.66989911, 7.944438109, 12.73206746, 9.968579163, 10.33456073, 0,
"9E", 13.2273732, 5.897944133, 12.93572598, 7.961684761, 7.603971928, 0,
"9F", 12.01597592, 7.742656418, 12.50113209, 9.234285958, 8.740734673, 3.127013175,
"9G", 11.7273935, 6.02662703, 13.45433503, 9.199789734, 8.083364299, 4.047046789,
"9H", 11.27384929, 8.291443316, 13.2031032, 9.177059498, 12.63156816, 0,
"10A", 12.15883653, 6.561980933, 12.30852365, 8.620689285, 8.590583764, 0,
"10B", 11.82490036, 5.667543281, 12.37734214, 8.901283426, 8.529130098, 0,
"10C", 12.62896778, 6.879147144, 12.85870416, 8.447478415, 9.643084661, 0,
"10D", 11.92409503, 6.217629093, 12.88171412, 7.662457613, 11.74790829, 0,
"10E", 12.22621752, 6.200638741, 12.49415335, 8.973083499, 9.185979115, 0,
"10F", 11.554951, 6.192786923, 12.15267055, 8.480753316, 8.78353641, 0,
"10G", 12.50196034, 6.025775368, 12.78712151, 8.386290523, 8.91718074, 0,
"10H", 12.22875566, 7.811085656, 12.3414332, 8.007302188, 8.559448922, 0,
"11A", 13.65923264, 5.099438908, 13.29766617, 7.937099378, 11.12705028, 0,
"11B", 12.45071573, 4.077092044, 12.87181408, 6.907933216, 7.000291564, 0,
"11C", 12.8786262, 6.903443025, 13.05641748, 7.994619685, 9.471106318, 0,
"11D", 12.65464126, 7.639517216, 12.8541468, 9.510818906, 9.530864184, 0,
"11E", 12.29092711, 4.599970653, 12.4461456, 10.00732116, 8.070479218, 2.813120291,
"11F", 11.64083914, 7.902414237, 12.77955034, 9.051583459, 6.43164937, 1.383961991,
"11G", 12.23623541, 6.563582946, 13.36986234, 9.706253321, 7.245004026, 0,
"11H", 11.94148574, 8.885978647, 13.167755, 9.835520714, 7.076415122, 0,
"12A", 11.76153901, 7.380108496, 12.35146494, 7.898313891, 8.839697031, 0,
"12B", 12.29721423, 7.375748057, 12.52955059, 8.019208437, 9.472604824, 0,
"12C", 12.205373, 6.665202752, 12.66651435, 7.460585438, 9.359881893, 0,
"12D", 12.36076381, 6.542737452, 12.5039474, 10.79811559, 8.373534187, 2.668531261,
"12E", 12.46852425, 8.129971394, 12.58057611, 9.363737086, 5.392030655, 0,
"12F", 12.2634949, 6.483940945, 12.4981267, 9.092231696, 8.844503392, 3.317648726,
"12G", 12.37979794, 6.41317005, 12.71403201, 10.57040416, 8.987126835, 2.899475144,
"12H", 11.93348787, 8.455203645, 12.36325413, 9.449850989, 9.385754427, 2.1759695
)
Metadata<- data.frame(stringsAsFactors=FALSE,
SampleID = c("1A", "1B", "1C", "1D", "1E", "1F", "1G", "1H", "2A",
"2B", "2C", "2D", "2E", "2F", "2G", "2H", "3A", "3B",
"3C", "3D", "3E", "3F", "3G", "3H", "4A", "4B", "4C",
"4D", "4E", "4F", "4G", "4H", "5A", "5B", "5C", "5D",
"5E", "5F", "5G", "5H", "6A", "6B", "6C", "6D", "6E", "6F",
"6G", "6H", "7A", "7B", "7C", "7D", "7E", "7F", "7G",
"7H", "8A", "8B", "8C", "8D", "8E", "8F", "8G", "8H",
"9A", "9B", "9C", "9D", "9E", "9F", "9G", "9H", "10A",
"10B", "10C", "10D", "10E", "10F", "10G", "10H", "11A",
"11B", "11C", "11D", "11E", "11F", "11G", "11H", "12A",
"12B", "12C", "12D", "12E", "12F", "12G", "12H"),
Month = c("July", "July", "July", "August", "August", "August",
"September", "September", "July", "July", "July",
"August", "August", "August", "September", "September",
"July", "July", "July", "August", "August", "August",
"September", "September", "July", "July", "July", "August",
"August", "August", "September", "September", "July",
"July", "July", "August", "August", "August", "September",
"September", "July", "July", "July", "August", "August",
"August", "September", "September", "July", "July",
"July", "August", "August", "September", "September",
"September", "July", "July", "July", "August", "August",
"September", "September", "September", "July", "July",
"July", "August", "August", "September", "September",
"September", "July", "July", "July", "August", "August",
"September", "September", "September", "July", "July", "July",
"August", "August", "September", "September",
"September", "July", "July", "August", "August", "August",
"September", "September", "September")
)
#Family <- t(Family)
attach(Family)
Family <- Family[,-1]
rownames(Family) <- SampleID
Metadata$Month <- factor(Metadata$Month,
levels = c("July", "August", "September"))
#remove the rare microbes. It keeps only the microbes that are present in at least 10% of the samples
dim(Family)
Family <- Family[,colMeans(Family) >=.1]
dim(Family)
library(vegan)
#NMDS
#autotransform=F evita la transformaci?n de datos
#sratmax=0.999999 para que el an?lisis no se detenga antes de tiempo cuando se tienen muchos datos
Family.mds <- metaMDS(Family, k=2, distance="bray", trymax=100, zerodist="add", autotransform=F)
stressplot(Family.mds)
#add the metadata information
##THIS HAS TO BE IN THE SAME ORDER IN THE METADATA FILE AND IN THE COUNTS FILE
colvec<- c("yellowgreen","turquoise4", "tomato2")
plot(Family.mds, type="n", xlim=c(-.5,.5), ylim=c(-0.6,0.6))
pl <-ordiellipse(Family.mds, Metadata$Month, kind="se", conf=0.95, lwd=2, col="gray30", label=T)
#plot with ggplot
data.scores <- as.data.frame(scores(Family.mds))
data.scores$site <- rownames(data.scores)
data.scores$grp <- Metadata$Month
head(data.scores)
library(ggplot2)
NMDS <- data.frame(MDS1 = Family.mds$points[,1], MDS2 = Family.mds$points[,2],group=Metadata$Month)
#NMDS.mean <- aggregate(NMDS[,1:2],list(group=group),mean)
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100)
{
theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))
}
#attach(NMDS.mean)
df_ell <- data.frame()
for(g in levels(NMDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$group==g,],
veganCovEllipse(pl[[g]]$cov,pl[[g]]$center,pl[[g]]$scale)))
,group=g))
}
library(dplyr)
data.scores <- data.scores %>%
rename(group = grp) %>%
mutate(group = as.character(group)) %>%
as_tibble()
df_ell <- df_ell %>%
mutate(group = as.character(group)) %>%
as_tibble()
#manual_colour_vals <- c("July" = "#C496A1", "August" = "turquoise4", "September" = "red")
manual_colour_vals <- c("July" ="#00AFBB", "August" ="#E7B800", "September" = "#FC4E07")
ggplot() +
geom_path(
data=df_ell,
aes(
x=NMDS1,
y=NMDS2,
colour=group),
size=1,
linetype=1,
alpha = 0.6)+
geom_polygon(
data=df_ell,
aes(
x=NMDS1,
y=NMDS2,
fill = group),
alpha = 0.2)+
geom_point(
data=data.scores,
aes(
x=NMDS1,
y=NMDS2,
colour=group,
shape = group),
size=3,
alpha = 0.6) +
theme_light()+
ylab("NMDS2")+
xlab("NMDS1")+
scale_colour_manual(values=manual_colour_vals) +
scale_fill_manual(values=manual_colour_vals)
#ALTO. EN ESTA PARTE CARLOS, ME HACE EL ANALISIS DE ADONIS ENTRE TODOS LOS NIVELES DE CADA CATEGORIA
#A MI ME GUSTARIA EL ANALISIS POR PARES, OSEA JUL VS AGOSTO, JUL VS SEPTIEMBRE Y SEP VS AGOSTO
dis <- vegdist(Family, method="bray")
adonis(dis ~ SamplingPoint*Depth*Month*Filter, Metadata, perm=999)
#ANOSIM #MARIBEL NO LO HIZO DIFERENCIA ATRIBUIBLE AL GRUPO
#provides a way to test statistically whether there is a significant difference between two or more groups of sampling units
anosim(Family, Metadata$SamplingPoint, permutations = 999, distance = "bray")
#the following example its what im trying to do... pairwise adonis
library(vegan)
library(dplyr)
dis = as.matrix((vegan::vegdist(Metadata[, -c(1:2)], "bray")))
site_combs <- combn(unique(Metadata$Month), 2)
df <- data.frame(July = Month_combs[1,], August = Month_combs[2,], adonis.p_value = NA)
for(i in 1:length(rownames(df))){
temp <- Metadata %>%
dplyr::filter(Month == df$July[i] | Site == df$August[i])
df$adonis.p_value[i] <- as.numeric(vegan::adonis(dis ~ Metadata$Month, strata = Metadata$Month, perm = 999)$factors$pvals)
}