Phyloseq microbiome

1 Data

Import raw data and assign sample key:

# extend map_corrected.txt with Diet and Flora
# setwd('~/DATA/Data_Laura_16S_2/core_diversity_e4753')
map_corrected <- read.csv("../map3_corrected.txt", sep = "\t", row.names = 1)
knitr::kable(map_corrected) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	BarcodeSequence	LinkerPrimerSequence	FileInput	Patient	Type	Group
P7.Nose	NA	NA	P7-Nose_merged.fasta	P7	Nose	P7-P15.Nose
P8.Nose	NA	NA	P8-Nose_merged.fasta	P8	Nose	P7-P15.Nose
P9.Nose	NA	NA	P9-Nose_merged.fasta	P9	Nose	P7-P15.Nose
P10.Nose	NA	NA	P10-Nose_merged.fasta	P10	Nose	P7-P15.Nose
P11.Nose	NA	NA	P11-Nose_merged.fasta	P11	Nose	P7-P15.Nose
P12.Nose	NA	NA	P12-Nose_merged.fasta	P12	Nose	P7-P15.Nose
P13.Nose	NA	NA	P13-Nose_merged.fasta	P13	Nose	P7-P15.Nose
P14.Nose	NA	NA	P14-Nose_merged.fasta	P14	Nose	P7-P15.Nose
P15.Nose	NA	NA	P15-Nose_merged.fasta	P15	Nose	P7-P15.Nose
P16.Foot	NA	NA	P16-Foot_merged.fasta	P16	Foot	P16-P20.Foot
P16.Nose	NA	NA	P16-Nose_merged.fasta	P16	Nose	P16-P20.Nose
P17.Foot	NA	NA	P17-Foot_merged.fasta	P17	Foot	P16-P20.Foot
P17.Nose	NA	NA	P17-Nose_merged.fasta	P17	Nose	P16-P20.Nose
P18.Foot	NA	NA	P18-Foot_merged.fasta	P18	Foot	P16-P20.Foot
P18.Nose	NA	NA	P18-Nose_merged.fasta	P18	Nose	P16-P20.Nose
P19.Foot	NA	NA	P19-Foot_merged.fasta	P19	Foot	P16-P20.Foot
P19.Nose	NA	NA	P19-Nose_merged.fasta	P19	Nose	P16-P20.Nose
P20.Foot	NA	NA	P20-Foot_merged.fasta	P20	Foot	P16-P20.Foot
P20.Nose	NA	NA	P20-Nose_merged.fasta	P20	Nose	P16-P20.Nose
AH.LH	NA	NA	AH-LH_merged.fasta	AH	LH	AH-XN.LH
AH.NLH	NA	NA	AH-NLH_merged.fasta	AH	NLH	AH-XN.NLH
AH.Nose	NA	NA	AH-Nose_merged.fasta	AH	Nose	AH-XN.Nose
AL.LH	NA	NA	AL-LH_merged.fasta	AL	LH	AH-XN.LH
AL.NLH	NA	NA	AL-NLH_merged.fasta	AL	NLH	AH-XN.NLH
AL.Nose	NA	NA	AL-Nose_merged.fasta	AL	Nose	AH-XN.Nose
CB.LH	NA	NA	CB-LH_merged.fasta	CB	LH	AH-XN.LH
CB.NLH	NA	NA	CB-NLH_merged.fasta	CB	NLH	AH-XN.NLH
CB.Nose	NA	NA	CB-Nose_merged.fasta	CB	Nose	AH-XN.Nose
HR.LH	NA	NA	HR-LH_merged.fasta	HR	LH	AH-XN.LH
HR.NLH	NA	NA	HR-NLH_merged.fasta	HR	NLH	AH-XN.NLH
HR.Nose	NA	NA	HR-Nose_merged.fasta	HR	Nose	AH-XN.Nose
KK.LH	NA	NA	KK-LH_merged.fasta	KK	LH	AH-XN.LH
KK.NLH	NA	NA	KK-NLH_merged.fasta	KK	NLH	AH-XN.NLH
KK.Nose	NA	NA	KK-Nose_merged.fasta	KK	Nose	AH-XN.Nose
MC.LH	NA	NA	MC-LH_merged.fasta	MC	LH	AH-XN.LH
MC.NLH	NA	NA	MC-NLH_merged.fasta	MC	NLH	AH-XN.NLH
MC.Nose	NA	NA	MC-Nose_merged.fasta	MC	Nose	AH-XN.Nose
MR.LH	NA	NA	MR-LH_merged.fasta	MR	LH	AH-XN.LH
MR.NLH	NA	NA	MR-NLH_merged.fasta	MR	NLH	AH-XN.NLH
MR.Nose	NA	NA	MR-Nose_merged.fasta	MR	Nose	AH-XN.Nose
PT2.LH	NA	NA	PT2-LH_merged.fasta	PT2	LH	AH-XN.LH
PT2.NLH	NA	NA	PT2-NLH_merged.fasta	PT2	NLH	AH-XN.NLH
PT2.Nose	NA	NA	PT2-Nose_merged.fasta	PT2	Nose	AH-XN.Nose
RP.LH	NA	NA	RP-LH_merged.fasta	RP	LH	AH-XN.LH
RP.NLH	NA	NA	RP-NLH_merged.fasta	RP	NLH	AH-XN.NLH
RP.Nose	NA	NA	RP-Nose_merged.fasta	RP	Nose	AH-XN.Nose
SA.LH	NA	NA	SA-LH_merged.fasta	SA	LH	AH-XN.LH
SA.NLH	NA	NA	SA-NLH_merged.fasta	SA	NLH	AH-XN.NLH
SA.Nose	NA	NA	SA-Nose_merged.fasta	SA	Nose	AH-XN.Nose
XN.LH	NA	NA	XN-LH_merged.fasta	XN	LH	AH-XN.LH
XN.NLH	NA	NA	XN-NLH_merged.fasta	XN	NLH	AH-XN.NLH
XN.Nose	NA	NA	XN-Nose_merged.fasta	XN	Nose	AH-XN.Nose

2 Prerequisites to be installed

R : https://pbil.univ-lyon1.fr/CRAN/
R studio : https://www.rstudio.com/products/rstudio/download/#download

install.packages("dplyr")     # To manipulate dataframes
install.packages("readxl")    # To read Excel files into R
install.packages("ggplot2")   # for high quality graphics
install.packages("heatmaply")
source("https://bioconductor.org/biocLite.R")
biocLite("phyloseq")

library("readxl")  # necessary to import the data from Excel file
library("ggplot2")  # graphics
library("picante")
library("microbiome")  # data analysis and visualisation
library("phyloseq")  # also the basis of data object. Data analysis and visualisation
library("ggpubr")  # publication quality figures, based on ggplot2
library("dplyr")  # data handling, filter and reformat data frames
library("RColorBrewer")  # nice color options
library("heatmaply")
library(vegan)
library(gplots)

3 Read the data and create phyloseq objects

Three tables are needed

OTU
Taxonomy
Samples

# Change your working directory to where the files are located
ps.ng.tax <- import_biom("./table_even33778.biom", "../clustering/rep_set.tre")
sample <- read.csv("../map3_corrected.txt", sep = "\t", row.names = 1)
SAM = sample_data(sample, errorIfNULL = T)
# rownames(SAM) <-
# c('P7.Nose','P8.Nose','P9.Nose','P10.Nose','P11.Nose','P12.Nose','P13.Nose','P14.Nose','P15.Nose',
# 'P16.Foot','P16.Nose','P17.Foot','P17.Nose','P18.Foot','P18.Nose','P19.Foot','P19.Nose','P20.Foot','P20.Nose','AH.LH','AH.NLH','AH.Nose','AL.LH','AL.NLH','AL.Nose','CB.LH','CB.NLH','CB.Nose','HR.LH','HR.NLH','HR.Nose','KK.LH','KK.NLH','KK.Nose','MC.LH','MC.NLH','MC.Nose','MR.LH','MR.NLH','MR.Nose','PT2.LH','PT2.NLH','PT2.Nose','RP.LH','RP.NLH','RP.Nose','SA.LH','SA.NLH','SA.Nose','XN.LH','XN.NLH','XN.Nose')
ps.ng.tax <- merge_phyloseq(ps.ng.tax, SAM)
print(ps.ng.tax)

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 15446 taxa and 52 samples ]
sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 15446 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 15446 tips and 15402 internal nodes ]

colnames(tax_table(ps.ng.tax)) <- c("Domain", "Phylum", "Class", "Order", "Family",
    "Genus", "Species")
saveRDS(ps.ng.tax, "./ps.ng.tax.rds")

Visualize data

sample_names(ps.ng.tax)

 [1] "P10.Nose" "RP.NLH"   "P14.Nose" "SA.LH"    "AH.NLH"   "P20.Nose"
 [7] "KK.Nose"  "AL.NLH"   "P16.Nose" "RP.LH"    "P13.Nose" "MR.LH"   
[13] "SA.NLH"   "MR.NLH"   "AL.Nose"  "MC.NLH"   "P15.Nose" "P8.Nose" 
[19] "MC.Nose"  "P17.Nose" "P16.Foot" "P7.Nose"  "P12.Nose" "P18.Nose"
[25] "HR.NLH"   "P11.Nose" "P9.Nose"  "XN.LH"    "P19.Nose" "XN.Nose" 
[31] "PT2.LH"   "SA.Nose"  "PT2.Nose" "CB.LH"    "P17.Foot" "HR.LH"   
[37] "P19.Foot" "AH.LH"    "AH.Nose"  "CB.NLH"   "P18.Foot" "RP.Nose" 
[43] "CB.Nose"  "HR.Nose"  "XN.NLH"   "MR.Nose"  "KK.NLH"   "AL.LH"   
[49] "KK.LH"    "MC.LH"    "P20.Foot" "PT2.NLH"

rank_names(ps.ng.tax)

[1] "Domain"  "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

sample_variables(ps.ng.tax)

[1] "BarcodeSequence"      "LinkerPrimerSequence" "FileInput"           
[4] "Patient"              "Type"                 "Group"

Normalize number of reads in each sample using median sequencing depth.

# RAREFACTION set.seed(9242) # This will help in reproducing the filtering and
# nomalisation.  ps.ng.tax <- rarefy_even_depth(ps.ng.tax, sample.size = 42369)
# total <- 42369

# NORMALIZE number of reads in each sample using median sequencing depth.
total = median(sample_sums(ps.ng.tax))
# > total [1] 42369
standf = function(x, t = total) round(t * (x/sum(x)))
ps.ng.tax = transform_sample_counts(ps.ng.tax, standf)
ps.ng.tax_rel <- microbiome::transform(ps.ng.tax, "compositional")

saveRDS(ps.ng.tax, "./ps.ng.tax.rds")
hmp.meta <- meta(ps.ng.tax)
hmp.meta$sam_name <- rownames(hmp.meta)

