longitudinal paper script
qcvp
2025-01-14
Loading Custom scripts
Before analyzing, run the following codes:
# OTU to matrix function
source("/Users/qcvp/R/custom script/OTU_to_matrix.R")
# ggrare function
source("/Users/qcvp/R/custom script/ggrare.R")Load context
Importing ASV table of Nem
asv_table <- utils::read.table("/Users/qcvp/R/outputs/dada2_total/asv_table/asv_table.tsv", sep="\t", row.names=1, header=T)Invert the ASV table in case sample IDs are on row not column
asv_inverted <- base::as.data.frame(t(asv_table))Load the ASV sequences (previously created with DADA2)
asv_seq <- Biostrings::readDNAStringSet(file = "/Users/qcvp/R/outputs/dada2_total/asv_table/asv.fasta",
format = "fasta",
nrec = -1L,
skip = 0L,
seek.first.rec = FALSE,
use.names = TRUE)
head(asv_seq, 1)## DNAStringSet object of length 1:
## width seq names
## [1] 370 TACGGAAGGTCCGGGCGTTATCC...TGGGGAGTACGCCGGCAACGGTG ASV1Load asv tax of Nem
asv_tax <- utils::read.table("/Users/qcvp/R/outputs/dada2_total/asv_table/taxonomy.tsv", header = TRUE, row.names = 1, sep="\t")
asv_tax <- base::as.matrix(asv_tax)
head(asv_tax, 1)## Kingdom Phylum Class Order Family
## ASV1 "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Prevotellaceae"
## Genus
## ASV1 "Prevotella_9"Import additional the metadata information of Nem
context <- utils::read.table("/Users/qcvp/R/outputs/dada2_total/asv_table/metadata.txt", header = TRUE, row.names = 1, sep="\t")
rownames(context) <- base::gsub("-",".",rownames(context))
head(context)## Sequencing Date Transect Station Biosample Biosample_type
## BD.S.1.M Metagenome Jul_2023 Bich_Dam BD Environment Sediment
## BD.S.2.M Metagenome Jul_2023 Bich_Dam BD Environment Sediment
## BD.S.3.M Metagenome Jul_2023 Bich_Dam BD Environment Sediment
## BD.S.4.B Barcoding Dec_2023 Bich_Dam BD Environment Sediment
## BD.S.5.B Barcoding Dec_2023 Bich_Dam BD Environment Sediment
## BD.S.6.B Barcoding Dec_2023 Bich_Dam BD Environment Sediment
## Biosample_cluster Host_type Host_subtype Host_genus Host_species
## BD.S.1.M Sediment Environment Sediment <NA> <NA>
## BD.S.2.M Sediment Environment Sediment <NA> <NA>
## BD.S.3.M Sediment Environment Sediment <NA> <NA>
## BD.S.4.B Sediment Environment Sediment <NA> <NA>
## BD.S.5.B Sediment Environment Sediment <NA> <NA>
## BD.S.6.B Sediment Environment Sediment <NA> <NA>
## Species_ID
## BD.S.1.M <NA>
## BD.S.2.M <NA>
## BD.S.3.M <NA>
## BD.S.4.B <NA>
## BD.S.5.B <NA>
## BD.S.6.B <NA>#Checking samples IDs uniformity
base::setdiff(x = row.names(asv_inverted), y = row.names(context)) ## character(0)Making raw phyloseq objects
physeq_Nem <- phyloseq::phyloseq(
phyloseq::otu_table(asv_inverted, taxa_are_rows = F), #if you have taxa on rows = T
phyloseq::sample_data(context),
phyloseq::refseq(asv_seq),
phyloseq::tax_table(asv_tax)
)
physeq_Nem## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111454 taxa and 412 samples ]
## sample_data() Sample Data: [ 412 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 111454 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111454 reference sequences ]Decontamination process
Quick checking library size of each Biomsaple types
library(ggplot2)
df <- base::as.data.frame(phyloseq::sample_data(physeq_Nem)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- phyloseq::sample_sums(physeq_Nem) # calculate sample size
df <- df[base::order(df$LibrarySize),] # Odering sample size
df$Index <- base::seq(base::nrow(df)) # Add STT
ggplot(data=df, aes(x=Index, y=LibrarySize, color=as.factor(Biosample_type))) + #plotting
geom_point() + geom_text(aes(label=Index), hjust=2, nudge_y = 0.05)
Identify contaminants ASVs with decontam function
#identify control samples
phyloseq::sample_data(physeq_Nem)$is.neg <- phyloseq::sample_data(physeq_Nem)$Biosample == "Control"
contam_df_prev05 <- decontam::isContaminant(physeq_Nem, method="prevalence", neg="is.neg", threshold=0.5)
table(contam_df_prev05$contaminant)##
## FALSE TRUE
## 111125 329#list of position ASVs which are contaminant
which(contam_df_prev05$contaminant) ## [1] 2 4 8 10 16 21 29 34 35 36 38 49
## [13] 53 70 73 102 107 128 130 138 142 168 171 172
## [25] 176 181 183 191 195 204 220 228 237 248 275 290
## [37] 301 306 309 311 334 350 353 363 376 380 381 382
## [49] 383 391 401 405 406 412 430 436 442 445 452 456
## [61] 468 471 481 490 509 520 521 522 558 570 577 581
## [73] 591 621 642 671 704 706 708 722 724 726 730 743
## [85] 744 747 761 762 795 798 799 802 803 816 820 824
## [97] 830 844 862 904 909 914 915 937 975 991 996 1000
## [109] 1017 1034 1056 1101 1104 1105 1116 1126 1127 1146 1151 1190
## [121] 1198 1215 1217 1224 1270 1280 1281 1284 1295 1301 1302 1332
## [133] 1347 1348 1349 1358 1367 1376 1391 1405 1462 1464 1472 1534
## [145] 1598 1607 1613 1640 1642 1666 1699 1711 1752 1758 1776 1798
## [157] 1816 1850 1859 1884 1915 1918 1919 1930 1952 1990 1991 2069
## [169] 2129 2137 2144 2178 2233 2274 2289 2321 2355 2360 2366 2390
## [181] 2397 2404 2413 2419 2427 2515 2530 2583 2595 2596 2601 2641
## [193] 2647 2694 2713 2763 2892 2924 2983 2995 3007 3041 3052 3079
## [205] 3178 3205 3213 3215 3228 3302 3414 3415 3485 3517 3520 3564
## [217] 3580 3598 3676 3747 3816 3831 3855 3869 3871 3881 3902 3922
## [229] 3967 3970 4028 4125 4134 4144 4150 4442 4453 4548 4626 4657
## [241] 4658 4710 4859 4905 4934 5066 5155 5163 5223 5244 5251 5327
## [253] 5333 5384 5403 5538 5574 5621 5650 5767 5932 5981 6077 6103
## [265] 6177 6265 6296 6337 6962 7029 7101 7144 7382 7548 7592 7645
## [277] 7647 7652 7808 7848 8222 8329 8389 8666 8947 8991 9080 9305
## [289] 9471 9497 9643 9739 9969 10144 10273 10372 10384 10513 10599 10706
## [301] 10956 11190 11359 11365 11902 12278 12435 12588 13964 14121 14642 15508
## [313] 16642 17814 19218 20331 21825 22059 24068 24298 24625 24837 26448 27917
## [325] 29044 30253 36085 46475 64190row_indices <- which(contam_df_prev05$contaminant) #get position of comtaminant
tax <- data.frame(physeq_Nem@tax_table) #get taxo table
Contam_taxonomy_table <- dplyr::tibble() # Create a blank data frame to store result
for (i in row_indices){
loc <- contam_df_prev05[i, 0] #extract comtanimant ASV
tax_key <- row.names(loc) #extract comtanimant ASV name
tax_value <- tax[tax_key, ] #extract taxo info
Contam_taxonomy_table <- base::rbind(Contam_taxonomy_table, tax_value) #adding contam ASV information
}
head(Contam_taxonomy_table)## Kingdom Phylum Class Order
## ASV2 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales
## ASV4 Bacteria Bacteroidota Bacteroidia Bacteroidales
## ASV8 Bacteria Firmicutes Bacilli Lactobacillales
## ASV10 Bacteria Firmicutes Bacilli Staphylococcales
## ASV16 Bacteria Bacteroidota Bacteroidia Bacteroidales
## ASV21 Bacteria Firmicutes Clostridia Clostridiales
## Family Genus
## ASV2 <NA> <NA>
## ASV4 Prevotellaceae Prevotella_9
## ASV8 Listeriaceae Listeria
## ASV10 Staphylococcaceae Staphylococcus
## ASV16 Prevotellaceae Prevotella_9
## ASV21 Clostridiaceae Clostridium sensu stricto 4Remove Control and MOCK samples
physeq_Nem_1 <- phyloseq::subset_samples(physeq_Nem, Biosample != "Control")
physeq_Nem_1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111454 taxa and 403 samples ]
## sample_data() Sample Data: [ 403 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 111454 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111454 reference sequences ]Remove comtaminant ASVs on list generated by decontam
physeq_Nem_2 <- phyloseq::subset_taxa(physeq_Nem_1, !taxa_names(physeq_Nem) %in% c("ASV2", "ASV8", "ASV52", "ASV97", "ASV128", "ASV318", "ASV356", "ASV471", "ASV708", "ASV914", "ASV1101", "ASV1127", "ASV1915", "ASV1952", "ASV2892", "ASV3517", "ASV4859", "ASV9739", "ASV26448", "ASV30253"))
physeq_Nem_2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111434 taxa and 403 samples ]
## sample_data() Sample Data: [ 403 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 111434 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111434 reference sequences ]Remove addtional genera
physeq_Nem_3 <- phyloseq::subset_taxa(physeq_Nem_2, !Genus %in% c("Cupriavidus", "Bradyrhizobium", "Listeria"))
physeq_Nem_3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111401 taxa and 403 samples ]
## sample_data() Sample Data: [ 403 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 111401 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111401 reference sequences ]Remove contaminant ASVs with kit
# Identify ASVs presence in both lists
contaminant_genera <- read.table("/Users/qcvp/R/Contaminant.txt", header = T)
contam_vec <- base::which(Contam_taxonomy_table[,6] %in% contaminant_genera$contaminant) #Identify comtaminant taxa in decontam function and contaminant kit
contam_ASV <- Contam_taxonomy_table[contam_vec,] #how this work?
contam_ASV_vec <- row.names(contam_ASV) #Extract ASV name
# Remove contam ASV present in contam kit
physeq_Nem_4 <- phyloseq::subset_taxa(physeq_Nem_3, !taxa_names(physeq_Nem) %in% contam_ASV_vec)
physeq_Nem_4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111338 taxa and 403 samples ]
## sample_data() Sample Data: [ 403 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 111338 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111338 reference sequences ]# Save contam list
writexl::write_xlsx(contam_ASV, "/Users/qcvp/R/outputs/Decomtaminated/contam_ASV_kit.xlsx")Removing samples with < 5000 reads
physeq_Nem_5 <- phyloseq::subset_samples(physeq_Nem_4, !sample_names(physeq_Nem_4) %in% c("Fi.4.3.GM", "Hu.35.SM","Sq.11.3.GB", "Fi.4.1.GB", "Hu.40.SM", "Oy.13.3.GM", "Fi.14.1.GM", "Hu.33.SM", "Hu.34.SM", "Hu.32.SM", "Fo.10.4.BM", "Fo.10.3.BM", "Fo.17.3.BM", "Fo.15.3.BM", "Fo.17.2.BM", "Oy.13.2.GB", "Fo.17.1.BM", "Fo.15.1.BM", "Fi.17.8.GM", "Fo.15.2.BM", "Oy.15.7.GB", "Lo.9.1.GB", "Lo.5.3.GB", "Lo.2.6.GB"))
physeq_Nem_5## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 111338 taxa and 379 samples ]
## sample_data() Sample Data: [ 379 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 111338 taxa by 6 taxonomic ranks ]
## refseq() DNAStringSet: [ 111338 reference sequences ]writexl::write_xlsx(as.data.frame(physeq_Nem_5@tax_table@.Data), "/Users/qcvp/R/outputs/Decomtaminated/ASV_after_decontam.xlsx")Checking data
max(rowSums(physeq_Nem_5@otu_table))## [1] 81231min(rowSums(physeq_Nem_5@otu_table))## [1] 6387mean(rowSums(physeq_Nem_5@otu_table))## [1] 43003.36Plotting bar plot of reads
library(ggplot2)
library(dplyr)
# Extract row sum of physeq_Nemt_2
OTU_asmatrix_2 <- otu.matrix(physeq_Nem_2, force.taxa.cols = TRUE)
sample_reads_2 <- rowSums(OTU_asmatrix_2)
# Extract row sum of physeq_Nemt_3
OTU_asmatrix_3 <- otu.matrix(physeq_Nem_3, force.taxa.cols = TRUE)
sample_reads_3 <- rowSums(OTU_asmatrix_3)
# Extract row sum of physeq_Nem_4
OTU_asmatrix_4 <- otu.matrix(physeq_Nem_4, force.taxa.cols = TRUE)
sample_reads_4 <- rowSums(OTU_asmatrix_4)
#Combine
sample_reads_total <- cbind(sample_reads_4, sample_reads_3, sample_reads_2)
# Load data
Abs_abun_arg <- sample_reads_total
# Reshape the data to long format
Abs_abun_arg_long <- reshape2::melt(Abs_abun_arg, id.vars = "Sample.name", variable.name = "Runs", value.name = "Reads")
Biosample_type <- physeq_Nem_4@sam_data$Biosample_type
Abs_abun_arg_long <- cbind( Abs_abun_arg_long, Biosample_type)
#Calculate errors
Abs_abun_summary <- Abs_abun_arg_long %>%
group_by(Biosample_type, Var2) %>%
summarise(
Mean = mean(Reads, na.rm = TRUE),
SD = sd(Reads, na.rm = TRUE))
Abs_abun_summary## # A tibble: 18 × 4
## # Groups: Biosample_type [6]
## Biosample_type Var2 Mean SD
## <chr> <fct> <dbl> <dbl>
## 1 Biofilm sample_reads_4 48132. 12695.
## 2 Biofilm sample_reads_3 48785. 12798.
## 3 Biofilm sample_reads_2 48793. 12809.
## 4 Gut sample_reads_4 41567. 17534.
## 5 Gut sample_reads_3 41936. 17616.
## 6 Gut sample_reads_2 41953. 17625.
## 7 Sediment sample_reads_4 31633. 12671.
## 8 Sediment sample_reads_3 31950. 12767.
## 9 Sediment sample_reads_2 31951. 12768.
## 10 Skin sample_reads_4 45718. 17601.
## 11 Skin sample_reads_3 46197. 17853.
## 12 Skin sample_reads_2 46206. 17859.
## 13 Water sample_reads_4 45229. 10217.
## 14 Water sample_reads_3 45538. 10277.
## 15 Water sample_reads_2 45539. 10277.
## 16 Wholebody sample_reads_4 20422. 20434.
## 17 Wholebody sample_reads_3 20523. 20539.
## 18 Wholebody sample_reads_2 20532. 20549.# Plot using ggplot2
# Plot for sediment data
ggplot(Abs_abun_summary, aes(x = Biosample_type, y = Mean, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = Mean-SD, ymax = Mean+SD), width = 0.2, position = position_dodge(0.9)) +
theme_minimal() +
theme(legend.text = element_text(size = 15),
legend.position = "top",
plot.title = element_text(size = 20, face = "bold"),
legend.title = element_text(size = 18),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size = 15),
strip.text = element_text(size = 15),
axis.text.x = element_text(angle = 45, colour = "black", vjust = 1, hjust = 1, size = 15)) +
labs(x = "biosample type", y = "number of read") +
facet_grid( scales = "free")
Calculate relative abundance
#Dataset with Relative abundance
physeq_Nem_rel <- phyloseq::transform_sample_counts(physeq_Nem_5, function(x) x / sum(x) )
rowSums(physeq_Nem_rel@otu_table@.Data)## BD.S.1.M BD.S.2.M BD.S.3.M BD.S.4.B BD.S.5.B
## 1 1 1 1 1
## BD.S.6.B BD.W.1.M BD.W.2.M BD.W.3.M BD.W.4.B
## 1 1 1 1 1
## BD.W.5.B BD.W.6.B Bds1b Bds2b Bds3b
## 1 1 1 1 1
## Bdw1b Bdw2b Bdw3b BioFilm.21.1 BioFilm.21.2
## 1 1 1 1 1
## BioFilm.21.3 CR.S.1.M CR.S.2.M CR.S.3.M CR.S.4.B
## 1 1 1 1 1
## CR.S.5.B CR.S.6.B CR.W.1.M CR.W.2.M CR.W.3.M
## 1 1 1 1 1
## CR.W.4.B CR.W.5.B CR.W.6.B F1s1b F1s2b
## 1 1 1 1 1
## F1s3b F1w1b F1w2b F1w3b F2s1b
## 1 1 1 1 1
## F2s2b F2s3b F2w1b F2w2b F2w3b
## 1 1 1 1 1
## F3s1b F3s2b F3s3b F3w1b F3w2b
## 1 1 1 1 1
## F3w3b F4s1b F4s2b F4s3b F4w1b
## 1 1 1 1 1
## F4w2b F4w3b Fi.10.3.GM Fi.10.4.GB Fi.10.5.GB
## 1 1 1 1 1
## Fi.10.6.GB Fi.12.10.GB Fi.12.11.GB Fi.12.13.GB Fi.12.4.GB
## 1 1 1 1 1
## Fi.12.5.GB Fi.12.6.GB Fi.12.7.GB Fi.12.8.GB Fi.12.9.GB
## 1 1 1 1 1
## Fi.13.1.GM Fi.13.2.GB Fi.13.3.GB Fi.14.2.GM Fi.14.5.GB
## 1 1 1 1 1
## Fi.14.6.GB Fi.18.1.GB Fi.18.2.GB Fi.18.3.GB Fi.19.1.GB
## 1 1 1 1 1
## Fi.20.1.GB Fi.20.2.GB Fi.3.6.GB Fi.3.7.GB Fi.4.4.GB
## 1 1 1 1 1
## Fi.5.3.GB Fi.5.4.GB Fi.5.5.GB Fi.6.4.GB Fi.6.5.GB
## 1 1 1 1 1
## Fi.6.6.GM Fi.6.7.GB Fi.6.8.GB Fi.6.9.GB Fi10.1gb
## 1 1 1 1 1
## Fi10.2gb Fi11.1gb Fi11.2gb Fi12.1gb Fi12.2gb
## 1 1 1 1 1
## Fi12.3gb Fi3.1gb Fi3.2gb Fi3.3gb Fi3.4gb
## 1 1 1 1 1
## Fi3.5gb Fi5.1gb Fi5.2gb Fi6.1gb Fi6.2gb
## 1 1 1 1 1
## Fi6.3gb Fi7.1gb Fi7.2gb Fi7.3gb Fi7.4gb
## 1 1 1 1 1
## Fi7.5gb Fo.17.5.BM Fo.4.1.BM Fo.4.2.BM Fo.4.3.BM
## 1 1 1 1 1
## Fo.4.4.BM Fo.4.5.BM HM.S.1.M HM.S.2.M HM.S.3.M
## 1 1 1 1 1
## HM.S.4.B HM.S.5.B HM.S.6.B HM.W.1.M HM.W.2.M
## 1 1 1 1 1
## HM.W.3.M HM.W.4.B HM.W.5.B HM.W.6.B Hu.18.GM
## 1 1 1 1 1
## Hu.18.SM Hu.19.GM Hu.19.SM Hu.20.GM Hu.20.SM
## 1 1 1 1 1
## Hu.21.GM Hu.21.SM Hu.22.GM Hu.22.SM Hu.23.GM
## 1 1 1 1 1
## Hu.23.SM Hu.25.GM Hu.25.SM Hu.26.GM Hu.26.SM
## 1 1 1 1 1
## Hu.27.GM Hu.27.SM Hu.28.GM Hu.28.SM Hu.29.GM
## 1 1 1 1 1
## Hu.29.SM Hu.30.GM Hu.30.SM Hu.31.GM Hu.31.SM
## 1 1 1 1 1
## Hu.32.GM Hu.33.GM Hu.34.GM Hu.35.GM Hu.36.GM
## 1 1 1 1 1
## Hu.36.SM Hu.37.GM Hu.37.SM Hu.38.GM Hu.38.SM
## 1 1 1 1 1
## Hu.39.GM Hu.39.SM Hu.41.GM Hu.41.SM Hu.43.GM
## 1 1 1 1 1
## Hu.43.SM Hu.44.GM Hu.44.SM Hu.45.GB Hu.45.SB
## 1 1 1 1 1
## Hu.46.GB Hu.46.SB Hu.47.GB Hu.47.SB Hu.48.GB
## 1 1 1 1 1
## Hu.48.SB Hu.49.GB Hu.49.SB Hu.50.GB Hu.50.SB
## 1 1 1 1 1
## Hu.51.GB Hu.51.SB Hu.52.GB Hu.52.SB Hu.53.GB
## 1 1 1 1 1
## Hu.53.SB Hu.54.GB Hu.54.SB Hu.55.GB Hu.55.SB
## 1 1 1 1 1
## Hu.56.GB Hu.56.SB Hu.57.GB Hu.57.SB Hu.58.GB
## 1 1 1 1 1
## Hu.58.SB Hu.59.GB Hu.59.SB Hu.60.GB Hu.60.SB
## 1 1 1 1 1
## Hu.61.GB Hu.61.SB Hu.63.GB Hu.63.SB Hu.64.GB
## 1 1 1 1 1
## Hu.64.SB Hu.66.GB Hu.66.SB Hu.67.GB Hu.67.SB
## 1 1 1 1 1
## Hu.68.GB Hu.68.SB Hu.69.GB Hu.69.SB Hu.70.GB
## 1 1 1 1 1
## Hu.70.SB Hu.71.GB Hu.71.SB Hu10.1.biof.6 Hu10.1gb
## 1 1 1 1 1
## Hu10.1sb Hu10.2.biof.6 Hu10.3.biof.6 Hu11.1.biof.7 Hu11.1gb
## 1 1 1 1 1
## Hu11.1sb Hu11.2.biof.7 Hu11.3.biof.7 Hu12.1.biof.8 Hu12.1gb
## 1 1 1 1 1
## Hu12.1sb Hu12.2.biof.8 Hu12.3.biof.8 Hu13.1.biof.9 Hu13.1gb
## 1 1 1 1 1
## Hu13.1sb Hu13.2.biof.9 Hu13.3.biof.9 Hu14.1gb Hu14.1sb
## 1 1 1 1 1
## Hu15.1.biof.11 Hu15.2.biof.11 Hu15.3.biof.11 Hu15.3gb Hu15.3sb
## 1 1 1 1 1
## Hu16.1.biof.12 Hu16.2.biof.12 Hu16.3.biof.12 Hu16.3gb Hu16.3sb
## 1 1 1 1 1
## Hu17.3gb Hu17.3sb Hu3.1gb Hu3.1sb Hu4.1gb
## 1 1 1 1 1
## Hu4.1sb Hu5.1.biof.3 Hu5.1gb Hu5.1sb Hu5.2.biof.3
## 1 1 1 1 1
## Hu5.3.biof.3 Hu6.1gb Hu6.1sb Hu7.1.biof.4 Hu7.1gb
## 1 1 1 1 1
## Hu7.1sb Hu7.2.biof.4 Hu7.3.biof.4 Hu8.1.biof.5 Hu8.1gb
## 1 1 1 1 1
## Hu8.1sb Hu8.2.biof.5 Hu8.3.biof.5 Hu9.1gb Hu9.1sb
## 1 1 1 1 1
## Lo.2.4.GB Lo.2.5.GB Lo.2.7.GB Lo.2.8.GB Lo.2.9.GB
## 1 1 1 1 1
## Lo.5.1.GM Lo.5.4.GB Lo.5.5.GB Lo.9.2.GB Lo2.1gb
## 1 1 1 1 1
## Lo2.2gb Lo2.3gb Lo4.1gb Lo4.2gb Lo8.1gb
## 1 1 1 1 1
## Lo8.2gb Lo8.3gb Lo8.4gb Lo8.5gb NP.S.1.M
## 1 1 1 1 1
## NP.S.2.M NP.S.3.M NP.S.4.B NP.S.5.B NP.S.6.B
## 1 1 1 1 1
## NP.W.1.M NP.W.2.M NP.W.3.M NP.W.4.B NP.W.5.B
## 1 1 1 1 1
## NP.W.6.B Oy.13.7.GB Oy.15.4.GM Oy.15.5.GM Oy.15.6.GB
## 1 1 1 1 1
## Oy.15.8.GB R1s1b R1s2b R1s3b R1w1b
## 1 1 1 1 1
## R1w2b R1w3b R2s1b R2s2b R2s3b
## 1 1 1 1 1
## R2w1b R2w2b R2w3b R3s1b R3s2b
## 1 1 1 1 1
## R3s3b R3w1b R3w2b R3w3b R4s1b
## 1 1 1 1 1
## R4s2b R4s3b R4w1b R4w2b R4w3b
## 1 1 1 1 1
## Rivs1b Rivs2b Rivs3b Rivw1b Rivw2b
## 1 1 1 1 1
## Rivw3b Sq.11.1.GM Sq.11.2.GM Sq.11.4.GB Sq.11.5.GB
## 1 1 1 1 1
## U1s1b U1s2b U1s3b U1w1b U1w2b
## 1 1 1 1 1
## U1w3b U2s1b U2s2b U2s3b U2w1b
## 1 1 1 1 1
## U2w2b U2w3b U3s1b U3s2b U3s3b
## 1 1 1 1 1
## U3w1b U3w2b U3w3b U4s1b U4s2b
## 1 1 1 1 1
## U4s3b U4w1b U4w2b U4w3b
## 1 1 1 1Rarefy at minimum sample size (for Alpha diversity analysis)
ps_Nem_rarefy <- phyloseq::rarefy_even_depth(physeq_Nem_5, rngseed = TRUE, replace = FALSE, sample.size = min(rowSums(physeq_Nem_5@otu_table)))
#Check rarefied curves
p01 <- ggrare(ps_Nem_rarefy, step = 100, color = "Biosample_type", se = TRUE)## rarefying sample BD.S.1.M
## rarefying sample BD.S.2.M
## rarefying sample BD.S.3.M
## rarefying sample BD.S.4.B
## rarefying sample BD.S.5.B
## rarefying sample BD.S.6.B
## rarefying sample BD.W.1.M
## rarefying sample BD.W.2.M
## rarefying sample BD.W.3.M
## rarefying sample BD.W.4.B
## rarefying sample BD.W.5.B
## rarefying sample BD.W.6.B
## rarefying sample Bds1b
## rarefying sample Bds2b
## rarefying sample Bds3b
## rarefying sample Bdw1b
## rarefying sample Bdw2b
## rarefying sample Bdw3b
## rarefying sample BioFilm.21.1
## rarefying sample BioFilm.21.2
## rarefying sample BioFilm.21.3
## rarefying sample CR.S.1.M
## rarefying sample CR.S.2.M
## rarefying sample CR.S.3.M
## rarefying sample CR.S.4.B
## rarefying sample CR.S.5.B
## rarefying sample CR.S.6.B
## rarefying sample CR.W.1.M
## rarefying sample CR.W.2.M
## rarefying sample CR.W.3.M
## rarefying sample CR.W.4.B
## rarefying sample CR.W.5.B
## rarefying sample CR.W.6.B
## rarefying sample F1s1b
## rarefying sample F1s2b
## rarefying sample F1s3b
## rarefying sample F1w1b
## rarefying sample F1w2b
## rarefying sample F1w3b
## rarefying sample F2s1b
## rarefying sample F2s2b
## rarefying sample F2s3b
## rarefying sample F2w1b
## rarefying sample F2w2b
## rarefying sample F2w3b
## rarefying sample F3s1b
## rarefying sample F3s2b
## rarefying sample F3s3b
## rarefying sample F3w1b
## rarefying sample F3w2b
## rarefying sample F3w3b
## rarefying sample F4s1b
## rarefying sample F4s2b
## rarefying sample F4s3b
## rarefying sample F4w1b
## rarefying sample F4w2b
## rarefying sample F4w3b
## rarefying sample Fi.10.3.GM
## rarefying sample Fi.10.4.GB
## rarefying sample Fi.10.5.GB
## rarefying sample Fi.10.6.GB
## rarefying sample Fi.12.10.GB
## rarefying sample Fi.12.11.GB
## rarefying sample Fi.12.13.GB
## rarefying sample Fi.12.4.GB
## rarefying sample Fi.12.5.GB
## rarefying sample Fi.12.6.GB
## rarefying sample Fi.12.7.GB
## rarefying sample Fi.12.8.GB
## rarefying sample Fi.12.9.GB
## rarefying sample Fi.13.1.GM
## rarefying sample Fi.13.2.GB
## rarefying sample Fi.13.3.GB
## rarefying sample Fi.14.2.GM
## rarefying sample Fi.14.5.GB
## rarefying sample Fi.14.6.GB
## rarefying sample Fi.18.1.GB
## rarefying sample Fi.18.2.GB
## rarefying sample Fi.18.3.GB
## rarefying sample Fi.19.1.GB
## rarefying sample Fi.20.1.GB
## rarefying sample Fi.20.2.GB
## rarefying sample Fi.3.6.GB
## rarefying sample Fi.3.7.GB
## rarefying sample Fi.4.4.GB
## rarefying sample Fi.5.3.GB
## rarefying sample Fi.5.4.GB
## rarefying sample Fi.5.5.GB
## rarefying sample Fi.6.4.GB
## rarefying sample Fi.6.5.GB
## rarefying sample Fi.6.6.GM
## rarefying sample Fi.6.7.GB
## rarefying sample Fi.6.8.GB
## rarefying sample Fi.6.9.GB
## rarefying sample Fi10.1gb
## rarefying sample Fi10.2gb
## rarefying sample Fi11.1gb
## rarefying sample Fi11.2gb
## rarefying sample Fi12.1gb
## rarefying sample Fi12.2gb
## rarefying sample Fi12.3gb
## rarefying sample Fi3.1gb
## rarefying sample Fi3.2gb
## rarefying sample Fi3.3gb
## rarefying sample Fi3.4gb
## rarefying sample Fi3.5gb
## rarefying sample Fi5.1gb## rarefying sample Fi5.2gb
## rarefying sample Fi6.1gb
## rarefying sample Fi6.2gb
## rarefying sample Fi6.3gb
## rarefying sample Fi7.1gb
## rarefying sample Fi7.2gb
## rarefying sample Fi7.3gb
## rarefying sample Fi7.4gb
## rarefying sample Fi7.5gb
## rarefying sample Fo.17.5.BM
## rarefying sample Fo.4.1.BM
## rarefying sample Fo.4.2.BM
## rarefying sample Fo.4.3.BM
## rarefying sample Fo.4.4.BM
## rarefying sample Fo.4.5.BM
## rarefying sample HM.S.1.M
## rarefying sample HM.S.2.M
## rarefying sample HM.S.3.M
## rarefying sample HM.S.4.B
## rarefying sample HM.S.5.B
## rarefying sample HM.S.6.B
## rarefying sample HM.W.1.M
## rarefying sample HM.W.2.M
## rarefying sample HM.W.3.M
## rarefying sample HM.W.4.B
## rarefying sample HM.W.5.B
## rarefying sample HM.W.6.B
## rarefying sample Hu.18.GM
## rarefying sample Hu.18.SM
## rarefying sample Hu.19.GM
## rarefying sample Hu.19.SM
## rarefying sample Hu.20.GM
## rarefying sample Hu.20.SM
## rarefying sample Hu.21.GM
## rarefying sample Hu.21.SM
## rarefying sample Hu.22.GM
## rarefying sample Hu.22.SM
## rarefying sample Hu.23.GM
## rarefying sample Hu.23.SM
## rarefying sample Hu.25.GM
## rarefying sample Hu.25.SM
## rarefying sample Hu.26.GM
## rarefying sample Hu.26.SM
## rarefying sample Hu.27.GM
## rarefying sample Hu.27.SM
## rarefying sample Hu.28.GM
## rarefying sample Hu.28.SM
## rarefying sample Hu.29.GM
## rarefying sample Hu.29.SM
## rarefying sample Hu.30.GM
## rarefying sample Hu.30.SM
## rarefying sample Hu.31.GM
## rarefying sample Hu.31.SM
## rarefying sample Hu.32.GM
## rarefying sample Hu.33.GM
## rarefying sample Hu.34.GM
## rarefying sample Hu.35.GM
## rarefying sample Hu.36.GM
## rarefying sample Hu.36.SM
## rarefying sample Hu.37.GM
## rarefying sample Hu.37.SM
## rarefying sample Hu.38.GM
## rarefying sample Hu.38.SM
## rarefying sample Hu.39.GM
## rarefying sample Hu.39.SM
## rarefying sample Hu.41.GM
## rarefying sample Hu.41.SM
## rarefying sample Hu.43.GM
## rarefying sample Hu.43.SM
## rarefying sample Hu.44.GM
## rarefying sample Hu.44.SM
## rarefying sample Hu.45.GB
## rarefying sample Hu.45.SB
## rarefying sample Hu.46.GB
## rarefying sample Hu.46.SB
## rarefying sample Hu.47.GB
## rarefying sample Hu.47.SB
## rarefying sample Hu.48.GB
## rarefying sample Hu.48.SB
## rarefying sample Hu.49.GB
## rarefying sample Hu.49.SB
## rarefying sample Hu.50.GB
## rarefying sample Hu.50.SB
## rarefying sample Hu.51.GB
## rarefying sample Hu.51.SB
## rarefying sample Hu.52.GB
## rarefying sample Hu.52.SB
## rarefying sample Hu.53.GB
## rarefying sample Hu.53.SB
## rarefying sample Hu.54.GB
## rarefying sample Hu.54.SB
## rarefying sample Hu.55.GB
## rarefying sample Hu.55.SB
## rarefying sample Hu.56.GB
## rarefying sample Hu.56.SB
## rarefying sample Hu.57.GB
## rarefying sample Hu.57.SB
## rarefying sample Hu.58.GB
## rarefying sample Hu.58.SB
## rarefying sample Hu.59.GB
## rarefying sample Hu.59.SB
## rarefying sample Hu.60.GB
## rarefying sample Hu.60.SB
## rarefying sample Hu.61.GB
## rarefying sample Hu.61.SB
## rarefying sample Hu.63.GB
## rarefying sample Hu.63.SB
## rarefying sample Hu.64.GB
## rarefying sample Hu.64.SB
## rarefying sample Hu.66.GB
## rarefying sample Hu.66.SB
## rarefying sample Hu.67.GB
## rarefying sample Hu.67.SB
## rarefying sample Hu.68.GB
## rarefying sample Hu.68.SB
## rarefying sample Hu.69.GB
## rarefying sample Hu.69.SB
## rarefying sample Hu.70.GB
## rarefying sample Hu.70.SB
## rarefying sample Hu.71.GB
## rarefying sample Hu.71.SB
## rarefying sample Hu10.1.biof.6
## rarefying sample Hu10.1gb
## rarefying sample Hu10.1sb
## rarefying sample Hu10.2.biof.6
## rarefying sample Hu10.3.biof.6
## rarefying sample Hu11.1.biof.7
## rarefying sample Hu11.1gb
## rarefying sample Hu11.1sb
## rarefying sample Hu11.2.biof.7
## rarefying sample Hu11.3.biof.7
## rarefying sample Hu12.1.biof.8
## rarefying sample Hu12.1gb
## rarefying sample Hu12.1sb
## rarefying sample Hu12.2.biof.8
## rarefying sample Hu12.3.biof.8
## rarefying sample Hu13.1.biof.9
## rarefying sample Hu13.1gb
## rarefying sample Hu13.1sb
## rarefying sample Hu13.2.biof.9
## rarefying sample Hu13.3.biof.9
## rarefying sample Hu14.1gb
## rarefying sample Hu14.1sb
## rarefying sample Hu15.1.biof.11
## rarefying sample Hu15.2.biof.11
## rarefying sample Hu15.3.biof.11
## rarefying sample Hu15.3gb
## rarefying sample Hu15.3sb
## rarefying sample Hu16.1.biof.12
## rarefying sample Hu16.2.biof.12
## rarefying sample Hu16.3.biof.12
## rarefying sample Hu16.3gb
## rarefying sample Hu16.3sb
## rarefying sample Hu17.3gb
## rarefying sample Hu17.3sb
## rarefying sample Hu3.1gb
## rarefying sample Hu3.1sb
## rarefying sample Hu4.1gb
## rarefying sample Hu4.1sb
## rarefying sample Hu5.1.biof.3
## rarefying sample Hu5.1gb
## rarefying sample Hu5.1sb
## rarefying sample Hu5.2.biof.3
## rarefying sample Hu5.3.biof.3
## rarefying sample Hu6.1gb
## rarefying sample Hu6.1sb
## rarefying sample Hu7.1.biof.4
## rarefying sample Hu7.1gb
## rarefying sample Hu7.1sb
## rarefying sample Hu7.2.biof.4
## rarefying sample Hu7.3.biof.4
## rarefying sample Hu8.1.biof.5
## rarefying sample Hu8.1gb
## rarefying sample Hu8.1sb
## rarefying sample Hu8.2.biof.5
## rarefying sample Hu8.3.biof.5
## rarefying sample Hu9.1gb
## rarefying sample Hu9.1sb
## rarefying sample Lo.2.4.GB
## rarefying sample Lo.2.5.GB
## rarefying sample Lo.2.7.GB
## rarefying sample Lo.2.8.GB## rarefying sample Lo.2.9.GB
## rarefying sample Lo.5.1.GM
## rarefying sample Lo.5.4.GB
## rarefying sample Lo.5.5.GB
## rarefying sample Lo.9.2.GB
## rarefying sample Lo2.1gb
## rarefying sample Lo2.2gb
## rarefying sample Lo2.3gb
## rarefying sample Lo4.1gb
## rarefying sample Lo4.2gb
## rarefying sample Lo8.1gb## rarefying sample Lo8.2gb
## rarefying sample Lo8.3gb
## rarefying sample Lo8.4gb
## rarefying sample Lo8.5gb
## rarefying sample NP.S.1.M
## rarefying sample NP.S.2.M
## rarefying sample NP.S.3.M
## rarefying sample NP.S.4.B
## rarefying sample NP.S.5.B
## rarefying sample NP.S.6.B
## rarefying sample NP.W.1.M
## rarefying sample NP.W.2.M
## rarefying sample NP.W.3.M
## rarefying sample NP.W.4.B
## rarefying sample NP.W.5.B
## rarefying sample NP.W.6.B
## rarefying sample Oy.13.7.GB
## rarefying sample Oy.15.4.GM
## rarefying sample Oy.15.5.GM
## rarefying sample Oy.15.6.GB
## rarefying sample Oy.15.8.GB
## rarefying sample R1s1b
## rarefying sample R1s2b
## rarefying sample R1s3b
## rarefying sample R1w1b
## rarefying sample R1w2b
## rarefying sample R1w3b
## rarefying sample R2s1b
## rarefying sample R2s2b
## rarefying sample R2s3b
## rarefying sample R2w1b
## rarefying sample R2w2b
## rarefying sample R2w3b
## rarefying sample R3s1b
## rarefying sample R3s2b
## rarefying sample R3s3b
## rarefying sample R3w1b
## rarefying sample R3w2b
## rarefying sample R3w3b
## rarefying sample R4s1b
## rarefying sample R4s2b
## rarefying sample R4s3b
## rarefying sample R4w1b
## rarefying sample R4w2b
## rarefying sample R4w3b
## rarefying sample Rivs1b
## rarefying sample Rivs2b
## rarefying sample Rivs3b
## rarefying sample Rivw1b
## rarefying sample Rivw2b
## rarefying sample Rivw3b
## rarefying sample Sq.11.1.GM
## rarefying sample Sq.11.2.GM
## rarefying sample Sq.11.4.GB
## rarefying sample Sq.11.5.GB
## rarefying sample U1s1b
## rarefying sample U1s2b
## rarefying sample U1s3b
## rarefying sample U1w1b
## rarefying sample U1w2b
## rarefying sample U1w3b
## rarefying sample U2s1b
## rarefying sample U2s2b
## rarefying sample U2s3b
## rarefying sample U2w1b
## rarefying sample U2w2b
## rarefying sample U2w3b
## rarefying sample U3s1b
## rarefying sample U3s2b
## rarefying sample U3s3b
## rarefying sample U3w1b
## rarefying sample U3w2b
## rarefying sample U3w3b
## rarefying sample U4s1b
## rarefying sample U4s2b
## rarefying sample U4s3b
## rarefying sample U4w1b
## rarefying sample U4w2b
## rarefying sample U4w3b
#check Zhang & Huang's coverage estimator
library(entropart)
Coverage_estimator_rarefied <- Coverage(ps_Nem_rarefy@otu_table, Estimator = "Best", Level = NULL, CheckArguments = TRUE)
Coverage_estimator_rarefied## ZhangHuang
## 0.9772663#Remove unnecessary samples
#ps_Nem_rarefy <- subset_samples(ps_Nem_rarefy, !(sample_names(ps_Nem_rarefy) %in% c("R1w2b", "JCA.U2S2B", "U1w2b")))
#Check sample size
dim(ps_Nem_rarefy@otu_table)## [1] 379 71157