Ecological index
Valerie De Anda
Instituto de Ecologia UNAM/ Estación Experimental de Aula Dei-CSIC, Zaragoza España http://valeriedeanda.com16/may./2018
`
source('http://bioconductor.org/biocLite.R')
biocLite('phyloseq')
install.packages("ggplot2")
install.packages("vegan")
install.packages("pvclust")library("ggplot2")
library("vegan")
library("pvclust")
library("phyloseq")1 Loading Metadata
metadata=read.table("../data/sample_data.txt", header=T)
c("1","1","1","2","2","2","3","3","3","4","4","4")->wl_order
cbind(wl_order,metadata)->metadata
metadata wl_order sample_ID site year samp_period water_level temp sp_comd
1 1 S1 A y2012 nov none NA NA
2 1 S2 B y2012 nov none NA NA
3 1 S3 C y2012 nov none NA NA
4 2 S4 A y2013 may low 20.00 8.360
5 2 S5 B y2013 may low 20.50 8.300
6 2 S6 C y2013 may low 20.30 8.300
7 3 S7 A y2013 oct high 26.50 6.530
8 3 S8 B y2013 oct high 26.16 6.400
9 3 S9 C y2013 oct high 25.42 6.478
10 4 S10 A y2014 may medium 24.99 11.590
11 4 S11 B y2014 may medium 25.45 11.580
12 4 S12 C y2014 may medium 25.22 11.600
salt pH season_year wl_order season
1 NA NA Autumn_12 1 Autumn
2 NA NA Autumn_12 1 Autumn
3 NA NA Autumn_12 1 Autumn
4 4.70 8.60 Spring_13 2 Spring
5 4.69 8.61 Spring_13 2 Spring
6 4.71 8.59 Spring_13 2 Spring
7 3.66 7.90 Autumn_13 3 Autumn
8 3.56 8.06 Autumn_13 3 Autumn
9 3.59 8.07 Autumn_13 3 Autumn
10 6.62 8.09 Spring_14 4 Spring
11 6.61 8.31 Spring_14 4 Spring
12 6.61 8.31 Spring_14 4 Spring2 Loading genera table
gen_otu<- as.matrix(read.table("../data/genus.tab", header=T, row.names=1))
genOTU=otu_table(gen_otu, taxa_are_rows=T)
mats2=phyloseq(genOTU)
data.frame(Sample_ID=metadata$sample_ID, Site=metadata$site, Year=metadata$year,
Period=metadata$samp_period, Water_level=metadata$water_level, Season=metadata$season_year,
wl_order=metadata$wl_order, Temperature=metadata$temp,Conductivity=metadata$sp_comd, Sal=metadata$salt,
pH=metadata$pH, row.names=sample_names((mats2)))-> o1
sampledata= sample_data(o1)
sampledata Sample_ID Site Year Period Water_level Season wl_order Temperature
S1 S1 A y2012 nov none Autumn_12 1 NA
S2 S2 B y2012 nov none Autumn_12 1 NA
S3 S3 C y2012 nov none Autumn_12 1 NA
S4 S4 A y2013 may low Spring_13 2 20.00
S5 S5 B y2013 may low Spring_13 2 20.50
S6 S6 C y2013 may low Spring_13 2 20.30
S7 S7 A y2013 oct high Autumn_13 3 26.50
S8 S8 B y2013 oct high Autumn_13 3 26.16
S9 S9 C y2013 oct high Autumn_13 3 25.42
S10 S10 A y2014 may medium Spring_14 4 24.99
S11 S11 B y2014 may medium Spring_14 4 25.45
S12 S12 C y2014 may medium Spring_14 4 25.22
Conductivity Sal pH
S1 NA NA NA
S2 NA NA NA
S3 NA NA NA
S4 8.360 4.70 8.60
S5 8.300 4.69 8.61
S6 8.300 4.71 8.59
S7 6.530 3.66 7.90
S8 6.400 3.56 8.06
S9 6.478 3.59 8.07
S10 11.590 6.62 8.09
S11 11.580 6.61 8.31
S12 11.600 6.61 8.313 Get the diversity index with Phyloseq
mats2=phyloseq(genOTU, sampledata)
estimate_richness(mats2)->div_alpha2
as.data.frame(t(gen_otu))->gen1
Pielou2<-vegan::diversity(gen1)/log(specnumber(gen1))
cbind (div_alpha2,Pielou2)->div_alpha2
write.table(div_alpha2, "../data/diversidad_genus.txt", sep="\t")
cbind (div_alpha2,Pielou2, metadata)->div_alpha2
div_alpha2 Observed Chao1 se.chao1 ACE se.ACE Shannon Simpson
S1 791 893.3529 22.24987 907.2252 15.07986 5.182797 0.9848085
S2 595 797.3239 41.70870 749.0884 13.54643 4.560236 0.9732375
S3 909 1082.7879 33.59760 1071.0955 16.51527 4.933919 0.9762544
S4 904 1078.6901 37.09365 1027.8185 16.05490 5.234497 0.9861490
S5 886 1053.0000 34.21966 1014.5400 15.83336 5.043441 0.9805729
S6 788 988.6800 40.76315 940.0626 15.41020 5.145337 0.9823072
S7 957 1103.3505 29.57780 1088.5515 16.52654 5.232292 0.9844902
S8 903 998.4130 21.66008 993.0575 15.65848 4.951102 0.9713006
S9 906 1087.8354 37.11671 1038.1193 15.94765 5.197987 0.9810412
S10 834 1002.7647 36.55029 942.0853 15.16222 5.159896 0.9815375
S11 774 933.0789 33.77470 896.7198 14.89700 4.786693 0.9678326
S12 893 1070.7051 36.58297 1025.1678 15.95616 5.195627 0.9822210
InvSimpson Fisher Pielou2 Pielou2 wl_order sample_ID site year
S1 65.82615 155.5618 0.7766470 0.7766470 1 S1 A y2012
S2 37.36565 126.3525 0.7138126 0.7138126 1 S2 B y2012
S3 42.11305 162.2370 0.7242615 0.7242615 1 S3 C y2012
S4 72.19675 162.7804 0.7690066 0.7690066 2 S4 A y2013
S5 51.47443 164.1531 0.7431341 0.7431341 2 S5 B y2013
S6 56.52022 160.9427 0.7714729 0.7714729 2 S6 C y2013
S7 64.47528 180.5816 0.7623021 0.7623021 3 S7 A y2013
S8 34.84397 165.0920 0.7274910 0.7274910 3 S8 B y2013
S9 52.74587 173.3913 0.7633951 0.7633951 3 S9 C y2013
S10 54.16386 162.2279 0.7671301 0.7671301 4 S10 A y2014
S11 31.08737 152.2825 0.7196333 0.7196333 4 S11 B y2014
S12 56.24616 168.2027 0.7646715 0.7646715 4 S12 C y2014
samp_period water_level temp sp_comd salt pH season_year wl_order
S1 nov none NA NA NA NA Autumn_12 1
S2 nov none NA NA NA NA Autumn_12 1
S3 nov none NA NA NA NA Autumn_12 1
S4 may low 20.00 8.360 4.70 8.60 Spring_13 2
S5 may low 20.50 8.300 4.69 8.61 Spring_13 2
S6 may low 20.30 8.300 4.71 8.59 Spring_13 2
S7 oct high 26.50 6.530 3.66 7.90 Autumn_13 3
S8 oct high 26.16 6.400 3.56 8.06 Autumn_13 3
S9 oct high 25.42 6.478 3.59 8.07 Autumn_13 3
S10 may medium 24.99 11.590 6.62 8.09 Spring_14 4
S11 may medium 25.45 11.580 6.61 8.31 Spring_14 4
S12 may medium 25.22 11.600 6.61 8.31 Spring_14 4
season
S1 Autumn
S2 Autumn
S3 Autumn
S4 Spring
S5 Spring
S6 Spring
S7 Autumn
S8 Autumn
S9 Autumn
S10 Spring
S11 Spring
S12 Spring4 Save output
mats2=phyloseq(genOTU, sampledata)
estimate_richness(mats2)->div_alpha2
as.data.frame(t(gen_otu))->gen1
Pielou2<-vegan::diversity(gen1)/log(specnumber(gen1))
cbind (div_alpha2,Pielou2)->div_alpha2
write.table(div_alpha2, "../data/alpharichness_genus.txt", sep="\t")
cbind (div_alpha2,Pielou2, metadata)->div_alpha2
div_alpha2 Observed Chao1 se.chao1 ACE se.ACE Shannon Simpson
S1 791 893.3529 22.24987 907.2252 15.07986 5.182797 0.9848085
S2 595 797.3239 41.70870 749.0884 13.54643 4.560236 0.9732375
S3 909 1082.7879 33.59760 1071.0955 16.51527 4.933919 0.9762544
S4 904 1078.6901 37.09365 1027.8185 16.05490 5.234497 0.9861490
S5 886 1053.0000 34.21966 1014.5400 15.83336 5.043441 0.9805729
S6 788 988.6800 40.76315 940.0626 15.41020 5.145337 0.9823072
S7 957 1103.3505 29.57780 1088.5515 16.52654 5.232292 0.9844902
S8 903 998.4130 21.66008 993.0575 15.65848 4.951102 0.9713006
S9 906 1087.8354 37.11671 1038.1193 15.94765 5.197987 0.9810412
S10 834 1002.7647 36.55029 942.0853 15.16222 5.159896 0.9815375
S11 774 933.0789 33.77470 896.7198 14.89700 4.786693 0.9678326
S12 893 1070.7051 36.58297 1025.1678 15.95616 5.195627 0.9822210
InvSimpson Fisher Pielou2 Pielou2 wl_order sample_ID site year
S1 65.82615 155.5618 0.7766470 0.7766470 1 S1 A y2012
S2 37.36565 126.3525 0.7138126 0.7138126 1 S2 B y2012
S3 42.11305 162.2370 0.7242615 0.7242615 1 S3 C y2012
S4 72.19675 162.7804 0.7690066 0.7690066 2 S4 A y2013
S5 51.47443 164.1531 0.7431341 0.7431341 2 S5 B y2013
S6 56.52022 160.9427 0.7714729 0.7714729 2 S6 C y2013
S7 64.47528 180.5816 0.7623021 0.7623021 3 S7 A y2013
S8 34.84397 165.0920 0.7274910 0.7274910 3 S8 B y2013
S9 52.74587 173.3913 0.7633951 0.7633951 3 S9 C y2013
S10 54.16386 162.2279 0.7671301 0.7671301 4 S10 A y2014
S11 31.08737 152.2825 0.7196333 0.7196333 4 S11 B y2014
S12 56.24616 168.2027 0.7646715 0.7646715 4 S12 C y2014
samp_period water_level temp sp_comd salt pH season_year wl_order
S1 nov none NA NA NA NA Autumn_12 1
S2 nov none NA NA NA NA Autumn_12 1
S3 nov none NA NA NA NA Autumn_12 1
S4 may low 20.00 8.360 4.70 8.60 Spring_13 2
S5 may low 20.50 8.300 4.69 8.61 Spring_13 2
S6 may low 20.30 8.300 4.71 8.59 Spring_13 2
S7 oct high 26.50 6.530 3.66 7.90 Autumn_13 3
S8 oct high 26.16 6.400 3.56 8.06 Autumn_13 3
S9 oct high 25.42 6.478 3.59 8.07 Autumn_13 3
S10 may medium 24.99 11.590 6.62 8.09 Spring_14 4
S11 may medium 25.45 11.580 6.61 8.31 Spring_14 4
S12 may medium 25.22 11.600 6.61 8.31 Spring_14 4
season
S1 Autumn
S2 Autumn
S3 Autumn
S4 Spring
S5 Spring
S6 Spring
S7 Autumn
S8 Autumn
S9 Autumn
S10 Spring
S11 Spring
S12 Spring5 Metabolic diversity
5.1 Create phyloseq object
pfam_tab<- as.matrix(read.table("../data/mats_abundances.sort.tab", header=T, row.names=1))
pfamOTU=otu_table(pfam_tab, taxa_are_rows=T)
pfam=phyloseq(pfamOTU)
cbind(wl_order,metadata)->metadata
data.frame(Sample_ID=metadata$sample_ID, Site=metadata$site, Year=metadata$year,
Period=metadata$samp_period, Water_level=metadata$water_level, Season=metadata$season_year,
wl_order=metadata$wl_order, Temperature=metadata$temp,Conductivity=metadata$sp_comd, Sal=metadata$salt,
pH=metadata$pH, row.names=sample_names((pfam)))-> o2
sampledata= sample_data(o2)
Sampledata= sample_data(o2)
pfam=phyloseq(pfamOTU, sampledata)6 Estimate richness
estimate_richness(pfam)->div_alpha_p
as.data.frame(t(pfam_tab))->sp2
Pielou<-vegan::diversity(sp2)/log(specnumber(sp2))
cbind (div_alpha_p,Pielou, metadata)->div_alpha_p
div_alpha_p Observed Chao1 se.chao1 ACE se.ACE Shannon Simpson
S1 6629 7916.856 92.49518 7786.558 44.12164 7.554437 0.9987244
S2 5167 6108.427 77.14891 5951.651 38.07480 7.564737 0.9988252
S3 6162 7151.772 80.13313 7034.084 41.61839 7.536980 0.9987023
S4 7090 8625.982 107.24925 8428.463 45.96354 7.568078 0.9987335
S5 6695 7978.183 94.12588 7803.488 44.01618 7.564408 0.9987472
S6 6123 7298.951 92.23206 7080.825 41.78794 7.485337 0.9986792
S7 7014 8649.187 114.58063 8350.412 45.63128 7.542102 0.9986914
S8 6378 7532.949 92.09878 7286.888 42.15693 7.491470 0.9986547
S9 6782 7972.291 89.60556 7793.669 43.75053 7.489023 0.9985754
S10 6399 7538.322 86.32907 7425.114 42.94990 7.517991 0.9986962
S11 5839 6744.615 74.79114 6683.076 40.48225 7.428460 0.9986031
S12 6490 7734.991 96.53433 7479.144 42.80138 7.475839 0.9986526
InvSimpson Fisher Pielou wl_order wl_order sample_ID site year
S1 783.9485 1332.195 0.8585359 1 1 S1 A y2012
S2 851.2295 1194.940 0.8847597 1 1 S2 B y2012
S3 770.6060 1164.599 0.8637227 1 1 S3 C y2012
S4 789.6073 1341.744 0.8535644 2 2 S4 A y2013
S5 798.1962 1296.143 0.8587023 2 2 S5 B y2013
S6 757.1067 1167.858 0.8584292 2 2 S6 C y2013
S7 764.1924 1315.046 0.8516699 3 3 S7 A y2013
S8 743.3185 1146.695 0.8551311 3 3 S8 B y2013
S9 701.9396 1241.477 0.8489004 3 3 S9 C y2013
S10 767.0120 1234.050 0.8578366 4 4 S10 A y2014
S11 715.8911 1079.173 0.8565718 4 4 S11 B y2014
S12 742.1443 1162.121 0.8516545 4 4 S12 C y2014
samp_period water_level temp sp_comd salt pH season_year wl_order
S1 nov none NA NA NA NA Autumn_12 1
S2 nov none NA NA NA NA Autumn_12 1
S3 nov none NA NA NA NA Autumn_12 1
S4 may low 20.00 8.360 4.70 8.60 Spring_13 2
S5 may low 20.50 8.300 4.69 8.61 Spring_13 2
S6 may low 20.30 8.300 4.71 8.59 Spring_13 2
S7 oct high 26.50 6.530 3.66 7.90 Autumn_13 3
S8 oct high 26.16 6.400 3.56 8.06 Autumn_13 3
S9 oct high 25.42 6.478 3.59 8.07 Autumn_13 3
S10 may medium 24.99 11.590 6.62 8.09 Spring_14 4
S11 may medium 25.45 11.580 6.61 8.31 Spring_14 4
S12 may medium 25.22 11.600 6.61 8.31 Spring_14 4
season
S1 Autumn
S2 Autumn
S3 Autumn
S4 Spring
S5 Spring
S6 Spring
S7 Autumn
S8 Autumn
S9 Autumn
S10 Spring
S11 Spring
S12 Springwrite.table(div_alpha_p, "../data/alpharichness_pfam.txt", sep="\t")7 Plot richness
# Plot diversity by Pfams
#pdf("../figures/alpharichnes.pdf")
par(mfrow=c(2,1))
pG=plot_richness(mats2, x="Site", color="Season")
pG+ geom_point(size=5, alpha=0.7)
pf=plot_richness(pfam, x="Site", color="Season")
pf + geom_point(size=5, alpha=0.7)
#dev.off()8 Rarefaction curves
8.0.1 Site A
Acum_1_A<-(gen_otu[,1])
Acum_2_A<-rowSums(gen_otu[,c(1,4)])
Acum_3_A<-rowSums(gen_otu[,c(1,4,7)])
Acum_4_A<-rowSums(gen_otu[,c(1,4,7,10)])
cbind(Acum_1_A, Acum_2_A, Acum_3_A, Acum_4_A)->acumulativa.A
acumulativa.A[1:3, 1:4] Acum_1_A Acum_2_A Acum_3_A Acum_4_A
Bacteria_unclassified 1438 3363 5301 7702
Desulfatitalea 187 246 1212 1711
Proteobacteria_unclassified 1089 2750 4341 58128.0.2 Site A
Acum_1_B<-(gen_otu[,2])
Acum_2_B<-rowSums(gen_otu[,c(2,5)])
Acum_3_B<-rowSums(gen_otu[,c(2,5,8)])
Acum_4_B<-rowSums(gen_otu[,c(2,5,8,11)])
cbind(Acum_1_B, Acum_2_B, Acum_3_B, Acum_4_B)->acumulativa.B8.0.3 Site B
Acum_1_C<-(gen_otu[,3])
Acum_2_C<-rowSums(gen_otu[,c(3,6)])
Acum_3_C<-rowSums(gen_otu[,c(3,6,9)])
Acum_4_C<-rowSums(gen_otu[,c(3,6,9,12)])
cbind(Acum_1_C, Acum_2_C, Acum_3_C, Acum_4_C)->acumulativa.C9 Plot graph
#pdf("../figures/curvas.rarefacción.pdf")
rarecurve(t(acumulativa.A), main="Rarefaction curve Site A",leyend=T, ylab="Genera", col=c("cornflowerblue","lightcoral","orange","slateblue"))
legend("right", legend=c("acum1 = NOV12", "acum2=+MAY13","acum3=+OCT13","acum4=+MAY14"),fill=c("cornflowerblue","lightcoral","orange","slateblue"),title="Season",
cex= .7, title.adj=0.3, bty="0")
rarecurve(t(acumulativa.B), main="Rarefaction curve Site B",leyend=T, ylab="Genera", col=c("cornflowerblue","lightcoral","orange","slateblue"))
legend("right", legend=c("acum1 = NOV12", "acum2=+MAY13","acum3=+OCT13","acum4=+MAY14"),fill=c("cornflowerblue","lightcoral","orange","slateblue"),title="Season",
cex= .7, title.adj=0.3, bty="0")
rarecurve(t(acumulativa.C), main="Rarefaction curve at genus level Site C",leyend=T, ylab="Genera", col=c("cornflowerblue","lightcoral","orange","slateblue"))
legend("right", legend=c("acum1 = NOV12", "acum2=+MAY13","acum3=+OCT13","acum4=+MAY14"),fill=c("cornflowerblue","lightcoral","orange","slateblue"),title="Season",
cex= .7, title.adj=0.3, bty="0")
#dev.off()