4 Heatmaps

# MOVE_FROM_ABOVE: The number of reads used for normalization is **`r
# sprintf('%.0f', total)`**.  A basic heatmap using the default parameters.
# plot_heatmap(ps.ng.tax, method = 'NMDS', distance = 'bray') NOTE that giving
# the correct OTU numbers in the text (1%, 0.5%, ...)!!!

We consider the most abundant OTUs for heatmaps. For example one can only take OTUs that represent at least 1% of reads in at least one sample. Remember we normalized all the sampples to median number of reads (total). We are left with only 168 OTUS which makes the reading much more easy.

# Custom function to plot a heatmap with the specified sample order
# plot_heatmap_custom <- function(ps, sample_order, method = 'NMDS', distance =
# 'bray') {
ps.ng.tax_abund <- phyloseq::filter_taxa(ps.ng.tax, function(x) sum(x > total * 0.01) >
    0, TRUE)
kable(otu_table(ps.ng.tax_abund)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	P10.Nose	RP.NLH	P14.Nose	SA.LH	AH.NLH	P20.Nose	KK.Nose	AL.NLH	P16.Nose	RP.LH	P13.Nose	MR.LH	SA.NLH	MR.NLH	AL.Nose	MC.NLH	P15.Nose	P8.Nose	MC.Nose	P17.Nose	P16.Foot	P7.Nose	P12.Nose	P18.Nose	HR.NLH	P11.Nose	P9.Nose	XN.LH	P19.Nose	XN.Nose	PT2.LH	SA.Nose	PT2.Nose	CB.LH	P17.Foot	HR.LH	P19.Foot	AH.LH	AH.Nose	CB.NLH	P18.Foot	RP.Nose	CB.Nose	HR.Nose	XN.NLH	MR.Nose	KK.NLH	AL.LH	KK.LH	MC.LH	P20.Foot	PT2.NLH
CP007016.16761.18305	287	0	0	353	0	0	0	0	0	0	0	0	119	61	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	6	0	0	0	8	0	0
AYSG01000002.292.2076	61	57	177	503	8	499	25	69	0	72	0	6	543	25	11	77	146	4	0	0	2	555	2	626	55	4	0	6	9	0	14	12	3	7	0	3	0	12	7	31	1	0	7	0	0	8	0	0	11	7	0	0
HM254595.1.1337	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	5	0	662	0
AB109770.615.2139	1108	509	15	198	8	1188	2	5	54	8	814	31	867	1	623	706	227	4	3989	1577	1976	5	195	232	35	114	437	10	494	330	62	348	273	30	99	33	31	5	2	7	3	21	221	73	6	246	5	0	7	83	34	0
New.ReferenceOTU156	0	123	0	38	0	0	7	521	0	13	0	29	227	126	0	401	0	0	0	0	0	0	0	0	99	0	0	3	0	0	22	0	8	1	0	27	0	7	1	6	0	0	9	0	0	8	0	0	9	53	0	0
FTOF01000024.348.1868	21	208	79	1	477	174	3	0	2	0	0	735	0	152	0	0	30	2	0	0	4	0	6	0	899	0	0	0	0	0	0	0	0	17	0	101	0	211	2	5	0	0	0	0	0	21	0	0	21	0	0	0
ACLJ02000003.1110917.1112424	34	0	0	0	7	0	2	0	4	14	0	0	0	0	0	0	0	0	1	0	1	0	0	0	1150	0	0	0	4	0	0	0	0	12	0	33	2	15	2	0	0	0	0	1	0	1	0	0	15	0	0	0
FJVG01000020.441229.442716	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	358	0	0	0	0	0	0	0	0	0	0	6	0	2	0	0	0	0	0	1	0	0	0	0	1	0	0	0
GQ006374.1.1344	7	15	0	0	50	57	1	0	0	0	21	130	0	51	0	0	1	1	0	0	1	1	3	4	6	0	0	0	1	0	0	0	0	5	47	8	1	15	0	1	19	0	0	4	2	6	410	0	86	0	317	0
CP009215.1123627.1125129	83	834	0	0	417	1142	28	0	864	18	239	765	0	395	0	19	20	0	87	98	9	14	31	227	346	0	0	5	7	1	21	0	1	76	556	59	12	94	2	67	38	0	19	7	65	13	4291	2	706	0	1875	4
CBYN010000177.386.1880	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	2	1	0	3	0	0	0	5865	0	0	0	0	1	0	0	0	0	0	241	0	1	0	3	0	0	0	0	0	0	1	0	0	0	2	0
CP017639.498785.500312	0	1	0	7	0	352	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	9	0	162	0	0	0	11	0	0	0	0	5	0	7	1	8	0	0	0	0	0	0	0	0	2	0	0	0	1	0
CR931997.108684.110210	0	0	0	0	10	0	7	0	3	0	0	8	7	49	0	0	0	0	0	0	85	0	0	0	1	0	0	0	2	0	0	0	0	33	730	8	2006	0	0	9	1	0	0	0	27	10	6	0	39	0	22	0
CP002857.605318.606853	0	0	0	0	0	0	0	50	0	0	0	0	0	0	3	0	0	0	0	0	415	0	0	0	12	0	0	0	0	0	0	0	0	0	0	0	223	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
LSCY01000011.94162.95670	31	84	42	28	2	0	0	16	2276	20	0	6	0	0	9	0	4	125	0	3	1	495	67	0	0	1	957	0	21	0	9	0	0	0	1	0	3	0	37	4	0	201	0	0	0	7	0	6	0	0	0	0
LBNF01000094.61.1632	0	105	0	0	0	0	0	124	0	21	0	0	0	0	0	0	0	0	0	0	0	4	0	0	0	0	0	1	5	0	346	0	30	0	0	0	0	0	49	0	0	0	0	0	0	0	0	0	0	0	0	0
LTIM01000071.135.1665	16	1	4	27	352	38	66	172	0	5	1	900	45	485	2	0	65	0	1	0	109	0	2	1	36	0	24	4	0	1	12	0	1	190	4015	15	4944	126	2	77	0	0	0	0	0	17	69	0	1025	21	1988	0
ADER01000013.440.1977	0	0	0	1066	0	0	0	57	0	0	0	0	457	0	6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	13	0	4	0	0	0	0	0	0	0	0	0	0	0	0	0	13	0	0	0	0
CP002280.1382146.1383660	0	0	0	171	0	424	48	0	0	16	0	0	61	2	0	0	203	0	0	0	0	0	5	0	9	73	0	0	1	0	0	0	0	0	0	14	0	0	0	0	0	0	0	0	43	6	0	4	0	0	0	0
ACSB01000014.3447.4963	1	0	249	112	0	305	52	0	0	66	1344	0	233	0	0	0	0	0	0	0	0	6	330	320	38	18	0	0	10	0	4	0	0	0	0	1206	0	0	0	0	0	0	0	188	0	0	0	0	0	0	0	1
New.ReferenceOTU98	1	8	0	0	4	0	6	11	7	1	0	0	4	14	0	3	2	0	0	0	54	1	0	26	2	1	0	0	2	0	100	15	4	0	0	1	0	1	0	2	0	2	2	0	5	0	0	0	25	0	0	578
CZJU01000044.3819.5330	34	234	0	7	74	3	37	314	430	122	0	21	59	320	0	79	13	0	0	0	1330	26	17	48	29	39	2	4	6	2	1360	273	98	20	0	48	29	50	4	19	16	125	44	0	34	24	0	16	269	0	0	31640
KY777463.1.1542	5	0	35	126	22	34	11	517	65	2362	8	54	153	536	2	1728	0	4	1	43	5	1	18	1	5433	2	8	4	7	15	34	2	1	0	4	2239	6	0	6	3	26	628	16	0	12	11	3	6	40	1342	3	3
ABYQ02000007.68757.70262	37	11	13	12	27	4	5	8	2	1	9	19	24	22	21	0	26	8	0	72	26	105	15	2	7	229	24	4	56	12	14	21	103	1	329	20	326	19	30	2	477	30	34	4	3	8	98	3	29	1	143	2
AJ439343.1.1503	0	0	2	0	4	1	946	0	3544	0	0	0	0	0	0	0	4	314	0	101	0	6712	0	0	0	2	4340	0	2232	420	0	0	1	11	0	1	2	0	0	54	0	1	12215	4	23	0	0	6	0	0	0	1
HK556006.10.1507	0	0	229	36	0	81	5	0	19	0	25	4	51	0	0	0	0	0	1	0	2676	0	0	0	14	1	1045	0	0	0	0	0	0	0	0	0	36	0	0	0	0	0	0	0	0	0	0	0	10	0	1	0
ATVG01000023.2522.4022	0	0	0	62	0	0	1	0	0	0	2	209	30	509	0	0	0	0	0	0	1	0	13	0	1	0	0	7	0	0	0	0	0	7	0	0	0	0	0	28	0	0	0	0	14	11	5	0	56	0	1211	0
AXMA01000012.53509.55015	2744	1053	4552	1745	728	5952	162	1545	225	208	5216	4343	2245	2222	2698	2	6545	15281	3342	3782	696	8757	3313	437	174	9889	5805	69	7098	14171	72	4045	6584	219	1814	2996	2623	531	2657	125	14296	7104	294	646	261	14285	648	256	1932	345	5196	6
CP009246.1261131.1262668	212	3	42	8	111	20	3	6	3	1	19	154	10	60	3	1	28	46	12	921	16	30	77	2	5	2210	177	2	24	58	57	66	1883	4	3188	11	1099	32	42	2	3217	18	0	1	22	57	684	1	50	1	708	0
LWPR01000046.100.1629	0	0	1	0	20	0	0	0	0	0	36	3	0	3	1	0	2	1	0	3	0	0	0	0	0	13	2	0	0	0	0	1	4	0	399	0	5	3	0	0	14	1	0	0	0	2	7	0	0	0	527	0
CP009312.1832900.1834429	221	228	819	349	135	31	0	233	189	28	220	677	78	1	674	15	67	5	2401	341	3	0	268	0	2	0	137	6	1232	1133	209	109	2438	18	3	32	6	92	5489	77	0	6870	164	8	19	1537	1	6	17	0	38	1
New.ReferenceOTU177	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	501	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
EF018956.1.1514	0	100	0	212	44	0	2	41	0	101	0	77	134	89	1	366	0	1	1	0	0	0	0	0	137	0	0	0	0	1	13	1	11	20	0	44	0	10	1	66	0	0	68	0	83	1	0	0	36	769	0	0
New.ReferenceOTU30	0	123	0	108	20	0	5	57	0	64	0	80	108	70	0	187	0	0	0	0	0	0	0	0	112	0	1	2	0	0	17	1	9	22	0	32	0	1	1	46	0	0	52	0	46	3	1	0	28	778	0	0
New.ReferenceOTU166	0	0	0	0	21	0	6	233	0	110	0	0	1	39	0	246	0	0	0	0	0	0	0	0	7	0	0	2	0	1	0	0	7	4	0	2	0	7	2	10	0	0	7	0	26	1	0	0	15	372	0	0
JX040386.1.1468	0	1053	0	655	17	1	7	311	0	103	0	132	438	220	2	507	0	1	0	0	0	0	0	0	72	0	0	0	0	2	58	0	33	35	0	8	0	9	4	118	0	0	30	0	32	7	0	0	29	676	0	0
GDAX01048398.486.2474	14	1287	4	314	13	0	25	6	0	1	0	1	0	54	628	0	0	0	0	7	0	29	2	0	0	0	0	2	0	0	0	0	0	0	3	0	0	1	1	16	0	0	0	0	90	2	0	12	0	39	0	0
GAOK01049512.6271.8188	2	1	0	1	13	0	799	902	1	4	0	0	0	0	0	0	0	0	0	2753	4	25	0	0	0	0	0	0	0	2	202	1	0	0	1	0	0	0	4	0	0	0	0	0	1375	1	7	0	4	1	1	0
AGYD01000015.142520.144019	184	0	0	8	0	176	0	0	30	0	41	0	0	0	0	0	0	2	230	0	4	0	22	0	0	0	0	0	41	0	0	0	0	0	0	1	3	0	0	0	0	0	37	3	0	1	0	0	1	0	1561	0
JF179313.1.1396	0	0	0	0	0	0	0	434	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
New.CleanUp.ReferenceOTU12322	1	397	0	0	11	0	6	0	0	9	0	0	28	0	0	150	0	0	0	0	0	0	0	0	0	0	0	0	0	0	4	0	0	1	0	0	0	0	0	0	0	0	8	0	0	3	0	0	0	0	0	0
KC494323.1.1465	0	0	0	0	0	0	0	654	11	0	0	0	1	0	0	0	0	0	0	0	2	1	0	0	0	0	0	0	0	0	0	0	0	0	0	14	0	0	0	0	0	0	0	0	0	0	0	18	0	0	0	0
APCS01000185.14.1495	0	0	0	0	0	0	0	419	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
AXNY01000032.20116.21623	449	0	1	30	0	21	1	0	0	2	14	0	18	0	0	0	0	0	0	0	0	1	4	2	0	0	0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0
FPLS01013939.17.1510	0	0	0	617	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
FJ820431.1.1504	2	2368	0	910	47	0	47	1865	1	454	1	606	1439	1232	4	2856	0	0	0	0	2	0	0	1	497	1	0	9	0	15	202	1	114	104	0	180	0	38	13	439	2	1	364	0	209	15	2	0	147	2013	0	1
DQ130022.1.1439	60	0	0	0	0	10	0	0	0	0	12	0	0	0	0	0	24	0	1	0	69	0	0	13	0	0	0	0	0	0	0	0	0	0	0	0	87	0	0	0	0	0	0	0	0	1	0	0	0	0	2724	0
LTXH01000090.666.2207	57	0	0	0	0	1	0	0	0	0	9	0	0	0	0	0	1	0	0	0	51	0	0	1	0	0	0	0	1	5	3	0	0	0	2	0	3	0	0	0	1	0	0	0	0	0	0	0	0	0	744	0
BABK02029359.1608.3968	0	399	0	0	12	0	4	0	0	0	0	1	0	25	66	323	0	0	0	0	0	18	0	0	0	0	48	0	0	0	0	1	0	0	3	0	0	0	0	0	1	0	0	0	24	0	1	0	0	43	0	0
FPLS01018017.9.1494	0	457	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
New.CleanUp.ReferenceOTU29046	0	0	0	0	0	0	0	0	0	36	0	0	0	351	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	28	0	0	0	0	0	0	1	0	5	0	0	0	0	0	0	0
FR667305.1.1373	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	470	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
FPLS01033006.16.1516	0	449	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
LN564828.1.1344	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	392	0	0
JF712662.1.1483	5	0	606	66	0	7	0	0	521	0	34	0	0	0	3	0	302	0	0	941	2	32	5	481	0	477	531	0	34	0	0	29	0	17	82	0	27	1	0	0	4	0	8	0	7	0	0	0	0	0	0	0
AY014171.1.1462	2	0	574	67	2	69	2	464	0	1	0	6	55	4	0	34	0	0	0	0	0	27	1	0	103	0	0	0	1	0	3	0	0	1	3	7	2	0	0	0	5	0	0	0	0	1	0	0	4	0	0	0
HQ863983.1.1492	0	0	162	34	0	1	0	134	0	1	0	0	491	172	0	20	0	0	0	0	0	0	0	0	4	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	10	0	3	0	0	0	0	0	0	0
KF037544.1.1453	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	730	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	5	0	0	0	0
GU454960.1.1486	2	0	105	563	11	94	96	247	1036	38	0	140	550	201	0	0	133	0	0	1	897	0	8	5	36	0	0	2	32	13	8	0	0	5	0	28	3	3	3	43	14	8	7	2	49	0	1	0	1	18	0	0
ATTI01000001.185818.187289	64	509	3190	48	18	588	4	833	2727	17	331	7	77	80	12	325	1675	7	155	4646	15	260	91	2892	119	2822	3446	0	48	1	10	124	24	9	326	21	81	1	1	7	22	29	4	1	12	1	1	1	14	883	1	0
AFUV01000004.339540.341080	3	0	102	3	0	2	0	0	0	0	754	0	2	0	1	1	0	0	0	0	0	0	683	0	7	0	0	1	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
DQ272508.1.1560	146	0	196	344	0	644	28	0	0	79	5	0	205	21	0	23	0	2	0	0	0	0	4	7	74	22	0	13	29	5	5	0	0	0	0	4	0	0	0	1	0	0	0	0	43	0	0	13	0	0	0	0
New.ReferenceOTU49	975	0	0	0	0	62	0	0	0	0	0	1	0	0	3	1	0	0	0	0	0	0	0	3	0	8	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
CCPS01000022.154.1916	1	175	93	43	9766	0	7911	0	0	0	16	47	52	20	251	1	12	6	8	6	0	2	0	50	89	134	32	5	2	0	167	0	2482	7683	14	16	1	597	0	10953	0	14	1793	1	20	131	3	22	12041	0	1	0
KM116008.1.1667	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	1962	0	0	0	0	0	0	1	0	0	0
KF625188.1.1765	0	1	0	10	0	0	0	0	0	0	0	431	0	0	7	0	1	2	121	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
KF625171.1.1580	0	0	0	7	0	0	0	0	0	0	0	54	4	0	0	0	0	1	0	1555	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	13	0	0	0	0	0	0	0
KF625183.1.1786	0	0	117	64	0	0	0	0	1	0	2	191	5	0	4019	58	230	4271	1551	8	1	1	0	53	0	0	0	1	0	0	8	0	0	0	0	11	0	0	0	2	2	238	3	0	0	0	0	0	1	0	0	0
AX067466.91472.92986	1	0	0	0	0	0	0	0	0	0	10351	0	5	5	5	0	0	0	0	0	0	1	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0
HAFP01051053.8.1513	0	504	0	59	3	0	5	722	0	93	0	47	141	126	0	392	0	1	0	0	0	0	0	1	134	0	0	1	0	4	11	0	8	7	0	33	0	7	1	25	0	0	28	0	11	4	1	0	13	417	0	0
EF517956.1.1666	3	100	2	40	8	16	1	100	0	59	56	253	108	10	0	760	0	0	0	1	2	1	4	0	16	0	0	6	1	1	1	11	0	5	7	8	0	1	0	23	0	0	0	0	23	0	0	1	5	309	0	0
AHBD01000036.106.1641	0	0	0	102	0	0	24	0	0	0	0	0	75	455	0	0	0	0	0	0	0	28	1	0	0	0	0	0	0	0	1	0	0	0	2	0	1	0	0	0	0	19	0	0	18	0	0	0	0	0	0	0
HK556428.10.1520	0	0	0	70	0	27	1	0	0	0	0	0	367	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
FM872796.1.1497	145	0	553	0	0	924	46	0	2310	5	0	46	0	0	0	0	0	0	0	0	0	0	106	0	0	5088	0	0	366	0	0	16	0	1	0	0	0	0	0	0	0	0	0	1	1	1100	1	0	0	0	0	0
FM873682.1.1495	0	0	0	0	0	0	101	0	423	0	0	38	0	0	0	946	0	0	3893	0	2	0	1	0	0	0	0	0	281	0	0	0	7	2	0	0	0	0	0	1	0	0	0	0	0	863	0	0	0	1940	0	0
HK556626.10.1530	354	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
ACQV01000002.9075.10620	0	59	395	78	0	236	2	0	4	160	0	12	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	22	14	1	0	0	0	2	0	0	0	0	0	31	2	0	0	0	0	0	0	0	0	0	0	0
FJ375378.1.1557	0	124	5	3	8	14	5	353	53	66	1	0	135	3	0	0	61	0	44	91	2	1	1	0	33	53	46	2	0	1	0	7	4	6	0	0	0	0	0	5	0	53	8	0	10	0	0	1	2	92	0	0
GAYR01026362.62.2014	27	31	181	214	5	26	55	2896	45	3	19	48	462	1474	13	295	22	1	0	0	31	36	43	155	666	0	264	5	3	2	27	325	0	6	10	204	12	7	6	15	1463	1362	8	0	72	1	0	36	11	1163	0	0
BCUL01000049.1202.2709	0	0	0	0	0	0	0	0	27	0	1	0	0	0	0	0	0	0	0	0	466	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
EU979056.1.1490	0	0	0	0	0	0	0	632	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	5	0	0	0	0	0	0
FLQT01000003.189534.191041	6	0	0	0	0	0	0	0	0	0	3	0	8	0	0	0	0	1	518	0	21	0	0	0	0	0	0	0	0	0	0	0	0	0	403	0	14	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
FXLP01000003.3248.4752	0	0	0	4	0	0	0	0	2	0	3	0	1	0	2	0	0	0	3	1	2	0	0	0	1	1	2090	0	0	1	0	0	1	0	0	0	0	0	0	0	0	5	0	0	0	1	0	0	0	0	0	0
AB559667.1.1517	815	288	471	222	1	470	0	0	29	0	1058	113	91	0	754	1	0	7	3807	1797	747	0	915	69	6	232	1020	10	1172	5	8	20	439	16	283	27	25	0	26	3	1	4792	183	78	0	1494	1	0	6	37	134	0
HM252077.1.1349	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	34	0	0	0	0	0	0	0	0	0	2	0	524	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0
JRMY01000047.88.1633	0	0	0	0	0	0	1	0	40	0	0	0	0	0	0	0	0	1	0	0	784	0	0	101	0	0	0	0	0	0	3	0	0	0	0	0	2	0	0	0	0	0	0	0	0	1	0	0	4	0	0	0
HE681314.1.1549	1432	0	0	0	0	0	8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
CP012071.1027746.1029299	696	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	24	0	0	0	2	6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
AM697408.1.1522	38	1	0	120	0	20	0	0	0	0	0	0	103	0	43	0	0	1	537	0	0	0	0	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	69	0	0
AY672076.1.1620	124	0	46	288	0	334	37	249	0	16	14	0	413	0	9	329	473	6	0	0	0	0	7	575	72	223	0	1	24	0	1	0	10	0	0	22	0	0	1	0	0	1	0	0	0	5	0	6	0	0	0	0
AY672087.1.1580	22	0	7	84	0	141	15	23	0	8	0	0	418	0	0	24	314	4	0	0	0	0	18	113	13	22	0	0	12	0	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0
JF010199.1.1375	167	41	18	75	233	61	95	29	26	78	103	120	121	98	34	30	49	360	16	39	114	213	60	51	141	43	30	22	97	640	103	59	267	163	99	91	216	44	213	157	197	68	276	53	88	429	157	6	198	40	52	8
FJ375394.1.1637	1	0	0	427	6	4	1	0	0	0	0	54	26	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	16	0	2	14	0	0	0	2	0	5	0	0	0	0	8	1	0	0	0	6	0	0
FJ375404.1.1601	0	0	0	0	0	0	0	0	0	9	0	1	0	0	0	190	0	0	0	0	0	0	0	0	30	0	0	0	0	0	0	11	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	618	0	0
ACWY01000141.90.1600	0	0	292	0	0	18	0	0	0	34	0	0	513	0	0	0	0	3	0	0	0	0	1	213	0	0	0	8	57	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
JNGW01000135.151.1667	438	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
CVNZ01002445.142.1661	49	0	0	101	0	0	0	0	3	10	0	0	110	0	0	0	0	0	671	0	0	0	0	0	0	0	0	0	10	1	0	1	9	8	0	0	0	0	0	0	0	0	0	0	0	3	0	0	0	497	0	0
HE681243.1.1520	378	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	9	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0
HM250376.1.1374	4	2	6	2	3	2	5	4	3	12	3	6	1	5	10	2	9	6	3	4	5	6	10	9	1	4	1	7	3	389	5	18	4	14	0	4	10	8	15	10	2	4	6	18	14	2	4	19	2	11	0	0
AGEF01000001.130585.132119	2	0	1	177	0	1	10846	3	5366	2	0	6	26	0	0	22	3	87	0	0	0	6	0	0	212	0	1	0	6750	1	0	5510	0	6	0	7171	1	0	0	0	0	0	26	2940	0	1	5	2	0	0	0	2
New.ReferenceOTU143	0	35	0	83	0	0	5	0	0	0	0	44	0	143	0	394	0	0	0	0	0	0	0	0	76	0	0	0	0	0	0	0	0	0	0	8	0	0	0	9	0	0	0	0	48	3	1	0	1	414	0	0
New.ReferenceOTU44	0	373	0	65	0	0	0	0	0	31	0	18	88	54	0	140	0	0	0	0	0	0	0	0	4	0	0	0	0	0	7	0	1	17	1	0	0	1	1	0	0	0	4	0	18	0	0	0	0	186	0	0
New.ReferenceOTU159	0	0	0	0	14	0	0	395	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
New.ReferenceOTU113	0	402	0	66	25	0	0	764	0	6	0	89	14	30	0	298	0	0	0	0	0	0	0	0	88	0	0	0	0	4	0	0	4	6	0	8	0	0	0	0	0	0	7	0	64	1	0	0	17	90	0	0
HK555703.10.1512	6	0	0	0	2	0	2	122	23	0	4	7	0	92	0	0	0	0	0	0	521	0	2	0	7	0	0	0	0	0	0	0	0	11	459	20	1617	0	0	0	0	0	0	0	0	3	1	1	0	0	0	0
JX079489.1.1518	3	391	44	0	10	134	0	4	40	6	1	0	0	0	3	9	12	2	18	74	0	3	6	0	0	10	14	0	15	16	105	6	46	3	0	0	1	9	16	3	0	88	2	0	0	10	0	0	0	77	0	0
KU978311.1.1503	0	0	0	0	10	0	5	0	5	0	0	0	0	8	7	0	0	1	0	0	4	4	0	36	1509	0	0	4	0	0	0	1	0	15	0	384	0	32	0	1	0	0	0	0	84	0	1	1	13	13	1	0
KF360241.1.1478	0	0	0	0	0	0	0	0	0	0	0	0	0	0	10	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	439	0	212	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3	0	0	0	0
AY534742.1.1519	0	0	0	0	0	0	0	0	37	0	0	0	0	0	0	0	0	0	0	0	918	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
EF517954.1.1705	14	172	211	2408	136	119	172	418	131	7	25	130	9	122	10	520	235	32	59	46	5540	6	18	15	150	201	112	59	130	35	88	17	6	15	332	209	371	22	5	66	6	2	29	7	18	68	11	6826	117	172	1	0
MRSV01000001.467980.469577	0	0	0	0	0	3	0	0	0	0	0	15	0	1	0	256	0	0	0	0	0	0	0	0	5	0	1	0	0	0	0	2	0	0	0	5	0	0	0	0	0	0	0	0	4	0	0	0	0	853	0	0
ADHG01000001.652476.654027	1	0	0	36	0	0	1	0	0	0	0	16	92	0	0	1468	0	0	0	0	0	0	1	11	0	0	0	0	0	6	0	2	1	199	1	0	0	0	0	231	0	0	2	0	0	4	0	0	0	1746	0	0
JF123767.1.1373	2	3	3	1	7	4	1	1	1	4	0	8	1	3	7	2	7	0	7	1	9	5	9	5	1	1	3	9	2	2	6	7	9	2	2	1	5	10	6	7	0	4	3	7	4	2	2122	9	55	1	0	0
FM873871.1.1482	2	0	0	0	3	43	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	4	0	0	0	1	0	0	0	0	0	0	0	495	0	75	0	0	0	0	0	0	0	0	2	0	0	0	0	161	0
AB853090.1.1522	1667	173	144	39	0	1310	1	0	52	1	52	73	23	0	0	0	0	0	1	120	0	0	917	0	126	192	90	1	658	1	1	0	548	2	20	46	4	2	63	0	43	44	0	29	0	76	0	0	0	0	61	0
CAJU020000019.27308.28816	320	6	55	0	11	107	0	0	0	0	171	0	0	20	173	0	0	0	1282	54	2	0	85	11	11	0	0	1	87	0	5	0	0	13	3	4	55	2	0	0	0	0	36	2	0	4	0	19	25	0	798	0
AM697027.1.1482	81	4	2	2	6	7	1	0	0	0	2	1	0	0	1	0	3	0	3	1	0	0	1	0	11	0	8	0	43	0	0	3	7	1	197	0	0	0	1	0	3632	0	19	0	0	10	1	0	69	0	26	1
AM697214.1.1482	3	0	0	0	0	94	0	0	0	0	9	0	0	0	0	0	0	2	0	0	2	0	0	0	1	0	0	0	0	0	0	0	0	8	306	0	12	0	1	1	4092	0	0	0	0	2	0	0	14	0	2983	0
DQ532162.1.1429	0	1	0	0	0	0	0	0	1	1	1	1	0	0	0	0	1	0	0	0	0	0	0	0	34	0	1	0	1	0	2	0	0	0	1103	19	0	0	0	0	0	0	0	0	0	1	0	0	0	0	11	0
KX504550.1.1489	549	6	8	0	5	632	0	3	1	0	434	2	2	0	12	0	13	0	850	27	0	1	4	2	0	0	26	0	543	0	2	53	1998	0	1	50	0	4	0	0	0	0	0	54	0	9	0	0	5	1	54	0
FM875546.1.1480	1563	395	592	60	0	1427	0	280	77	2	178	108	14	5	1062	0	877	1	483	1685	1	1	396	102	8	56	2695	0	453	2	6	374	3	0	5	7	1	1	12	0	39	0	0	44	7	664	0	6	0	3	2	0
AM697289.1.1481	7	0	1	0	0	14	0	1	0	0	2	0	1	47	1	0	2	1	0	3	0	0	3	0	1	1	1	0	2	0	0	1	1	1	3	0	0	0	0	0	6	0	0	0	0	2	0	0	10	0	1289	0
KJ569766.1.1491	0	0	0	0	0	4	0	1	0	0	0	0	0	0	0	0	0	0	2	0	0	0	1	0	881	0	0	0	0	0	0	0	0	0	4	22	0	1	0	0	2	0	2	0	0	1	0	0	0	0	2	0
EF419348.1.1393	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	472	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
KF625179.1.1776	0	50	0	34	4	0	0	446	0	8	0	15	57	42	0	22	0	0	0	0	0	0	0	0	11	0	0	0	0	0	3	0	1	0	0	2	0	0	0	0	0	0	0	0	1	1	0	0	1	0	1	0
CVJT01000011.50.2173	9826	10183	12307	5049	15738	9363	7612	5156	9962	21385	8302	16282	6618	12829	19681	7467	17449	9784	7964	10067	8401	12742	21811	20561	6260	8286	7270	30547	8308	12980	25586	17062	12679	20455	13186	13338	12506	28377	21933	15618	4053	7759	13930	26203	25461	9195	21333	24077	10956	6063	1929	33
New.ReferenceOTU77	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	569	0	173	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
AGZC01000026.457.1994	8	0	0	5	5	198	1	0	40	0	1	0	45	0	0	0	0	3	2	0	2	0	0	68	0	0	0	0	7	0	0	0	0	0	991	0	1524	0	0	11	0	0	0	0	0	0	2	0	3	0	3631	0
CRKN01000054.4139.5702	225	419	765	1527	15	1482	130	514	0	1084	1220	30	1334	7	69	10	862	6	2	45	0	19	308	1666	139	150	40	10	112	6	0	1	3	2	0	55	0	3	1	46	2	12	22	1	18	29	0	45	2	103	0	0

# Calculate the relative abundance for each sample
ps.ng.tax_abund_rel <- transform_sample_counts(ps.ng.tax_abund, function(x) x/sum(x))

datamat_ = as.data.frame(otu_table(ps.ng.tax_abund))
datamat <- datamat_[c("P7.Nose", "P8.Nose", "P9.Nose", "P10.Nose", "P11.Nose", "P12.Nose",
    "P13.Nose", "P14.Nose", "P15.Nose", "P16.Foot", "P17.Foot", "P18.Foot", "P19.Foot",
    "P20.Foot", "P16.Nose", "P17.Nose", "P18.Nose", "P19.Nose", "P20.Nose", "AH.LH",
    "AL.LH", "CB.LH", "HR.LH", "KK.LH", "MC.LH", "MR.LH", "PT2.LH", "RP.LH", "SA.LH",
    "XN.LH", "AH.NLH", "AL.NLH", "CB.NLH", "HR.NLH", "KK.NLH", "MC.NLH", "MR.NLH",
    "PT2.NLH", "RP.NLH", "SA.NLH", "XN.NLH", "AH.Nose", "AL.Nose", "CB.Nose", "HR.Nose",
    "KK.Nose", "MC.Nose", "MR.Nose", "PT2.Nose", "RP.Nose", "SA.Nose", "XN.Nose")]
hr <- hclust(as.dist(1 - cor(t(datamat), method = "pearson")), method = "complete")
hc <- hclust(as.dist(1 - cor(datamat, method = "spearman")), method = "complete")
mycl = cutree(hr, h = max(hr$height)/1.08)
mycol = c("YELLOW", "DARKBLUE", "DARKORANGE", "DARKMAGENTA", "DARKCYAN", "DARKRED",
    "MAROON", "DARKGREEN", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN", "LIGHTGREEN",
    "BLUE", "ORANGE", "CYAN", "RED", "GREEN")

mycol = mycol[as.vector(mycl)]
sampleCols <- rep("GREY", ncol(datamat))

sampleCols[colnames(datamat) == "P7.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P8.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P9.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P10.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P11.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P12.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P13.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P14.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P15.Nose"] <- "#1f78b4"

sampleCols[colnames(datamat) == "P16.Foot"] <- "#b2df8a"
sampleCols[colnames(datamat) == "P17.Foot"] <- "#b2df8a"
sampleCols[colnames(datamat) == "P18.Foot"] <- "#b2df8a"
sampleCols[colnames(datamat) == "P19.Foot"] <- "#b2df8a"
sampleCols[colnames(datamat) == "P20.Foot"] <- "#b2df8a"
sampleCols[colnames(datamat) == "P16.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P17.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P18.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P19.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "P20.Nose"] <- "#1f78b4"

# sampleCols[colnames(datamat)=='P16.Nose'] <- '#33a02c'
# sampleCols[colnames(datamat)=='P16.Nose'] <- '#33a02c'

# sampleCols[colnames(datamat)=='7'] <- '#a6cee3'
sampleCols[colnames(datamat) == "AH.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "AL.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "CB.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "HR.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "KK.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "MC.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "MR.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "PT2.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "RP.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "SA.LH"] <- "#fb9a99"
sampleCols[colnames(datamat) == "XN.LH"] <- "#fb9a99"

sampleCols[colnames(datamat) == "AH.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "AL.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "CB.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "HR.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "KK.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "MC.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "MR.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "PT2.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "RP.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "SA.NLH"] <- "#e31a1c"
sampleCols[colnames(datamat) == "XN.NLH"] <- "#e31a1c"

sampleCols[colnames(datamat) == "AH.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "AL.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "CB.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "HR.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "KK.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "MC.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "MR.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "PT2.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "RP.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "SA.Nose"] <- "#1f78b4"
sampleCols[colnames(datamat) == "XN.Nose"] <- "#1f78b4"

# sampleCols[colnames(datamat)=='AH.Nose'] <- '#cab2d6'
# sampleCols[colnames(datamat)=='P16.Nose'] <- '#6a3d9a'

# bluered(75) color_pattern <- colorRampPalette(c('blue', 'white', 'red'))(100)
library(RColorBrewer)
custom_palette <- colorRampPalette(brewer.pal(9, "Blues"))
heatmap_colors <- custom_palette(100)
# colors <- heatmap_color_default(100)
png("figures/heatmap.png", width = 1200, height = 2400)
# par(mar=c(2, 2, 2, 2)) , lwid=1 lhei=c(0.7, 10)) # Adjust height of color
# keys keysize=0.3,
heatmap.2(as.matrix(datamat), Rowv = as.dendrogram(hr), Colv = NA, dendrogram = "row",
    scale = "row", trace = "none", col = heatmap_colors, cexRow = 1.2, cexCol = 1.5,
    RowSideColors = mycol, ColSideColors = sampleCols, srtCol = 60, labRow = row.names(datamat),
    key = TRUE, margins = c(10, 15), lhei = c(0.7, 15), lwid = c(1, 8))
dev.off()

png 
  2

knitr::include_graphics("./figures/heatmap.png")

Heatmap

5 Taxonomic summary at the Phylum Level

5.1 Bar plots in P16-P20 samples

library(ggplot2)
geom.text.size = 6
theme.size = 8  #(14/5) * geom.text.size
# ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in%
# c('D_1__Actinobacteria', 'D_1__Bacteroidetes', 'D_1__Firmicutes',
# 'D_1__Proteobacteria', 'D_1__Verrucomicrobia', NA))
ps.ng.tax_most = phyloseq::filter_taxa(ps.ng.tax_rel, function(x) mean(x) > 0.001,
    TRUE)
# CONSOLE(OPTIONAL): for sampleid in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
# 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
# 70 71 72 73; do echo
# 'otu_table(ps.ng.tax_most)[,${sampleid}]=otu_table(ps.ng.tax_most)[,${sampleid}]/sum(otu_table(ps.ng.tax_most)[,${sampleid}])'
# done OR
ps.ng.tax_most_ = transform_sample_counts(ps.ng.tax_most, function(x) x/sum(x))

## ---- Phylum ---- aes(color='Phylum', fill='Phylum') --> aes()
## ggplot(data=data, aes(x=Sample, y=Abundance, fill=Phylum))

ps_subset_2_type <- merge_samples(ps_subset_2, "Type")
ps_subset_2_type_ = transform_sample_counts(ps_subset_2_type, function(x) x/sum(x))
ps_subset_2_copied <- data.table::copy(ps_subset_2)
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Phylum') +
# geom_bar(aes(color=Phylum, fill=Phylum), stat='identity', position='stack')
# plot_bar(ps_subset_2_type_, fill='Phylum') + geom_bar(aes(), stat='identity',
# position='stack') + scale_fill_manual(values = c('darkblue',
# 'darkgoldenrod1', 'darkseagreen', 'darkorchid', 'darkolivegreen1',
# 'lightskyblue', 'darkgreen', 'deeppink', 'khaki2', 'firebrick', 'brown1',
# 'darkorange1', 'cyan1', 'royalblue4', 'darksalmon', 'darkblue','royalblue4',
# 'dodgerblue3', 'steelblue1', 'lightskyblue', 'darkseagreen',
# 'darkgoldenrod1', 'darkseagreen', 'darkorchid', 'darkolivegreen1', 'brown1',
# 'darkorange1', 'cyan1', 'darkgrey')) + theme(axis.text = element_text(size =
# theme.size, colour='black'))


# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P16.Foot", "P17.Foot", "P18.Foot", "P19.Foot", "P20.Foot",
    "P16.Nose", "P17.Nose", "P18.Nose", "P19.Nose", "P20.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Phylum'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Phylum), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

Regroup together Foot vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Phylum') +
# geom_bar(aes(color=Phylum, fill=Phylum), stat='identity', position='stack')
plot_bar(ps_subset_type_, fill = "Phylum") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset_2, x = "Phylum", fill = "Phylum", facet_grid = Patient ~ Type) +
    geom_bar(aes(), stat = "identity", position = "stack") + scale_fill_manual(values = c("darkblue",
    "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue",
    "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1",
    "royalblue4", "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1",
    "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid",
    "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black"), axis.text.x = element_blank(), axis.ticks = element_blank()) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

# samples_to_keep <-
# c('P7.Nose','P8.Nose','P9.Nose','P10.Nose','P11.Nose','P12.Nose','P13.Nose','P14.Nose','P15.Nose')

5.2 Bar plots in AH-XN samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("AH.LH", "AL.LH", "CB.LH", "HR.LH", "KK.LH", "MC.LH", "MR.LH",
    "PT2.LH", "RP.LH", "SA.LH", "XN.LH", "AH.NLH", "AL.NLH", "CB.NLH", "HR.NLH",
    "KK.NLH", "MC.NLH", "MR.NLH", "PT2.NLH", "RP.NLH", "SA.NLH", "XN.NLH", "AH.Nose",
    "AL.Nose", "CB.Nose", "HR.Nose", "KK.Nose", "MC.Nose", "MR.Nose", "PT2.Nose",
    "RP.Nose", "SA.Nose", "XN.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Phylum'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Phylum), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 9,
    angle = 60, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

Regroup together LH vs NLH vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Phylum') +
# geom_bar(aes(color=Phylum, fill=Phylum), stat='identity', position='stack')
plot_bar(ps_subset_type_, fill = "Phylum") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset, x = "Phylum", fill = "Phylum", facet_grid = Patient ~ Type) +
    geom_bar(aes(), stat = "identity", position = "stack") + scale_fill_manual(values = c("darkblue",
    "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue",
    "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1",
    "royalblue4", "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1",
    "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid",
    "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black"), axis.text.x = element_blank(), axis.ticks = element_blank()) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

5.3 Bar plots in P7-P15 samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P7.Nose", "P8.Nose", "P9.Nose", "P10.Nose", "P11.Nose", "P12.Nose",
    "P13.Nose", "P14.Nose", "P15.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Phylum'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Phylum), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 9,
    angle = 60, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))

6 Taxonomic summary at the Class Level

6.1 Bar plots in P7-P15 samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P7.Nose", "P8.Nose", "P9.Nose", "P10.Nose", "P11.Nose", "P12.Nose",
    "P13.Nose", "P14.Nose", "P15.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Class), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 9,
    angle = 60, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

# aes(color='Class', fill='Class') --> aes() ggplot(data=data, aes(x=Sample,
# y=Abundance, fill=Class))

ps_subset_2_type <- merge_samples(ps_subset_2, "Type")
ps_subset_2_type_ = transform_sample_counts(ps_subset_2_type, function(x) x/sum(x))
ps_subset_2_copied <- data.table::copy(ps_subset_2)
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Class') +
# geom_bar(aes(color=Class, fill=Class), stat='identity', position='stack')
# plot_bar(ps_subset_2_type_, fill='Class') + geom_bar(aes(), stat='identity',
# position='stack') + scale_fill_manual(values = c('darkblue',
# 'darkgoldenrod1', 'darkseagreen', 'darkorchid', 'darkolivegreen1',
# 'lightskyblue', 'darkgreen', 'deeppink', 'khaki2', 'firebrick', 'brown1',
# 'darkorange1', 'cyan1', 'royalblue4', 'darksalmon', 'darkblue','royalblue4',
# 'dodgerblue3', 'steelblue1', 'lightskyblue', 'darkseagreen',
# 'darkgoldenrod1', 'darkseagreen', 'darkorchid', 'darkolivegreen1', 'brown1',
# 'darkorange1', 'cyan1', 'darkgrey')) + theme(axis.text = element_text(size =
# theme.size, colour='black'))


# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P16.Foot", "P17.Foot", "P18.Foot", "P19.Foot", "P20.Foot",
    "P16.Nose", "P17.Nose", "P18.Nose", "P19.Nose", "P20.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Class), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Regroup together Foot vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Class') +
# geom_bar(aes(color=Class, fill=Class), stat='identity', position='stack')
plot_bar(ps_subset_type_, fill = "Class") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset_2, x = "Class", fill = "Class", facet_grid = Patient ~ Type) +
    geom_bar(aes(), stat = "identity", position = "stack") + scale_fill_manual(values = c("darkblue",
    "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue",
    "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1",
    "royalblue4", "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1",
    "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid",
    "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black"), axis.text.x = element_blank(), axis.ticks = element_blank()) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

# samples_to_keep <-
# c('P7.Nose','P8.Nose','P9.Nose','P10.Nose','P11.Nose','P12.Nose','P13.Nose','P14.Nose','P15.Nose')

6.2 Bar plots in AH-XN samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("AH.LH", "AL.LH", "CB.LH", "HR.LH", "KK.LH", "MC.LH", "MR.LH",
    "PT2.LH", "RP.LH", "SA.LH", "XN.LH", "AH.NLH", "AL.NLH", "CB.NLH", "HR.NLH",
    "KK.NLH", "MC.NLH", "MR.NLH", "PT2.NLH", "RP.NLH", "SA.NLH", "XN.NLH", "AH.Nose",
    "AL.Nose", "CB.Nose", "HR.Nose", "KK.Nose", "MC.Nose", "MR.Nose", "PT2.Nose",
    "RP.Nose", "SA.Nose", "XN.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Class), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 9,
    angle = 60, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Regroup together LH vs NLH vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
# plot_bar(ps.ng.tax_most_SampleType_, fill = 'Class') +
# geom_bar(aes(color=Class, fill=Class), stat='identity', position='stack')
plot_bar(ps_subset_type_, fill = "Class") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 10,
    angle = 15, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset, x = "Class", fill = "Class", facet_grid = Patient ~ Type) + geom_bar(aes(),
    stat = "identity", position = "stack") + scale_fill_manual(values = c("darkblue",
    "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue",
    "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1",
    "royalblue4", "darksalmon", "darkblue", "royalblue4", "dodgerblue3", "steelblue1",
    "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid",
    "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5,
    colour = "black"), axis.text.x = element_blank(), axis.ticks = element_blank()) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

6.3 Bar plots in P7-P15 samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P7.Nose", "P8.Nose", "P9.Nose", "P10.Nose", "P11.Nose", "P12.Nose",
    "P13.Nose", "P14.Nose", "P15.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps.ng.tax_most_, rownames(sample_data(ps.ng.tax_most_)) %in%
    samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = Class), stat = "identity") +
    scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick",
        "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue",
        "royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1",
        "darkorange1", "cyan1", "darkgrey")) + theme(axis.text.x = element_text(size = 9,
    angle = 60, hjust = 1), axis.text.y = element_text(size = 9, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

7 Alpha diversity

Plot Chao1 richness estimator, Observed OTUs, Shannon index, and Phylogenetic diversity. Regroup together samples from the same group.

hmp.div_qiime <- read.csv("adiv_even.txt", sep = "\t")
colnames(hmp.div_qiime) <- c("sam_name", "chao1", "observed_otus", "shannon", "PD_whole_tree")
row.names(hmp.div_qiime) <- hmp.div_qiime$sam_name
div.df <- merge(hmp.div_qiime, hmp.meta, by = "sam_name")
div.df2 <- div.df[, c("Group", "chao1", "shannon", "observed_otus", "PD_whole_tree")]
colnames(div.df2) <- c("Group", "Chao-1", "Shannon", "OTU", "Phylogenetic Diversity")
# colnames(div.df2)
options(max.print = 999999)
# 27 H47 830.5000 5.008482 319 10.60177 FITTING4: if occuring 'Computation
# failed in `stat_signif()`:not enough 'y' observations' means: the patient H47
# contains only one sample, it should be removed for the statistical p-values
# calculations.  delete H47(1) div.df2 <- div.df2[-c(3), ] div.df2 <-
# div.df2[-c(55,54, 45,40,39,27,26,25,1), ] knitr::kable(div.df2) %>%
# kable_styling(bootstrap_options = c('striped', 'hover', 'condensed',
# 'responsive'))

# https://uc-r.github.io/t_test We can perform the test with t.test and
# transform our data and we can also perform the nonparametric test with the
# wilcox.test function.
stat.test.Shannon <- compare_means(Shannon ~ Group, data = div.df2, method = "t.test")
knitr::kable(stat.test.Shannon) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

.y.	group1	group2	p	p.adj	p.format	p.signif	method
Shannon	AH-XN.LH	AH-XN.NLH	0.2985783	1.00	0.299	ns	T-test
Shannon	AH-XN.LH	AH-XN.Nose	0.4783590	1.00	0.478	ns	T-test
Shannon	AH-XN.LH	P7-P15.Nose	0.8404760	1.00	0.840	ns	T-test
Shannon	AH-XN.LH	P16-P20.Foot	0.3047879	1.00	0.305	ns	T-test
Shannon	AH-XN.LH	P16-P20.Nose	0.4767965	1.00	0.477	ns	T-test
Shannon	AH-XN.NLH	AH-XN.Nose	0.0584535	0.76	0.058	ns	T-test
Shannon	AH-XN.NLH	P7-P15.Nose	0.2788795	1.00	0.279	ns	T-test
Shannon	AH-XN.NLH	P16-P20.Foot	0.8243016	1.00	0.824	ns	T-test
Shannon	AH-XN.NLH	P16-P20.Nose	0.5220542	1.00	0.522	ns	T-test
Shannon	AH-XN.Nose	P7-P15.Nose	0.1206787	1.00	0.121	ns	T-test
Shannon	AH-XN.Nose	P16-P20.Foot	0.0314331	0.44	0.031		T-test
Shannon	AH-XN.Nose	P16-P20.Nose	0.0194719	0.29	0.019		T-test
Shannon	P7-P15.Nose	P16-P20.Foot	0.2420622	1.00	0.242	ns	T-test
Shannon	P7-P15.Nose	P16-P20.Nose	0.4146842	1.00	0.415	ns	T-test
Shannon	P16-P20.Foot	P16-P20.Nose	0.5675750	1.00	0.568	ns	T-test

div_df_melt <- reshape2::melt(div.df2)
# head(div_df_melt)

# https://plot.ly/r/box-plots/#horizontal-boxplot
# http://www.sthda.com/english/wiki/print.php?id=177
# https://rpkgs.datanovia.com/ggpubr/reference/as_ggplot.html
# http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/82-ggplot2-easy-way-to-change-graphical-parameters/
# https://plot.ly/r/box-plots/#horizontal-boxplot library('gridExtra')
# par(mfrow=c(4,1))
p <- ggboxplot(div_df_melt, x = "Group", y = "value", facet.by = "variable", scales = "free",
    width = 0.5, fill = "gray", legend = "right")
# ggpar(p, xlab = FALSE, ylab = FALSE)
lev <- levels(factor(div_df_melt$Group))  # get the variables
# FITTING4: delete H47(1) in lev lev <- lev[-c(3)] make a pairwise list that we
# want to compare. my_stat_compare_means
# https://stackoverflow.com/questions/47839988/indicating-significance-with-ggplot2-in-a-boxplot-with-multiple-groups
L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i])  #%>% filter(p.signif != 'ns')
my_stat_compare_means <- function(mapping = NULL, data = NULL, method = NULL, paired = FALSE,
    method.args = list(), ref.group = NULL, comparisons = NULL, hide.ns = FALSE,
    label.sep = ", ", label = NULL, label.x.npc = "left", label.y.npc = "top", label.x = NULL,
    label.y = NULL, tip.length = 0.03, symnum.args = list(), geom = "text", position = "identity",
    na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, ...) {
    if (!is.null(comparisons)) {
        method.info <- ggpubr:::.method_info(method)
        method <- method.info$method
        method.args <- ggpubr:::.add_item(method.args, paired = paired)
        if (method == "wilcox.test")
            method.args$exact <- FALSE
        pms <- list(...)
        size <- ifelse(is.null(pms$size), 0.3, pms$size)
        color <- ifelse(is.null(pms$color), "black", pms$color)
        map_signif_level <- FALSE
        if (is.null(label))
            label <- "p.format"
        if (ggpubr:::.is_p.signif_in_mapping(mapping) | (label %in% "p.signif")) {
            if (ggpubr:::.is_empty(symnum.args)) {
                map_signif_level <- c(`****` = 1e-04, `***` = 0.001, `**` = 0.01,
                  `*` = 0.05, ns = 1)
            } else {
                map_signif_level <- symnum.args
            }
            if (hide.ns)
                names(map_signif_level)[5] <- " "
        }
        step_increase <- ifelse(is.null(label.y), 0.12, 0)
        ggsignif::geom_signif(comparisons = comparisons, y_position = label.y, test = method,
            test.args = method.args, step_increase = step_increase, size = size,
            color = color, map_signif_level = map_signif_level, tip_length = tip.length,
            data = data)
    } else {
        mapping <- ggpubr:::.update_mapping(mapping, label)
        layer(stat = StatCompareMeans, data = data, mapping = mapping, geom = geom,
            position = position, show.legend = show.legend, inherit.aes = inherit.aes,
            params = list(label.x.npc = label.x.npc, label.y.npc = label.y.npc, label.x = label.x,
                label.y = label.y, label.sep = label.sep, method = method, method.args = method.args,
                paired = paired, ref.group = ref.group, symnum.args = symnum.args,
                hide.ns = hide.ns, na.rm = na.rm, ...))
    }
}

