Spatio-temporal dynamics, composition, diversity and ecological processes shaping fungal assemblages in forest soils along the transect of Iztaccíhuatl-Popocatépetl (PNIP) and La Malinche (PNLM) National Parks
Stephanie Hereira-Pacheco
Universidad Autónoma de Tlaxcalasehereirap@uatx.mx
08 - 01 - 2025
Preprocessing sequences
IMPORTING TO QIIME2
- Make directories for each library
mkdir ITS1
mkdir ITS2
mkdir ITS3
mkdir ITS4
mkdir ITS5
mkdir ITS6
mkdir ITS7
mkdir ITS8
mkdir ITS9
- Loop to rename forward seqs
for i in ITS* ;do mv ITS*/ITS_*1.fastq.gz ITS*/forward.fastq.gz ; done- Loop to rename reverse seqs
for i in ITS* ;do mv ITS*/ITS_*2.fastq.gz ITS*/reverse.fastq.gz ; done- Loop to import files to QIIME2
for i in ITS*; do qiime tools import --type MultiplexedPairedEndBarcodeInSequence --input-path $i --output-path $i.qza ; done
DEMULTIPLEXING WITH CUTADAPT
qiime cutadapt demux-paired\
--i-seqs its1.qza\
--m-forward-barcodes-file ../../maps_its/its1.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its1_demux
qiime cutadapt demux-paired\
--i-seqs its2.qza\
--m-forward-barcodes-file ../../maps_its/its2.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its2_demux
qiime cutadapt demux-paired\
--i-seqs its3.qza\
--m-forward-barcodes-file ../../maps_its/its3.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its3_demux
qiime cutadapt demux-paired\
--i-seqs its4.qza --m-forward-barcodes-file ../../maps_its/its4.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its4_demux
qiime cutadapt demux-paired --i-seqs its5.qza\
--m-forward-barcodes-file ../../maps_its/its5.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its5_demux
qiime cutadapt demux-paired --i-seqs its6.qza\
--m-forward-barcodes-file ../../maps_its/its6.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its6_demux
qiime cutadapt demux-paired --i-seqs its7.qza\
--m-forward-barcodes-file ../../maps_its/its7.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its7_demux
qiime cutadapt demux-paired\
--i-seqs its8.qza\
--m-forward-barcodes-file ../../maps_its/its8.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its8_demux
qiime cutadapt demux-paired\
--i-seqs its9.qza\
--m-forward-barcodes-file ../../maps_its/its9.txt\
--m-forward-barcodes-column BarcodeSequence\
--output-dir its9_demux
EXPORTING DATA
for i in its*_demux; do qiime tools export --input-path $i/per_sample_sequences.qza --output-path exported_$i; doneITSXPRESS STANDALONE
- Move all data to one directory and run this bash script:
#!/bin/bash
# Set the path to the itsxpress executable
ITSXPRESS_EXECUTABLE="itsxpress"
# Loop through the files
for forward_file in *_R1_001.fastq.gz; do
# Extract the sample name from the forward file
sample_name=$(basename "$forward_file" _R1_001.fastq.gz)
# Construct the reverse file name
reverse_file="${sample_name}_R2_001.fastq.gz"
# Run itsxpress command
$ITSXPRESS_EXECUTABLE --fastq "$forward_file" --fastq2 "$reverse_file" --region ITS2 \
--taxa Fungi --log "${sample_name}_logfile.txt" --outfile "${sample_name}_trimmed_reads.fastq.gz" --threads 2
# Optionally, you can print a message indicating the completion of each iteration
echo "Processing $forward_file and $reverse_file"
doneIMPORTING STANDALONE DATA
qiime tools import\
--type 'SampleData[SequencesWithQuality]'\
--input-path manifest.txt\
--output-path single-end-demux-standalone.qza\
--input-format SingleEndFastqManifestPhred33V2OTU’s FOR MERGED SEQUENCES
For qiime2-itsxpress data, we have to export and run:
for i in *.fastq.gz; do reformat in=$i out=clean_$i minconsecutivebases=1; doneThis script is from (BBMap)[https://github.com/BioInfoTools/BBMap/blob/master/sh/reformat.sh]
Then import again and run all. This is due to an error of ITSxpress that generates sequences with 0 lentgh.
qiime tools import\
--type 'SampleData[SequencesWithQuality]'\
--input-path manifest_reformat.txt\
--output-path single-end-demux-qiime2-reformat.qza\
--input-format SingleEndFastqManifestPhred33V2- Filter by q-score
for i in *;do qiime quality-filter q-score --i-demux $i --output-dir filterqscore_$i; done
- Derreplication
for i in filterqscore_*; do qiime vsearch dereplicate-sequences --i-sequences $i/filtered_sequences.qza --output-dir derep_$i;done- Clustering de novo
for i in derep_* ; do qiime vsearch cluster-features-de-novo --i-sequences $i/dereplicated_sequences.qza --i-table $i/dereplicated_table.qza --p-perc-identity 0.97 --p-threads 4 --output-dir cluster97_$i; done- Chimera checking and filter from table
for i in cluster97_*; do qiime vsearch uchime-denovo --i-sequences $i/clustered_sequences.qza --i-table $i/clustered_table.qza --output-dir chimera97_$i;done
qiime feature-table filter-features\
--i-table cluster97_derep_filterqscore_single-end-demux-qiime2-reformat.qza/clustered_table.qza\
--m-metadata-file chimera97_cluster_derep_filterqscore_single-end-demux-qiime2-reformat.qza/nonchimeras.qza \
--o-filtered-table chimera97_cluster_derep_filterqscore_single-end-demux-qiime2-reformat.qza/table97-nonchimeras-qiime2.qza
- Filtering singletons from table and seqs
for i in chimera* ; do qiime feature-table filter-features --i-table $i/table*.qza --p-min-frequency 2 --o-filtered-table $i/filtered_$i ; donefor i in chimera*/; do
data_file=$(find $i -name 'nonchimeras*' -type f)
table_file=$(find $i -name 'filtered-table*' -type f)
qiime feature-table filter-seqs \
--i-data $data_file \
--i-table $table_file \
--o-filtered-data ${data_file%.qza}-filtered.qza
done
qiime feature-table summarize \
--i-table filtered-table.qza \
--o-visualization visualization.qzv
qiime diversity alpha-rarefaction \
--i-table filtered-table.qza \
--p-max-depth 41834 \
--p-min-depth 1 \
--p-steps 10 \
--p-iterations 10 \
--o-visualization visualization-rar.qzv
qiime feature-table rarefy \
--i-table filtered-table-2.qza \
--p-sampling-depth 2313 \
--p-no-with-replacement \
--o-rarefied-table rarefied-table.qza
Phycochemical characteristics - Both seasons
- Loading packages and data
library(reshape2)
library(picante)
library(readxl)
library(tidyverse)
library(pgirmess)
library(car)
library(ggplot2)
library(ggpubr)
dry<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_secas_all") %>% mutate(Season="Dry")
wet<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_lluvias_all") %>% mutate(Season="Wet")
# select characteristics in both datasets
drys<- dry %>% dplyr::select( Poligono,Sitio,Transecto,pH,MO,N,P,Season)
wets<- wet %>% dplyr::select(Poligono,Sitio,Transecto,pH,MO,N,P,Season)
env_season<- rbind(drys,wets) %>% mutate(sea=c(rep("D",36), rep("R",36))) %>% unite("SampleID", c("Poligono", "Sitio", "Transecto", "sea"),remove = F, sep = "")
env_season$Poligono<- factor(env_season$Poligono, levels = c(1,2,3,4,5,6),
labels = c("P1", "P2", "P3", "P4", "P5", "P6"))
env_season$Season<- factor(env_season$Season)
env_season<-env_season %>% unite("interact", c("Poligono", "Season"), remove = F)
df <-env_season %>% column_to_rownames(var="SampleID") %>% dplyr::select(pH,MO,N,P)
dfs=data.frame(scale(df, center = T, scale = T)) #standardize environmental data
df$Poligono<-NA
df$Poligono<-env_season$Poligono
df$Season<-env_season$Season
dfs$Poligono<-NA
dfs$Poligono<-env_season$Poligono
#melted_df=melt(dfs, id="Poligono")
#melted_df=melt(df, id=c("Poligono", "Season"))
#mean_table=Rmisc::summarySE(melted_df, measurevar = "value", groupvars=c("variable","Poligono", "Season"), na.rm = T)
# checking collinearity
cor.table(na.omit(dfs[,1:4]))## $r
## pH MO N P
## pH 1.00000000 -0.2026648 -0.02146947 -0.1934269
## MO -0.20266479 1.0000000 0.17447231 0.1834227
## N -0.02146947 0.1744723 1.00000000 0.3432520
## P -0.19342691 0.1834227 0.34325203 1.0000000
##
## $df
## [1] 70
##
## $P
## pH MO N P
## pH 0.00000000 0.08775927 0.857932381 0.103530901
## MO 0.08775927 0.00000000 0.142701322 0.123007768
## N 0.85793238 0.14270132 0.000000000 0.003158629
## P 0.10353090 0.12300777 0.003158629 0.000000000##
## Shapiro-Wilk normality test
##
## data: env_season$pH
## W = 0.98364, p-value = 0.4738##
## Shapiro-Wilk normality test
##
## data: env_season$MO
## W = 0.95713, p-value = 0.01526##
## Shapiro-Wilk normality test
##
## data: env_season$N
## W = 0.87959, p-value = 5.103e-06##
## Shapiro-Wilk normality test
##
## data: env_season$P
## W = 0.88619, p-value = 8.936e-06k.ph.s<-round(kruskal.test(env_season$pH~env_season$Season)$p.value, digits = 3)
k.ph.p<-round(kruskal.test(env_season$pH~env_season$Poligono)$p.value, digits = 3)
leveneTest(env_season$pH~env_season$Season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3122 0.2559
## 70## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 1.4212 0.2282
## 66## [1] 0.005## [1] 0.001k.mo.s<-round(kruskal.test(env_season$MO~env_season$Season)$p.value, 3)
k.mo.p<-kruskal.test(env_season$MO~env_season$Poligono)$p.value
leveneTest(env_season$MO~env_season$Season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0155 0.9013
## 70## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.9802 0.4366
## 66## [1] 0.093## [1] 2.182728e-05k.n.s<-round(kruskal.test(env_season$N~env_season$Season)$p.value,3)
k.n.p<-round(kruskal.test(env_season$N~env_season$Poligono)$p.value,3)
leveneTest(env_season$N~env_season$Season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 4.6553 0.03439 *
## 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 1.4895 0.2052
## 66## [1] 0.001## [1] 0.006k.p.s<-round(kruskal.test(env_season$P~env_season$Season)$p.value,3)
k.p.p<-kruskal.test(env_season$P~env_season$Poligono)$p.value
leveneTest(env_season$P~env_season$Season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.1499 0.08028 .
## 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 3.0447 0.0156 *
## 66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] 0.222## [1] 7.84677e-08all.kw<- data.frame(pH=c(k.ph.s,k.ph.p),
OM=c(k.mo.s, k.mo.p),
N=c(k.n.s, k.n.p),
P=c(k.p.s, k.p.p))
rownames(all.kw)<-c("Season", "Poligon")
#write.csv(all.kw,"oldies/Data/all.kw.csv")#linear models and normality and heterogenity test
model1<- lm(pH ~ Poligono*Season, data = env_season)
car::leveneTest(pH ~ Poligono*Season, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 0.8665 0.5767
## 60
model2<- lm(MO ~ Poligono*Season, data = env_season)
car::leveneTest(MO ~ Poligono*Season, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 0.9248 0.5232
## 60
model3<- lm(N ~ Poligono*Season, data = env_season)
car::leveneTest(N ~ Poligono*Season, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.4854 0.161
## 60
model4<- lm(P ~ Poligono*Season, data = env_season)
car::leveneTest(P ~ Poligono*Season, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 0.7239 0.7113
## 60
#linear models for test with interaction
library(tidyverse)
lm.ph<-lm(pH~Season*Poligono ,data = env_season)
lm.ph.perm<-summary(aov(lm.ph))
lm.ph.perm<- lm.ph.perm[[1]][5] %>%as.data.frame() %>% slice(1:3) %>% mutate_if(is.numeric, ~round(.,digits = 3))
lm.ph.perm## Pr(>F)
## Season 0.002
## Poligono 0.000
## Season:Poligono 0.007lm.mo<-lm(MO~Season*Poligono ,data = env_season)
lm.mo.perm<-summary(aov(lm.mo))
lm.mo.perm<-lm.mo.perm[[1]][5] %>%as.data.frame() %>% slice(1:3) %>% mutate_if(is.numeric, ~round(.,digits = 3))
lm.mo.perm## Pr(>F)
## Season 0.010
## Poligono 0.000
## Season:Poligono 0.224lm.n<-glm(N~Season*Poligono ,data = env_season, family=Gamma(link = "identity"))
lm.n.perm<-summary(aov(lm.n))
lm.n.perm<- lm.n.perm[[1]][5] %>%as.data.frame() %>% slice(1:3) %>% mutate_if(is.numeric, ~round(.,digits = 3))
lm.n.perm## Pr(>F)
## Season 0.000
## Poligono 0.001
## Season:Poligono 0.016lm.p<-lm(P~Season*Poligono ,data = env_season,family=Gamma(link = "identity"))
lm.p.perm<-summary(aov(lm.p))
lm.p.perm<- lm.p.perm[[1]][5] %>%as.data.frame() %>% slice(1:3) %>% mutate_if(is.numeric, ~round(.,digits = 3))
lm.p.perm## Pr(>F)
## Season 0.034
## Poligono 0.000
## Season:Poligono 0.010all.lm<- cbind(lm.ph.perm, lm.mo.perm, lm.n.perm, lm.p.perm)
colnames(all.lm)<-c("pH", "MO", "N", "P")
write.csv(all.lm,"../Data/all.lm.csv")# Plot the temporal trends of the selected environmental variables over time
df$Season<-NA
df$Season<-env_season$Season
df$Season = factor(df$Season, levels = c("Dry", "Wet"), labels = c("Dry", "Rainy"))
melted_df=melt(df, id=c("Poligono", "Season"))
mean=Rmisc::summarySE(melted_df, measurevar = "value", groupvars=c("Poligono","variable", "Season"), na.rm = T)
mean$variable=factor(mean$variable, levels = c("pH", "MO", "N", "P"))
conservation_status <- c(
pH = "pH",
MO = "OM (%)",
N = "N (mg/Kg)",
P= "P (mg/Kg)")
library(plyr)
mean$conservation2 <- plyr::revalue(mean$variable, conservation_status)
melted_df$conservation2 <- plyr::revalue(melted_df$variable, conservation_status)
labels <- all.lm %>% t() %>%
as.data.frame() %>%
mutate(letter= c("A", "B", "C", "D")) %>%
rownames_to_column(var = "variable") %>%
mutate(`Poligono `=case_when(
`Poligono `==0 ~ "0.000",
TRUE ~ as.character(`Poligono `))) %>%
mutate(`Season `=case_when(
`Season `==0 ~ "0.000" ,
TRUE ~ as.character(`Season `))) %>%
mutate(labels = paste0(letter, " Season = ",`Season `, ", Polygon = ", `Poligono `, ", Season:Polygon = ", signif(`Season:Poligono`, digits = 3)))
#labels <- labels %>% rownames_to_column(var = "conservation2")
labels$conservation2<-factor(labels$variable, levels = c("pH", "MO", "N", "P"), labels =c("pH", "OM (%)", "N (mg/Kg)", "P (mg/Kg)") )
labels<- labels %>% mutate(Poligono="P6") %>% mutate(value=c(7.5,25,130,80))
fq<-ggline(melted_df, x="Poligono", y="value", shape = "Season",
linetype="Season", facet.by = "conservation2", add = "mean_sd")+
facet_grid(conservation2~., scales = "free_y")+labs(x="Polygon",y="Mean±SD")+
theme_grey()+
theme(legend.position="right")+
theme(axis.text = element_text(size=12),
axis.title = element_text(size=14),
strip.text = element_text(size=12),
plot.title = element_text(size=8, face="bold") )+
geom_text(data = labels, aes(label=labels),
fontface="bold",
hjust = 1, vjust = 1)+ theme(
strip.background = element_rect(
color="black", fill="black", size=1.5, linetype="solid"),
strip.text = element_text(colour = "white", face = "bold"))
fq
Phycochemical characteristics - Only dry season
- Loading packages and data
dry<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_secas_all") %>% mutate(Season="Dry")
# select characteristics in both datasets
drys<- dry %>% dplyr::select( Poligono,Sitio,Transecto,Season,K, Ca, Mg, Fe, Cu, Mn, WHC, EC=CONDUC, Clay=ARCILLA, Silt=LIMO)
env_season<- drys %>% mutate(sea=c(rep("D",36))) %>% unite("SampleID", c("Poligono", "Sitio", "Transecto", "sea"),remove = F, sep = "")
env_season$Poligono<- factor(env_season$Poligono, levels = c(1,2,3,4,5,6),
labels = c("P1", "P2", "P3", "P4", "P5", "P6"))
env_season$Season<- factor(env_season$Season)
env_season<-env_season %>% unite("interact", c("Poligono", "Season"), remove = F)
vars= c("K", "Ca", "Mg","Mn", "EC", "Fe", "Cu", "WHC", "Clay", "Silt")
df <-env_season %>% column_to_rownames(var="SampleID") %>% dplyr::select(vars)
dfs=data.frame(scale(df, center = T, scale = T)) #standardize environmental data
df$Poligono<-NA
df$Poligono<-env_season$Poligono
df$Season<-env_season$Season
dfs$Poligono<-NA
dfs$Poligono<-env_season$Poligono
#melted_df=melt(dfs, id="Poligono")
#melted_df=melt(df, id=c("Poligono", "Season"))
#mean_table=Rmisc::summarySE(melted_df, measurevar = "value", groupvars=c("variable","Poligono", "Season"), na.rm = T)
# checking collinearity
cor.table(na.omit(dfs[,1:4]))## $r
## K Ca Mg Mn
## K 1.0000000 0.4698482 0.6787645 0.1380205
## Ca 0.4698482 1.0000000 0.6240184 0.3059809
## Mg 0.6787645 0.6240184 1.0000000 0.6417935
## Mn 0.1380205 0.3059809 0.6417935 1.0000000
##
## $df
## [1] 34
##
## $P
## K Ca Mg Mn
## K 0.000000e+00 3.836723e-03 5.357233e-06 4.221192e-01
## Ca 3.836723e-03 0.000000e+00 4.768783e-05 6.953911e-02
## Mg 5.357233e-06 4.768783e-05 0.000000e+00 2.457224e-05
## Mn 4.221192e-01 6.953911e-02 2.457224e-05 0.000000e+00#linear models and normality and heterogenity test
model1<- lm(K ~ Poligono, data = env_season)
car::leveneTest(K ~ Poligono, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.9335 0.4734
## 30
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 1.7016 0.1647
## 30
model3<- lm(Mg ~ Poligono, data = env_season)
car::leveneTest(Mg ~ Poligono*Season, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.2841 0.9181
## 30
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.1672 0.9727
## 30
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.6226 0.04407 *
## 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.5517 0.7358
## 30
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.7102 0.6205
## 30
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.2602 0.07388 .
## 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.9074 0.02951 *
## 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model10<- lm(Silt ~ Poligono, data = env_season)
car::leveneTest(Silt ~ Poligono, data = env_season)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.8875 0.5017
## 30
# Plot the temporal trends of the selected environmental variables over time
melted_df=melt(df, id=c("Poligono", "Season"))
mean=Rmisc::summarySE(melted_df, measurevar = "value", groupvars=c("Poligono","variable", "Season"), na.rm = T)
mean$variable=factor(mean$variable, levels = c("K", "Ca", "Mg","Mn", "EC", "Fe", "Cu", "WHC", "Clay", "Silt"))
conservation_status <- c(
K = "K (mg/Kg)",
Ca = "Ca (mg/Kg)",
Mg = "Mg (mg/Kg)",
Mn = "Mn (mg/Kg)",
EC = "EC (dS/m)",
Fe= "Fe (mg/Kg)",
Cu = "Cu (mg/Kg)",
WHC= "WHC (g/Kg) ",
Clay = "Clay (%)",
Silt= "Silt (%)"
)
library(plyr)
mean$conservation2 <- plyr::revalue(mean$variable, conservation_status)
melted_df$conservation2 <- plyr::revalue(melted_df$variable, conservation_status)
fq<-ggline(melted_df, x="Poligono", y="value",
linetype="Season", facet.by = "conservation2", add = "mean_se")+
facet_wrap(conservation2~., scales = "free_y", ncol = 5)+
labs(x="Polygon",y="Mean±SD")+
theme_grey()+
theme(legend.position="right")+
theme(axis.text = element_text(size=12),
axis.title = element_text(size=14),
strip.text = element_text(size=12),
plot.title = element_text(size=8, face="bold") )+ theme(
strip.background = element_rect(
color="black", fill="black", size=1.5, linetype="solid"),
strip.text = element_text(colour = "white", face = "bold"),
legend.position = "none")+
stat_compare_means(method = "anova", label = "p.format", label.x = 2.5)
fq
Alpha diversity
- Load files and packages
library(qiime2R)
library(tidyverse)
library(hillR)
library(ggpubr)
table_otus<- read_qza("../Data/table_SMOQ_filt.qza")$data
table_asvs<- read_qza("../Data/table_SMDR.qza")$data
meta<- read.delim("../Data/its_map.txt")- Calculate hill numbers and join
hill0<- hill_taxa(comm = table_otus, q = 0, MARGIN = 2)
hill1<- hill_taxa(comm = table_otus, q = 1, MARGIN = 2)
hill2<- hill_taxa(comm = table_otus, q = 2, MARGIN = 2)
hill0a<- hill_taxa(comm = table_asvs, q = 0, MARGIN = 2)
hill1a<- hill_taxa(comm = table_asvs, q = 1, MARGIN = 2)
hill2a<- hill_taxa(comm = table_asvs, q = 2, MARGIN = 2)
q0<- as.data.frame(hill0) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill0") %>%
mutate(order="q0") %>% rownames_to_column(var = "SampleID")
q1<- as.data.frame(hill1) %>% #colMeans() %>%
round() %>% as.data.frame()%>% dplyr::rename(q="hill1") %>%
mutate(order="q1")%>% rownames_to_column(var = "SampleID")
q2<- as.data.frame(hill2) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill2") %>%
mutate(order="q2")%>% rownames_to_column(var = "SampleID")
q0a<- as.data.frame(hill0a) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill0a") %>%
mutate(order="q0") %>% rownames_to_column(var = "SampleID")
q1a<- as.data.frame(hill1a) %>% #colMeans() %>%
round() %>% as.data.frame()%>% dplyr::rename(q="hill1a") %>%
mutate(order="q1")%>% rownames_to_column(var = "SampleID")
q2a<- as.data.frame(hill2a) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill2a") %>%
mutate(order="q2")%>% rownames_to_column(var = "SampleID")
qs<- q0 %>% full_join(q1) %>% full_join(q2) %>% pivot_longer(
cols = c(-SampleID,-order), values_to = "val", names_to = "ids") %>% inner_join(meta) %>% unite("interact",c("Poligono", "Season"),remove = F)
qsa<- q0a %>% full_join(q1a) %>% full_join(q2a) %>% pivot_longer(
cols = c(-SampleID,-order), values_to = "val", names_to = "ids") %>% inner_join(meta)%>% unite("interact",c("Poligono", "Season"),remove = F)
#qs<-qsa
#write.csv(qs, "../Data/qs.csv", row.names = F)#install.packages("ggpubr")
library(ggpubr)
depth<- read.delim("../Data/depth.txt")
q0_data<- q0 %>% dplyr::rename("q0"="q") %>% dplyr::select(-order)
q1_data<- q1 %>% dplyr::rename("q1"="q") %>% dplyr::select(-order)
q2_data<- q2 %>% dplyr::rename("q2"="q") %>% dplyr::select(-order)
q_data<- q0_data %>% inner_join(q1_data) %>% inner_join(q2_data) %>%
inner_join(depth, by = c("SampleID"="Sample.ID"))
q0_vs_depth_fltr <- ggscatter(q_data ,
x = "Frequency",
y = "q0",
xlab= "Sequencing depth (number of reads)",
ylab = "q=0",
#ylab="Alpha diversity q=0 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#B03A2E", fill = "#566573"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+theme_grey()+
theme(legend.title = element_blank(), legend.position = "none")+ labs(y=expression(paste(italic("q"), "=0", " (number of total OTUs/OTUs)")))
q1_vs_depth_fltr <- ggscatter(q_data,
x = "Frequency",
y = "q1",
xlab= "Sequencing depth (number of reads)",
ylab = "q=1",
#ylab="Alpha diversity q=1 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#B03A2E", fill = "#566573"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+theme_grey()+
theme(legend.title = element_blank(), legend.position = "none")+ labs(y=expression(paste(italic("q"), "=1", " (number of frequent OTUs/OTUs)")))
q2_vs_depth_fltr <- ggscatter(q_data,
x = "Frequency",
y = "q2",
xlab= "Sequencing depth (number of reads)",
ylab = "q=2",
#ylab="Alpha diversity q=2 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#B03A2E", fill = "#566573"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+theme_grey()+
theme(legend.title = element_blank(), legend.position = "none")+ labs(y=expression(paste(italic("q"), "=2", " (number of dominant OTUs/OTUs)")))
library(ggplot2)
library(cowplot)
#save plots
#Combine plot
title_corr_plot <- ggdraw() + draw_label("Alpha diversity depth correlation with samples")
correlation_plot_q012_fltr <- plot_grid(q0_vs_depth_fltr,q1_vs_depth_fltr,q2_vs_depth_fltr, labels = c("A", "B", "C"), nrow = 1, rel_heights = c(1, 1, 1))
correlation_plot_q012_fltr
ggsave('../Plots/alpha_corr.png', width = 11, height = 6, dpi = 300, plot =correlation_plot_q012_fltr)lm.q0<-qs %>% filter(order=="q0")
model0<- lm(val ~ Poligono*Season, data = lm.q0)
car::leveneTest(val ~ Poligono*Season, data = lm.q0)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 0.8236 0.6171
## 60
#shapiro.test(model0$residuals)
a0<- aov(model0)
a00<-summary(a0)
pvalues0<-a00[[1]][[5]]
modelinteract0<- aov(lm(val ~ interact, data = lm.q0))
agricolae::HSD.test(modelinteract0, "interact", console = TRUE)##
## Study: modelinteract0 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 38637.89
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 556.6667 166.31015 6 80.24742 304 814 503.75 559.0 603.75
## P1_Rainy 671.1667 181.69471 6 80.24742 403 929 585.25 669.0 767.00
## P2_Dry 720.5000 168.89020 6 80.24742 499 979 639.00 685.0 809.75
## P2_Rainy 744.6667 199.63834 6 80.24742 546 1037 576.00 714.0 873.75
## P3_Dry 969.3333 289.83283 6 80.24742 704 1527 820.50 912.5 955.75
## P3_Rainy 619.0000 130.62925 6 80.24742 475 822 519.00 613.5 680.25
## P4_Dry 1351.5000 291.65099 6 80.24742 1061 1859 1166.25 1277.5 1448.00
## P4_Rainy 827.5000 137.61214 6 80.24742 672 1039 720.50 825.5 894.50
## P5_Dry 705.6667 212.39837 6 80.24742 447 930 521.75 746.0 875.00
## P5_Rainy 731.8333 196.62494 6 80.24742 510 973 581.00 718.5 881.50
## P6_Dry 565.8333 96.61763 6 80.24742 477 751 513.50 548.5 563.25
## P6_Rainy 719.1667 190.86479 6 80.24742 461 1036 654.25 684.5 774.75
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 385.8678
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 1351.5000 a
## P3_Dry 969.3333 ab
## P4_Rainy 827.5000 bc
## P2_Rainy 744.6667 bc
## P5_Rainy 731.8333 bc
## P2_Dry 720.5000 bc
## P6_Rainy 719.1667 bc
## P5_Dry 705.6667 bc
## P1_Rainy 671.1667 bc
## P3_Rainy 619.0000 bc
## P6_Dry 565.8333 c
## P1_Dry 556.6667 c##
## Study: modelinteract0 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 38637.89
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 556.6667 166.31015 6 80.24742 304 814 503.75 559.0 603.75
## P1_Rainy 671.1667 181.69471 6 80.24742 403 929 585.25 669.0 767.00
## P2_Dry 720.5000 168.89020 6 80.24742 499 979 639.00 685.0 809.75
## P2_Rainy 744.6667 199.63834 6 80.24742 546 1037 576.00 714.0 873.75
## P3_Dry 969.3333 289.83283 6 80.24742 704 1527 820.50 912.5 955.75
## P3_Rainy 619.0000 130.62925 6 80.24742 475 822 519.00 613.5 680.25
## P4_Dry 1351.5000 291.65099 6 80.24742 1061 1859 1166.25 1277.5 1448.00
## P4_Rainy 827.5000 137.61214 6 80.24742 672 1039 720.50 825.5 894.50
## P5_Dry 705.6667 212.39837 6 80.24742 447 930 521.75 746.0 875.00
## P5_Rainy 731.8333 196.62494 6 80.24742 510 973 581.00 718.5 881.50
## P6_Dry 565.8333 96.61763 6 80.24742 477 751 513.50 548.5 563.25
## P6_Rainy 719.1667 190.86479 6 80.24742 461 1036 654.25 684.5 774.75
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 385.8678
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 1351.5000 a
## P3_Dry 969.3333 ab
## P4_Rainy 827.5000 bc
## P2_Rainy 744.6667 bc
## P5_Rainy 731.8333 bc
## P2_Dry 720.5000 bc
## P6_Rainy 719.1667 bc
## P5_Dry 705.6667 bc
## P1_Rainy 671.1667 bc
## P3_Rainy 619.0000 bc
## P6_Dry 565.8333 c
## P1_Dry 556.6667 clm.q1<-qs %>% filter(order=="q1")
model1<- lm(val ~ Poligono*Season, data = lm.q1)
car::leveneTest(val ~ Poligono*Season, data = lm.q1)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.5605 0.1344
## 60
#shapiro.test(model1$residuals)
a1<- aov(model1)
a11<-summary(a1)
pvalues1<-a11[[1]][[5]]
modelinteract1<- aov(lm(val ~ interact, data = lm.q1))
agricolae::HSD.test(modelinteract1, "interact", console = TRUE)##
## Study: modelinteract1 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 1137.844
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 88.33333 13.76469 6 13.77101 70 102 77.50 91.0 100.00
## P1_Rainy 57.33333 37.97719 6 13.77101 22 115 28.75 46.0 80.50
## P2_Dry 100.16667 32.99949 6 13.77101 63 146 72.50 101.0 120.50
## P2_Rainy 119.33333 39.60135 6 13.77101 72 172 90.50 116.0 147.50
## P3_Dry 137.83333 46.48190 6 13.77101 75 202 106.25 143.5 162.75
## P3_Rainy 95.00000 42.70831 6 13.77101 40 146 61.25 102.5 124.25
## P4_Dry 202.16667 53.30447 6 13.77101 125 278 174.25 206.0 227.25
## P4_Rainy 120.00000 15.41428 6 13.77101 102 143 108.75 118.5 129.00
## P5_Dry 117.33333 22.29499 6 13.77101 78 144 114.00 119.0 128.50
## P5_Rainy 100.50000 26.62142 6 13.77101 69 137 80.00 99.0 118.75
## P6_Dry 108.66667 26.40959 6 13.77101 62 135 100.50 116.0 124.75
## P6_Rainy 101.00000 19.95996 6 13.77101 75 129 88.00 100.0 113.50
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 66.21759
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 202.16667 a
## P3_Dry 137.83333 ab
## P4_Rainy 120.00000 bc
## P2_Rainy 119.33333 bc
## P5_Dry 117.33333 bc
## P6_Dry 108.66667 bc
## P6_Rainy 101.00000 bc
## P5_Rainy 100.50000 bc
## P2_Dry 100.16667 bc
## P3_Rainy 95.00000 bc
## P1_Dry 88.33333 bc
## P1_Rainy 57.33333 c##
## Study: modelinteract1 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 1137.844
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 88.33333 13.76469 6 13.77101 70 102 77.50 91.0 100.00
## P1_Rainy 57.33333 37.97719 6 13.77101 22 115 28.75 46.0 80.50
## P2_Dry 100.16667 32.99949 6 13.77101 63 146 72.50 101.0 120.50
## P2_Rainy 119.33333 39.60135 6 13.77101 72 172 90.50 116.0 147.50
## P3_Dry 137.83333 46.48190 6 13.77101 75 202 106.25 143.5 162.75
## P3_Rainy 95.00000 42.70831 6 13.77101 40 146 61.25 102.5 124.25
## P4_Dry 202.16667 53.30447 6 13.77101 125 278 174.25 206.0 227.25
## P4_Rainy 120.00000 15.41428 6 13.77101 102 143 108.75 118.5 129.00
## P5_Dry 117.33333 22.29499 6 13.77101 78 144 114.00 119.0 128.50
## P5_Rainy 100.50000 26.62142 6 13.77101 69 137 80.00 99.0 118.75
## P6_Dry 108.66667 26.40959 6 13.77101 62 135 100.50 116.0 124.75
## P6_Rainy 101.00000 19.95996 6 13.77101 75 129 88.00 100.0 113.50
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 66.21759
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 202.16667 a
## P3_Dry 137.83333 ab
## P4_Rainy 120.00000 bc
## P2_Rainy 119.33333 bc
## P5_Dry 117.33333 bc
## P6_Dry 108.66667 bc
## P6_Rainy 101.00000 bc
## P5_Rainy 100.50000 bc
## P2_Dry 100.16667 bc
## P3_Rainy 95.00000 bc
## P1_Dry 88.33333 bc
## P1_Rainy 57.33333 clm.q2<-qs %>% filter(order=="q2")
model2<- lm(val ~ Poligono*Season, data = lm.q2)
car::leveneTest(val ~ Poligono*Season, data = lm.q2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.3283 0.2317
## 60
#shapiro.test(model2$residuals)
a2<- aov(model2)
a22<-summary(a2)
pvalues2<-a22[[1]][[5]]
modelinteract2<- aov(lm(val ~ interact, data = lm.q2))
let2<-agricolae::HSD.test(modelinteract2, "interact", console = TRUE)$groups%>%
mutate(order="q2")##
## Study: modelinteract2 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 273.2417
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 32.50000 7.422937 6 6.748354 26 46 27.00 31.5 33.75
## P1_Rainy 18.16667 14.048725 6 6.748354 5 41 8.50 13.0 25.75
## P2_Dry 29.00000 11.696153 6 6.748354 14 42 20.25 30.0 38.25
## P2_Rainy 42.00000 19.636700 6 6.748354 17 74 30.75 41.5 48.50
## P3_Dry 44.16667 20.517472 6 6.748354 22 69 26.00 45.0 59.50
## P3_Rainy 31.66667 18.532854 6 6.748354 10 57 16.25 35.5 40.50
## P4_Dry 60.83333 27.629091 6 6.748354 19 87 43.50 67.5 83.25
## P4_Rainy 41.16667 11.771434 6 6.748354 26 60 36.50 38.0 46.25
## P5_Dry 40.16667 17.814788 6 6.748354 20 65 26.25 38.5 52.25
## P5_Rainy 32.83333 13.585532 6 6.748354 19 56 23.75 29.5 38.25
## P6_Dry 38.66667 14.264174 6 6.748354 20 59 30.75 36.5 47.50
## P6_Rainy 34.66667 11.500725 6 6.748354 19 48 26.25 39.0 40.50
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 32.4493
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 60.83333 a
## P3_Dry 44.16667 ab
## P2_Rainy 42.00000 ab
## P4_Rainy 41.16667 ab
## P5_Dry 40.16667 ab
## P6_Dry 38.66667 ab
## P6_Rainy 34.66667 ab
## P5_Rainy 32.83333 ab
## P1_Dry 32.50000 ab
## P3_Rainy 31.66667 ab
## P2_Dry 29.00000 ab
## P1_Rainy 18.16667 b#letters_int<-rbind(let0, let1, let2)#%>% rownames_to_column(var="interact") %>%
#dplyr::rename(mean=val)library(ggh4x)
library(rcartocolor)
#colors<-viridis::turbo(6, alpha = 1, begin = 0, end = 1, direction = 1)
colors = carto_pal(6, "Safe")
my_colors=carto_pal(2, "Earth")
my_colors## [1] "#A16928" "#2887a1"strip <- strip_themed(background_x = elem_list_rect(fill = colors))
labels <- let0 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=1600)
a <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q0") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues0[1], digits = 3), "; ",
"Season =",signif(pvalues0[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues0[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))
labels <- let1 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=250)
b <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q1") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues1[1], digits = 3), "; ",
"Season =",signif(pvalues1[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues1[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))+
theme(strip.text.x = element_blank())
#strip.background = element_rect(color="black", size=1.5, linetype="solid" ))
labels <- let2 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=80)
c <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q2") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_text(size = 14, colour = "black"),
# axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues2[1], digits = 3), "; ",
"Season =",signif(pvalues2[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues2[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))+
theme(strip.text.x = element_blank())library(cowplot)
alpha<- plot_grid(a + theme(legend.position = "none", axis.title.x = element_blank())+ylab(""),
b+theme(legend.position = "none", axis.title.x = element_blank(),
axis.title.y = element_text(size = 15)),
c+ theme(legend.position = "none")+ylab(""),
labels = c("A", "B", "C"),
ncol = 1)
alpha
Alpha diversity - rarified
- Load files and packages
library(qiime2R)
library(tidyverse)
library(hillR)
library(ggpubr)
table_otus<- read_qza("../Data/table_SMOQ_rar.qza")$data
table_asvs<- read_qza("../Data/table_SMDR.qza")$data
meta<- read.delim("../Data/its_map.txt")- Calculate hill numbers and join
hill0<- hill_taxa(comm = table_otus, q = 0, MARGIN = 2)
hill1<- hill_taxa(comm = table_otus, q = 1, MARGIN = 2)
hill2<- hill_taxa(comm = table_otus, q = 2, MARGIN = 2)
hill0a<- hill_taxa(comm = table_asvs, q = 0, MARGIN = 2)
hill1a<- hill_taxa(comm = table_asvs, q = 1, MARGIN = 2)
hill2a<- hill_taxa(comm = table_asvs, q = 2, MARGIN = 2)
q0<- as.data.frame(hill0) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill0") %>%
mutate(order="q0") %>% rownames_to_column(var = "SampleID")
q1<- as.data.frame(hill1) %>% #colMeans() %>%
round() %>% as.data.frame()%>% dplyr::rename(q="hill1") %>%
mutate(order="q1")%>% rownames_to_column(var = "SampleID")
q2<- as.data.frame(hill2) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill2") %>%
mutate(order="q2")%>% rownames_to_column(var = "SampleID")
q0a<- as.data.frame(hill0a) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill0a") %>%
mutate(order="q0") %>% rownames_to_column(var = "SampleID")
q1a<- as.data.frame(hill1a) %>% #colMeans() %>%
round() %>% as.data.frame()%>% dplyr::rename(q="hill1a") %>%
mutate(order="q1")%>% rownames_to_column(var = "SampleID")
q2a<- as.data.frame(hill2a) %>% #colMeans() %>%
round() %>% as.data.frame() %>% dplyr::rename(q="hill2a") %>%
mutate(order="q2")%>% rownames_to_column(var = "SampleID")
qs<- q0 %>% full_join(q1) %>% full_join(q2) %>% pivot_longer(
cols = c(-SampleID,-order), values_to = "val", names_to = "ids") %>% inner_join(meta) %>% unite("interact",c("Poligono", "Season"),remove = F)
qsa<- q0a %>% full_join(q1a) %>% full_join(q2a) %>% pivot_longer(
cols = c(-SampleID,-order), values_to = "val", names_to = "ids") %>% inner_join(meta)%>% unite("interact",c("Poligono", "Season"),remove = F)
#qs<-qsa
#write.csv(qs, "../Data/qs.csv", row.names = F)lm.q0<-qs %>% filter(order=="q0")
model0<- lm(val ~ Poligono*Season, data = lm.q0)
car::leveneTest(val ~ Poligono*Season, data = lm.q0)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.3825 0.2048
## 60
#shapiro.test(model0$residuals)
a0<- aov(model0)
a00<-summary(a0)
pvalues0<-a00[[1]][[5]]
modelinteract0<- aov(lm(val ~ interact, data = lm.q0))
agricolae::HSD.test(modelinteract0, "interact", console = TRUE)##
## Study: modelinteract0 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 2120.25
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 285.0000 27.89982 6 18.79827 260 333 269.25 272.0 296.50
## P1_Rainy 240.0000 62.04192 6 18.79827 170 326 191.25 233.0 283.75
## P2_Dry 322.0000 50.49752 6 18.79827 242 377 306.25 316.5 362.00
## P2_Rainy 328.0000 50.76219 6 18.79827 281 403 286.00 314.5 362.50
## P3_Dry 382.8333 77.04652 6 18.79827 280 514 350.75 382.5 394.00
## P3_Rainy 294.5000 60.46735 6 18.79827 211 358 252.75 298.5 347.25
## P4_Dry 452.8333 30.38695 6 18.79827 411 502 438.50 454.0 460.50
## P4_Rainy 339.3333 25.40604 6 18.79827 290 363 341.00 346.0 350.25
## P5_Dry 333.1667 38.38446 6 18.79827 269 364 314.25 350.0 360.25
## P5_Rainy 310.5000 31.72854 6 18.79827 279 366 294.25 298.0 321.25
## P6_Dry 305.1667 39.13779 6 18.79827 257 371 284.25 301.0 317.00
## P6_Rainy 307.3333 22.70389 6 18.79827 277 341 296.00 303.5 320.00
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 90.39105
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 452.8333 a
## P3_Dry 382.8333 ab
## P4_Rainy 339.3333 bc
## P5_Dry 333.1667 bc
## P2_Rainy 328.0000 bcd
## P2_Dry 322.0000 bcd
## P5_Rainy 310.5000 bcd
## P6_Rainy 307.3333 bcd
## P6_Dry 305.1667 bcd
## P3_Rainy 294.5000 bcd
## P1_Dry 285.0000 cd
## P1_Rainy 240.0000 d##
## Study: modelinteract0 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 2120.25
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 285.0000 27.89982 6 18.79827 260 333 269.25 272.0 296.50
## P1_Rainy 240.0000 62.04192 6 18.79827 170 326 191.25 233.0 283.75
## P2_Dry 322.0000 50.49752 6 18.79827 242 377 306.25 316.5 362.00
## P2_Rainy 328.0000 50.76219 6 18.79827 281 403 286.00 314.5 362.50
## P3_Dry 382.8333 77.04652 6 18.79827 280 514 350.75 382.5 394.00
## P3_Rainy 294.5000 60.46735 6 18.79827 211 358 252.75 298.5 347.25
## P4_Dry 452.8333 30.38695 6 18.79827 411 502 438.50 454.0 460.50
## P4_Rainy 339.3333 25.40604 6 18.79827 290 363 341.00 346.0 350.25
## P5_Dry 333.1667 38.38446 6 18.79827 269 364 314.25 350.0 360.25
## P5_Rainy 310.5000 31.72854 6 18.79827 279 366 294.25 298.0 321.25
## P6_Dry 305.1667 39.13779 6 18.79827 257 371 284.25 301.0 317.00
## P6_Rainy 307.3333 22.70389 6 18.79827 277 341 296.00 303.5 320.00
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 90.39105
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 452.8333 a
## P3_Dry 382.8333 ab
## P4_Rainy 339.3333 bc
## P5_Dry 333.1667 bc
## P2_Rainy 328.0000 bcd
## P2_Dry 322.0000 bcd
## P5_Rainy 310.5000 bcd
## P6_Rainy 307.3333 bcd
## P6_Dry 305.1667 bcd
## P3_Rainy 294.5000 bcd
## P1_Dry 285.0000 cd
## P1_Rainy 240.0000 dlm.q1<-qs %>% filter(order=="q1")
model1<- lm(val ~ Poligono*Season, data = lm.q1)
car::leveneTest(val ~ Poligono*Season, data = lm.q1)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.3542 0.2185
## 60
#shapiro.test(model1$residuals)
a1<- aov(model1)
a11<-summary(a1)
pvalues1<-a11[[1]][[5]]
modelinteract1<- aov(lm(val ~ interact, data = lm.q1))
agricolae::HSD.test(modelinteract1, "interact", console = TRUE)##
## Study: modelinteract1 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 661.55
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 75.50000 11.04083 6 10.5004 60 86 67.25 78.5 84.50
## P1_Rainy 46.33333 29.64906 6 10.5004 17 92 24.00 39.5 63.25
## P2_Dry 78.33333 23.21781 6 10.5004 52 111 58.25 81.5 90.50
## P2_Rainy 98.66667 33.18232 6 10.5004 56 145 77.25 95.5 120.50
## P3_Dry 108.83333 36.42755 6 10.5004 63 159 83.75 108.5 131.00
## P3_Rainy 76.66667 33.04341 6 10.5004 33 116 51.00 86.5 95.75
## P4_Dry 144.50000 34.83533 6 10.5004 92 192 126.50 149.5 161.25
## P4_Rainy 95.66667 13.95230 6 10.5004 81 115 85.00 93.0 105.50
## P5_Dry 97.83333 19.02016 6 10.5004 66 120 91.00 101.0 108.75
## P5_Rainy 81.66667 20.53939 6 10.5004 54 109 68.50 81.0 95.75
## P6_Dry 86.16667 19.82339 6 10.5004 56 110 75.00 89.0 99.25
## P6_Rainy 81.50000 17.14351 6 10.5004 62 107 68.25 81.5 90.25
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 50.49091
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 144.50000 a
## P3_Dry 108.83333 ab
## P2_Rainy 98.66667 ab
## P5_Dry 97.83333 ab
## P4_Rainy 95.66667 abc
## P6_Dry 86.16667 bc
## P5_Rainy 81.66667 bc
## P6_Rainy 81.50000 bc
## P2_Dry 78.33333 bc
## P3_Rainy 76.66667 bc
## P1_Dry 75.50000 bc
## P1_Rainy 46.33333 c##
## Study: modelinteract1 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 661.55
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 75.50000 11.04083 6 10.5004 60 86 67.25 78.5 84.50
## P1_Rainy 46.33333 29.64906 6 10.5004 17 92 24.00 39.5 63.25
## P2_Dry 78.33333 23.21781 6 10.5004 52 111 58.25 81.5 90.50
## P2_Rainy 98.66667 33.18232 6 10.5004 56 145 77.25 95.5 120.50
## P3_Dry 108.83333 36.42755 6 10.5004 63 159 83.75 108.5 131.00
## P3_Rainy 76.66667 33.04341 6 10.5004 33 116 51.00 86.5 95.75
## P4_Dry 144.50000 34.83533 6 10.5004 92 192 126.50 149.5 161.25
## P4_Rainy 95.66667 13.95230 6 10.5004 81 115 85.00 93.0 105.50
## P5_Dry 97.83333 19.02016 6 10.5004 66 120 91.00 101.0 108.75
## P5_Rainy 81.66667 20.53939 6 10.5004 54 109 68.50 81.0 95.75
## P6_Dry 86.16667 19.82339 6 10.5004 56 110 75.00 89.0 99.25
## P6_Rainy 81.50000 17.14351 6 10.5004 62 107 68.25 81.5 90.25
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 50.49091
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 144.50000 a
## P3_Dry 108.83333 ab
## P2_Rainy 98.66667 ab
## P5_Dry 97.83333 ab
## P4_Rainy 95.66667 abc
## P6_Dry 86.16667 bc
## P5_Rainy 81.66667 bc
## P6_Rainy 81.50000 bc
## P2_Dry 78.33333 bc
## P3_Rainy 76.66667 bc
## P1_Dry 75.50000 bc
## P1_Rainy 46.33333 clm.q2<-qs %>% filter(order=="q2")
model2<- lm(val ~ Poligono*Season, data = lm.q2)
car::leveneTest(val ~ Poligono*Season, data = lm.q2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.2133 0.2983
## 60
#shapiro.test(model2$residuals)
a2<- aov(model2)
a22<-summary(a2)
pvalues2<-a22[[1]][[5]]
modelinteract2<- aov(lm(val ~ interact, data = lm.q2))
let2<-agricolae::HSD.test(modelinteract2, "interact", console = TRUE)$groups%>%
mutate(order="q2")##
## Study: modelinteract2 ~ "interact"
##
## HSD Test for val
##
## Mean Square Error: 226.9083
##
## interact, means
##
## val std r se Min Max Q25 Q50 Q75
## P1_Dry 31.00000 8.221922 6 6.149639 22 45 25.25 31.0 33.00
## P1_Rainy 17.33333 12.769756 6 6.149639 5 38 8.50 13.0 24.25
## P2_Dry 26.66667 10.013324 6 6.149639 13 38 19.50 28.0 34.25
## P2_Rainy 40.33333 18.672618 6 6.149639 15 69 29.50 41.5 47.50
## P3_Dry 42.16667 19.813295 6 6.149639 22 65 24.75 41.5 58.25
## P3_Rainy 29.66667 16.693312 6 6.149639 10 52 15.50 33.5 38.00
## P4_Dry 54.66667 23.457763 6 6.149639 18 79 41.00 61.0 71.25
## P4_Rainy 39.16667 10.980285 6 6.149639 25 56 33.25 37.5 44.75
## P5_Dry 38.00000 15.283979 6 6.149639 20 57 25.75 38.0 49.50
## P5_Rainy 31.66667 13.559744 6 6.149639 17 55 23.50 28.5 36.50
## P6_Dry 35.33333 12.596296 6 6.149639 21 55 27.25 32.5 42.25
## P6_Rainy 32.50000 11.184811 6 6.149639 19 46 22.75 35.5 39.25
##
## Alpha: 0.05 ; DF Error: 60
## Critical Value of Studentized Range: 4.808477
##
## Minimun Significant Difference: 29.57039
##
## Treatments with the same letter are not significantly different.
##
## val groups
## P4_Dry 54.66667 a
## P3_Dry 42.16667 ab
## P2_Rainy 40.33333 ab
## P4_Rainy 39.16667 ab
## P5_Dry 38.00000 ab
## P6_Dry 35.33333 ab
## P6_Rainy 32.50000 ab
## P5_Rainy 31.66667 ab
## P1_Dry 31.00000 ab
## P3_Rainy 29.66667 ab
## P2_Dry 26.66667 ab
## P1_Rainy 17.33333 b#letters_int<-rbind(let0, let1, let2)#%>% rownames_to_column(var="interact") %>%
#dplyr::rename(mean=val)library(ggh4x)
library(rcartocolor)
#colors<-viridis::turbo(6, alpha = 1, begin = 0, end = 1, direction = 1)
colors = carto_pal(6, "Safe")
my_colors=carto_pal(2, "Earth")
my_colors## [1] "#A16928" "#2887a1"strip <- strip_themed(background_x = elem_list_rect(fill = colors))
labels <- let0 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=600)
a <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q0") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues0[1], digits = 3), "; ",
"Season =",signif(pvalues0[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues0[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))
labels <- let1 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=180)
b <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q1") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues1[1], digits = 3), "; ",
"Season =",signif(pvalues1[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues1[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))+
theme(strip.text.x = element_blank())
#strip.background = element_rect(color="black", size=1.5, linetype="solid" ))
labels <- let2 %>% rownames_to_column(var = "interact") %>%
separate("interact", c("Poligono", "Season"), remove = F)
labels<- labels %>% mutate(val=70)
c <- qs %>%
mutate(
orders = case_when(
order == "q0" ~ "q=0",
order == "q1" ~ "q=1",
order == "q2" ~ "q=2" )) %>%filter(
order=="q2") %>%
ggbarplot(x = "Season", y = "val", fill="Season", facet.by = "Poligono", add=c("mean_se", "jitter"))+
facet_grid2(vars(orders), vars(Poligono), scales = "free", strip = strip)+
theme(axis.text.x = element_blank() ,
axis.title.x = element_blank(),
axis.ticks.x = element_blank())+
ylab("Effective number of OTUs")+
stat_compare_means(aes(group = Season), method = "t.test",label = "p.format")+
scale_fill_manual(values = c("#A16928" ,"#2887a1"))+theme_grey()+
theme(axis.text.x = element_text(size = 14, colour = "black"),
# axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face="bold.italic", colour = "white"),
axis.title = element_text(size = 18))+
geom_label(data = labels, aes(label=groups),
hjust = 1, vjust = 1)+
ggtitle(paste0("Polygon =",signif(pvalues2[1], digits = 3), "; ",
"Season =",signif(pvalues2[2], digits = 3), "; ",
"Polygon*Season =",signif(pvalues2[3], digits = 3)))+ theme(
strip.background.y = element_rect(
color="grey", fill="black", size=1, linetype="solid"),
strip.background.x = element_rect(colour = "black", linetype = "solid", size = 1))+
theme(strip.text.x = element_blank())library(cowplot)
alpha<- plot_grid(a + theme(legend.position = "none", axis.title.x = element_blank())+ylab(""),
b+theme(legend.position = "none", axis.title.x = element_blank(),
axis.title.y = element_text(size = 15)),
c+ theme(legend.position = "none")+ylab(""),
labels = c("A", "B", "C"),
ncol = 1)
alpha
Heatmap
- Load data and packages and format
library(tidyverse)
library(qiime2R)
library(ComplexHeatmap)
library(circlize)
library(viridis)
library(ggh4x)
tab <- read_qza("../Data/table_SMOQ.qza")$data %>% as.data.frame() %>%
rownames_to_column(var = "Feature.ID")
taxa<- qiime2R::read_qza("../Data/taxonomy_SMOQ_ver10.qza")$data %>% as.data.frame()
#tab <- read_qza("../Data/table_SMDR.qza")$data %>% as.data.frame() %>%
#rownames_to_column(var = "Feature.ID")
#taxa<- qiime2R::read_qza("../Data/taxonomy_SMDR.qza")$data %>% as.data.frame()
tab_tax<- tab %>% inner_join(taxa)
tab_taxa<- tab_tax %>% dplyr::select(OTUID=Feature.ID, `121D`:`221D`, taxonomy=Taxon)
#tab_taxa<- tab_tax %>% dplyr::select(OTUID=Feature.ID, `111D`:`623R`, taxonomy=Taxon)
#guilds
#guild1<- read.delim("../Data/otus/table_otus.guilds_matched.txt", check.names = F) %>%
# filter(!"Trophic Mode"=="na") %>%
#mutate(Guild = str_replace_all(Guild, "\\|", ""))
guild1<- read.delim("../Data/asvs/table_asvs.guilds_matched.txt", check.names = F) %>%
filter(!"Trophic Mode"=="na") %>%
mutate(Guild = str_replace_all(Guild, "\\|", ""))
guild<- guild1 %>% group_by(Genus=Taxon, Guild) %>% dplyr::count()
trophic<- guild1 %>% group_by(Genus=Taxon, `Trophic Mode`) %>% dplyr::count()
relabunda<- function(x){(as.data.frame(t(t(x)/colSums(x)))*100)}
metadata<- read.delim("../Data/its_map.txt")
vector_order<- c("111D" , "112D", "113D" ,"121D", "122D" ,"123D", "111R" ,"112R", "113R", "121R", "122R" ,"123R",
"211D" ,"212D" ,"213D" ,"221D" ,"222D" ,"223D", "211R" ,"212R" ,"213R", "221R", "222R" ,"223R",
"311D" ,"312D" ,"313D" ,"321D", "323D" ,"322D","311R" ,"312R" ,"313R" ,"321R", "322R" ,"323R" ,
"411D" ,"412D" ,"413D" ,"421D" ,"422D" ,"423D", "411R", "412R" ,"413R" ,"421R" ,"422R" ,"423R",
"511D" ,"512D" ,"513D" ,"521D" ,"522D", "523D", "511R", "512R" ,"513R" ,"521R", "522R" ,"523R",
"611D" ,"612D" ,"613D", "621D" ,"622D", "623D", "611R", "612R","613R", "621R" ,"622R", "623R")
#ordering
#taxones_color<- read_csv("../Data/taxones_color.csv") %>% dplyr::rename(Genus=Taxon)
table_genus <- tab_taxa %>% separate(taxonomy,
c(
"Kingdom",
"Phylum",
"Order",
"Class",
"Family",
"Genus",
"Species"
),
sep = ";") %>% mutate_at(c("Genus"), ~ str_replace(., "g__", "")) %>%
dplyr::mutate(Genus = stringr::str_trim(Genus, side = "both")) %>% mutate_if(is.character, ~
replace_na(., "Unassigned")) %>% group_by(Genus) %>% summarise_if(is.numeric, sum) %>%
column_to_rownames(var = "Genus") %>% t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
inner_join(metadata) %>% group_by(Poligono, Season) %>%
summarise_if(is.numeric, sum) %>%
unite("ids", Poligono:Season, sep = ".") %>% column_to_rownames(var = "ids") %>%
t() %>% as.data.frame() %>% mutate(all = rowSums(.)) %>% dplyr::arrange(-all) %>% relabunda(.) %>%
as.data.frame() %>% rownames_to_column(var = "Genus") %>% filter(!Genus ==
"unidentified" ,
!Genus == "Unassigned") %>%
filter(!grepl('Incertae_sedis', Genus)) %>% slice(c(1:50)) %>% pivot_longer(.,
cols = -Genus,
names_to = "SampleID",
values_to = "relab") %>% filter(!SampleID ==
"all")
#cols <-
# table_genus %>% inner_join(taxones_color) %>% arrange(Genus)
#col <- as.character(cols$color)
#names(col) <- as.character(cols$Genus)
barplot_genus<- table_genus %>% dplyr::select(
Genus, SampleID, relab) %>% pivot_wider(
names_from = Genus, values_from = relab) %>% column_to_rownames(
var = "SampleID") %>% t() %>% as.data.frame()
data_fungi_order=barplot_genus%>%
rownames_to_column(var="Genus")
merge_data<- data_fungi_order %>% mutate(Genus=case_when(
Genus=="Leptosphaeria maculans species complex" ~ "Leptosphaeria",
TRUE ~ as.character(Genus)))%>% column_to_rownames(
var = "Genus") %>% mutate(proms=rowMeans(.)) %>% dplyr::arrange(-proms) %>% slice(
1:80)%>% dplyr::select(-proms) %>% mutate_all(., funs(R = case_when(
. <= 0.001 ~ 0,
. > 0.001 & . <= 0.005 ~ 1,
. > 0.005 & . <= 0.01 ~ 2,
. > 0.01 & . <= 0.10 ~ 3,
. > 0.10 & . <= 0.20 ~ 4,
. > 0.20 & . <= 1.00 ~ 5,
. > 1.00 & . <= 2.00 ~ 6,
. > 2.00 & . <= 5.00 ~ 7,
. > 5.00 & . <= 10.00 ~ 8,
. > 10.00 & . <= 25.00 ~ 9,
. > 25.00 & . <= 50.00 ~ 10,
. > 50.00 & . <= 75.00 ~ 11,
. > 75.00 ~ 12))) %>%select_at(
vars(contains("_R"))) %>% select_all(~str_replace(., "_R", ""))- Constructing data and annotations for heatmap
#annotatio seccion y origen (columnas)
annotation_columns<- data.frame(id=colnames(merge_data))
rownames(annotation_columns) <- colnames(heatmap)
#set.seed(123)
#split = rep(1:6, each = 12)
split= c(rep("1D", 1),rep("1R", 1),
rep("2D", 1),rep("2R", 1),
rep("3D", 1),rep("3R", 1),
rep("4D", 1),rep("4R", 1),
rep("5D", 1),rep("5R", 1),
rep("6D", 1),rep("6R", 1) )
split=c("P1", "P2", "P3", "P4", "P5", "P6")
Pol<- c(1,1,
2,2,
3,3,
4,4,
5,5,
6,6)
#pal<- viridis(6, option = "H")
library(rcartocolor)
pal=rcartocolor::carto_pal(n = 6, name = "Safe")
ha = HeatmapAnnotation("Pol" = anno_block(gp = gpar(
fill = pal, col="black"),
labels = c("P1", "P2", "P3", "P4","P5", "P6" ),
labels_gp = gpar(col = "white", fontsize = 11, fontface= "bold")))
cols_ho<- list("Season" = c("Dry"="#A16928","Rainy"="#2887a1"))
ho = HeatmapAnnotation("Season" = c(rep("Dry", 1), rep("Rainy", "1"),
rep("Dry", 1), rep("Rainy", "1"),
rep("Dry", 1), rep("Rainy", "1"),
rep("Dry", 1), rep("Rainy", "1"),
rep("Dry", 1), rep("Rainy", "1"),
rep("Dry", 1), rep("Rainy", "1")),
which = "col", col = cols_ho,
annotation_name_gp = gpar(fontsize=8,
fontface="bold"),
show_legend = F, gp = gpar(
col = "white", fontize=12),
simple_anno_size = unit(0.25, "cm"),
show_annotation_name = T)
#annotatio row (filas)
guilds <- guild %>%
mutate(Guilds = case_when(
Guild == "Dung Saprotroph-Undefined Saprotroph-Wood Saprotroph" ~ "Dung Sap.-Undef Sap.-Wood Sap.",
Guild == "Ectomycorrhizal-Fungal Parasite" ~ "Ectomycorrhizal-Fung.Parasit",
Guild == "Animal Pathogen-Endophyte-Endosymbiont-Epiphyte-Soil Saprotroph-Undefined Saprotroph" ~ "Animal P-Endoph/symb-Undef. Sap.",
Guild == "Animal Pathogen-Endophyte-Fungal Parasite-Undefined Saprotroph" ~ "Animal P-Endoph/symb-Undef. Sap.",
Guild == "Animal Pathogen-Endophyte-Fungal Parasite-Plant Pathogen-Wood Saprotroph" ~ "Animal/Plant P-Endoph-Fungal Para-Wood. Sap.",
Guild == "Bryophyte Parasite-Dung Saprotroph-Ectomycorrhizal-Fungal Parasite-Leaf Saprotroph-Plant Parasite-Undefined Saprotroph-Wood Saprotroph" ~ "Bryophyte Para-Dung Sap.-Ecto-Para-Undef. Sap.",
Guild== "Animal Parasite-Dung Saprotroph-Endophyte-Fungal Parasite-Plant Pathogen-Undefined Saprotroph-Wood Saprotroph"~ "Animal Par-Pathogen-Undef. Sap",
Guild == "Animal Parasite-Animal Pathogen-Undefined Saprotroph" ~ "Animal Par-Pathogen-Undef. Sap",
Guild=="Animal Parasite-Animal Pathogen-Endophyte-Plant Saprotroph-Undefined Saprotroph" ~ "Animal Par-Pathogen-Undef. Sap",
Guild == "Dung Saprotroph-Ectomycorrhizal-Soil Saprotroph-Wood Saprotroph" ~ "Dung Sap.-Ectomycorrhizal-Soil Sap.-Wood Sap.",
Guild=="Ectomycorrhizal-Plant Saprotroph-Wood Saprotroph" ~"Ectomycorrhizal-Undef. Sap.",
Guild == "Endophyte-Plant Pathogen-Wood Saprotroph" ~ "Endophyte-Pathogen-Saprotroph",
Guild=="Endophyte-Plant Saprotroph-Undefined Saprotroph" ~ "Endophyte-Pathogen-Saprotroph",
Guild == "Animal Pathogen-Endophyte-Epiphyte-Fungal Parasite-Plant Pathogen-Wood Saprotroph" ~ "Animal-Plant P-Endophyte-Saprotroph",
Guild == "Endophyte-Litter Saprotroph-Soil Saprotroph-Undefined Saprotroph" ~ "Endophyte-Undefined Sap.",
Guild=="Endophyte-Ericoid Mycorrhizal-Undefined Saprotroph" ~ "Endophyte-Undefined Sap.",
Guild == "Animal Endosymbiont-Animal Pathogen-Plant Pathogen-Undefined Saprotroph" ~ "Animal-Plant endosymb-P-Undef. Sap.",
Guild == "Endophyte-Plant Pathogen-Undefined Saprotroph" ~ "Endophyte-Plant Pathogen-Undef. Sap.",
Guild=="Ectomycorrhizal-Plant Saprotroph"~"Ectomycorrhizal-Undef. Sap.",
Guild=="Ectomycorrhizal-Endophyte-Plant Pathogen-Plant Saprotroph" ~ "Ectomycorrhizal-Undef. Sap.",
Guild == "Ectomycorrhizal-Undefined Saprotroph" ~ "Ectomycorrhizal-Undef. Sap.",
Guild== "Ectomycorrhizal-Ericoid Mycorrhizal-Plant Pathogen-Plant Saprotroph-Undefined Saprotroph-Wood Saprotroph" ~ "Ectomycorrhizal-Undef. Sap.",
Guild == "Animal Pathogen-Fungal Parasite-Undefined Saprotroph" ~ "Animal P-Fungal Parasite-Undef. Sap.",
Guild == "Animal Pathogen-Endophyte-Plant Pathogen-Wood Saprotroph" ~ "Animal P-Fungal Parasite-Undef. Sap.",
Guild == "Fungal Parasite-Plant Pathogen-Plant Saprotroph" ~ "Fungal Parasite-Plant Path-Plant Sap.",
Guild == "Leaf Saprotroph-Plant Pathogen-Undefined Saprotroph-Wood Saprotroph" ~ "Leaf Sap.-Plant P/Undef Sap.",
Guild == "Dung Saprotroph-Ectomycorrhizal-Litter Saprotroph-Undefined Saprotroph" ~ "Dung Sap.-Ectomycorrhizal-Soil Sap.-Wood Sap.",
Guild == "Endophyte-Epiphyte-Fungal Parasite-Insect Parasite" ~ "Endophyte-Epiphyte-Fungal/Insect Para",
Guild == "Animal Pathogen-Endophyte-Plant Pathogen-Soil Saprotroph-Wood Saprotroph" ~ "Animal Pathogen-Endoph-Plant P-Soil/Wood Sap",
Guild == "Ectomycorrhizal-Orchid Mycorrhizal-Root Associated Biotroph" ~ "Ectomycorrhizal-Mycorrhizal-Biotroph",
Guild == "Clavicipitaceous Endophyte-Plant Pathogen" ~ "Clavicipitaceous Endophyte-Plant Pat.",
Guild=="Plant Pathogen-Plant Saprotroph-Undefined Saprotroph-Wood Saprotroph" ~ "Plant Pathogen-Undefined Saprotroph",
Guild=="Animal Pathogen-Dung Saprotroph"~"Animal Pathogen-Soil Saprotroph",
Guild=="Animal Pathogen-Dung Saprotroph-Endophyte-Epiphyte-Plant Pathogen-Plant Saprotroph-Wood Saprotroph"~"Animal Pathogen-Soil Saprotroph",
TRUE ~ as.character(Guild)
))
cols_guild <- list('Guilds' = c(
"Unassigned"= "#85929e",
#"Animal Pathogen-Endophyte-Plant Pathogen-Soil Saprotroph-Wood Saprotroph" = "#FCDE9C",
"Animal Pathogen-Endoph-Plant P-Soil/Wood Sap" = "#FAAC7B",
"Animal Pathogen-Endophyte-Epiphyte-Fungal Parasite-Plant Pathogen-Wood Saprotroph" = "#F28170",
"Animal Pathogen-Soil Saprotroph" = "#FCDE9C",
"Animal Endosymbiont-Animal Pathogen-Plant Pathogen-Undefined Saprotroph" = "#E04572",
"Animal Pathogen" = "#D23377",
"Animal Pathogen-Endophyte-Endosymbiont-Epiphyte-Soil Saprotroph-Undefined Saprotroph" = "#B02378",
"Animal P-Endoph/symb-Undef. Sap." = "#F28170",
"Animal Pathogen-Undefined Saprotroph" = "#FCDE9C" ,
"Plant Pathogen" = "#FDE0C5",
"Fungal Parasite" = "#FACEAA",
"Animal-Plant P-Endophyte-Saprotroph" = "#F8BB92",
"Animal-Plant endosymb-P-Undef. Sap." = "#F6A77C",
"Animal Par-Pathogen-Undef. Sap" = "#F39369",
"Animal P-Fungal Parasite-Undef. Sap." = "#F17D58",
"Fungal Parasite-Plant Pathogen-Plant Saprotroph" = "#EE654A",
"Fungal Parasite-Plant Path-Plant Sap." = "#EB4A40",
"Plant Pathogen-Undefined Saprotroph" = "#EB4A40" ,
"Ectomycorrhizal-Orchid Mycorrhizal-Root Associated Biotroph"="#b4ff68",
"Ectomycorrhizal-Mycorrhizal-Biotroph"="#b4ff68",
"Ectomycorrhizal" = "#7FC97F",
"Ectomycorrhizal-Fung.Parasit"="#58d68d",
"Ectomycorrhizal-Undefined Saprotroph"="#edf2a3",
"Ectomycorrhizal-Undef. Sap."="#008000",
"Ericoid Mycorrhizal"="#a2ab16",
"Arbuscular Mycorrhizal"="#008000",
#"Ectomycorrhizal-Orchid"="#0000FF", #no sale
"Ectomycorrhizal-Orchid Mycorrhizal-Root Associated Biotroph" = "#B0F2BC",
"Ectomycorrhizal-Mycorrhizal-Biotroph" = "#8EE9AE",
"Ectomycorrhizal" = "#70DEA7",
"Ectomycorrhizal-Fungal Parasite" = "#57D0A3",
"Ectomycorrhizal-Undefined Saprotroph" = "#43BEA3",
"Ectomycorrhizal-Undef. Sap." = "#34AAA2",
"Ericoid Mycorrhizal" = "#2A949E",
"Arbuscular Mycorrhizal" = "#257D98",
#"Epiphyte"="#02f0a5",
#"Lichenized"="#DAF7A6",
"Endophyte" = "#D1EEEA",
"Endophyte-Plant Pathogen" = "#ADDDDB",
"Endophyte-Plant Pathogen-Undefined Saprotroph" = "#8FCACD",
"Endophyte-Plant Pathogen-Undef. Sap." = "#74B5BF",
"Endophyte-Pathogen-Saprotroph" = "#5D9FB0",
"Endophyte-Plant Pathogen-Wood Saprotroph" = "#49879F",
"Endophyte-Litter Saprotroph-Soil Saprotroph-Undefined Saprotroph" = "#386E8B",
"Endophyte-Undefined Sap." = "#2A5674",
"Clavicipitaceous Endophyte-Plant Pat."="#8FCACD",
"Dung Saprotroph" = "#EDE5CF",
"Soil Saprotroph" = "#E1C7A8",
"Undefined Saprotroph" = "#D6A68B",
"Wood Saprotroph" = "#C88677",
"Plant Saprotroph-Undefined Saprotroph" = "#B56769",
"Dung Saprotroph-Wood Saprotroph" = "#9B4B5D",
"Dung Sap.-Endo.-Epiph.-Wood Sap." = "#7A3350",
"Dung Sap.-Undef Sap.-Wood Sap." = "#541F3F",
"Undefined Saprotroph-Undefined Symbiotroph" = "#EDE5CF" ,
"Undefined Saprotroph-Wood Saprotroph"= "#B56769"
))
cols_tro <- list('Trophic' = c(
"Pathotroph"= "#EE654A",
"Pathotroph-Saprotroph"= "#541F3F",
"Pathotroph-Symbiotroph"="#EB4A40",
"Saprotroph"= "#C88677",
"Saprotroph-Symbiotroph"= "#EDE5CF",
"Symbiotroph"= "#008000",
"Unassigned"= "#85929e",
"Pathotroph-Saprotroph-Symbiotroph"= "#E04572"))
table_phylum <- tab_taxa %>% separate(
taxonomy, c("Kingdom","Phylum","Order","Class","Family","Genus","Species"), sep = ";" ) %>% mutate_at(
c("Phylum"), ~str_replace(., "p__", "") )%>% mutate_at(
c("Genus"), ~str_replace(., "g__", ""))%>% group_by(Genus, Phylum) %>% dplyr::count()
annotation_rows<- merge_data %>% rownames_to_column(
var = "Genus") %>% left_join(
guilds) %>% left_join(table_phylum, by = "Genus") %>% dplyr::select_if(
is.character) %>% replace(
is.na(.), "Unassigned") %>% left_join(trophic) %>% column_to_rownames(
var = "Genus") %>% dplyr::rename(Trophic=`Trophic Mode`)
rcartocolor::carto_pal(n = 4, name = "ag_Sunset")## [1] "#4B2991" "#C0369D" "#FA7876" "#EDD9A3"cols_phy<-list("Phylum"=c(
"Basidiomycota" ="#4B2991",
"Ascomycota"="#C0369D" ,
"Mortierellomycota"="#FA7876" ,
"Mucoromycota"="#EDD9A3"))
annphyl= HeatmapAnnotation("Phylum" = annotation_rows$Phylum,
which = "row", col = cols_phy,
show_legend = T,
show_annotation_name = T,
annotation_name_gp =gpar(
fontsize = 7, fontface="bold"),
annotation_legend_param = list(
title_gp = gpar(fontsize = 7,
fontface="bold"),
labels_gp = gpar(fontsize = 7),
direction ="horizontal"
#,
#grid_width = unit(0.3, "cm"),
#grid_height = unit(0.1, "cm")
),
simple_anno_size = unit(0.3, "cm"),
gp = gpar(col = "white"))
annguild = HeatmapAnnotation("Guilds" = annotation_rows$Guilds,
which = "row", col = cols_guild,
show_legend = T,
show_annotation_name = T,
annotation_name_gp =gpar(
fontsize = 7, fontface="bold"),
annotation_legend_param = list(
title_gp = gpar(fontsize = 7,
fontface="bold"),
labels_gp = gpar(fontsize = 7),
direction ="horizontal"#,
# grid_width = unit(0.45, "cm"),
# grid_height = unit(0.1, "cm")
),
simple_anno_size = unit(0.3, "cm"),
gp = gpar(col = "white"))
anntro = HeatmapAnnotation("Trophic" = annotation_rows$Trophic,
which = "row", col = cols_tro,
show_legend = T,
show_annotation_name = T,
annotation_name_gp =gpar(
fontsize = 7, fontface="bold"),
annotation_legend_param = list(
title_gp = gpar(fontsize = 7,
fontface="bold"),
labels_gp = gpar(fontsize =7),
direction ="horizontal"#,
# grid_width = unit(0.45, "cm"),
#grid_height = unit(0.1, "cm")
),
simple_anno_size = unit(0.3, "cm"),
gp = gpar(col = "white"))
my_palette <- carto_pal(8, "ag_Sunset")
my_palette <- rev(my_palette)
heats<-ComplexHeatmap::Heatmap(
merge_data,
col = my_palette,
row_dend_width = unit(0.4, "cm"),
width = ncol(merge_data)*unit(6, "mm"),
#height = nrow(merge_data)*unit(2.4, "mm"),
heatmap_legend_param = list(direction = "horizontal",
title = "Relative abund(%)",
grid_height = unit(0.2, "cm"),
legend_height = unit(1, "cm"),
labels_gp = gpar(fontsize = 7),
title_gp = gpar(fontsize = 6,
fontface="bold"),
at = c(0,1,2,3,5,8,10, 50,100),
break_dist = 3),
rect_gp = gpar(col = "white"),
cluster_columns = F, cluster_rows = F,
show_heatmap_legend = T, top_annotation = c(ha,ho),
right_annotation = c(anntro, annguild, annphyl),
#column_order = sort(colnames(merge_data)),
column_split = c(1,1,2,2,3,3,4,4,5,5,6,6) ,
column_title = NULL,
show_column_names = T,
row_names_gp = gpar(fontsize=8, fontface="italic"),
column_names_rot=0,
column_labels = gt_render(rep(c("D", "R"), 6)),
column_names_centered = TRUE,
column_names_gp = gpar(fontsize=9, fontface="bold"))
lgd1 = Legend(at = c("Dry", "Rainy"),
title = "Season", nrow = 1,
title_position = "leftcenter",
legend_gp = gpar(fill = cols_ho$Season),
labels_gp = gpar(fontsize = 8),
title_gp = gpar(fontsize = 9, fontface="bold"))
draw(heats, heatmap_legend_side = "right",
annotation_legend_side = "top",
merge_legend=F,
annotation_legend_list = list(lgd1))
heatm<-grid.grabExpr(draw(heats, heatmap_legend_side = "right",
annotation_legend_side = "top",
merge_legend=F,
annotation_legend_list = list(lgd1)))
heatm## gTree[GRID.gTree.3973]ALDEx2
- Load data and packages and format
library(tidyverse)
library(qiime2R)
library(ComplexHeatmap)
library(circlize)
library(viridis)
library(ggh4x)
tab <- read_qza("../Data/table_SMOQ_filt.qza")$data %>% as.data.frame() %>%
rownames_to_column(var = "Feature.ID")
taxa<- qiime2R::read_qza("../Data/taxonomy_SMOQ_ver10.qza")$data %>% as.data.frame()
metadata<- read.delim("../Data/its_map.txt")
tab_tax<- tab %>% inner_join(taxa)
tab_taxa<- tab_tax %>% dplyr::select(OTUID=Feature.ID, `121D`:`221D`, taxonomy=Taxon)
vector_order<- c("111D" , "112D", "113D" ,"121D", "122D" ,"123D", "111R" ,"112R", "113R", "121R", "122R" ,"123R",
"211D" ,"212D" ,"213D" ,"221D" ,"222D" ,"223D", "211R" ,"212R" ,"213R", "221R", "222R" ,"223R",
"311D" ,"312D" ,"313D" ,"321D", "323D" ,"322D","311R" ,"312R" ,"313R" ,"321R", "322R" ,"323R" ,
"411D" ,"412D" ,"413D" ,"421D" ,"422D" ,"423D", "411R", "412R" ,"413R" ,"421R" ,"422R" ,"423R",
"511D" ,"512D" ,"513D" ,"521D" ,"522D", "523D", "511R", "512R" ,"513R" ,"521R", "522R" ,"523R",
"611D" ,"612D" ,"613D", "621D" ,"622D", "623D", "611R", "612R","613R", "621R" ,"622R", "623R")
#ordering
tab_taxa_order = tab_taxa %>% dplyr::select(OTUID, vector_order, taxonomy)
tab_order_metadata = tab_taxa_order %>% column_to_rownames(var = "OTUID") %>% t() %>%
as.data.frame() %>% rownames_to_column(var = "SampleID") %>% inner_join(metadata)
#season
tab_season = tab_order_metadata %>% arrange(Season) %>%dplyr::select(-overhang:-LIB) %>% column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_if(is.character, as.numeric)
vector_season = c(rep("Dry", 36), rep("Rainy", 36))
#polygon
tab_polygon = tab_order_metadata %>% arrange(Poligono)%>% dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric)
vector_polygon = c(rep("P1", 12),rep("P2", 12), rep("P3", 12),
rep("P4", 12),rep("P5", 12),rep("P6", 12))
#season-polygon
tab_polygon_dry = tab_order_metadata %>% filter(Season=="Dry") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_polygon_dry = c(rep("P1", 6),rep("P2", 6), rep("P3", 6),
rep("P4", 6),rep("P5", 6),rep("P6", 6))
tab_polygon_rainy = tab_order_metadata %>% filter(Season=="Rainy") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_polygon_rainy = c(rep("P1", 6),rep("P2", 6), rep("P3", 6),
rep("P4", 6),rep("P5", 6),rep("P6", 6))
#polygon-season
tab_p1 = tab_order_metadata %>% filter(Poligono=="P1") %>% arrange(Season)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
tab_p2 = tab_order_metadata %>% filter(Poligono=="P2") %>% arrange(Season)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
tab_p3 = tab_order_metadata %>% filter(Poligono=="P3") %>% arrange(Season)%>%dplyr::select(-overhang:-LIB) %>% column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric)
tab_p4 = tab_order_metadata %>% filter(Poligono=="P4") %>% arrange(Season)%>%dplyr::select(-overhang:-LIB) %>% column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric)
tab_p5 = tab_order_metadata %>% filter(Poligono=="P5") %>% arrange(Season) %>%dplyr::select(-overhang:-LIB)%>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
tab_p6 = tab_order_metadata %>% filter(Poligono=="P6") %>% arrange(Season)%>%dplyr::select(-overhang:-LIB)%>% column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p = c(rep("Dry", 6), rep("Rainy", 6))#polygon-paired
tab_p12 = tab_order_metadata %>% filter(Poligono=="P1"|Poligono=="P2") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p12 = c(rep("P1", 12), rep("P2", 12))
tab_p13 = tab_order_metadata %>% filter(Poligono=="P1"|Poligono=="P3") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p13 = c(rep("P1", 12), rep("P3", 12))
tab_p14 = tab_order_metadata %>% filter(Poligono=="P1"|Poligono=="P4") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p14 = c(rep("P1", 12), rep("P4", 12))
tab_p15 = tab_order_metadata %>% filter(Poligono=="P1"|Poligono=="P5") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p15 = c(rep("P1", 12), rep("P5", 12))
tab_p16 = tab_order_metadata %>% filter(Poligono=="P1"|Poligono=="P6") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p16 = c(rep("P1", 12), rep("P6", 12))
tab_p23 = tab_order_metadata %>% filter(Poligono=="P2"|Poligono=="P3") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p23 = c(rep("P2", 12), rep("P3", 12))
tab_p24 = tab_order_metadata %>% filter(Poligono=="P2"|Poligono=="P4") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p24 = c(rep("P2", 12), rep("P4", 12))
tab_p25 = tab_order_metadata %>% filter(Poligono=="P2"|Poligono=="P5") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p25 = c(rep("P2", 12), rep("P5", 12))
tab_p26 = tab_order_metadata %>% filter(Poligono=="P2"|Poligono=="P6") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p26 = c(rep("P2", 12), rep("P6", 12))
tab_p34 = tab_order_metadata %>% filter(Poligono=="P3"|Poligono=="P4") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p34 = c(rep("P3", 12), rep("P4", 12))
tab_p35 = tab_order_metadata %>% filter(Poligono=="P3"|Poligono=="P5") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p35 = c(rep("P3", 12), rep("P5", 12))
tab_p36 = tab_order_metadata %>% filter(Poligono=="P3"|Poligono=="P6") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p36 = c(rep("P3", 12), rep("P6", 12))
tab_p45 = tab_order_metadata %>% filter(Poligono=="P4"|Poligono=="P5") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p45 = c(rep("P4", 12), rep("P5", 12))
tab_p46 = tab_order_metadata %>% filter(Poligono=="P4"|Poligono=="P6") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p46 = c(rep("P4", 12), rep("P6", 12))
tab_p56 = tab_order_metadata %>% filter(Poligono=="P5"|Poligono=="P6") %>% arrange(Poligono)%>%dplyr::select(-overhang:-LIB) %>%
column_to_rownames(var = "SampleID") %>% t() %>% as.data.frame() %>% mutate_all(as.numeric) %>%
filter(rowSums(.) != 0)
vector_p56 = c(rep("P5", 12), rep("P6", 12))library(ALDEx2)
aldex_p1 = aldex(tab_p1, conditions = vector_p, mc.samples = 128,denom = "all")
aldex_p2 = aldex(tab_p2, conditions = vector_p, mc.samples = 128,denom = "all")
aldex_p3 = aldex(tab_p3, conditions = vector_p, mc.samples = 128,denom = "all")
aldex_p4 = aldex(tab_p4, conditions = vector_p, mc.samples = 128,denom = "all")
aldex_p5 = aldex(tab_p5, conditions = vector_p, mc.samples = 128,denom = "all")
aldex_p6 = aldex(tab_p6, conditions = vector_p, mc.samples = 128,denom = "all")p1= plot_volcanop(aldex_p1, taxa)
p2= plot_volcanop(aldex_p2, taxa)
p3= plot_volcanop(aldex_p3, taxa)
p4= plot_volcanop(aldex_p4, taxa)
p5= plot_volcanop(aldex_p5, taxa)
p6= plot_volcanop(aldex_p6, taxa)library(cowplot)
all=plot_grid(p1 +xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
p2+xlab("")+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
p3 +xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
p4+ xlab("")+ylab("")+theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
p5+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
p6+ ylab("")+theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
ncol = 2, nrow = 3, align = "hv", labels = "AUTO")
all
aldex_polp = function(x, conds){
aldex = aldex(x, conditions = conds, mc.samples = 128,denom = "all")
#aldex_filt = aldex %>% filter(wi.ep <0.05)
return(aldex)
}
#x = aldex_polp(x=tab_p4, conds = vector_p)
aldex_polp2 = function(x, conds){
aldex = aldex(x, conditions = conds, mc.samples = 128,denom = "all")
aldex_filt = aldex %>% filter(wi.eBH <0.05)
return(aldex_filt)
}
library(dplyr)
library(purrr)
# Listar las tablas y sus correspondientes vectores
datos_list <- list(
tab_p12 = tab_p12,
tab_p13 = tab_p13,
tab_p14 = tab_p14,
tab_p15 = tab_p15,
tab_p16 = tab_p16,
tab_p23 = tab_p23,
tab_p24 = tab_p24,
tab_p25 = tab_p25,
tab_p26 = tab_p26,
tab_p34 = tab_p34,
tab_p35 = tab_p35,
tab_p36 = tab_p36,
tab_p45 = tab_p45,
tab_p46 = tab_p46,
tab_p56 = tab_p56
)
# Listar los vectores de condiciones
vectores_list <- list(
vector_p12 = vector_p12,
vector_p13 = vector_p13,
vector_p14 = vector_p14,
vector_p15 = vector_p15,
vector_p16 = vector_p16,
vector_p23 = vector_p23,
vector_p24 = vector_p24,
vector_p25 = vector_p25,
vector_p26 = vector_p26,
vector_p34 = vector_p34,
vector_p35 = vector_p35,
vector_p36 = vector_p36,
vector_p45 = vector_p45,
vector_p46 = vector_p46,
vector_p56 = vector_p56
)# Aplicar la función aldex_polp a cada conjunto de datos y sus correspondientes vectores
resultados <- map2(datos_list, vectores_list, aldex_polp)
#resultados2 <- map2(datos_list, vectores_list, aldex_polp2)library(cowplot)
pp1= plot_volcano(x = resultados[[1]], taxa = taxa)
pp2= plot_volcano(x = resultados[[2]], taxa = taxa)
pp3= plot_volcano(x = resultados[[3]], taxa = taxa)
pp4= plot_volcano(x = resultados[[4]], taxa = taxa)
pp5= plot_volcano(x = resultados[[5]], taxa = taxa)
pp6= plot_volcano(x = resultados[[6]], taxa = taxa)
pp7= plot_volcano(x = resultados[[7]], taxa = taxa)
pp8= plot_volcano(x = resultados[[8]], taxa = taxa)
pp9= plot_volcano(x = resultados[[9]], taxa = taxa)
pp10= plot_volcano(x = resultados[[10]], taxa = taxa)
pp11= plot_volcano(x = resultados[[11]], taxa = taxa)
pp12= plot_volcano(x = resultados[[12]], taxa = taxa)
pp13= plot_volcano(x = resultados[[13]], taxa = taxa)
pp14= plot_volcano(x = resultados[[14]], taxa = taxa)
pp15= plot_volcano(x = resultados[[15]], taxa = taxa)
pps=plot_grid(pp1+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp2+ylab("")+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp3+ylab("")+xlab("")+ theme(plot.margin = unit(c(0.01, 0.1, 0.01, 0.01), "cm")),
pp4+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp5+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm"))+xlab(""), pp6+ylab("")+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
pp7+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp8+ylab("")+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp9+xlab("")+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),
pp10+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),pp11+ylab("")+xlab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp12+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm"))+xlab(""),
pp13+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")),pp14+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), pp15+ylab("")+ theme(plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm")), ncol = 3, nrow = 5, align = "hv", labels = "AUTO")
pps
NMDS and perMANOVA - BETA
library(qiime2R)
library(tidyverse)
library(ALDEx2)
library(vegan)
library(rcartocolor)
#taxonomys
taxa<- read_qza("../Data/taxonomy_SMOQ_ver10.qza")$data
table<- read_qza("../Data/table_SMOQ_filt.qza")$data
set.seed(124)
metadata<- read.delim("../Data/its_map.txt")
#betadisper test
aldex.clr.transform <- aldex.clr(table, mc.samples = 2, denom="all",
verbose = FALSE, useMC=FALSE)
aldex.clr.transform.data<- t(getMonteCarloSample(aldex.clr.transform,1) )
mat<- dist(aldex.clr.transform.data, method = "euclidean")#nmds using aitchison
nmds1=metaMDS(t(aldex.clr.transform@reads), distance = "aitchison",trymax = 500, k = 2)nmds11= nmds1 %>%
vegan:::scores.metaMDS(display = "sites") %>%
as_tibble(., rownames = "sample")%>%
left_join(metadata, by = join_by(sample == SampleID)) %>%
group_by(Poligono, Season) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup()
set.seed(124)
tab<-aldex.clr.transform@reads %>% t() %>% as.data.frame() %>% rownames_to_column(
var = "SampleID") %>% inner_join(metadata)
permdisp<-betadisper(mat, tab$Poligono)
permdisp2<- permutest( permdisp, permutations = 999)
permdisp2##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 5 83.13 16.625 0.7328 999 0.606
## Residuals 66 1497.36 22.687perm<-adonis2(t(aldex.clr.transform@reads) ~Poligono*Season, data = tab, method =
"aitchison", permutations =999, by = "terms")
perm## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = t(aldex.clr.transform@reads) ~ Poligono * Season, data = tab, permutations = 999, method = "aitchison", by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Poligono 5 46659 0.21538 3.7362 0.001 ***
## Season 1 4273 0.01972 1.7106 0.011 *
## Poligono:Season 5 15845 0.07314 1.2687 0.008 **
## Residual 60 149864 0.69176
## Total 71 216640 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1titles<- paste0("adonis:", " Polygon, F=", round(perm$F[1], digits = 2), ", p-value=", round(perm$`Pr(>F)`[1], digits=3), ";", " Season, F=", round(perm$F[2], digits = 2), ", p-value=", round(perm$`Pr(>F)`[2], digits=3), ";", " \nPolygon*Season, F=", round(perm$F[3], digits = 2), ", p-value=", round(perm$`Pr(>F)`[3], digits=3))
source("../Scripts/gg_ordiplots2.R")
metadata1<- as.data.frame(nmds11) %>% dplyr::rename(SampleID=sample) %>%
inner_join(metadata)
y<-gg_ordiplot2(nmds1, groups = metadata1$Poligono, groups2 = metadata1$Season,ellipse = F,
spiders = T, hull = F, plot = F, scaling = 3)
xlab <- y$plot$labels$x
ylab <- y$plot$labels$y
z <- ggplot() + geom_point(
data = y$df_ord %>% rownames_to_column(var = "SampleID") %>%
inner_join(metadata1),
aes(
x = x,
y = y,
color = Group,
shape = Season
),
size = 4
) +
xlab("NMDS1") +
ylab("NMDS2") +
geom_segment(
data = y$df_spiders,
aes(
x = cntr.x,
xend = x,
y = cntr.y,
yend = y,
color = Group
),
show.legend = FALSE
) +
annotate(
"text", label = paste("Stress", round(nmds1$stress, digits = 3)),
x = 25, y = -15, size = 5, colour = "black"
)+
guides(color = guide_legend(title = "Sites")) + theme_linedraw() +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw() +
scale_color_carto_d(name = "Polygons: ", palette = "Safe") +
theme_grey() +
theme(
axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
panel.border = element_rect(
colour = "black",
fill = NA,
size = 1
)
) +
ggtitle(titles) + theme(
panel.background = element_rect(
fill = "#F2F4F4",
colour = "#F2F4F4",
size = 0.5,
linetype = "solid"
),
panel.grid.major = element_line(
size = 0.5,
linetype = 'solid',
colour = "white"
),
panel.grid.minor = element_line(
size = 0.25,
linetype = 'solid',
colour = "white"
)
) +
theme(aspect.ratio = 7 / 10)
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank())
z
Beta partition
#beta
library(tidyverse)
library(betapart)
library(qiime2R)
library(vegan)
library(colorRamps)
table<- read_qza("../Data/table_SMOQ_rar.qza")$data
metadata<- read.delim("../Data/its_map.txt")
#convrt data
presabs<- function(x){x[x>0]=1
return(t(x))}
tables_presabs<- presabs(table)
table_presabs<- tables_presabs %>% as.data.frame() %>% rownames_to_column(var = "SampleID") %>%
arrange(SampleID) %>% column_to_rownames(var = "SampleID")
#betapart function with jaccard
x<- beta.pair(table_presabs, index.family = "jaccard")
#extract each part
jac<- x$beta.jac
jtu<- x$beta.jtu
jne<- x$beta.jne
#jac<- x$beta.sor
#jtu<- x$beta.sim
#jne<- x$beta.sne
env1<- table_presabs%>% as.data.frame() %>% rownames_to_column(var="SampleID") %>% inner_join(metadata)
jacs<- betadisper(jac,factor(env1$Poligono))
jtus<- betadisper(jtu,factor(env1$Poligono))
jnes<- betadisper(jne,factor(env1$Poligono))
library(tidyverse)
library(ggordiplots)
function_plot_beta <- function(x, env) {
y <- gg_ordiplot(x, groups = env$Poligono, hull = FALSE, spiders = TRUE,
ellipse = FALSE, plot = FALSE, label = TRUE)
xlabs <- y$plot$labels$x
ylabs <- y$plot$labels$y
z <- ggplot() +
geom_point(data = y$df_ord %>% rownames_to_column(var = "SampleID") %>%
inner_join(env),
aes(x = x, y = y, color = Group, shape = Season), size = 3) +
xlab(xlabs) +
ylab(ylabs) +
geom_segment(data = y$df_spiders, aes(x = cntr.x, xend = x, y = cntr.y, yend = y, color = Group),
show.legend = FALSE)
a <- z +
geom_label(data = y$df_mean.ord, aes(x = x, y = y, label = Group)) +
guides(color = guide_legend(title = "Polygon")) +
theme_linedraw() +
geom_vline(xintercept = 0, linetype = 2) +
geom_hline(yintercept = 0, linetype = 2) +
scale_color_carto_d(name = "Polygons: ", palette = "Safe") +
theme_gray() +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
panel.background = element_rect(fill = "#F2F4F4", colour = "#F2F4F4", size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "white"),
panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "white"),
panel.border = element_rect(colour = "black", fill = NA, size = 1),
aspect.ratio = 7 / 10)
return(a)
}
#jaccard
plot_jac<- function_plot_beta(jacs, env1)
plot_jac=plot_jac+ guides(colour = guide_legend(nrow = 1, title = "Polygons"),
shape = guide_legend(nrow = 1) )+theme(legend.position = "top")
#turnover
plot_turn<- function_plot_beta(jtus, env1)
#nestedness
plot_nes<-function_plot_beta(jnes, env1)
heads=betapart::beta.multi(table_presabs, index.family = "jaccard")
library(cowplot)
leg<- get_legend(plot_jac)
a<-plot_grid(plot_jac+theme(legend.position = "none")+ylab("DIM2")+xlab("")+theme(aspect.ratio =6/10)+ggtitle(paste("Global Jaccard dissimilarity =", round(heads$beta.JAC[1], digits = 3))),
plot_turn+theme(legend.position = "none")+ylab("DIM2")+xlab("")+theme(aspect.ratio =6/10)+ggtitle(paste("Global Turnover component =", round(heads$beta.JTU[1], digits = 3))),
plot_nes+theme(legend.position = "none")+ylab("DIM2")+xlab("DIM1")+theme(aspect.ratio =6/10)+ggtitle(paste("Global Nestedness component =", round(heads$beta.JNE[1], digits = 3))),
ncol = 1, align = "hv", labels = c("A", "B", "C"), label_x = 0.1)
b= plot_grid(leg,a, ncol = 1, rel_heights = c(0.2,1))
b
#beta
library(tidyverse)
library(betapart)
library(qiime2R)
library(vegan)
library(colorRamps)
library(ggpubr)
library(ggordiplots)
table<- read_qza("../Data/table_SMOQ_rar.qza")$data
metadata<- read.delim("../Data/its_map.txt")
#convrt data
presabs<- function(x){x[x>0]=1
return(t(x))}
tables_presabs<- presabs(table)
table_presabs<- tables_presabs %>% as.data.frame() %>% rownames_to_column(var = "SampleID") %>%
arrange(SampleID) %>% column_to_rownames(var = "SampleID")#betapart season
nam<- colnames(table)
nam2<- str_extract(nam, pattern = "^\\d+")
nam3<- unique(nam2)
drys<- paste0(nam3,"D")
rainys<- paste0(nam3,"R")
tables_dry<- table_presabs %>% as.data.frame() %>%rownames_to_column(var="SampleID") %>%
filter(SampleID %in% drys) %>% mutate(SampleID = str_extract(SampleID, "^\\d+")) %>%
column_to_rownames(var = "SampleID")
tables_rainy<- table_presabs %>% as.data.frame() %>%rownames_to_column(var="SampleID") %>%
filter(SampleID %in% rainys) %>% mutate(SampleID = str_extract(SampleID, "^\\d+")) %>%
column_to_rownames(var = "SampleID")
beta.t <-beta.temp(tables_dry, tables_rainy,index.family="jaccard")
nams<-paste0(rownames(beta.t), "D")
table_jac<-rbind(beta.t[[3]])
colnames(table_jac)<-nams
table_turn<-rbind(beta.t[[1]])
colnames(table_turn)<-nams
table_nes<-rbind(beta.t[[2]])
colnames(table_nes)<-nams
table_jacs<- table_jac %>%t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata) %>%
dplyr::rename(beta=V1)## Joining with `by = join_by(SampleID)`table_turns<- table_turn %>%t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata) %>%
dplyr::rename(beta=V1)## Joining with `by = join_by(SampleID)`table_ness<- table_nes %>%t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata) %>%
dplyr::rename(beta=V1)## Joining with `by = join_by(SampleID)`model<-aov(beta~Poligono, data=table_turns)
out <- agricolae::HSD.test(model,"Poligono", group=TRUE,console=FALSE)
g<- data.frame(out$groups) %>% rownames_to_column(var = "Poligono") %>% dplyr::select(-beta)
data<- data.frame(Poligono=c("P1", "P2", "P3", "P4", "P5", "P6"),
beta=rep(0.85,6)) %>% inner_join(g, by = "Poligono")
library(ggpubr)
t<-ggboxplot(data = table_turns, x = "Poligono", y="beta", fill="Poligono")+
scale_fill_carto_d(name = "Polygons ", palette = "Safe") +
theme_grey() +theme(axis.text = element_text(colour = "black", size = 14),
axis.title = element_text(colour = "black", size = 16),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=1,))+
ylab("Temporal turnonver of Jaccard dissimilarity")+
xlab("Polygon")+
stat_compare_means(method = "anova", label.y = 0.86, size=5)+
geom_label(data = data,aes(Poligono, beta, label=groups), size=5)+
stat_summary(fun.y=mean, geom="point", shape=23, size=4, color="black", fill="white")## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Computation failed in `stat_summary()`.
## Caused by error in `call_f()`:
## ! Can't convert `fun`, a data frame, to a function.
## Warning: Computation failed in `stat_summary()`.
## Caused by error in `call_f()`:
## ! Can't convert `fun`, a data frame, to a function.Spatial turnover
# TURNOVER SPATIAL
library(qiime2R)
library(tidyverse)
library(ggpubr)
library(hilldiv)
library(rcartocolor)
table<- read_qza("../Data/table_SMOQ_rar.qza")$data
metadata<- read.delim("../Data/its_map.txt")
#convert data to incidence
presabs<- function(x){x[x>0]=1
return(t(x))}
tables_presabs<- presabs(table)
table_presabs<- tables_presabs %>% as.data.frame() %>% rownames_to_column(var = "SampleID") %>% arrange(SampleID) %>% column_to_rownames(var = "SampleID")library(betapart)
library(vegan)
#betapart function with jaccard
x<- beta.pair(table_presabs, index.family = "jaccard")
jac<- x$beta.jac
jtu<- x$beta.jtu
jne<- x$beta.jne## Run 0 stress 0.1636985
## Run 1 stress 0.1637304
## ... Procrustes: rmse 0.001647559 max resid 0.01081117
## Run 2 stress 0.1636985
## ... New best solution
## ... Procrustes: rmse 8.712242e-06 max resid 4.832512e-05
## ... Similar to previous best
## Run 3 stress 0.1672362
## Run 4 stress 0.1648359
## Run 5 stress 0.1648359
## Run 6 stress 0.1659383
## Run 7 stress 0.1659383
## Run 8 stress 0.1636985
## ... Procrustes: rmse 1.439475e-05 max resid 6.688429e-05
## ... Similar to previous best
## Run 9 stress 0.1636985
## ... Procrustes: rmse 1.208017e-05 max resid 5.491042e-05
## ... Similar to previous best
## Run 10 stress 0.1659383
## Run 11 stress 0.1636985
## ... Procrustes: rmse 3.077747e-06 max resid 1.092183e-05
## ... Similar to previous best
## Run 12 stress 0.1636986
## ... Procrustes: rmse 2.808939e-05 max resid 0.0001515109
## ... Similar to previous best
## Run 13 stress 0.1636985
## ... Procrustes: rmse 9.653192e-06 max resid 3.407932e-05
## ... Similar to previous best
## Run 14 stress 0.1636986
## ... Procrustes: rmse 1.734427e-05 max resid 7.781983e-05
## ... Similar to previous best
## Run 15 stress 0.1648358
## Run 16 stress 0.1636985
## ... Procrustes: rmse 8.902241e-06 max resid 5.040696e-05
## ... Similar to previous best
## Run 17 stress 0.1672411
## Run 18 stress 0.1648447
## Run 19 stress 0.1648447
## Run 20 stress 0.167241
## *** Best solution repeated 8 timesnmds_data <- data.frame(
NMDS1 = resultado$points[, 1],
NMDS2 = resultado$points[, 2] )
source("../Scripts/gg_ordiplots4.R")
metadata1 <- metadata[match(colnames(as.matrix(jtu)), metadata$SampleID),]
nmds_data_metadata<- nmds_data %>% rownames_to_column(var = "SampleID") %>%
inner_join(metadata) %>% column_to_rownames(var = "SampleID")# %>% filter(Season=="Dry")
nmds_data_metadata_dry<- nmds_data %>% rownames_to_column(var = "SampleID") %>%
inner_join(metadata) %>% column_to_rownames(var = "SampleID") %>% filter(Season=="Dry")
nmds_data_metadata_rainy<- nmds_data %>% rownames_to_column(var = "SampleID") %>%
inner_join(metadata) %>% column_to_rownames(var = "SampleID") %>% filter(Season=="Rainy")
y1<-gg_ordiplot4(nmds_data_metadata, groups = nmds_data_metadata$Poligono, groups2 = nmds_data_metadata$Season, spiders = T, plot = F, scaling = 3)
y2<-gg_ordiplot4(nmds_data_metadata_dry, groups = nmds_data_metadata_dry$Poligono, groups2 = nmds_data_metadata_dry$Season, spiders = T, plot = F, scaling = 3)
y3<-gg_ordiplot4(nmds_data_metadata_rainy, groups = nmds_data_metadata_rainy$Poligono, groups2 = nmds_data_metadata_rainy$Season, spiders = T, plot = F, scaling = 3)
#xlab <- y1$plot$labels$x
#ylab <- y1$plot$labels$y
z1<-ggplot()+ geom_point(data = y1$df_ord %>% rownames_to_column(var="SampleID") %>%
inner_join(metadata1),
aes(x = x*-1, y = y*-1, color = Group, shape=Season), size = 3) +
xlab("NMDS1") +
ylab("NMDS2")+
geom_segment(data = y1$df_spiders,
aes(x = cntr.x*-1, xend = x*-1, y = cntr.y*-1, yend = y*-1, color = Group),
show.legend = FALSE)+
guides(
color=guide_legend(title="Sites"))+theme_linedraw() +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw()+
#scale_fill_viridis_d(option ="turbo", name="Poligono")+#color of points
scale_color_carto_d(name = "Polygons", palette = "Safe") +
theme_grey()+#color of points
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(color = "black")) +
scale_shape_manual(values = c(16,17))
z2<-ggplot()+ geom_point(data = y2$df_ord %>% rownames_to_column(var="SampleID") %>%
inner_join(metadata1),
aes(x = x*-1, y = y*-1, color = Group, shape=Season), size = 3) +
xlab("NMDS1") +
ylab("NMDS2")+
geom_segment(data = y2$df_spiders,
aes(x = cntr.x*-1, xend = x*-1, y = cntr.y*-1, yend = y*-1, color = Group),
show.legend = FALSE)+
guides(
color=guide_legend(title="Sites"))+theme_linedraw() +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw()+
#scale_fill_viridis_d(option ="turbo", name="Poligono")+#color of points
scale_color_carto_d(name = "Polygons", palette = "Safe") +
theme_grey()+#color of points
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(color = "black")) +
scale_shape_manual(values = c(16,17))
z3<-ggplot()+ geom_point(data = y3$df_ord %>% rownames_to_column(var="SampleID") %>%
inner_join(metadata1),
aes(x = x*-1, y = y*-1, color = Group, shape=Season), size = 3) +
xlab("NMDS1") +
ylab("NMDS2")+
geom_segment(data = y3$df_spiders,
aes(x = cntr.x*-1, xend = x*-1, y = cntr.y*-1, yend = y*-1, color = Group),
show.legend = FALSE)+
guides(
color=guide_legend(title="Sites"))+theme_linedraw() +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw()+
#scale_fill_viridis_d(option ="turbo", name="Poligono")+#color of points
scale_color_carto_d(name = "Polygons", palette = "Safe") +
theme_grey()+#color of points
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(color = "black")) +
scale_shape_manual(values = c(17,16))#dendograms
library(usedist)
library(ggdendro)
source("../Scripts/ggdendro_funct.R")
dry_ids=metadata1 %>% filter(Season=="Dry") %>% dplyr::select(SampleID)
rainy_ids=metadata1 %>% filter(Season=="Rainy") %>% dplyr::select(SampleID)
jtu2 <- dist_subset(jtu, dry_ids$SampleID)
jtu3 <- dist_subset(jtu, rainy_ids$SampleID)
dendrogram <- hclust(as.dist(jtu), method='ward.D2')
dendrogram2 <- hclust(as.dist(jtu2), method='ward.D2')
dendrogram3 <- hclust(as.dist(jtu3), method='ward.D2')
hcdata <- dendro_data_k(dendrogram, 2)
hcdata2 <- dendro_data_k(dendrogram2, 2)
hcdata3 <- dendro_data_k(dendrogram3, 2)
clusters <- dendrogram
#clusters$labels
clusters=clusters$labels[clusters$order]
table_all = table[,match(clusters, colnames(table))]
table_all=table_all %>% t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata)
clusters2 <- dendrogram2
#clusters2$labels
clusters2=clusters2$labels[clusters2$order]
table_dry = table[,match(clusters2, colnames(table))]
table_dry=table_dry %>% t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata) %>% filter(Season=="Dry")
clusters3 <- dendrogram3
#clusters3$labels
clusters3=clusters3$labels[clusters3$order]
table_rainy = table[,match(clusters3, colnames(table))]
table_rainy=table_rainy %>% t() %>% as.data.frame() %>%
rownames_to_column(var = "SampleID") %>% inner_join(metadata)%>% filter(Season=="Rainy")
cols <- c("#2E86C1", "#800080","#30123BFF", "#3E9BFEFF", "#46F884FF", "#E1DD37FF", "#F05B12FF", "#7A0403FF")
p <- plot_ggdendro(hcdata,
direction = "tb",
scale.color = cols,
label.size = 2,
branch.size = 0.5,
expand.y = 0.2)
p=p + theme_void() + expand_limits(x = c(-1, 32))+
geom_point(data = table_all,
aes(x = match(table_all$SampleID, hcdata$labels$label),
y = -0.5,
fill = as.factor(Poligono),
shape= as.factor(Season)),
size = 4,
show.legend = FALSE) +
scale_fill_carto_d(name = "Polygons", palette = "Safe") +
scale_shape_manual(values = c(21,24))
p2 <- plot_ggdendro(hcdata2,
direction = "tb",
scale.color = cols,
label.size = 3,
branch.size = 0.5,
expand.y = 0.2)
p2=p2 + theme_void() + expand_limits(x = c(-1, 32))+
geom_point(data = table_dry,
aes(x = match(table_dry$SampleID, hcdata2$labels$label),
y = -0.5,
fill = as.factor(Poligono),
shape= as.factor(Season)),
size = 4,
show.legend = FALSE) +
scale_fill_carto_d(name = "Polygons", palette = "Safe") +
scale_shape_manual(values = c(21,24))
p3 <- plot_ggdendro(hcdata3,
direction = "tb",
scale.color = c("#800080","#2E86C1"),
label.size = 3,
branch.size = 0.5,
expand.y = 0.2)
p3=p3 + theme_void() + expand_limits(x = c(-1, 32))+
geom_point(data = table_rainy,
aes(x = match(table_rainy$SampleID, hcdata3$labels$label),
y = -0.5,
fill = as.factor(Poligono),
shape= as.factor(Season)),
size = 4,
show.legend = FALSE) +
scale_fill_carto_d(name = "Polygons", palette = "Safe") +
scale_shape_manual(values = c(24,21))z1a = z1+ guides(colour = guide_legend(nrow = 1, title = "Polygon"),
shape = guide_legend(nrow = 1) )+
theme(legend.position = "top", legend.box = "horizontal")
leg<- get_legend(z1a)
library(cowplot)
dendos=plot_grid(p2+theme(legend.position = "none"),
p3+theme(legend.position = "none"),
p+theme(legend.position = "none"),
rel_widths = c(1,1,1.5),labels = c("A", "B", "C"),
nrow = 1)
ndmss=plot_grid( z2+theme(legend.position = "none"),
z3+theme(legend.position = "none"),
z1+theme(legend.position = "none") ,
nrow = 1, rel_widths = c(1,1,1.5))
all = plot_grid(dendos,ndmss, nrow = 2, rel_heights = c(0.3,1))
all2 = plot_grid(leg, all, nrow = 2, rel_heights = c(0.1,1))
all2
Incidence-based null approach
# Incidence-based Raup-Crick beta-diversity
library(qiime2R)
library(tidyverse)
library(vegan)
library(readxl)
library(reshape2)
library(ggpubr)
library(ggh4x)
library(rcartocolor)
abund_table=read_qza("../Data/table_SMOQ_rar.qza")$data #
#abund_table=read_qza("../Data/table_SMOQ.qza")$data #
dry<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_secas_all") %>% mutate(Season="Dry")
wet<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_lluvias_all") %>% mutate(Season="Wet")
drys<- dry %>% dplyr::select( Poligono,Sitio,Transecto,pH,MO,N,P,Season)
wets<- wet %>% dplyr::select(Poligono,Sitio,Transecto,pH,MO,N,P,Season)
metadata<- rbind(drys,wets) %>% mutate(sea=c(rep("D",36), rep("R",36))) %>%
unite("SampleID", c("Poligono", "Sitio", "Transecto", "sea"),remove = F, sep = "") %>% column_to_rownames(var = "SampleID")
meta_table=metadata
meta_table$Poligono<- as.factor(meta_table$Poligono)
tp1d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="1"),]
tp2d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="2"),]
tp3d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="3"),]
tp4d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="4"),]
tp5d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="5"),]
tp6d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="6"),]
tp1r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="1"),]
tp2r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="2"),]
tp3r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="3"),]
tp4r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="4"),]
tp5r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="5"),]
tp6r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="6"),]
tp1di=row.names(tp1d)
tp2di=row.names(tp2d)
tp3di=row.names(tp3d)
tp4di=row.names(tp4d)
tp5di=row.names(tp5d)
tp6di=row.names(tp6d)
tp1ri=row.names(tp1r)
tp2ri=row.names(tp2r)
tp3ri=row.names(tp3r)
tp4ri=row.names(tp4r)
tp5ri=row.names(tp5r)
tp6ri=row.names(tp6r)
abund_table=t(abund_table)
row.names(abund_table)## [1] "121D" "122D" "122R" "123R" "211R" "212R" "213R" "221R" "222R" "223R"
## [11] "522R" "621R" "111R" "113D" "121R" "123D" "211D" "212D" "222D" "223D"
## [21] "311R" "312D" "313D" "313R" "321D" "322D" "323D" "411D" "411R" "412D"
## [31] "413D" "421D" "421R" "422D" "422R" "423D" "512D" "513D" "513R" "521D"
## [41] "521R" "613D" "621D" "622D" "622R" "623D" "623R" "112R" "423R" "511D"
## [51] "511R" "512R" "522D" "523D" "611D" "611R" "612D" "612R" "613R" "111D"
## [61] "112D" "113R" "413R" "523R" "312R" "322R" "323R" "412R" "213D" "311D"
## [71] "321R" "221D"abund_table1d=subset(abund_table, row.names(abund_table) %in% tp1di)
abund_table2d=subset(abund_table, row.names(abund_table) %in% tp2di)
abund_table3d=subset(abund_table, row.names(abund_table) %in% tp3di)
abund_table4d=subset(abund_table, row.names(abund_table) %in% tp4di)
abund_table5d=subset(abund_table, row.names(abund_table) %in% tp5di)
abund_table6d=subset(abund_table, row.names(abund_table) %in% tp6di)
abund_table1r=subset(abund_table, row.names(abund_table) %in% tp1ri)
abund_table2r=subset(abund_table, row.names(abund_table) %in% tp2ri)
abund_table3r=subset(abund_table, row.names(abund_table) %in% tp3ri)
abund_table4r=subset(abund_table, row.names(abund_table) %in% tp4ri)
abund_table5r=subset(abund_table, row.names(abund_table) %in% tp5ri)
abund_table6r=subset(abund_table, row.names(abund_table) %in% tp6ri)
abund_table1d = abund_table1d[,which(colSums(abund_table1d) != 0)]
abund_table2d = abund_table2d[,which(colSums(abund_table2d) != 0)]
abund_table3d = abund_table3d[,which(colSums(abund_table3d) != 0)]
abund_table4d = abund_table4d[,which(colSums(abund_table4d) != 0)]
abund_table5d = abund_table5d[,which(colSums(abund_table5d) != 0)]
abund_table6d = abund_table6d[,which(colSums(abund_table6d) != 0)]
abund_table1r = abund_table1r[,which(colSums(abund_table1r) != 0)]
abund_table2r = abund_table2r[,which(colSums(abund_table2r) != 0)]
abund_table3r = abund_table3r[,which(colSums(abund_table3r) != 0)]
abund_table4r = abund_table4r[,which(colSums(abund_table4r) != 0)]
abund_table5r = abund_table5r[,which(colSums(abund_table5r) != 0)]
abund_table6r = abund_table6r[,which(colSums(abund_table6r) != 0)]rc1d=QsNullModels::raup_crick(abund_table1d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc2d=QsNullModels::raup_crick(abund_table2d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc3d=QsNullModels::raup_crick(abund_table3d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc4d=QsNullModels::raup_crick(abund_table4d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc5d=QsNullModels::raup_crick(abund_table5d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc6d=QsNullModels::raup_crick(abund_table6d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc1r=QsNullModels::raup_crick(abund_table1r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc2r=QsNullModels::raup_crick(abund_table2r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc3r=QsNullModels::raup_crick(abund_table3r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc4r=QsNullModels::raup_crick(abund_table4r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc5r=QsNullModels::raup_crick(abund_table5r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc6r=QsNullModels::raup_crick(abund_table6r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rcc1d=subset(melt(as.matrix(rc1d)),value!=0)
rcc2d=subset(melt(as.matrix(rc2d)),value!=0)
rcc3d=subset(melt(as.matrix(rc3d)),value!=0)
rcc4d=subset(melt(as.matrix(rc4d)),value!=0)
rcc5d=subset(melt(as.matrix(rc5d)),value!=0)
rcc6d=subset(melt(as.matrix(rc6d)),value!=0)
rcc1r=subset(melt(as.matrix(rc1r)),value!=0)
rcc2r=subset(melt(as.matrix(rc2r)),value!=0)
rcc3r=subset(melt(as.matrix(rc3r)),value!=0)
rcc4r=subset(melt(as.matrix(rc4r)),value!=0)
rcc5r=subset(melt(as.matrix(rc5r)),value!=0)
rcc6r=subset(melt(as.matrix(rc6r)),value!=0)
rcc1d$factor="P1-dry"
rcc2d$factor="P2-dry"
rcc3d$factor="P3-dry"
rcc4d$factor="P4-dry"
rcc5d$factor="P5-dry"
rcc6d$factor="P6-dry"
rcc1r$factor="P1-wet"
rcc2r$factor="P2-wet"
rcc3r$factor="P3-wet"
rcc4r$factor="P4-wet"
rcc5r$factor="P5-wet"
rcc6r$factor="P6-wet"
rcc_dmp=rbind(rcc1r, rcc2r,rcc3r,rcc4r,rcc5r,rcc6r,
rcc1d, rcc2d, rcc3d, rcc4d, rcc5d, rcc6d)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
summary_se3=summarySE(rcc_dmp, measurevar = "value", groupvars = "factor")
# Let's save the summary_se object
#saveRDS(betadiv, "../Data/betadiv.RDS")betadiv=rcc_dmp %>% separate(factor, c("Pol", "Season"))
#colors<-viridis::turbo(6, alpha = 1, begin = 0, end = 1, direction = 1)
colors= rcartocolor::carto_pal(n = 6, name = "Safe")
strip <- strip_themed(background_x = elem_list_rect(fill = colors))
betadiv$Season = factor(betadiv$Season, levels = c("dry", "wet"), labels = c("Dry", "Rainy"))
plot<-ggerrorplot(betadiv, x="Season", y="value", desc_stat = "mean_se",
color = "Season",
facet.by = "Pol", size = 1,
error.plot = "errorbar", # Change error plot type
add = "mean")+facet_grid2(.~Pol, scales = "free", strip = strip)+
geom_hline(yintercept = 0.95, linetype=2)+
geom_hline(yintercept = -0.95, linetype=2)+theme_linedraw()+
theme(#legend.position="none",
#panel.border = element_blank(),
panel.spacing.x = unit(0.3,"line"),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"))+
stat_compare_means(aes(group = Season), label = "p.signif")+
ylab("Incidence-based beta diversity - Raup-Crick metric (±SE)")+
theme(axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank())+
scale_y_continuous(breaks = c(-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1),
sec.axis = sec_axis(~.,
breaks = c(-0.95,0.95),
name=" ‘stochastic’ range (−0.95 < βRC < + 0.95)",
labels = c("Deterministic \n more similar",
"Deterministic \n dissimilar" )))+
scale_color_manual(values = c("#A16928" ,"#2887a1"),labels=c("Dry", "Rainy"))+theme_grey()+
theme(axis.text.x = element_text(
size = 10, colour = "black" ),
#axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face = "bold"),
panel.border = element_rect(colour = "black", fill=NA, size=0.3),
axis.title = element_text(size = 14))+ theme(
strip.background.y = element_rect(
color="black", fill="black", size=1.5, linetype="solid"),
strip.background = element_rect(color="black", fill="black", size=0.7, linetype="solid"),
strip.text.y = element_text(colour = "white", face = "bold"))+
xlab("Season,Polygon")
plot
QPS
library(reshape2)
library(picante)
library(readxl)
library(tidyverse)
library(pgirmess)
library(car)
library(ggplot2)
library(ggpubr)
library(ggvegan)
library(ggordiplots)
library(qiime2R)
library(picante)
library(rcartocolor)dry<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_secas_all") %>% mutate(Season="Dry")
wet<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_lluvias_all") %>% mutate(Season="Wet")
# select characteristics in both datasets
drys<- dry %>% dplyr::select( Poligono,Sitio,Transecto,pH,MO,N,P,Season)
wets<- wet %>% dplyr::select(Poligono,Sitio,Transecto,pH,MO,N,P,Season)
env_season<- rbind(drys,wets) %>% mutate(sea=c(rep("D",36), rep("R",36))) %>% unite("SampleID", c("Poligono", "Sitio", "Transecto", "sea"),remove = F, sep = "")
env_season$Poligono<- factor(env_season$Poligono, levels = c(1,2,3,4,5,6),
labels = c("P1", "P2", "P3", "P4", "P5", "P6"))
env_season$Season<- factor(env_season$Season)
env_season<-env_season %>% unite("interact", c("Poligono", "Season"), remove = F)
df <-env_season %>% column_to_rownames(var="SampleID") %>% dplyr::select(pH,MO,N,P)
df_ph = df %>% dplyr::select(pH)
dist_ph = dist(df_ph)
df_mo = df %>% dplyr::select(MO)
dist_mo = dist(df_mo)
df_n = df %>% dplyr::select(N)
dist_n = dist(df_n)
df_p = df %>% dplyr::select(P)
dist_p = dist(df_p)
abund_table=read_qza("../Data/table_SMOQ_rar.qza")$data #tree = read_qza("../Data/tree_SMOQ_rar/rooted_tree.qza")$data
phydist <- cophenetic(tree)
comdis = comdist(t(abund_table), phydist, abundance.weighted = TRUE)
data(phylocom)
comdist(phylocom$sample, cophenetic(phylocom$phylo), abundance.weighted=TRUE)
# Qantitative Process Estimate (QPE) from Stegen et al. 2013
# This approach requires that phylogenetic distances (PD) among taxa reflect differences in
# the ecological niches they inhabit, thus, carry a phylogenetic signal.
# The presence of phylogenetic signals was tested using Mantel correlograms, as described in Stegen et al. 2013
# The necessary R script was provided by Jianjun Wang and used in the study of Langenheder et al. 2017 FEMS Microbiology Ecology
# After performing phylogenetic signal tests, let's use mantel.correlog function and plot the results
# The following script was run separately on 16S and 18S dataset!
library(readr)
library(vegan)
library(picante)
library(ape)
library(dplyr)
phylo.sig.correlog = mantel.correlog(comdis,dist_ph,
nperm=999,
cutoff=FALSE,
n.class=50,
mult="bonferroni");
phylo.sig.correlog[["mantel.res"]]
write.csv(phylo.sig.correlog$mantel.res, "ph_signal.csv", quote=F);
phylo.sig.correlog = mantel.correlog(comdis, dist_mo,
nperm=999,
cutoff=FALSE,
n.class=50,
mult="bonferroni");
phylo.sig.correlog[["mantel.res"]]
write.csv(phylo.sig.correlog$mantel.res, "mo_signal.csv", quote=F);
phylo.sig.correlog = mantel.correlog(comdis, dist_p,
nperm=999,
cutoff=FALSE,
n.class=50,
mult="bonferroni");
phylo.sig.correlog[["mantel.res"]]
write.csv(phylo.sig.correlog$mantel.res, "p_signal.csv", quote=F);
phylo.sig.correlog = mantel.correlog(comdis, dist_n,
nperm=999,
cutoff=FALSE,
n.class=50,
mult="bonferroni");
phylo.sig.correlog[["mantel.res"]]
write.csv(phylo.sig.correlog$mantel.res, "n_signal.csv", quote=F);# Plotting the Mantel correlograms
# Figure S5 and S6
# 16S dataset (Figure S5)
#pdf("phylo_signals.pdf")
#png(#"phylo_signals.png", width = 1600, height = 1200, res = 200)
par(mfrow=c(2,2))
phylo.sig<-read.csv("../Data/ph_signal.csv", row.names = 1, header=T, sep=",")
plot(phylo.sig[,c(1,3)],
xlab="Phylogenetic Distance Class", ylab="Mantel Test Statistic") +lines(phylo.sig[,c(1,3)]) + points(phylo.sig[,c(1,3)], pch=21, bg="white", col="black", pty=4)+
points(phylo.sig[,c(1,3)][phylo.sig[,4] < 0.05,], pch=21, bg="black", col="black", pty=4) + abline(h=0, lty=2, col="red") +
title("pH")## integer(0)phylo.sig<-read.csv("../Data/mo_signal.csv", row.names = 1, header=T, sep=",")
plot(phylo.sig[,c(1,3)],
xlab="Phylogenetic Distance Class", ylab="Mantel Test Statistic") +lines(phylo.sig[,c(1,3)]) + points(phylo.sig[,c(1,3)], pch=21, bg="white", col="black", pty=4)+
points(phylo.sig[,c(1,3)][phylo.sig[,4] < 0.05,], pch=21, bg="black", col="black", pty=4) + abline(h=0, lty=2, col="red") +
title("OM")## integer(0)phylo.sig<-read.csv("../Data/n_signal.csv", row.names = 1, header=T, sep=",")
phylo.sig = drop_na(phylo.sig)
plot(phylo.sig[,c(1,3)],
xlab="Phylogenetic Distance Class", ylab="Mantel Test Statistic") +lines(phylo.sig[,c(1,3)]) + points(phylo.sig[,c(1,3)], pch=21, bg="white", col="black", pty=4)+
points(phylo.sig[,c(1,3)][phylo.sig[,4] < 0.05,], pch=21, bg="black", col="black", pty=4) + abline(h=0, lty=2, col="red") +
title("Total nitrogen")## integer(0)phylo.sig<-read.csv("../Data/p_signal.csv", row.names = 1, header=T, sep=",")
plot(phylo.sig[,c(1,3)],
xlab="Phylogenetic Distance Class", ylab="Mantel Test Statistic") +lines(phylo.sig[,c(1,3)]) + points(phylo.sig[,c(1,3)], pch=21, bg="white", col="black", pty=4)+
points(phylo.sig[,c(1,3)][phylo.sig[,4] < 0.05,], pch=21, bg="black", col="black", pty=4) + abline(h=0, lty=2, col="red") +
title("Total phosphorus")
## integer(0)# After I evaluated the existence of phylogenetic signal, QPE analyses can be performed
# First step of the QPE approach (it can be time-consuming process, so I run it on UPPMAX server)
dry<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_secas_all") %>% mutate(Season="Dry")
wet<- read_excel("../Data/Fisicoquimicos.xlsx", sheet = "metadata_lluvias_all") %>% mutate(Season="Wet")
drys<- dry %>% dplyr::select( Poligono,Sitio,Transecto,pH,MO,N,P,Season)
wets<- wet %>% dplyr::select(Poligono,Sitio,Transecto,pH,MO,N,P,Season)
metadata<- rbind(drys,wets) %>% mutate(sea=c(rep("D",36), rep("R",36))) %>%
unite("SampleID", c("Poligono", "Sitio", "Transecto", "sea"),remove = F, sep = "") %>% column_to_rownames(var = "SampleID")
meta_table=metadata
meta_table$Poligono<- as.factor(meta_table$Poligono)
tp1d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="1"),]
tp2d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="2"),]
tp3d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="3"),]
tp4d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="4"),]
tp5d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="5"),]
tp6d<- meta_table[which(meta_table$Season=="Dry" & meta_table$Poligono=="6"),]
tp1r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="1"),]
tp2r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="2"),]
tp3r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="3"),]
tp4r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="4"),]
tp5r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="5"),]
tp6r<- meta_table[which(meta_table$Season=="Wet" & meta_table$Poligono=="6"),]
tp1di=row.names(tp1d)
tp2di=row.names(tp2d)
tp3di=row.names(tp3d)
tp4di=row.names(tp4d)
tp5di=row.names(tp5d)
tp6di=row.names(tp6d)
tp1ri=row.names(tp1r)
tp2ri=row.names(tp2r)
tp3ri=row.names(tp3r)
tp4ri=row.names(tp4r)
tp5ri=row.names(tp5r)
tp6ri=row.names(tp6r)
library(picante)
library(qiime2R)
abund_table=read_qza("../Data/table_SMOQ_rar.qza")$data #
name= colnames(abund_table)
tree = read_qza("../Data/tree_SMOQ_rar/rooted_tree.qza")$data
m=match.phylo.data(tree, data.frame(abund_table, check.names = F, check.rows = F))
tree=m$phy
abund_table=m$data
abund_table=t(abund_table)
row.names(abund_table)
table1d=subset(abund_table, row.names(abund_table) %in% tp1di)
table2d=subset(abund_table, row.names(abund_table) %in% tp2di)
table3d=subset(abund_table, row.names(abund_table) %in% tp3di)
table4d=subset(abund_table, row.names(abund_table) %in% tp4di)
table5d=subset(abund_table, row.names(abund_table) %in% tp5di)
table6d=subset(abund_table, row.names(abund_table) %in% tp6di)
table1r=subset(abund_table, row.names(abund_table) %in% tp1ri)
table2r=subset(abund_table, row.names(abund_table) %in% tp2ri)
table3r=subset(abund_table, row.names(abund_table) %in% tp3ri)
table4r=subset(abund_table, row.names(abund_table) %in% tp4ri)
table5r=subset(abund_table, row.names(abund_table) %in% tp5ri)
table6r=subset(abund_table, row.names(abund_table) %in% tp6ri)
m=match.phylo.data(tree, data.frame(t(table1d), check.names = F, check.rows = F))
abund_table1d=data.frame(m$data, check.names = F)
abund_table1d = abund_table1d[which(rowSums(abund_table1d) != 0),]
m=match.phylo.data(tree, abund_table1d)
OTU_tree1d=m$phy
m=match.phylo.data(tree, data.frame(t(table2d), check.names = F, check.rows = F))
abund_table2d=data.frame(m$data, check.names = F)
abund_table2d = abund_table2d[which(rowSums(abund_table2d) != 0),]
m=match.phylo.data(tree, abund_table2d)
OTU_tree2d=m$phy
m=match.phylo.data(tree, data.frame(t(table3d), check.names = F, check.rows = F))
abund_table3d=data.frame(m$data, check.names = F)
abund_table3d = abund_table3d[which(rowSums(abund_table3d) != 0),]
m=match.phylo.data(tree, abund_table3d)
OTU_tree3d=m$phy
m=match.phylo.data(tree, data.frame(t(table4d), check.names = F, check.rows = F))
abund_table4d=data.frame(m$data, check.names = F)
abund_table4d = abund_table4d[which(rowSums(abund_table4d) != 0),]
m=match.phylo.data(tree, abund_table4d)
OTU_tree4d=m$phy
m=match.phylo.data(tree, data.frame(t(table5d), check.names = F, check.rows = F))
abund_table5d=data.frame(m$data, check.names = F)
abund_table5d = abund_table5d[which(rowSums(abund_table5d) != 0),]
m=match.phylo.data(tree, abund_table5d)
OTU_tree5d=m$phy
m=match.phylo.data(tree, data.frame(t(table6d), check.names = F, check.rows = F))
abund_table6d=data.frame(m$data, check.names = F)
abund_table6d = abund_table6d[which(rowSums(abund_table6d) != 0),]
m=match.phylo.data(tree, abund_table6d)
OTU_tree6d=m$phy
m=match.phylo.data(tree, data.frame(t(table1r), check.names = F, check.rows = F))
abund_table1r=data.frame(m$data, check.names = F)
abund_table1r = abund_table1r[which(rowSums(abund_table1r) != 0),]
m=match.phylo.data(tree, abund_table1r)
OTU_tree1r=m$phy
m=match.phylo.data(tree, data.frame(t(table2r), check.names = F, check.rows = F))
abund_table2r=data.frame(m$data, check.names = F)
abund_table2r = abund_table2r[which(rowSums(abund_table2d) != 0),]
m=match.phylo.data(tree, abund_table2r)
OTU_tree2r=m$phy
m=match.phylo.data(tree, data.frame(t(table3r), check.names = F, check.rows = F))
abund_table3r=data.frame(m$data, check.names = F)
abund_table3r = abund_table3r[which(rowSums(abund_table3r) != 0),]
m=match.phylo.data(tree, abund_table3r)
OTU_tree3r=m$phy
m=match.phylo.data(tree, data.frame(t(table4r), check.names = F, check.rows = F))
abund_table4r=data.frame(m$data, check.names = F)
abund_table4r = abund_table4r[which(rowSums(abund_table4r) != 0),]
m=match.phylo.data(tree, abund_table4r)
OTU_tree4r=m$phy
m=match.phylo.data(tree, data.frame(t(table5r), check.names = F, check.rows = F))
abund_table5r=data.frame(m$data, check.names = F)
abund_table5r = abund_table5r[which(rowSums(abund_table5r) != 0),]
m=match.phylo.data(tree, abund_table5r)
OTU_tree5r=m$phy
m=match.phylo.data(tree, data.frame(t(table6r), check.names = F, check.rows = F))
abund_table6r=data.frame(m$data, check.names = F)
abund_table6r = abund_table6r[which(rowSums(abund_table6r) != 0),]
m=match.phylo.data(tree, abund_table6r)
OTU_tree6r=m$phy
#1d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table1d),cophenetic(OTU_tree1d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table1d),colnames(beta.mntd.weighted));
identical(colnames(abund_table1d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table1d),ncol(abund_table1d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table1d),taxaShuffle(cophenetic(OTU_tree1d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table1d),ncol=ncol(abund_table1d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table1d)-1)) {
for (rows in (columns+1):ncol(abund_table1d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table1d);
colnames(weighted.bNTI) = colnames(abund_table1d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI1d.csv",quote=F)
#2d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table2d),cophenetic(OTU_tree2d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table2d),colnames(beta.mntd.weighted));
identical(colnames(abund_table2d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table2d),ncol(abund_table2d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table2d),taxaShuffle(cophenetic(OTU_tree2d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table2d),ncol=ncol(abund_table2d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table2d)-1)) {
for (rows in (columns+1):ncol(abund_table2d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table2d);
colnames(weighted.bNTI) = colnames(abund_table2d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI2d.csv",quote=F)
#3d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table3d),cophenetic(OTU_tree3d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table3d),colnames(beta.mntd.weighted));
identical(colnames(abund_table3d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table3d),ncol(abund_table3d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table3d),taxaShuffle(cophenetic(OTU_tree3d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table3d),ncol=ncol(abund_table3d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table3d)-1)) {
for (rows in (columns+1):ncol(abund_table3d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table3d);
colnames(weighted.bNTI) = colnames(abund_table3d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI3d.csv",quote=F)
#4d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table4d),cophenetic(OTU_tree4d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table4d),colnames(beta.mntd.weighted));
identical(colnames(abund_table4d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table4d),ncol(abund_table4d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table4d),taxaShuffle(cophenetic(OTU_tree4d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table4d),ncol=ncol(abund_table4d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table4d)-1)) {
for (rows in (columns+1):ncol(abund_table4d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table4d);
colnames(weighted.bNTI) = colnames(abund_table4d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI4d.csv",quote=F)
#5d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table5d),cophenetic(OTU_tree5d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table5d),colnames(beta.mntd.weighted));
identical(colnames(abund_table5d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table5d),ncol(abund_table5d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table5d),taxaShuffle(cophenetic(OTU_tree5d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table5d),ncol=ncol(abund_table5d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table5d)-1)) {
for (rows in (columns+1):ncol(abund_table5d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table5d);
colnames(weighted.bNTI) = colnames(abund_table5d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI5d.csv",quote=F)
#6d
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table6d),cophenetic(OTU_tree6d),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table6d),colnames(beta.mntd.weighted));
identical(colnames(abund_table6d),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table6d),ncol(abund_table6d),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table6d),taxaShuffle(cophenetic(OTU_tree6d)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table6d),ncol=ncol(abund_table6d));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table6d)-1)) {
for (rows in (columns+1):ncol(abund_table6d)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table6d);
colnames(weighted.bNTI) = colnames(abund_table6d);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI6d.csv",quote=F)
#1r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table1r),cophenetic(OTU_tree1r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table1r),colnames(beta.mntd.weighted));
identical(colnames(abund_table1r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table1r),ncol(abund_table1r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table1r),taxaShuffle(cophenetic(OTU_tree1r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table1r),ncol=ncol(abund_table1r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table1r)-1)) {
for (rows in (columns+1):ncol(abund_table1r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table1r);
colnames(weighted.bNTI) = colnames(abund_table1r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI1r.csv",quote=F)
#2r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table2r),cophenetic(OTU_tree2r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table2r),colnames(beta.mntd.weighted));
identical(colnames(abund_table2r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table2r),ncol(abund_table2r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table2r),taxaShuffle(cophenetic(OTU_tree2r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table2r),ncol=ncol(abund_table2r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table2r)-1)) {
for (rows in (columns+1):ncol(abund_table2r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table2r);
colnames(weighted.bNTI) = colnames(abund_table2r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI2r.csv",quote=F)
#3r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table3r),cophenetic(OTU_tree3r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table3r),colnames(beta.mntd.weighted));
identical(colnames(abund_table3r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table3r),ncol(abund_table3r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table3r),taxaShuffle(cophenetic(OTU_tree3r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table3r),ncol=ncol(abund_table3r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table3r)-1)) {
for (rows in (columns+1):ncol(abund_table3r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table3r);
colnames(weighted.bNTI) = colnames(abund_table3r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI3r.csv",quote=F)
#4r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table4r),cophenetic(OTU_tree4r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table4r),colnames(beta.mntd.weighted));
identical(colnames(abund_table4r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table4r),ncol(abund_table4r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table4r),taxaShuffle(cophenetic(OTU_tree4r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table4r),ncol=ncol(abund_table4r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table4r)-1)) {
for (rows in (columns+1):ncol(abund_table4r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table4r);
colnames(weighted.bNTI) = colnames(abund_table4r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI4r.csv",quote=F)
#5r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table5r),cophenetic(OTU_tree5r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table5r),colnames(beta.mntd.weighted));
identical(colnames(abund_table5r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table5r),ncol(abund_table5r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table5r),taxaShuffle(cophenetic(OTU_tree5r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table5r),ncol=ncol(abund_table5r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table5r)-1)) {
for (rows in (columns+1):ncol(abund_table5r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table5r);
colnames(weighted.bNTI) = colnames(abund_table5r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI5r.csv",quote=F)
#6r
beta.mntd.weighted = as.matrix(comdistnt(t(abund_table6r),cophenetic(OTU_tree6r),abundance.weighted=T));
dim(beta.mntd.weighted);
identical(colnames(abund_table6r),colnames(beta.mntd.weighted));
identical(colnames(abund_table6r),rownames(beta.mntd.weighted));
beta.reps = 999;
rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(abund_table6r),ncol(abund_table6r),beta.reps));
dim(rand.weighted.bMNTD.comp);
for (rep in 1:beta.reps) {
rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(abund_table6r),taxaShuffle(cophenetic(OTU_tree6r)),abundance.weighted=T,exclude.conspecifics = F));
print(c(date(),rep));
}
weighted.bNTI = matrix(c(NA),nrow=ncol(abund_table6r),ncol=ncol(abund_table6r));
dim(weighted.bNTI);
for (columns in 1:(ncol(abund_table6r)-1)) {
for (rows in (columns+1):ncol(abund_table6r)) {
rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
rm("rand.vals");
};
};
rownames(weighted.bNTI) = colnames(abund_table6r);
colnames(weighted.bNTI) = colnames(abund_table6r);
weighted.bNTI;
write.csv(weighted.bNTI, "./weighted_bNTI6r.csv",quote=F)# For 16S dataset
comm=read_qza("../Data/table_SMOQ_rar.qza")$data #
Tree = read_qza("../Data/tree_SMOQ_rar/rooted_tree.qza")$data
mc=match.phylo.data(Tree, data.frame(comm, check.names = F, check.rows = F))
comm=mc$data
comm=data.frame(t(comm))
abund_table1d=subset(abund_table, row.names(abund_table) %in% tp1di)
abund_table2d=subset(abund_table, row.names(abund_table) %in% tp2di)
abund_table3d=subset(abund_table, row.names(abund_table) %in% tp3di)
abund_table4d=subset(abund_table, row.names(abund_table) %in% tp4di)
abund_table5d=subset(abund_table, row.names(abund_table) %in% tp5di)
abund_table6d=subset(abund_table, row.names(abund_table) %in% tp6di)
abund_table1r=subset(abund_table, row.names(abund_table) %in% tp1ri)
abund_table2r=subset(abund_table, row.names(abund_table) %in% tp2ri)
abund_table3r=subset(abund_table, row.names(abund_table) %in% tp3ri)
abund_table4r=subset(abund_table, row.names(abund_table) %in% tp4ri)
abund_table5r=subset(abund_table, row.names(abund_table) %in% tp5ri)
abund_table6r=subset(abund_table, row.names(abund_table) %in% tp6ri)
abund_table1d = abund_table1d[,which(colSums(abund_table1d) != 0)]
abund_table2d = abund_table2d[,which(colSums(abund_table2d) != 0)]
abund_table3d = abund_table3d[,which(colSums(abund_table3d) != 0)]
abund_table4d = abund_table4d[,which(colSums(abund_table4d) != 0)]
abund_table5d = abund_table5d[,which(colSums(abund_table5d) != 0)]
abund_table6d = abund_table6d[,which(colSums(abund_table6d) != 0)]
abund_table1r = abund_table1r[,which(colSums(abund_table1r) != 0)]
abund_table2r = abund_table2r[,which(colSums(abund_table2r) != 0)]
abund_table3r = abund_table3r[,which(colSums(abund_table3r) != 0)]
abund_table4r = abund_table4r[,which(colSums(abund_table4r) != 0)]
abund_table5r = abund_table5r[,which(colSums(abund_table5r) != 0)]
abund_table6r = abund_table6r[,which(colSums(abund_table6r) != 0)]
source("../Scripts/raup_crick_abundance.R")
library(phyloseq)
library(picante)
library(ape)
rc1da=raup_crick_abundance2(abund_table1d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc2da=raup_crick_abundance2(abund_table2d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc3da=raup_crick_abundance2(abund_table3d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc4da=raup_crick_abundance2(abund_table4d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc5da=raup_crick_abundance2(abund_table5d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc6da=raup_crick_abundance2(abund_table6d, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc1ra=raup_crick_abundance2(abund_table1r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc2ra=raup_crick_abundance2(abund_table2r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc3ra=raup_crick_abundance2(abund_table3r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc4ra=raup_crick_abundance2(abund_table4r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc5ra=raup_crick_abundance2(abund_table5r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
rc6ra=raup_crick_abundance2(abund_table6r, plot_names_in_col1 = F, reps = 999, as.distance.matrix = T, set_all_species_equal = F)
write.csv(as.matrix(rc1da), "./rc1da_abundance.csv")
write.csv(as.matrix(rc2da), "./rc2da_abundance.csv")
write.csv(as.matrix(rc3da), "./rc3da_abundance.csv")
write.csv(as.matrix(rc4da), "./rc4da_abundance.csv")
write.csv(as.matrix(rc5da), "./rc5da_abundance.csv")
write.csv(as.matrix(rc6da), "./rc6da_abundance.csv")
write.csv(as.matrix(rc1ra), "./rc1ra_abundance.csv")
write.csv(as.matrix(rc2ra), "./rc2ra_abundance.csv")
write.csv(as.matrix(rc3ra), "./rc3ra_abundance.csv")
write.csv(as.matrix(rc4ra), "./rc4ra_abundance.csv")
write.csv(as.matrix(rc5ra), "./rc5ra_abundance.csv")
write.csv(as.matrix(rc6ra), "./rc6ra_abundance.csv")library(reshape2)
#1d
select_1d=read.csv("../Data/weighted_bNTI1d.csv", check.names = F, row.names = 1)
head(select_1d)
sel_1d = melt(as.matrix(select_1d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_1d) <- c("Row", "Column", "Value")
sel_1d = sel_1d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#2d
select_2d=read.csv("../Data/weighted_bNTI2d.csv", check.names = F, row.names = 1)
head(select_2d)
sel_2d = melt(as.matrix(select_2d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_2d) <- c("Row", "Column", "Value")
sel_2d = sel_2d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#3d
select_3d=read.csv("../Data/weighted_bNTI3d.csv", check.names = F, row.names = 1)
head(select_3d)
sel_3d = melt(as.matrix(select_3d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_3d) <- c("Row", "Column", "Value")
sel_3d = sel_3d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#4d
select_4d=read.csv("../Data/weighted_bNTI4d.csv", check.names = F, row.names = 1)
head(select_4d)
sel_4d = melt(as.matrix(select_4d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_4d) <- c("Row", "Column", "Value")
sel_4d = sel_4d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#5d
select_5d=read.csv("../Data/weighted_bNTI5d.csv", check.names = F, row.names = 1)
head(select_5d)
sel_5d = melt(as.matrix(select_5d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_5d) <- c("Row", "Column", "Value")
sel_5d = sel_5d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#6d
select_6d=read.csv("../Data/weighted_bNTI6d.csv", check.names = F, row.names = 1)
head(select_6d)
sel_6d = melt(as.matrix(select_6d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_6d) <- c("Row", "Column", "Value")
sel_6d = sel_6d %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#1r
select_1r=read.csv("../Data/weighted_bNTI1r.csv", check.names = F, row.names = 1)
head(select_1r)
sel_1r = melt(as.matrix(select_1r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_1r) <- c("Row", "Column", "Value")
sel_1r = sel_1r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#2r
select_2r=read.csv("../Data/weighted_bNTI2r.csv", check.names = F, row.names = 1)
head(select_2r)
sel_2r = melt(as.matrix(select_2r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_2r) <- c("Row", "Column", "Value")
sel_2r = sel_2r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#3r
select_3r=read.csv("../Data/weighted_bNTI3r.csv", check.names = F, row.names = 1)
head(select_3r)
sel_3r = melt(as.matrix(select_3r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_3r) <- c("Row", "Column", "Value")
sel_3r = sel_3r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#4r
select_4r=read.csv("../Data/weighted_bNTI4r.csv", check.names = F, row.names = 1)
head(select_4r)
sel_4r = melt(as.matrix(select_4r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_4r) <- c("Row", "Column", "Value")
sel_4r = sel_4r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#5r
select_5r=read.csv("../Data/weighted_bNTI5r.csv", check.names = F, row.names = 1)
head(select_5r)
sel_5r = melt(as.matrix(select_5r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_5r) <- c("Row", "Column", "Value")
sel_5r = sel_5r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
#6r
select_6r=read.csv("../Data/weighted_bNTI6r.csv", check.names = F, row.names = 1)
head(select_6r)
sel_6r = melt(as.matrix(select_6r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(sel_6r) <- c("Row", "Column", "Value")
sel_6r = sel_6r %>% mutate(
selection = ifelse(
abs(Value) > 2,
ifelse(Value > 2, "variable", "homogeneous"),
"no selection"
),
percentage_variable = ifelse(
abs(Value) > 2 & selection == "variable",
100,
0
),
percentage_homog = ifelse(
abs(Value) > 2 & selection == "homogeneous",
100,
0
)
)
sel1d = sel_1d %>% mutate(pol="1_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel2d = sel_2d %>% mutate(pol="2_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel3d = sel_3d %>% mutate(pol="3_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel4d = sel_4d %>% mutate(pol="4_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel5d = sel_5d %>% mutate(pol="5_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel6d = sel_6d %>% mutate(pol="6_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel1r = sel_1r %>% mutate(pol="1_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel2r = sel_2r %>% mutate(pol="2_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel3r = sel_3r %>% mutate(pol="3_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel4r = sel_4r %>% mutate(pol="4_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel5r = sel_5r %>% mutate(pol="5_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
sel6r = sel_6r %>% mutate(pol="6_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
selection = rbind(sel1d, sel2d, sel3d, sel4d, sel5d, sel6d,
sel1r, sel2r, sel3r, sel4r, sel5r, sel6r)
selection = selection %>% mutate(selection=percentage_variable+percentage_homog)rc1da = read.csv("../Data/rc1da_abundance.csv", check.names = F)
#1d
nosel_1d= as.matrix(rc1da)
lower_1d <- lower.tri(nosel_1d)
nosel_1d[!lower_1d] <- NA
nosel_1d = melt(as.matrix(nosel_1d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_1d) <- c("Row", "Column", "Value")
nosel_1d = nosel_1d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#2d
rc2da = read.csv("../Data/rc2da_abundance.csv", check.names = F)
nosel_2d= as.matrix(rc2da)
lower_2d <- lower.tri(nosel_2d)
nosel_2d[!lower_2d] <- NA
nosel_2d = melt(as.matrix(nosel_2d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_2d) <- c("Row", "Column", "Value")
nosel_2d = nosel_2d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#3d
rc3da = read.csv("../Data/rc3da_abundance.csv", check.names = F)
nosel_3d= as.matrix(rc3da)
lower_3d <- lower.tri(nosel_3d)
nosel_3d[!lower_3d] <- NA
nosel_3d = melt(as.matrix(nosel_3d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_3d) <- c("Row", "Column", "Value")
nosel_3d = nosel_3d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#4d
rc4da = read.csv("../Data/rc4da_abundance.csv", check.names = F)
nosel_4d= as.matrix(rc4da)
lower_4d <- lower.tri(nosel_4d)
nosel_4d[!lower_4d] <- NA
nosel_4d = melt(as.matrix(nosel_4d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_4d) <- c("Row", "Column", "Value")
nosel_4d = nosel_4d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#5d
rc5da = read.csv("../Data/rc5da_abundance.csv", check.names = F)
nosel_5d= as.matrix(rc5da)
lower_5d <- lower.tri(nosel_5d)
nosel_5d[!lower_5d] <- NA
nosel_5d = melt(as.matrix(nosel_5d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_5d) <- c("Row", "Column", "Value")
nosel_5d = nosel_5d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#6d
rc6da = read.csv("../Data/rc6da_abundance.csv", check.names = F)
nosel_6d= as.matrix(rc6da)
lower_6d <- lower.tri(nosel_6d)
nosel_6d[!lower_6d] <- NA
nosel_6d = melt(as.matrix(nosel_6d), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_6d) <- c("Row", "Column", "Value")
nosel_6d = nosel_6d %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#1r
rc1ra = read.csv("../Data/rc1ra_abundance.csv", check.names = F)
nosel_1r= as.matrix(rc1ra)
lower_1r <- lower.tri(nosel_1r)
nosel_1r[!lower_1r] <- NA
nosel_1r = melt(as.matrix(nosel_1r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_1r) <- c("Row", "Column", "Value")
nosel_1r = nosel_1r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#2r
rc2ra = read.csv("../Data/rc1ra_abundance.csv", check.names = F)
nosel_2r= as.matrix(rc2ra)
lower_2r <- lower.tri(nosel_2r)
nosel_2r[!lower_2r] <- NA
nosel_2r = melt(as.matrix(nosel_2r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_2r) <- c("Row", "Column", "Value")
nosel_2r = nosel_2r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#3r
rc3ra = read.csv("../Data/rc3ra_abundance.csv", check.names = F)
nosel_3r= as.matrix(rc3ra)
lower_3r <- lower.tri(nosel_3r)
nosel_3r[!lower_3r] <- NA
nosel_3r = melt(as.matrix(nosel_3r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_3r) <- c("Row", "Column", "Value")
nosel_3r = nosel_3r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#4r
rc4ra = read.csv("../Data/rc4ra_abundance.csv", check.names = F)
nosel_4r= as.matrix(rc4ra)
lower_4r <- lower.tri(nosel_4r)
nosel_4r[!lower_4r] <- NA
nosel_4r = melt(as.matrix(nosel_4r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_4r) <- c("Row", "Column", "Value")
nosel_4r = nosel_4r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#5r
rc5ra = read.csv("../Data/rc5ra_abundance.csv", check.names = F)
nosel_5r= as.matrix(rc5ra)
lower_5r <- lower.tri(nosel_5r)
nosel_5r[!lower_5r] <- NA
nosel_5r = melt(as.matrix(nosel_5r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_5r) <- c("Row", "Column", "Value")
nosel_5r = nosel_5r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
#6r
rc6ra = read.csv("../Data/rc6ra_abundance.csv", check.names = F)
nosel_6r= as.matrix(rc6ra)
lower_6r <- lower.tri(nosel_6r)
nosel_6r[!lower_6r] <- NA
nosel_6r = melt(as.matrix(nosel_6r), varnames = c("Row", "Column"), value.name = "Value", na.rm = TRUE)
colnames(nosel_6r) <- c("Row", "Column", "Value")
nosel_6r = nosel_6r %>%
mutate(
# Categorizar los valores
category = case_when(
Value > 0.95 ~ "histcont_displimit",
Value <= 0.95 & Value >= -0.95 ~ "eco_drift",
Value < -0.95 ~ "homog_disp"
),
# Agregar columnas de porcentaje basadas en la nueva categoría
percentage_histocont_displimit = ifelse(category == "histcont_displimit", 100, 0),
percentage_eco_drift = ifelse(category == "eco_drift", 100, 0),
percentage_homog_disp = ifelse(category == "homog_disp", 100, 0)
)
nosel1d = nosel_1d %>% mutate(pol="1_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel2d = nosel_2d %>% mutate(pol="2_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel3d = nosel_3d %>% mutate(pol="3_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel4d = nosel_4d %>% mutate(pol="4_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel5d = nosel_5d %>% mutate(pol="5_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel6d = nosel_6d %>% mutate(pol="6_d") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel1r = nosel_1r %>% mutate(pol="1_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel2r = nosel_2r %>% mutate(pol="2_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel3r = nosel_3r %>% mutate(pol="3_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel4r = nosel_4r %>% mutate(pol="4_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel5r = nosel_5r %>% mutate(pol="5_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
nosel6r = nosel_6r %>% mutate(pol="6_r") %>% group_by(pol) %>% summarise_if(is.numeric, mean)
noselection = rbind(nosel1d, nosel2d, nosel3d, nosel4d, nosel5d, nosel6d,
nosel1r, nosel2r, nosel3r, nosel4r, nosel5r, nosel6r)
noselection = noselection %>% mutate(noselection= percentage_histocont_displimit+percentage_eco_drift+percentage_homog_disp)all = selection %>% full_join(noselection, by = "pol")
all= all %>% mutate(noselection_ok = noselection-selection) %>%
mutate(percentage_histocont_displimit_ok = (percentage_histocont_displimit/100)*noselection_ok,
percentage_eco_drift_ok = (percentage_eco_drift/100)*noselection_ok,
percentage_homog_disp_ok = (percentage_homog_disp/100)*noselection_ok) %>%
mutate(total=percentage_variable+percentage_homog+percentage_eco_drift_ok+percentage_homog_disp_ok+percentage_histocont_displimit_ok )
data_plot = all %>% dplyr::select(
pol,
Selection_variable=percentage_variable,
Selection_Homogeneous=percentage_homog,
Ecological_drift=percentage_eco_drift_ok,
Homogenizing_dispersal=percentage_homog_disp_ok,
"Historical_cont/Dispersal_limit"=percentage_histocont_displimit_ok ) %>% pivot_longer(cols = -pol, names_to = "Process", values_to = "Proportion") %>% separate("pol", into = c("Polygon", "Season"), remove = F)
data_plot$Season = factor(data_plot$Season, levels = c("d", "r"), labels = c("Dry", "Rainy"))
data_plot$Polygon = factor(data_plot$Polygon, levels=c("1","2","3","4","5","6"), labels=c("P1", "P2", "P3", "P4", "P5", "P6"))
library(ggh4x)
library(ggpubr)
bseason=ggbarplot(data = data_plot, x = "Season", y = "Proportion", fill = "Process", add = "mean",facet.by = "Season", size = 1)+facet_grid2(.~Season, scales = "free")+
scale_fill_manual(values = c("#008000", "#2ecc71", "#808000", "#800000", "#FF0000"))+ theme_grey()+
theme(legend.position="right")+
theme(axis.text = element_text(size=12),
axis.title = element_text(size=14),
strip.text = element_text(size=12),
axis.text.x = element_blank(),
plot.title = element_text(size=8, face="bold") )+ theme(
strip.background = element_rect(
color="black", fill="black", size=1.5, linetype="solid"),
strip.text = element_text(colour = "white", face = "bold"),
legend.position = "none")+ylab("Proportion (%)")
colors<-viridis::turbo(6, alpha = 1, begin = 0, end = 1, direction = 1)
combos = data.frame(
Season=c(rep("Dry", 6), rep("Rainy", "6")),
Polygon=rep(paste0("P", 1:6),2)
)
combos = data.frame(
Polygon=c("P1", "P1", "P2", "P2", "P3", "P3", "P4", "P4", "P5", "P5", "P6", "P6"),
Season=c("Dry", "Rainy","Dry", "Rainy","Dry", "Rainy","Dry", "Rainy","Dry", "Rainy","Dry", "Rainy")
)
strip_background <- strip_nested(
text_x = elem_list_text(colour = "white", face = "bold"),
background_x =
elem_list_rect(
fill = c(
# level2 colors
case_match(
unique(combos$Polygon),
"P1" ~ "#88CCEE" ,
"P2" ~ "#CC6677",
"P3" ~ "#DDCC77",
"P4" ~ "#117733",
"P5"~ "#332288" ,
"P6" ~ "#888888",
.default = "blue"),
# level1 colors
case_match(
combos$Season,
"Dry" ~ "#A16928",
"Rainy" ~ "#2887a1",
.default = "grey"
)
)
)
)
p=ggbarplot(data = data_plot, x = "Season", y = "Proportion", fill = "Process", add = "mean", facet.by = "Polygon")+
#scale_fill_manual(values = c("#008000", "#2ecc71", "#808000","#800000","#FF0000"))+
scale_fill_carto_d(name = "Process", palette = "Geyser") +
theme_gray()+
facet_nested(.~Polygon, scales = "free", strip = strip_background)+
ylab("Proportion (%)")+
theme_grey()+
theme(#axis.text.x = element_blank(),
#axis.ticks.x = element_blank(),
strip.text.x = element_text(size = 14, face = "bold", color="white"),
strip.text.y = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 14))+ theme(
strip.background = element_rect(color="black", fill="black", size=0.7,linetype="solid"),
panel.border = element_rect(colour = "black", fill=NA, size=0.3),
axis.text.x = element_text(size = 10, colour = "black" ),
strip.text.y = element_text(colour = "white", face = "bold"))+
xlab("Season, Polygon")
p
Test dispersal limitation
- Loading packages
Calculating distances between points
coords<- read_csv("../Data/coord.csv")
coords2=read.csv("../Data/coord_lluvi.csv")
coords_mod=coords %>% mutate(SampleID=paste0(pol, Sitio, Transecto, "D")) %>% dplyr::select(
SampleID, Longitude, Latitude)
coords_mod2=coords2 %>% mutate(SampleID=paste0(pol, Sitio, Transecto, "R")) %>% dplyr::select(
SampleID, Longitude, Latitude)
#distance<- distm(coords_mat)/1000
#colnames(distance)<- rownames(coords_mat)
#rownames(distance)<- rownames(coords_mat)
#distance_complete<- distance
#distance[upper.tri(distance)] <- NA
#distance_dm<-melt(as.matrix(distance), varnames = c(
# "SampleID.x", "SampleID.y")) %>% drop_na() %>% filter(!value==0)
#BRAY-euclidean
#dry
otu_dry=abund_table[,match(coords_mod$SampleID, colnames(abund_table))]
otu_dist_dry = vegdist(t(otu_dry), method = "bray")
spatial_dist_dry=vegdist(coords_mod %>% column_to_rownames(var = "SampleID"), "euc")
mantel(otu_dist_dry,spatial_dist_dry)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = otu_dist_dry, ydis = spatial_dist_dry)
##
## Mantel statistic r: 0.4443
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0673 0.0853 0.1040 0.1249
## Permutation: free
## Number of permutations: 999#rainy
otu_rainy=abund_table[,match(coords_mod2$SampleID, colnames(abund_table))]
otu_dist_rainy = vegdist(t(otu_rainy), method = "bray")
spatial_dist_rainy=vegdist(coords_mod2 %>% column_to_rownames(var = "SampleID"), "euc")
mantel(otu_dist_rainy,spatial_dist_rainy)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = otu_dist_rainy, ydis = spatial_dist_rainy)
##
## Mantel statistic r: 0.4107
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0561 0.0757 0.1013 0.1323
## Permutation: free
## Number of permutations: 999#BRAY-distance in km
distance_dry<- distm(coords_mod %>% column_to_rownames(var = "SampleID"))/1000
colnames(distance_dry)<- coords_mod$SampleID
rownames(distance_dry)<- coords_mod$SampleID
distance_complete_dry<- distance_dry
distance_dry[upper.tri(distance_dry)] <- NA
distance_dm_dry<-melt(as.matrix(distance_dry), varnames = c(
"SampleID.x", "SampleID.y")) %>% drop_na() %>% filter(!value==0) %>%
dplyr::rename(spatial=value)
distance_otu_dry<-melt(as.matrix(otu_dist_dry), varnames = c(
"SampleID.x", "SampleID.y")) %>% drop_na() %>% filter(!value==0)%>%
dplyr::rename(comm=value)
mant1=mantel(otu_dist_dry,distance_complete_dry, method = "spearman")
dat1= distance_dm_dry %>% inner_join(distance_otu_dry)
fit1 <- lm(spatial ~ comm, data = dat1)
pdry=dat1 %>% ggplot(aes(spatial,comm))+ geom_point() +
stat_smooth(method = "lm", col = "red")+
labs(title = paste(" Slope of lm=",signif(fit1$coef[[2]], 5),
" r-mantel =",signif(mant1$statistic, 3),
"p-value =", signif(mant1$signif, 3)))+
theme_grey() +
theme(
axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 18),
legend.text = element_text(size = 10),
legend.title = element_text(size = 18),
legend.position = "right",
legend.box = "vertical",
panel.border = element_rect(
colour = "black",
fill = NA,
size = 1),
plot.title = element_text(size=18),
aspect.ratio = 7 / 10)+xlab("Spatial Distance (km)")+
ylab("Bray-curtis dissimilarity")
#BRAY-distance in km
distance_rainy<- distm(coords_mod2 %>% column_to_rownames(var = "SampleID"))/1000
colnames(distance_rainy)<- coords_mod2$SampleID
rownames(distance_rainy)<- coords_mod2$SampleID
distance_complete_rainy<- distance_rainy
distance_rainy[upper.tri(distance_rainy)] <- NA
distance_dm_rainy<-melt(as.matrix(distance_rainy), varnames = c(
"SampleID.x", "SampleID.y")) %>% drop_na() %>% filter(!value==0) %>%
dplyr::rename(spatial=value)
distance_otu_rainy<-melt(as.matrix(otu_dist_rainy), varnames = c(
"SampleID.x", "SampleID.y")) %>% drop_na() %>% filter(!value==0)%>%
dplyr::rename(comm=value)
mant2=mantel(otu_dist_rainy,distance_complete_rainy, method = "spearman")
dat2= distance_dm_rainy %>% inner_join(distance_otu_rainy) %>% filter(!comm < 0.3)
fit2 <- lm(spatial ~ comm, data = dat2)
prainy=dat2 %>% ggplot(aes(spatial,comm))+ geom_point() +
stat_smooth(method = "lm", col = "red")+
labs(title = paste(" Slope of lm=",signif(fit2$coef[[2]], 5),
" r-mantel =",signif(mant2$statistic, 3),
"p-value =", signif(mant2$signif, 3))) +
theme_grey() +
theme(
axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 18),
legend.text = element_text(size = 10),
legend.title = element_text(size = 18),
legend.position = "right",
legend.box = "vertical",
panel.border = element_rect(
colour = "black",
fill = NA,
size = 1),
plot.title = element_text(size=18),
aspect.ratio = 7 / 10)+xlab("Spatial Distance (km)")+
ylab("Bray-curtis dissimilarity")
library(cowplot)
plot<- plot_grid(pdry,prainy+ylab(""),
labels = c("A", "B"),
label_size = 20,
ncol = 2)
plot