Hello, I hope you are great. I am having a problem knowing how to represent in my PCOA my microorganisms greater than 0.01 of the total sample (all microorganisms), but on a monthly basis.
The only thing that I achieved and that will be observed in my script, are all the microorganisms that were greater than 0.01 during the whole study. I would first like to obtain those that were greater than 0.01 in the month of July, then in the month of August and finally September, to later gather that information, in my PCOA. Because in a barplot analysis I conducted earlier there are differences in the microorganism communities in the three months, and in the PCOA it only tells me that there is a difference in July and September, because in this analysis I did not establish that I only want to represent the communities greater than 0.01 (in the one of barplots if I established it).
library(tidyverse)
Family <- data.frame (tibble::tribble(
~Index, ~Cyanobacteriaceae, ~Microcystaceae, ~Nostocaceae, ~Phormidiaceae, ~Un.Nostocales, ~Nodosilineaceae, ~Cyanobacteriaceae,
"1A", 1, 0, 0, 0, 0, 1, 1,
"1C", 1, 0, 0, 2, 0, 1, 1,
"1D", 0, 0, 1, 0, 0, 1, 0,
"1E", 1, 0, 1, 0, 0, 0, 0,
"1F", 1, 0, 0, 0, 0, 0, 17,
"1G", 1, 0, 0, 0, 0, 0, 2,
"1H", 1, 112, 1, 583, 0, 1, 0,
"2A", 0, 0, 1, 0, 0, 1, 2,
"2B", 1, 0, 0, 35, 0, 1, 29,
"2C", 0, 1, 0, 0, 0, 1, 0,
"2D", 3, 0, 0, 0, 0, 0, 0,
"2E", 0, 0, 1, 0, 0, 0, 0,
"2F", 0, 3, 0, 1, 0, 0, 1,
"2H", 0, 0, 9, 18, 0, 0, 17,
"3A", 0, 0, 44, 88, 0, 0, 0,
"3B", 0, 0, 572, 60, 0, 38, 0,
"4A", 0, 0, 0, 0, 0, 0, 27,
"4B", 0, 0, 0, 0, 0, 0, 12,
"4E", 0, 0, 0, 0, 0, 0, 8,
"4G", 0, 0, 0, 0, 0, 0, 8,
"4H", 0, 3, 16, 28, 0, 0, 14,
"5E", 0, 0, 0, 17, 0, 0, 0,
"5F", 0, 2, 0, 2, 0, 0, 10,
"5G", 0, 2, 0, 28, 0, 0, 0,
"6A", 0, 0, 37, 0, 0, 0, 0,
"6F", 0, 0, 0, 0, 0, 0, 4,
"6H", 0, 55, 0, 123, 0, 0, 33,
"7D", 0, 0, 0, 128, 0, 0, 2,
"7F", 0, 2, 0, 269, 0, 0, 3,
"7G", 0, 0, 0, 0, 0, 0, 5,
"8A", 0, 0, 0, 0, 0, 0, 16,
"8B", 0, 0, 0, 0, 0, 0, 3,
"8C", 0, 0, 0, 29, 0, 0, 0,
"8D", 0, 3, 0, 0, 0, 0, 0,
"8E", 0, 0, 0, 102, 0, 16, 2,
"8G", 0, 93, 45, 439, 0, 0, 34,
"8H", 0, 35, 0, 24, 0, 0, 17,
"9A", 0, 0, 1, 0, 0, 1, 0,
"9B", 0, 0, 1, 0, 0, 1, 0,
"9D", 0, 0, 1, 0, 0, 0, 0,
"9E", 0, 0, 0, 125, 0, 0, 0,
"9F", 0, 2, 8, 0, 0, 0, 31,
"9G", 0, 0, 0, 0, 0, 0, 20,
"9H", 0, 340, 3670, 2668, 33, 38, 0,
"10A", 0, 0, 1, 0, 0, 0, 32,
"10B", 0, 0, 0, 1, 0, 0, 14,
"10D", 0, 0, 0, 47, 0, 0, 0,
"10E", 0, 1, 0, 0, 0, 0, 0,
"10F", 0, 0, 0, 0, 0, 0, 4,
"10G", 0, 1, 0, 40, 0, 0, 0,
"10H", 0, 1, 0, 0, 1, 0, 0,
"11A", 0, 0, 1, 0, 0, 1, 0,
"11B", 0, 0, 1, 0, 1, 0, 0,
"11C", 0, 0, 1, 0, 0, 0, 0,
"11D", 0, 0, 1, 0, 0, 0, 0,
"11E", 0, 0, 1, 0, 0, 0, 0,
"11F", 0, 0, 0, 1, 0, 0, 0,
"11G", 0, 0, 0, 1, 0, 0, 0,
"11H", 0, 0, 0, 1, 0, 0, 0,
"12D", 0, 20, 77, 145, 0, 0, 22,
"12F", 0, 0, 0, 0, 0, 0, 2,
"12G", 0, 68, 2, 296, 0, 0, 0,
"12H", 0, 148, 280, 216, 0, 0, 28
)
)
t(Family)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> Index "1A" "1C" "1D" "1E" "1F" "1G" "1H"
#> Cyanobacteriaceae "1" "1" "0" "1" "1" "1" "1"
#> Microcystaceae " 0" " 0" " 0" " 0" " 0" " 0" "112"
#> Nostocaceae " 0" " 0" " 1" " 1" " 0" " 0" " 1"
#> Phormidiaceae " 0" " 2" " 0" " 0" " 0" " 0" " 583"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 1" " 1" " 1" " 0" " 0" " 0" " 1"
#> Cyanobacteriaceae.1 " 1" " 1" " 0" " 0" "17" " 2" " 0"
#> [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> Index "2A" "2B" "2C" "2D" "2E" "2F" "2H"
#> Cyanobacteriaceae "0" "1" "0" "3" "0" "0" "0"
#> Microcystaceae " 0" " 0" " 1" " 0" " 0" " 3" " 0"
#> Nostocaceae " 1" " 0" " 0" " 0" " 1" " 0" " 9"
#> Phormidiaceae " 0" " 35" " 0" " 0" " 0" " 1" " 18"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 1" " 1" " 1" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 " 2" "29" " 0" " 0" " 0" " 1" "17"
#> [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> Index "3A" "3B" "4A" "4B" "4E" "4G" "4H"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 0" " 0" " 0" " 0" " 0" " 0" " 3"
#> Nostocaceae " 44" " 572" " 0" " 0" " 0" " 0" " 16"
#> Phormidiaceae " 88" " 60" " 0" " 0" " 0" " 0" " 28"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" "38" " 0" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 " 0" " 0" "27" "12" " 8" " 8" "14"
#> [,22] [,23] [,24] [,25] [,26] [,27] [,28]
#> Index "5E" "5F" "5G" "6A" "6F" "6H" "7D"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 0" " 2" " 2" " 0" " 0" " 55" " 0"
#> Nostocaceae " 0" " 0" " 0" " 37" " 0" " 0" " 0"
#> Phormidiaceae " 17" " 2" " 28" " 0" " 0" " 123" " 128"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 " 0" "10" " 0" " 0" " 4" "33" " 2"
#> [,29] [,30] [,31] [,32] [,33] [,34] [,35]
#> Index "7F" "7G" "8A" "8B" "8C" "8D" "8E"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 2" " 0" " 0" " 0" " 0" " 3" " 0"
#> Nostocaceae " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Phormidiaceae " 269" " 0" " 0" " 0" " 29" " 0" " 102"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" " 0" " 0" " 0" " 0" " 0" "16"
#> Cyanobacteriaceae.1 " 3" " 5" "16" " 3" " 0" " 0" " 2"
#> [,36] [,37] [,38] [,39] [,40] [,41] [,42]
#> Index "8G" "8H" "9A" "9B" "9D" "9E" "9F"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 93" " 35" " 0" " 0" " 0" " 0" " 2"
#> Nostocaceae " 45" " 0" " 1" " 1" " 1" " 0" " 8"
#> Phormidiaceae " 439" " 24" " 0" " 0" " 0" " 125" " 0"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" " 0" " 1" " 1" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 "34" "17" " 0" " 0" " 0" " 0" "31"
#> [,43] [,44] [,45] [,46] [,47] [,48] [,49]
#> Index "9G" "9H" "10A" "10B" "10D" "10E" "10F"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 0" "340" " 0" " 0" " 0" " 1" " 0"
#> Nostocaceae " 0" "3670" " 1" " 0" " 0" " 0" " 0"
#> Phormidiaceae " 0" "2668" " 0" " 1" " 47" " 0" " 0"
#> Un.Nostocales " 0" "33" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" "38" " 0" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 "20" " 0" "32" "14" " 0" " 0" " 4"
#> [,50] [,51] [,52] [,53] [,54] [,55] [,56]
#> Index "10G" "10H" "11A" "11B" "11C" "11D" "11E"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 1" " 1" " 0" " 0" " 0" " 0" " 0"
#> Nostocaceae " 0" " 0" " 1" " 1" " 1" " 1" " 1"
#> Phormidiaceae " 40" " 0" " 0" " 0" " 0" " 0" " 0"
#> Un.Nostocales " 0" " 1" " 0" " 1" " 0" " 0" " 0"
#> Nodosilineaceae " 0" " 0" " 1" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> [,57] [,58] [,59] [,60] [,61] [,62] [,63]
#> Index "11F" "11G" "11H" "12D" "12F" "12G" "12H"
#> Cyanobacteriaceae "0" "0" "0" "0" "0" "0" "0"
#> Microcystaceae " 0" " 0" " 0" " 20" " 0" " 68" "148"
#> Nostocaceae " 0" " 0" " 0" " 77" " 0" " 2" " 280"
#> Phormidiaceae " 1" " 1" " 1" " 145" " 0" " 296" " 216"
#> Un.Nostocales " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Nodosilineaceae " 0" " 0" " 0" " 0" " 0" " 0" " 0"
#> Cyanobacteriaceae.1 " 0" " 0" " 0" "22" " 2" " 0" "28"
#NORMALIZATION OF RAW READ TABLE LOG2
#OTUS en filas, columnas son las muestras, es al reves que para el resto de pasos
data<- Family
replicates <- as.data.frame(colnames(data)[-1])
colnames(replicates) <- "replicates"
attach(Family)
rwnames <- Index
data <- as.matrix(data[,-1])
rownames(data) <- rwnames
data[is.na(data)] <- 0
library(DESeq2)
#> Loading required package: S4Vectors
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: parallel
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:parallel':
#>
#> clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
#> clusterExport, clusterMap, parApply, parCapply, parLapply,
#> parLapplyLB, parRapply, parSapply, parSapplyLB
#> The following objects are masked from 'package:dplyr':
#>
#> combine, intersect, setdiff, union
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, append, as.data.frame, basename, cbind,
#> colMeans, colnames, colSums, dirname, do.call, duplicated,
#> eval, evalq, Filter, Find, get, grep, grepl, intersect,
#> is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
#> paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
#> Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
#> table, tapply, union, unique, unsplit, which, which.max,
#> which.min
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:dplyr':
#>
#> first, rename
#> The following object is masked from 'package:tidyr':
#>
#> expand
#> The following object is masked from 'package:base':
#>
#> expand.grid
#> Loading required package: IRanges
#>
#> Attaching package: 'IRanges'
#> The following objects are masked from 'package:dplyr':
#>
#> collapse, desc, slice
#> The following object is masked from 'package:purrr':
#>
#> reduce
#> The following object is masked from 'package:grDevices':
#>
#> windows
#> Loading required package: GenomicRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: SummarizedExperiment
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#> Loading required package: DelayedArray
#> Loading required package: matrixStats
#>
#> Attaching package: 'matrixStats'
#> The following objects are masked from 'package:Biobase':
#>
#> anyMissing, rowMedians
#> The following object is masked from 'package:dplyr':
#>
#> count
#> Loading required package: BiocParallel
#>
#> Attaching package: 'DelayedArray'
#> The following objects are masked from 'package:matrixStats':
#>
#> colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
#> The following object is masked from 'package:purrr':
#>
#> simplify
#> The following objects are masked from 'package:base':
#>
#> aperm, apply
dds <- DESeqDataSetFromMatrix(data, replicates, ~ replicates)
#> converting counts to integer mode
cts <- counts(dds)
geoMeans <- apply(cts, 1, function(row) if (all(row == 0)) 0 else exp(mean(log(row[row != 0]))))
dds <- estimateSizeFactors(dds, geoMeans=geoMeans)
norm <- counts(dds, normalized=TRUE)
logcounts <- log2( counts(dds, normalized=TRUE) + 1 )# in case you want a log transformed normalized count
write.csv(logcounts, file = "~/RSTUDIO/Datos_cianobacterias/ciano_rdp_normalizados_pruebaa.csv")
#log counts are like variance stabilized counts (https://support.bioconductor.org/p/76712/)
library(readxl)
Family <- read_excel("~/RSTUDIO/Datos_cianobacterias/ciano_rdp_normalizados_pruebaa.xlsx")
attach(Family)
#> The following objects are masked from Family (pos = 17):
#>
#> Cyanobacteriaceae, Cyanobacteriaceae.1, Microcystaceae,
#> Nodosilineaceae, Nostocaceae, Phormidiaceae, Un.Nostocales
Family <- Family[,-1]
rownames(Family) <- SampleID
#> Warning: Setting row names on a tibble is deprecated.
#remove the rare microbes. It keeps only the microbes that are present in at least 10% of the samples
dim(Family)
#> [1] 64 7
Family <- Family[,colMeans(Family) >=.1]
dim(Family)
#> [1] 64 7
#Family
Family <- data.frame(Family)
Family_counts <- colSums(Family)
Counts <- unname(Family_counts)
Family_counts <- data.frame(Family_counts)
Family_counts <- t(Family_counts)
total <- sum(Counts)
rel_ab <- Family_counts/total
Others <- rel_ab[,colMeans(rel_ab)<.01]
Others <- sum(Others)
rel_ab <- rel_ab[,colMeans(rel_ab)>=.01]
rel_ab <- data.frame(t(rel_ab), Others)
rel_ab_P <- t(rel_ab)
abundance <- c("abundance")
rel_ab_P <- data.frame(rel_ab_P)
write.csv(rel_ab_P, file = "~/RSTUDIO/Datos_cianobacterias/ciano_rdp_family.csv")
Famiy <- read.csv("~/RSTUDIO/Datos_cianobacterias/ciano_rdp_family.csv")
#(Transposing)
library(readxl)
Family <- read_excel("~/RSTUDIO/Datos_cianobacterias/ciano_rdp_normalizados_transpuesto.xlsx")
Metadata <-data.frame (tibble::tribble(
~SampleID, ~SamplingPoint, ~Depth, ~Month, ~Filter, ~SampleorReplica,
"1A", "CEA1", "80 cm", "July", "20-25 µm", "Sample",
"1C", "CEA4", "Interstitial", "July", "20-25 µm", "Sample",
"1D", "CEA1", "80 cm", "August", "20-25 µm", "Sample",
"1E", "CEA3", "80 cm", "August", "20-25 µm", "Sample",
"1F", "CEA5", "80 cm", "August", "20-25 µm", "Sample",
"1G", "CEA2", "80 cm", "September", "20-25 µm", "Sample",
"1H", "CEA4", "80 cm", "September", "20-25 µm", "Sample",
"2A", "CEA1", "80 cm", "July", "0.45 µm", "Sample",
"2B", "CEA2", "Interstitial", "July", "20-25 µm", "Sample",
"2C", "CEA4", "Interstitial", "July", "0.45 µm", "Sample",
"2D", "CEA1", "80 cm", "August", "0.45 µm", "Sample",
"2E", "CEA3", "80 cm", "August", "0.45 µm", "Sample",
"2F", "CEA5", "80 cm", "August", "0.45 µm", "Sample",
"2H", "CEA4", "80 cm", "September", "0.45 µm", "Sample",
"3A", "CEA1", "Interstitial", "July", "20-25 µm", "Sample",
"3B", "CEA2", "Interstitial", "July", "0.45 µm", "Sample",
"4A", "CEA1", "Interstitial", "July", "0.45 µm", "Sample",
"4B", "CEA3", "80 cm", "July", "20-25 µm", "Sample",
"4E", "CEA3", "Interstitial", "August", "0.45 µm", "Sample",
"4G", "CEA2", "Interstitial", "September", "0.45 µm", "Sample",
"4H", "CEA4", "Interstitial", "September", "0.45 µm", "Sample",
"5E", "CEA3", "80 cm", "August", "0.45 µm", "Replica",
"5F", "CEA5", "80 cm", "August", "0.45 µm", "Replica",
"5G", "CEA2", "80 cm", "September", "0.45 µm", "Replica",
"6A", "CEA1", "80 cm", "July", "0.45 µm", "Replica",
"6F", "CEA5", "Interstitial", "August", "0.45 µm", "Replica",
"6H", "CEA4", "Interstitial", "September", "0.45 µm", "Replica",
"7D", "CEA2", "80 cm", "August", "20-25 µm", "Sample",
"7F", "CEA1", "80 cm", "September", "20-25 µm", "Sample",
"7G", "CEA3", "80 cm", "September", "20-25 µm", "Sample",
"8A", "CEA1", "Interstitial", "July", "0.45 µm", "Replica",
"8B", "CEA3", "80 cm", "July", "0.45 µm", "Replica",
"8C", "CEA5", "Interstitial", "July", "20-25 µm", "Sample",
"8D", "CEA2", "80 cm", "August", "0.45 µm", "Sample",
"8E", "CEA4", "80 cm", "August", "0.45 µm", "Sample",
"8G", "CEA3", "80 cm", "September", "0.45 µm", "Sample",
"8H", "CEA5", "80 cm", "September", "0.45 µm", "Sample",
"9A", "CEA2", "80 cm", "July", "20-25 µm", "Sample",
"9B", "CEA3", "Interstitial", "July", "20-25 µm", "Replica",
"9D", "CEA2", "Interstitial", "August", "20-25 µm", "Sample",
"9E", "CEA4", "Interstitial", "August", "20-25 µm", "Sample",
"9F", "CEA1", "Interstitial", "September", "20-25 µm", "Sample",
"9G", "CEA3", "Interstitial", "September", "20-25 µm", "Sample",
"9H", "CEA5", "Interstitial", "September", "20-25 µm", "Sample",
"10A", "CEA2", "80 cm", "July", "0.45 µm", "Replica",
"10B", "CEA3", "Interstitial", "July", "0.45 µm", "Replica",
"10D", "CEA2", "Interstitial", "August", "0.45 µm", "Sample",
"10E", "CEA4", "Interstitial", "August", "0.45 µm", "Sample",
"10F", "CEA1", "Interstitial", "September", "0.45 µm", "Sample",
"10G", "CEA3", "Interstitial", "September", "0.45 µm", "Sample",
"10H", "CEA5", "Interstitial", "September", "0.45 µm", "Sample",
"11A", "CEA2", "Interstitial", "July", "20-25 µm", "Replica",
"11B", "CEA4", "80 cm", "July", "20-25 µm", "Sample",
"11C", "CEA5", "Interstitial", "July", "20-25 µm", "Replica",
"11D", "CEA2", "80 cm", "August", "0.45 µm", "Replica",
"11E", "CEA4", "80 cm", "August", "0.45 µm", "Replica",
"11F", "CEA1", "80 cm", "September", "0.45 µm", "Replica",
"11G", "CEA3", "80 cm", "September", "0.45 µm", "Replica",
"11H", "CEA5", "80 cm", "September", "0.45 µm", "Replica",
"12D", "CEA2", "Interstitial", "August", "0.45 µm", "Replica",
"12F", "CEA1", "Interstitial", "September", "0.45 µm", "Replica",
"12G", "CEA3", "Interstitial", "September", "0.45 µm", "Replica",
"12H", "CEA5", "Interstitial", "September", "0.45 µm", "Replica"
)
)
rownames(Metadata) <- SampleID
#> Error in `.rowNamesDF<-`(x, value = value): invalid 'row.names' length
Metadata$Month <- factor(Metadata$Month,
levels = c("July", "August", "September"))
library(vegan)
#PCoA
#create a distance matrix
dis.bray<-vegdist(Family, method="bray")
#> Error in vegdist(Family, method = "bray"): input data must be numeric
Family.cmd <- cmdscale(dis.bray, k=(NROW(Family)-1), eig=TRUE, add= TRUE)
#> Error in cmdscale(dis.bray, k = (NROW(Family) - 1), eig = TRUE, add = TRUE): objeto 'dis.bray' no encontrado
Family.cmd$eig
#> Error in eval(expr, envir, enclos): objeto 'Family.cmd' no encontrado
eig2<-eigenvals(Family.cmd)
#> Error in eigenvals(Family.cmd): objeto 'Family.cmd' no encontrado
eig2
#> Error in eval(expr, envir, enclos): objeto 'eig2' no encontrado
eig2/sum(eig2) #the first two are the axes values of axes 1 and 2 (https://stat.ethz.ch/pipermail/r-sig-ecology/2011-June/002183.html)
#> Error in eval(expr, envir, enclos): objeto 'eig2' no encontrado
#0.16225839 0.11075892
x<-Family.cmd$points[,1]
#> Error in eval(expr, envir, enclos): objeto 'Family.cmd' no encontrado
y<-Family.cmd$points[,2]
#> Error in eval(expr, envir, enclos): objeto 'Family.cmd' no encontrado
plot(x,y, type="n", xlab = "PCoA1 (%)", ylab = "PCoA2 (%)", asp = 1, axes = TRUE)
#> Error in plot(x, y, type = "n", xlab = "PCoA1 (%)", ylab = "PCoA2 (%)", : objeto 'x' no encontrado
text(x,y, rownames(Family.cmd$points), cex =0.7)
#> Error in text(x, y, rownames(Family.cmd$points), cex = 0.7): objeto 'x' no encontrado
Family.pco <- data.frame(x,y)
#> Error in data.frame(x, y): objeto 'x' no encontrado
ordiplot(Family.pco, choices = c(1, 2), display = "sites", type = "none", xlab = "PCoA 1",ylim = c(-0.4,0.4), ylab = "PCoA 2")
#> Error in ordiplot(Family.pco, choices = c(1, 2), display = "sites", type = "none", : objeto 'Family.pco' no encontrado
pl <-ordiellipse(Family.pco, Metadata$Month, kind="se", conf=0.95, lwd=2, col="gray30", label=T)
#> Error in scores(ord, display = display, ...): objeto 'Family.pco' no encontrado
#with(Metadata, points(Clase.pco$x, Clase.pco$y, display="sites", col=colvec[Crop], pch=21, bg=colvec[Crop]))
#with(Metadata, legend("topright", legend=levels(Crop), bty="n", col=colvec, pch=21, pt.bg=colvec))
#plot with ggplot
data.scores <- as.data.frame(scores(Family.cmd))
#> Error in scores(Family.cmd): objeto 'Family.cmd' no encontrado
data.scores$site <- rownames(data.scores)
#> Error in rownames(data.scores): objeto 'data.scores' no encontrado
data.scores$grp <- Metadata$Month
#> Error in data.scores$grp <- Metadata$Month: objeto 'data.scores' no encontrado
head(data.scores)
#> Error in head(data.scores): objeto 'data.scores' no encontrado
library(ggplot2)
NMDS <- data.frame(Dim1 = Family.mds$points[,1], Dim2 = Family.mds$points[,2],group=Metadata$Month)
#> Error in data.frame(Dim1 = Family.mds$points[, 1], Dim2 = Family.mds$points[, : objeto 'Family.mds' no encontrado
#NMDS.mean <- aggregate(NMDS[,1:2],list(group=group),mean)
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100)
{
theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))
}
#attach(NMDS.mean)
df_ell <- data.frame()
for(g in levels(NMDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$group==g,],
veganCovEllipse(pl[[g]]$cov,pl[[g]]$center,pl[[g]]$scale)))
,group=g))
}
#> Error in levels(NMDS$group): objeto 'NMDS' no encontrado
ggplot() +
geom_path(data=df_ell, aes(x=x, y=y,colour=group), size=1, linetype=1)+
geom_point(data=data.scores,aes(x=Dim1,y=Dim2,colour=grp),size=1.5) +
# add the point markers
theme_light()+
#geom_text(data=data.scores,aes(x=NMDS1,y=NMDS2,label=site),size=6,vjust=0) +
scale_colour_manual(values=c("July" = "deeppink", "August" = "turquoise4", "September"="darkorchid"))+
ylab("Family PCoA2 (22.00%)")+
xlab("Family PCoA1 (24.10%)")
#> Error in fortify(data): objeto 'data.scores' no encontrado
<sup>Created on 2020-02-22 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)</sup>