p3 <- p + stat_compare_means(method = "t.test", comparisons = list(c("P16-P20.Foot",
    "P16-P20.Nose"), c("AH-XN.LH", "AH-XN.Nose"), c("AH-XN.NLH", "AH-XN.Nose")),
    label = "p.signif", symnum.args <- list(cutpoints = c(0, 1e-04, 0.001, 0.01,
        0.05, 1), symbols = c("****", "***", "**", "*", "ns")), ) + theme(axis.text.x = element_text(angle = 30,
    hjust = 1))  # add this line to rotate x-axis text
print(p3)

# https://stackoverflow.com/questions/20500706/saving-multiple-ggplots-from-ls-into-one-and-separate-files-in-r
# FITTING3: mkdir figures
ggsave("./figures/alpha_diversity_Group2.png", device = "png", height = 10, width = 12)
ggsave("./figures/alpha_diversity_Group2.svg", device = "svg", height = 10, width = 12)

8 ST summary

## Change your working directory to where the files are located
count_data <- read.table("../count_table_.txt", header = TRUE, row.names = 1, sep = ",")
OTU = otu_table(as.matrix(count_data), taxa_are_rows = TRUE)

sample_data <- read.csv("../map3_corrected.txt", sep = "\t", row.names = 1)
SAM = sample_data(sample_data, errorIfNULL = TRUE)
# rownames(SAM) <- c('P7.Nose', 'P8.Nose', 'P9.Nose', 'P10.Nose', 'P11.Nose',
# 'P12.Nose', 'P13.Nose', 'P14.Nose', 'P15.Nose', 'P16.Foot', 'P16.Nose',
# 'P17.Foot', 'P17.Nose', 'P18.Foot', 'P18.Nose', 'P19.Foot', 'P19.Nose',
# 'P20.Foot', 'P20.Nose', 'AH.LH', 'AH.NLH', 'AH.Nose', 'AL.LH', 'AL.NLH',
# 'AL.Nose', 'CB.LH', 'CB.NLH', 'CB.Nose', 'HR.LH', 'HR.NLH', 'HR.Nose',
# 'KK.LH', 'KK.NLH', 'KK.Nose', 'MC.LH', 'MC.NLH', 'MC.Nose', 'MR.LH',
# 'MR.NLH', 'MR.Nose', 'PT2.LH', 'PT2.NLH', 'PT2.Nose', 'RP.LH', 'RP.NLH',
# 'RP.Nose', 'SA.LH', 'SA.NLH', 'SA.Nose', 'XN.LH', 'XN.NLH', 'XN.Nose')

# tree_file <- read_tree('../clustering/rep_set.tre')
ps <- phyloseq(OTU, SAM)  #, tree_file)
saveRDS(ps, "./ST_ps.rds")

Visualize data

sample_names(ps)

 [1] "P7.Nose"  "P8.Nose"  "P9.Nose"  "P10.Nose" "P11.Nose" "P12.Nose"
 [7] "P13.Nose" "P14.Nose" "P15.Nose" "P16.Foot" "P16.Nose" "P17.Foot"
[13] "P17.Nose" "P18.Foot" "P18.Nose" "P19.Foot" "P19.Nose" "P20.Foot"
[19] "P20.Nose" "AH.LH"    "AH.NLH"   "AH.Nose"  "AL.LH"    "AL.NLH"  
[25] "AL.Nose"  "CB.LH"    "CB.NLH"   "CB.Nose"  "HR.LH"    "HR.NLH"  
[31] "HR.Nose"  "KK.LH"    "KK.NLH"   "KK.Nose"  "MC.LH"    "MC.NLH"  
[37] "MC.Nose"  "MR.LH"    "MR.NLH"   "MR.Nose"  "PT2.LH"   "PT2.NLH" 
[43] "PT2.Nose" "RP.LH"    "RP.NLH"   "RP.Nose"  "SA.LH"    "SA.NLH"  
[49] "SA.Nose"  "XN.LH"    "XN.NLH"   "XN.Nose"

sample_variables(ps)

[1] "BarcodeSequence"      "LinkerPrimerSequence" "FileInput"           
[4] "Patient"              "Type"                 "Group"

Normalize number of reads in each sample using median sequencing depth.

# NORMALIZE number of reads in each sample using median sequencing depth.
total = median(sample_sums(ps))
# > total [1] 42369
standf = function(x, t = total) round(t * (x/sum(x)))
ps = transform_sample_counts(ps, standf)
ps_rel <- microbiome::transform(ps, "compositional")

saveRDS(ps_rel, "./ps.rds")
hmp.meta <- meta(ps)
hmp.meta$sam_name <- rownames(hmp.meta)

knitr::kable(otu_table(ps_rel)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	P7.Nose	P8.Nose	P9.Nose	P10.Nose	P11.Nose	P12.Nose	P13.Nose	P14.Nose	P15.Nose	P16.Foot	P16.Nose	P17.Foot	P17.Nose	P18.Foot	P18.Nose	P19.Foot	P19.Nose	P20.Foot	P20.Nose	AH.LH	AH.NLH	AH.Nose	AL.LH	AL.NLH	AL.Nose	CB.LH	CB.NLH	CB.Nose	HR.LH	HR.NLH	HR.Nose	KK.LH	KK.NLH	KK.Nose	MC.LH	MC.NLH	MC.Nose	MR.LH	MR.NLH	MR.Nose	PT2.LH	PT2.NLH	PT2.Nose	RP.LH	RP.NLH	RP.Nose	SA.LH	SA.NLH	SA.Nose	XN.LH	XN.NLH	XN.Nose
ST215	0.057120387	0.028166963	0.331595351	0.027787487	0.263751948	0.0366567634	0.075938306	0.326998136	0.0548751183	0.0062770342	0.27594594	0.064053830	0.256921486	0.0056508695	0.3643301341	0.030336571	0.176062799	0.000000000	0.01367549	0.0035244600	0.1479588121	0.2284087268	0.02041916	0.0950307302	0.020855786	0.7842766013	0.981532726	0.878477178	0.1402793353	0.075150738	0.038273888	0.1690600332	0.002213743	0.6278399161	0.2016399	0.15442087	0.011301574	0.063863564	0.061313954	0.230382733	0.007908304	0.0212344528	0.08436858	0.0110539883	0.1157700038	0.8075500277	0.1706887261	0.384180709	0.000000000	0.028545629	0.012758149	0.00000000
ST130	0.007209228	0.000000000	0.041478547	0.037210183	0.094825452	0.0484679017	0.327099809	0.053565187	0.0116944586	0.0081266475	0.14386415	0.025822519	0.000000000	0.0000000000	0.1571734413	0.010325794	0.097942123	0.531130303	0.09342732	0.0000000000	0.0000000000	0.0000000000	0.04458135	0.2384433894	0.011025021	0.0000000000	0.000000000	0.000000000	0.0186272082	0.061897411	0.041914860	0.3114823338	0.003757537	0.0000000000	0.0000000	0.38691872	0.004514804	0.154728962	0.501871459	0.522530001	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.5901285991	0.000000000	0.040371676	0.000000000	0.000000000	0.00000000
ST278	0.000000000	0.684920334	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.1680939066	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST200	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.309336906	0.0040051266	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0026943186	0.0040051849	0.4580699659	0.00000000	0.0000000000	0.000000000	0.0000000000	0.002315691	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.3985902	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.007733535	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.406769294	0.000000000	0.000000000	0.000000000	0.00000000
ST5	0.770935889	0.000000000	0.000000000	0.000000000	0.000000000	0.2836420832	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.003174927	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.3897497888	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.556960182	0.0003932306	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST59	0.000000000	0.050144185	0.003815793	0.521539671	0.001267058	0.0000000000	0.169208139	0.001368985	0.6486273939	0.0004514804	0.03603157	0.010005680	0.038113686	0.0000000000	0.0011068552	0.004383729	0.077902218	0.001514644	0.47549626	0.0000000000	0.0000000000	0.0001747692	0.02468650	0.0000000000	0.145757479	0.0196760945	0.000000000	0.001631179	0.0005388637	0.003058460	0.000000000	0.0000000000	0.007718971	0.0000000000	0.0000000	0.12125890	0.214292996	0.003655588	0.001470952	0.006393452	0.000000000	0.0000000000	0.12108414	0.0000000000	0.0000000000	0.0000000000	0.1284534611	0.025617873	0.047377006	0.469968833	0.940054178	0.30314726
ST83	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0004514870	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0005971193	0.0031458907	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.003393435	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.1912647093	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0003495383	0.0004806152	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST114	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0015728995	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST297	0.005301331	0.008753022	0.223937549	0.000000000	0.039278796	0.0379529302	0.108518664	0.044025982	0.0198354329	0.0040633238	0.01489907	0.000000000	0.022384690	0.0004223588	0.0857812796	0.007864498	0.154042206	0.000000000	0.01128701	0.0019515605	0.0015146881	0.0000000000	0.63158125	0.6657831115	0.438568641	0.1794005418	0.002883691	0.059115668	0.1146760264	0.430369054	0.018714592	0.1210713349	0.018758556	0.0000000000	0.0000000	0.00000000	0.573278185	0.054542542	0.082023797	0.059652802	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0731835195	0.051483332	0.443971920	0.000000000	0.000000000	0.00000000
ST384	0.005257639	0.001995281	0.000000000	0.007471164	0.156008331	0.0000000000	0.045484336	0.000000000	0.0743173378	0.0137774347	0.00000000	0.036672929	0.200821403	0.0000000000	0.2178320202	0.024947934	0.141779415	0.003102107	0.16068334	0.0000000000	0.0000000000	0.2462206169	0.00000000	0.0001165128	0.000305846	0.0007136407	0.000000000	0.000000000	0.0081703392	0.001106871	0.013733743	0.0000000000	0.000000000	0.0000000000	0.0000000	0.15012452	0.028370447	0.002228307	0.000000000	0.002140860	0.000000000	0.0000000000	0.03019093	0.0000000000	0.0000000000	0.0000000000	0.0003204055	0.001325313	0.000000000	0.000000000	0.000000000	0.00000000
ST14	0.000000000	0.000000000	0.218257551	0.000000000	0.000000000	0.0000000000	0.000000000	0.208027496	0.0000000000	0.0000000000	0.40955696	0.000000000	0.000000000	0.0046896391	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0014126968	0.0002330289	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.001878768	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.218257551	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST89	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.3367557963	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST210	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0007864383	0.2077131	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
Novel	0.000451487	0.015772917	0.025501733	0.014286963	0.192825831	0.0093353140	0.004136264	0.020884306	0.0106313260	0.9632261917	0.01102502	0.497764379	0.048643374	0.6077451866	0.0521532703	0.901941366	0.044201389	0.004121579	0.04881814	0.0055051483	0.0013544807	0.0010340509	0.15183292	0.0006262562	0.160307011	0.0019807171	0.000000000	0.000000000	0.6199991262	0.327604789	0.836782547	0.0000000000	0.000000000	0.0013252942	0.1920567	0.06237712	0.004005068	0.650214092	0.328386467	0.159137248	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0006553843	0.0000000000	0.0010049080	0.000000000	0.000000000	0.001004923	0.000000000	0.00000000
ST328	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.355908012	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST331	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.8264516975	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.856441700	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.3911188139	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST73	0.003874050	0.003145845	0.013093123	0.000000000	0.015801815	0.0755854596	0.000000000	0.088867529	0.0323891357	0.0006262470	0.04054644	0.001762281	0.000000000	0.0004369229	0.0525610591	0.003422513	0.019195200	0.000000000	0.02902582	0.0055488400	0.0010486302	0.0403716757	0.00000000	0.0000000000	0.000000000	0.0000000000	0.001194256	0.000000000	0.0000000000	0.000000000	0.000000000	0.0009175381	0.000000000	0.0004077828	0.0000000	0.00000000	0.111180694	0.025370656	0.022311871	0.006058488	0.072223355	0.0000000000	0.03338042	0.0000000000	0.0812967872	0.0803501209	0.0335260621	0.056041827	0.280635577	0.478750983	0.035434447	0.56440878
ST2	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.025371026	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.040953643	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.9777755123	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0249919897	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST88	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0046313153	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST100	0.000597128	0.002534153	0.000000000	0.003160317	0.020491385	0.0004369102	0.003961492	0.001893277	0.0004951576	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0009320886	0.000000000	0.013427901	0.001354441	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.08770627	0.0000000000	0.014170866	0.0067577408	0.000000000	0.000000000	0.0386816772	0.010355073	0.045366500	0.0000000000	0.111109493	0.0000000000	0.0000000	0.00000000	0.000000000	0.042512598	0.000000000	0.009029477	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.064692775	0.000000000	0.000000000	0.000000000	0.00000000
ST10	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0012379482	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST290	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0102387090	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.5026943579	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST87	0.140354199	0.000000000	0.000000000	0.000000000	0.000000000	0.0198939765	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.388913971	0.02739467	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0024030059	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST23	0.000000000	0.000000000	0.000000000	0.000000000	0.145245620	0.0000000000	0.000000000	0.000000000	0.1065462754	0.0000000000	0.00000000	0.004456678	0.002330222	0.3709038478	0.0607605260	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.053566747	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0249191693	0.00000000	0.0009466525	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.003888614	0.00000000
ST218	0.003670152	0.000000000	0.059785616	0.000000000	0.000000000	0.0362052895	0.148628770	0.155714785	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0003204055	0.0000000000	0.0072674842	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.040808599	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.11056901	0.049123982	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0512219277	0.0000000000	0.000000000	0.178381055	0.000000000	0.000000000	0.00000000
ST329	0.000000000	0.033948909	0.000000000	0.002286497	0.044711125	0.0162676220	0.082929174	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.018729155	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0014564097	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.014054353	0.004442049	0.10177242
ST19	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0112577187	0.000000000	0.000000000	0.0000000000	0.0000000000	0.04195916	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.009655855	0.000000000	0.00000000	0.0000000000	0.0024759325	0.0022719991	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0555320915	0.034618858	0.005213871	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST225	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.01433086	0.001325313	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
ST170	0.000000000	0.000000000	0.000000000	0.011432483	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.00000000	0.000000000	0.000000000	0.0000000000	0.0000000000	0.000000000	0.012641452	0.000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.00000000	0.0000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.0000000000	0.00000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.000000000	0.000000000	0.000000000	0.00000000
Unclassified	0.005228511	0.170618392	0.082534735	0.374825236	0.025792639	0.0875422346	0.008724021	0.098654317	0.0405883638	0.0034516406	0.02617168	0.003553691	0.080494590	0.0044566136	0.0073693255	0.016777595	0.234420285	0.066688027	0.14019195	0.0006699387	0.0015729453	0.0147243017	0.03919255	0.0000000000	0.209009350	0.0071946637	0.006801433	0.019967376	0.0034953323	0.002271999	0.000000000	0.0077189712	0.000000000	0.0078789468	0.0000000	0.00000000	0.002606935	0.002883691	0.002621499	0.004674939	0.136917072	0.0346479858	0.73097593	0.9817951444	0.8019282864	0.0603973086	0.0026943186	0.009888878	0.009262765	0.007675279	0.003422563	0.03067154

8.1 Bar plots in P16-P20 samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P16.Foot", "P17.Foot", "P18.Foot", "P19.Foot", "P20.Foot",
    "P16.Nose", "P17.Nose", "P18.Nose", "P19.Nose", "P20.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps_rel, rownames(sample_data(ps_rel)) %in% samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = OTU), stat = "identity") +
    scale_fill_manual(values = c("cornflowerblue", "darkgreen", "seagreen3", "tan",
        "red", "navyblue", "gold", "green", "orange", "pink", "purple", "magenta",
        "brown", "darksalmon", "chocolate4", "darkkhaki", "lightcyan3", "maroon",
        "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet", "darkgoldenrod",
        "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) + labs(fill = "") +
    theme(axis.text.x = element_text(size = 9, angle = 60, hjust = 1), axis.text.y = element_text(size = 9,
        colour = "black"), legend.position = "bottom", legend.text = element_text(size = 8)) +
    guides(fill = guide_legend(nrow = 3))

Regroup together Foot vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
plot_bar(ps_subset_type_, fill = "OTU") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("cornflowerblue", "darkgreen", "seagreen3", "tan",
        "red", "navyblue", "gold", "green", "orange", "pink", "purple", "magenta",
        "brown", "darksalmon", "chocolate4", "darkkhaki", "lightcyan3", "maroon",
        "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet", "darkgoldenrod",
        "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) + labs(fill = "") +
    theme(axis.text.x = element_text(size = 10, angle = 15, hjust = 1), axis.text.y = element_text(size = 9,
        colour = "black"), legend.position = "bottom", legend.text = element_text(size = 8)) +
    guides(fill = guide_legend(nrow = 3))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset, x = "OTU", fill = "OTU", facet_grid = Patient ~ Type) + geom_bar(aes(),
    stat = "identity", position = "stack") + scale_fill_manual(values = c("cornflowerblue",
    "darkgreen", "seagreen3", "tan", "red", "navyblue", "gold", "green", "orange",
    "pink", "purple", "magenta", "brown", "darksalmon", "chocolate4", "darkkhaki",
    "lightcyan3", "maroon", "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet",
    "darkgoldenrod", "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) +
    theme(axis.text = element_text(size = 5, colour = "black"), axis.text.x = element_blank(),
        axis.ticks = element_blank()) + labs(fill = "") + theme(legend.position = "bottom",
    legend.text = element_text(size = 8)) + guides(fill = guide_legend(nrow = 3))

8.2 Bar plots in AH-XN samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("AH.LH", "AL.LH", "CB.LH", "HR.LH", "KK.LH", "MC.LH", "MR.LH",
    "PT2.LH", "RP.LH", "SA.LH", "XN.LH", "AH.NLH", "AL.NLH", "CB.NLH", "HR.NLH",
    "KK.NLH", "MC.NLH", "MR.NLH", "PT2.NLH", "RP.NLH", "SA.NLH", "XN.NLH", "AH.Nose",
    "AL.Nose", "CB.Nose", "HR.Nose", "KK.Nose", "MC.Nose", "MR.Nose", "PT2.Nose",
    "RP.Nose", "SA.Nose", "XN.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps_rel, rownames(sample_data(ps_rel)) %in% samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = OTU), stat = "identity") +
    scale_fill_manual(values = c("cornflowerblue", "darkgreen", "seagreen3", "tan",
        "red", "navyblue", "gold", "green", "orange", "pink", "purple", "magenta",
        "brown", "darksalmon", "chocolate4", "darkkhaki", "lightcyan3", "maroon",
        "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet", "darkgoldenrod",
        "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) + labs(fill = "") +
    theme(axis.text.x = element_text(size = 9, angle = 60, hjust = 1), axis.text.y = element_text(size = 9,
        colour = "black"), legend.position = "bottom", legend.text = element_text(size = 8)) +
    guides(fill = guide_legend(nrow = 3))

Regroup together LH vs NLH vs Nose samples and normalize number of reads in each group using median sequencing depth.

ps_subset_type <- merge_samples(ps_subset, "Type")
ps_subset_type_ = transform_sample_counts(ps_subset_type, function(x) x/sum(x))
plot_bar(ps_subset_type_, fill = "OTU") + geom_bar(aes(), stat = "identity", position = "stack") +
    scale_fill_manual(values = c("cornflowerblue", "darkgreen", "seagreen3", "tan",
        "red", "navyblue", "gold", "green", "orange", "pink", "purple", "magenta",
        "brown", "darksalmon", "chocolate4", "darkkhaki", "lightcyan3", "maroon",
        "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet", "darkgoldenrod",
        "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) + labs(fill = "") +
    theme(axis.text.x = element_text(size = 10, angle = 15, hjust = 1), axis.text.y = element_text(size = 9,
        colour = "black"), legend.position = "bottom", legend.text = element_text(size = 8)) +
    guides(fill = guide_legend(nrow = 3))

Do separate panels Patient and Type.

#-- If existing replicates, to be processed as follows --
plot_bar(ps_subset, x = "OTU", fill = "OTU", facet_grid = Patient ~ Type) + geom_bar(aes(),
    stat = "identity", position = "stack") + scale_fill_manual(values = c("cornflowerblue",
    "darkgreen", "seagreen3", "tan", "red", "navyblue", "gold", "green", "orange",
    "pink", "purple", "magenta", "brown", "darksalmon", "chocolate4", "darkkhaki",
    "lightcyan3", "maroon", "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet",
    "darkgoldenrod", "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) +
    theme(axis.text = element_text(size = 5, colour = "black"), axis.text.x = element_blank(),
        axis.ticks = element_blank()) + labs(fill = "") + theme(legend.position = "bottom",
    legend.text = element_text(size = 8)) + guides(fill = guide_legend(nrow = 3))

8.3 Bar plots in P7-P15 samples

# Assuming you have a vector of sample IDs that you are interested in.
samples_to_keep <- c("P7.Nose", "P8.Nose", "P9.Nose", "P10.Nose", "P11.Nose", "P12.Nose",
    "P13.Nose", "P14.Nose", "P15.Nose")

# Subsetting the phyloseq object.
ps_subset <- subset_samples(ps_rel, rownames(sample_data(ps_rel)) %in% samples_to_keep)
# Convert to a data frame for manipulation with dplyr and ggplot
ps_data <- psmelt(ps_subset)
# Assuming 'SampleID' is the column with your sample names
ps_data$Sample <- factor(ps_data$Sample, levels = samples_to_keep)

# Now plot using ggplot, with fill determined by 'Class'
ggplot(ps_data, aes(x = Sample, y = Abundance)) + geom_bar(aes(fill = OTU), stat = "identity") +
    scale_fill_manual(values = c("cornflowerblue", "darkgreen", "seagreen3", "tan",
        "red", "navyblue", "gold", "green", "orange", "pink", "purple", "magenta",
        "brown", "darksalmon", "chocolate4", "darkkhaki", "lightcyan3", "maroon",
        "lightgreen", "blue", "cyan", "skyblue2", "azure3", "blueviolet", "darkgoldenrod",
        "tomato", "mediumpurple4", "indianred", "orange", "darkgrey")) + labs(fill = "") +
    theme(axis.text.x = element_text(size = 9, angle = 60, hjust = 1), axis.text.y = element_text(size = 9,
        colour = "black"), legend.position = "bottom", legend.text = element_text(size = 8)) +
    guides(fill = guide_legend(nrow = 3))