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("../map_corrected.txt", sep = "\t", row.names = 1)
knitr::kable(map_corrected) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

BarcodeSequence	LinkerPrimerSequence	FileInput	Group	Sex_age	pre_post_stroke	Conc	Vol_50ng	Vol_PCR	Description
NA	NA	6341_S1_R1.fastq.gz_merged.fasta	Group1	f.aged	post	19.2	2.60	2.60	PCR1
NA	NA	6340_S2_R1.fastq.gz_merged.fasta	Group1	f.aged	post	16.2	3.09	3.09	PCR1
NA	NA	6342_S3_R1.fastq.gz_merged.fasta	Group1	f.aged	post	6.3	7.94	5.00	PCR1
NA	NA	8129_S4_R1.fastq.gz_merged.fasta	Group1	f.aged	post	28.3	1.77	1.77	PCR1
NA	NA	8130_S5_R1.fastq.gz_merged.fasta	Group1	f.aged	post	7.3	6.85	5.00	PCR1
NA	NA	8128_S6_R1.fastq.gz_merged.fasta	Group1	f.aged	post	11.1	4.50	4.50	PCR1
NA	NA	9989_S7_R1.fastq.gz_merged.fasta	Group1	f.aged	post	19.7	2.54	2.54	PCR1
NA	NA	6341_S8_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	18.0	2.78	2.78	PCR1
NA	NA	6340_S9_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	57.9	0.86	0.86	PCR2
NA	NA	6342_S10_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	17.0	2.94	2.94	PCR2
NA	NA	8129_S11_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	14.2	3.52	3.52	PCR2
NA	NA	8130_S12_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	19.5	2.56	2.56	PCR2
NA	NA	8128_S13_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	35.1	1.42	1.42	PCR2
NA	NA	9989_S14_R1.fastq.gz_merged.fasta	Group2	f.aged	pre	7.9	6.33	5.00	PCR2
NA	NA	16880_S15_R1.fastq.gz_merged.fasta	Group3	f.young	post	8.2	6.10	5.00	PCR2
NA	NA	16681_S16_R1.fastq.gz_merged.fasta	Group3	f.young	post	7.2	6.94	5.00	PCR2
NA	NA	16685_S17_R1.fastq.gz_merged.fasta	Group3	f.young	post	13.1	3.82	3.82	PCR3
NA	NA	16686_S18_R1.fastq.gz_merged.fasta	Group3	f.young	post	13.2	3.79	3.79	PCR3
NA	NA	16819_S19_R1.fastq.gz_merged.fasta	Group3	f.young	post	16.3	3.07	3.07	PCR3
NA	NA	21909_S20_R1.fastq.gz_merged.fasta	Group3	f.young	post	2.6	19.23	5.00	PCR3
NA	NA	16880_S21_R1.fastq.gz_merged.fasta	Group4	f.young	pre	17.0	2.94	2.94	PCR3
NA	NA	16681_S22_R1.fastq.gz_merged.fasta	Group4	f.young	pre	17.6	2.84	2.84	PCR3
NA	NA	16685_S23_R1.fastq.gz_merged.fasta	Group4	f.young	pre	20.2	2.48	2.48	PCR3
NA	NA	16686_S24_R1.fastq.gz_merged.fasta	Group4	f.young	pre	29.0	1.72	1.72	PCR3
NA	NA	16819_S25_R1.fastq.gz_merged.fasta	Group4	f.young	pre	27.0	1.85	1.85	PCR4
NA	NA	16684_S26_R1.fastq.gz_merged.fasta	Group4	f.young	pre	4.3	11.63	5.00	PCR4
NA	NA	21908_S27_R1.fastq.gz_merged.fasta	Group4	f.young	pre	18.0	2.78	2.78	PCR4
NA	NA	21909_S28_R1.fastq.gz_merged.fasta	Group4	f.young	pre	9.4	5.32	5.00	PCR4
NA	NA	4896_S29_R1.fastq.gz_merged.fasta	Group5	m.aged	post	25.0	2.00	2.00	PCR4
NA	NA	4897_S30_R1.fastq.gz_merged.fasta	Group5	m.aged	post	25.2	1.98	1.98	PCR4
NA	NA	4900_S31_R1.fastq.gz_merged.fasta	Group5	m.aged	post	15.7	3.18	3.18	PCR4
NA	NA	9976_S32_R1.fastq.gz_merged.fasta	Group5	m.aged	post	34.0	1.47	1.47	PCR4
NA	NA	4896_S33_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	44.0	1.14	1.14	PCR5
NA	NA	4897_S34_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	41.0	1.22	1.22	PCR5
NA	NA	4900_S35_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	28.5	1.75	1.75	PCR5
NA	NA	4898_S36_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	96.0	0.52	0.52	PCR5
NA	NA	5114_S37_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	21.9	2.28	2.28	PCR5
NA	NA	9975_S38_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	15.7	3.18	3.18	PCR5
NA	NA	9976_S39_R1.fastq.gz_merged.fasta	Group6	m.aged	pre	6.5	7.69	5.00	PCR5
NA	NA	16888_S40_R1.fastq.gz_merged.fasta	Group7	m.young	post	38.0	1.32	1.32	PCR5
NA	NA	16625_S41_R1.fastq.gz_merged.fasta	Group7	m.young	post	7.8	6.41	5.00	PCR6
NA	NA	16824_S42_R1.fastq.gz_merged.fasta	Group7	m.young	post	42.1	1.19	1.19	PCR6
NA	NA	16826_S43_R1.fastq.gz_merged.fasta	Group7	m.young	post	18.1	2.76	2.76	PCR6
NA	NA	16827_S44_R1.fastq.gz_merged.fasta	Group7	m.young	post	9.6	5.21	5.00	PCR6
NA	NA	21911_S45_R1.fastq.gz_merged.fasta	Group7	m.young	post	29.6	1.69	1.69	PCR6
NA	NA	21914_S46_R1.fastq.gz_merged.fasta	Group7	m.young	post	62.3	0.80	0.80	PCR6
NA	NA	16888_S47_R1.fastq.gz_merged.fasta	Group8	m.young	pre	13.0	3.85	3.85	PCR6
NA	NA	16625_S48_R1.fastq.gz_merged.fasta	Group8	m.young	pre	43.1	1.16	1.16	PCR6
NA	NA	16824_S49_R1.fastq.gz_merged.fasta	Group8	m.young	pre	13.2	3.79	3.79	PCR7
NA	NA	16826_S50_R1.fastq.gz_merged.fasta	Group8	m.young	pre	32.2	1.55	1.55	PCR7
NA	NA	5115_S51_R1.fastq.gz_merged.fasta	Group8	m.young	pre	33.2	1.51	1.51	PCR7
NA	NA	16827_S52_R1.fastq.gz_merged.fasta	Group8	m.young	pre	12.8	3.91	3.91	PCR7
NA	NA	16691_S53_R1.fastq.gz_merged.fasta	Group8	m.young	pre	20.7	2.42	2.42	PCR7
NA	NA	21911_S54_R1.fastq.gz_merged.fasta	Group8	m.young	pre	9.8	5.10	5.00	PCR7
NA	NA	21914_S55_R1.fastq.gz_merged.fasta	Group8	m.young	pre	6.8	7.35	5.00	PCR7

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
source("https://bioconductor.org/biocLite.R")
biocLite("phyloseq")

# library('phyloseq') library('ggplot2') # graphics library('readxl') #
# necessary to import the data from Excel file library('dplyr') # filter and
# reformat data frames
library("ggplot2")
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  
library("RColorBrewer")  # nice color options

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_even42434.biom", "../clustering/rep_set.tre")
sample <- read.csv("../map_corrected.txt", sep = "\t", row.names = 1)
# MODIFIED > rownames(SAM) [1] 'sa1' 'sa2' 'sa3' 'sa4' 'sa5' 'sa6' 'sa7' 'sa8'
# 'sa9' 'sa10' [11] 'sa11' 'sa12' 'sa13' 'sa14' 'sa15' 'sa16' 'sa17' 'sa18'
# 'sa19' 'sa20' [21] 'sa21' 'sa22' 'sa23' 'sa24' 'sa25' 'sa26' 'sa27' 'sa28'
# 'sa29' 'sa30' [31] 'sa31' 'sa32' 'sa33' 'sa34' 'sa35' 'sa36' 'sa37' 'sa38'
# 'sa39' 'sa40' [41] 'sa41' 'sa42' 'sa43' 'sa44' 'sa45' 'sa46' 'sa47' 'sa48'
# 'sa49' 'sa50' [51] 'sa51' 'sa52' 'sa53' 'sa54' 'sa55'
SAM = sample_data(sample, errorIfNULL = T)
rownames(SAM) <- c("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")
ps.ng.tax <- merge_phyloseq(ps.ng.tax, SAM)
print(ps.ng.tax)

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

# otu_table(ps.ng.tax) 17 53 49 26 36 10 51 8 14 29 20 31 24 34 2 35 27 21 13
# 23 38 9 32 37 48 52 33 25 12 39 22 18 6 11 50 55 1 44 43 46 16 47 28 5 41 15
# 30 42 7 45 4 19 54 40 ID01 ID03 ID11 ID12 ID02 ID16 ID15 ID04 ID05 ID08 ID10
# ID09 ID07 ID14 ID18 ID06 ID17 > colnames(sample_data(ps.ng.tax)) [1]
# 'BarcodeSequence' 'LinkerPrimerSequence' 'FileInput' [4] 'Group' 'Sex_age'
# 'pre_post_stroke' [7] 'Conc' 'Vol_50ng' 'Vol_PCR' [10] 'Description'
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] "17" "53" "49" "26" "36" "10" "51" "8"  "14" "29" "20" "31" "24" "34" "2" 
[16] "35" "27" "21" "13" "23" "38" "9"  "32" "37" "48" "52" "33" "25" "12" "39"
[31] "22" "18" "6"  "11" "50" "55" "1"  "44" "43" "46" "16" "47" "28" "5"  "41"
[46] "15" "30" "42" "7"  "45" "4"  "19" "54" "40"

rank_names(ps.ng.tax)

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

sample_variables(ps.ng.tax)

 [1] "BarcodeSequence"      "LinkerPrimerSequence" "FileInput"           
 [4] "Group"                "Sex_age"              "pre_post_stroke"     
 [7] "Conc"                 "Vol_50ng"             "Vol_PCR"             
[10] "Description"

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 = 42434)
# total <- 42434

# NORMALIZE number of reads in each sample using median sequencing depth.
total = median(sample_sums(ps.ng.tax))
# > total [1] 42434
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 166 OTUS which makes the reading much more easy.

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"))

	17	53	49	26	36	10	51	8	14	29	20	31	24	34	2	35	27	21	13	23	38	9	32	37	48	52	33	25	12	39	22	18	6	11	50	55	1	44	43	46	16	47	28	5	41	15	30	42	7	45	4	19	54	40
EU505095.1.1391	0	0	443	0	0	19	0	0	0	0	0	0	0	0	0	0	0	0	9	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	69	0	0	0	0	0	0
HK240365.1.1492	3	13	4455	8	3713	2	6824	5	8144	4375	1	3778	5	5516	2	4897	5	4	33	11	46	6	11	8747	9	1043	1548	3	2	6	6043	2	1094	40	2022	1	2	9440	7671	3	11535	2	2	3	8	0	4566	7614	6920	10	1179	2	5	16
EU791193.1.1502	2189	335	210	199	52	415	339	731	273	166	57	571	89	93	1538	92	168	520	240	408	185	456	526	97	1224	262	333	251	389	267	262	774	661	1220	196	7	1296	206	217	622	111	318	1	770	1682	267	464	274	750	0	1531	702	0	125
EF406536.1.1499	368	102	159	101	179	285	600	647	381	622	298	865	264	669	135	454	192	137	353	370	659	194	1246	350	615	833	636	234	345	537	377	450	194	452	188	192	684	226	193	247	342	222	251	14	33	99	1410	490	1334	379	892	133	127	183
EU791223.1.1512	159	43	60	36	101	96	179	234	139	211	101	256	101	318	47	205	67	66	122	192	283	65	405	160	245	301	303	51	150	248	92	203	70	168	46	71	255	54	64	85	120	61	118	9	11	36	431	139	462	116	292	28	39	104
JX198570.1.1508	0	0	0	0	4	378	19	61	8	4	4	30	0	59	269	5	0	0	2	0	173	16	79	7	0	0	1	0	58	32	0	0	2	432	0	0	119	0	1	0	0	0	1	143	16	0	28	0	1	0	3	0	0	56
CVUG01000013.3917.5427	33	8	12	8	10	37	10	0	23	26	71	15	5	69	27	40	105	0	95	2	205	60	271	15	0	47	58	247	43	187	30	12	17	195	42	17	35	32	29	53	81	7	65	119	7	121	14	29	86	56	200	585	23	9
HL966144.1.1491	0	7	3	11	70	99	9	0	11	0	26	28	0	94	0	22	22	5	0	0	0	398	427	3	23	12	0	33	0	54	5	0	0	0	10	4	6	2	1	18	9	6	48	0	10	30	51	0	52	49	0	167	8	2
AB606239.1.1494	722	27	24	38	120	97	64	5	14	238	146	365	45	193	40	720	107	0	7	40	24	179	98	86	6	32	313	92	84	14	14	200	10	1	33	91	403	23	9	199	44	51	45	138	34	445	196	2	118	631	10	177	178	18
AB622814.1.1518	39	55	34	5	11	18	14	1	9	7	109	12	16	14	3	14	260	4	6	2	29	77	36	13	16	60	5	153	84	39	29	4	0	9	92	2	10	18	5	46	94	24	61	103	98	104	8	3	22	100	13	522	66	36
FJ880076.1.1490	123	14	6	0	181	11	32	34	9	0	2	21	0	37	78	180	4	15	2	14	24	3	85	54	11	1	12	6	7	820	42	0	11	10	35	0	11	2	13	2	1	0	1	25	22	0	3	19	27	1	11	1	1	7
JQ084417.1.1495	64	18	3	10	1	9	1	1	5	163	23	9	31	32	0	22	26	1	9	16	27	14	19	2	32	12	546	24	24	0	3	26	12	7	20	3	3	8	13	75	17	17	30	18	33	51	11	12	3	28	15	117	17	163
JQ084693.1.1491	143	76	71	42	429	215	75	17	41	131	130	1	331	178	197	0	114	24	214	33	194	72	147	13	186	50	163	634	258	45	42	136	21	92	49	52	174	42	13	196	37	23	204	167	194	26	90	43	22	23	88	227	74	106
AB606325.1.1494	5	88	16	0	30	14	61	0	0	20	232	24	14	86	37	96	231	12	88	0	38	15	28	11	6	52	99	330	242	8	62	8	123	199	222	8	61	57	108	54	70	22	91	533	10	198	47	394	20	111	259	376	46	28
ASTC01000013.23144.24660	629	217	33	29	257	88	316	2	7	70	54	26	15	58	148	10	29	52	21	29	32	131	183	46	30	24	47	45	31	26	14	27	28	15	43	15	390	30	2	141	8	104	30	17	95	53	89	8	21	41	21	12	68	1025
JQ084374.1.1490	3	0	0	0	1325	56	19	2	4	18	6	1	3	127	136	0	9	6	0	12	18	235	2	9	80	0	106	32	61	3	1	9	0	4	0	6	39	0	0	16	0	0	7	29	6	0	19	0	0	1	1	0	8	7
EF602810.1.1476	220	4	1	198	1	172	5	35	48	9	4	1	144	8	73	1	5	17	690	36	3	14	2	0	51	3	82	41	66	0	0	20	1	3	2	2	27	1	4	45	1	0	5	4	9	11	7	3	33	1	2	0	4	3
HL953771.1.1511	316	161	37	39	47	57	38	2	53	221	301	37	230	55	33	34	453	256	103	46	126	95	1414	14	264	154	274	459	94	398	30	68	42	79	129	49	174	49	23	277	42	66	93	59	783	125	54	44	411	152	129	0	125	1373
AB969468.1.1539	2	4	81	2	2829	299	1152	2	16	2611	28	1128	3	2719	807	3943	126	34	18	1	2	4	1287	702	391	0	622	17	48	1414	80	4	4006	1717	115	150	312	4	2	34	1	29	16	109	90	5	903	383	180	1	1483	26	3	5
ASSR01000026.1887.3427	532	185	177	1987	41	2	2	1	2	1	52	0	37	2	3	3	746	830	0	1072	0	6	5	2	91	407	0	5	2083	4	1	5	5	6	853	84	3	7	1	131	0	87	592	514	101	3	9	17	0	13	8	0	390	5
New.ReferenceOTU515	1379	2776	690	6	128	8	3	27	3	12	2	224	6613	416	921	477	10	6	291	576	789	4342	151	3	7	2916	179	3708	2	1	656	380	24	201	1033	1	32	1337	1182	2	286	1454	5	2	2	1792	106	1	4	1	109	1	9	1
EF406647.1.1503	106	186	52	0	14	3	0	2	1	4	0	12	550	29	76	45	1	0	35	36	76	321	16	0	1	228	21	273	1	2	45	30	3	15	95	1	3	86	60	1	24	88	0	1	0	120	7	0	2	0	9	0	0	0
ASSX01000038.243.1764	327	9	5	0	0	0	903	0	32	1	15	0	0	0	0	1	4	6	82	813	0	0	0	660	56	0	0	0	23	0	0	0	0	465	0	39	0	3	1	20	0	0	82	7	8	0	0	0	2	4	248	2	0	2
HK693730.8.1517	1	4	2	0	0	460	0	12	0	0	0	1	13	1	10	0	0	1	3	0	5	27	1	0	0	6	0	12	3	0	2	2	0	3	1	0	7	4	0	0	2	2	0	4	1	3	4	0	0	0	3	0	0	0
New.ReferenceOTU186	3	39	5	14	8	1050	44	40	0	8	7	10	1	62	52	1	83	26	52	0	175	50	8	36	20	2	0	26	58	3	9	5	11	8	39	4	23	1	4	22	1	6	5	100	10	62	49	12	1	5	13	3	16	4
AB606268.1.1497	0	186	3	0	56	508	203	26	0	2	34	30	0	257	10	0	466	31	0	1	844	24	8	192	8	0	1	4	5	19	42	2	0	1	58	2	10	0	0	1	0	9	32	3	25	346	125	6	2	0	0	0	0	2
New.ReferenceOTU16	0	110	1	0	31	284	154	16	0	0	20	19	0	166	3	0	258	25	1	0	569	12	5	98	6	0	1	2	6	13	12	2	0	2	36	0	8	0	0	0	0	9	9	2	25	248	101	6	0	0	1	0	0	0
JPNB01000002.737568.739079	2	217	4	4	43	664	323	19	0	2	46	51	4	266	15	0	728	48	3	1	830	16	23	158	4	2	3	15	18	36	44	2	0	6	83	0	20	1	5	19	3	29	30	7	35	459	233	13	0	5	1	0	13	6
New.ReferenceOTU336	342	162	50	110	110	303	112	38	59	339	280	390	265	233	128	173	209	32	566	90	86	176	267	33	46	315	3888	214	1105	22	78	236	150	290	145	77	290	89	149	336	26	148	139	416	55	373	140	109	104	131	413	191	72	59
New.ReferenceOTU471	0	5	18	14	1	69	0	9	0	3	2	2	0	4	79	0	1	3	73	1	1	15	7	5	11	9	0	3	168	1	2	6	15	17	0	0	8	1	12	6	3	0	3	467	2	2	3	23	4	9	32	1	6	1
ACTP02000009.47.1566	6	3	0	636	0	1	0	3	0	0	11	0	5	3	4	0	0	1	256	3	3	0	3	1	2	1	0	4	58	0	0	19	0	0	4	1	2	0	0	6	0	0	0	21	6	0	1	1	1	8	0	0	44	0
AB627556.1.1524	1	277	133	0	143	54	58	195	272	57	1349	13	0	2	49	6	655	1	0	1	333	9	172	4	0	131	1	0	2	23	103	0	2	1	342	1	24	178	2	4	0	338	129	0	0	398	18	8	4	0	0	2	2	0
New.ReferenceOTU153	116	45	14	1	37	23	6	13	27	7	181	4	16	0	3	1	91	8	95	3	68	0	36	3	2	22	0	79	0	3	11	18	1	25	48	165	2	21	0	254	0	24	625	2	10	69	9	0	1	1	11	6	528	1
FJ879443.1.1488	0	0	0	655	0	0	0	0	0	0	0	0	0	0	0	0	0	25	0	11	0	0	0	0	49	0	0	0	0	0	0	1	0	0	0	4	0	0	0	4	0	0	0	0	17	0	0	0	0	0	0	0	0	0
JQ083995.1.1490	6	0	0	1160	0	0	0	0	0	0	0	0	0	0	0	0	0	54	0	18	0	0	0	0	78	0	0	0	0	0	0	1	0	0	1	1	0	0	0	2	0	0	0	0	20	0	0	0	0	0	0	0	0	0
AB606324.1.1498	20	52	66	12	1	12	7	2	24	159	10	99	85	15	6	6	71	5	158	0	5	14	29	5	22	29	104	123	71	9	32	36	5	3	26	77	116	33	34	96	10	89	14	84	36	146	70	105	44	129	10	532	34	13
AB606364.1.1499	2	7	10	0	0	13	3	2	0	0	0	24	0	0	28	1	1	44	0	0	9	0	10	0	14	0	0	107	0	1	3	0	0	1	0	1	0	2	4	15	1	0	6	0	838	0	10	1	2	49	0	7	55	35
AB626913.1.1497	0	0	0	0	0	0	0	0	0	0	30	0	1	0	0	0	11	0	0	0	0	0	813	0	0	0	0	0	0	32	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
AB626943.1.1505	0	11	11	3	11	14	77	44	157	308	0	11	174	96	45	47	39	2	20	0	7	69	5	2	19	5	77	13	20	0	0	159	4	14	26	8	85	7	1	6	0	50	8	50	78	23	177	43	566	14	11	241	9	50
New.ReferenceOTU371	1	0	0	4	0	1	0	0	0	0	0	0	2	0	1	0	0	4	0	0	1824	0	29	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	12
AATF01001259.61.1571	6	17	32	506	0	16	66	59	0	35	14	4	6	5	26	7	33	189	0	129	4	145	28	8	4	11	6	30	37	1	54	32	2	40	92	22	304	0	20	14	1	8	21	59	60	7	5	11	9	56	23	7	39	43
AB606370.1.1495	385	10	25	747	2	715	3	187	85	287	32	385	754	66	648	12	93	33	1763	237	2	170	31	8	195	33	51	20	92	13	41	256	25	155	88	24	457	10	57	151	3	38	41	132	274	208	135	67	163	5	293	200	56	113
JN012616.1.1459	0	10	118	0	0	1	0	0	1	0	12	0	0	0	0	0	0	420	0	0	0	0	0	0	501	0	0	0	0	0	130	0	0	0	164	6	0	4	0	85	0	0	0	0	133	0	0	62	0	0	0	0	0	230
EF603872.1.1478	594	2	41	1	9	0	0	0	0	1	8	0	0	0	0	0	3	6	0	1108	0	0	0	5	0	0	0	0	591	1	0	0	0	0	2	1	0	0	0	29	0	0	61	171	0	0	0	2	0	0	0	0	9	0
AQFU01000029.6100.7620	73	12	13	96	14	46	38	206	92	11	37	71	26	242	82	8	59	12	176	50	149	175	115	59	14	82	525	41	319	30	30	13	68	108	18	2	26	22	9	8	29	24	57	565	7	94	127	32	79	157	106	77	20	119
AB969484.1.1523	0	80	82	24	734	39	27	4	6	67	101	36	60	133	12	21	71	17	83	0	155	34	98	3	39	1	107	115	22	13	100	51	16	86	41	10	46	10	2	109	0	43	13	35	53	96	35	177	9	31	23	83	8	50
EF096325.1.1399	104	264	77	91	0	17	19	89	19	48	2	5	0	0	236	1	10	440	667	285	21	319	128	21	252	27	1	48	550	3	82	2	0	20	52	105	105	0	10	107	1	3	110	727	903	1	95	2	74	572	76	15	197	61
DQ777894.1.1513	347	54	32	14	81	250	33	1110	229	112	92	27	376	242	467	24	194	147	679	58	238	2554	458	37	174	178	985	423	63	50	54	232	479	611	22	9	392	193	76	38	168	39	43	21	199	1023	152	152	404	205	560	228	43	1039
EU622695.1.1497	943	2055	574	0	3088	2	1432	2214	950	1178	388	655	327	1149	4347	2361	2680	2994	2137	17	1056	0	636	1424	1734	3645	1124	2825	2207	768	3164	1296	3272	3679	2108	36	0	3610	5342	1569	1828	3683	12	5694	3001	2631	366	174	150	14	2723	396	118	129
AB702763.1.1498	0	0	0	0	0	927	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	309	0	0	0	0	0	1	0	0	0	1	0	0	0	0	429	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
EU622762.1.1520	47	164	0	0	0	7	0	282	3	1	0	0	4	0	188	0	110	1157	1	2	42	0	121	0	20	67	0	7	0	299	300	8	0	0	292	1	1	596	4	0	14	440	0	0	284	480	1	1	23	0	0	42	0	94
JQ084406.1.1485	0	440	162	1	2	0	1	0	0	0	601	0	0	0	0	0	363	0	449	1	627	1	181	1	0	105	0	549	475	294	214	0	259	622	177	436	0	318	523	131	169	378	932	529	0	228	0	267	1	543	42	30	717	0
New.ReferenceOTU511	0	131	1	0	0	0	0	0	0	0	468	1	1	0	0	0	3	0	0	0	0	0	1	0	1	191	0	179	0	0	109	0	71	1	23	8	0	367	1390	18	745	2	11	13	0	119	0	4	0	5	27	466	41	0
AB021157.1.1508	0	752	339	0	0	999	527	455	372	0	642	1	2	1	2131	1	107	0	17	1	311	245	4170	418	1	1077	0	490	3	1144	692	0	1581	995	286	3	1848	96	403	542	1195	891	229	2	0	1406	0	3212	2376	1902	5309	515	251	0
BAIW01000188.173.1680	1	15	48	0	327	36	32	121	71	1817	37	2079	0	410	150	600	4	2	21	3	140	20	144	60	0	31	329	63	2	38	7	0	835	104	5	16	1547	0	22	104	36	23	42	2	1	35	1474	4	963	81	73	539	27	0
CP015401.3482423.3483958	522	98	120	20	74	121	161	265	117	76	36	191	376	67	330	75	83	99	128	115	75	210	84	98	176	109	135	99	74	156	101	1242	435	112	45	71	215	69	259	391	161	271	83	372	1529	348	92	48	63	67	16	156	258	198
LT622246.247875.249393	0	55	146	1	938	248	118	627	421	7905	161	11080	0	1374	868	1953	51	1	149	9	487	157	705	220	1	115	1153	370	37	106	51	2	4352	590	34	54	9319	10	157	796	48	137	194	16	2	150	7669	38	4938	449	392	2634	134	3
AB559603.1.1522	0	92	112	0	15	101	160	23	66	149	50	234	1	66	257	20	8	184	145	0	53	4	376	159	234	209	89	17	5	411	231	1	100	122	44	118	123	167	308	271	97	166	32	57	111	163	339	87	560	980	306	119	103	161
AB627719.1.1500	0	278	0	0	1	141	99	0	199	262	0	57	0	42	104	5	6	26	111	0	32	27	143	38	0	0	54	1	0	20	1	0	58	484	0	0	0	0	1	0	4	0	0	0	0	0	57	0	5	0	335	0	2	155
AB606293.1.1476	0	305	1	0	1	0	0	0	0	0	197	0	0	0	0	0	111	44	0	0	0	0	0	0	0	18	0	90	0	0	30	0	0	0	249	0	0	651	0	382	3	0	256	0	0	0	0	42	0	0	0	11	0	72
EF406551.1.1491	665	16	3	243	20	4	14	3	12	20	10	52	134	17	42	3	18	3	8	113	27	38	27	2	98	12	33	26	27	0	8	119	2	0	42	4	161	6	2	21	17	4	37	14	47	30	41	41	62	55	11	101	14	11
JQ084431.1.1502	13	3	2	10	249	1	40	9	135	17	62	17	6	5	134	126	93	52	114	5	3	2	7	70	0	33	149	31	15	86	125	4	355	173	15	3	2	6	5	12	24	14	40	147	116	36	0	4	93	202	129	110	20	667
JX218602.1.1508	1	1	75	0	0	442	1	74	203	0	109	0	0	0	59	3	3	0	377	0	0	0	0	0	0	35	0	0	7	1	26	0	1229	545	0	0	1	13	2	8	0	335	197	63	0	69	1	43	134	556	825	142	479	0
EU510524.1.1382	709	0	0	0	1	331	0	125	6	3	0	4	0	21	1	53	0	0	0	197	14	0	36	279	0	0	0	0	0	1640	0	0	0	0	1	0	1	0	0	0	0	1	0	0	0	0	49	0	2	0	0	0	0	0
HL953748.1.1486	0	1	0	0	143	0	124	49	51	0	0	18	0	0	0	18	14	25	0	0	0	0	21	97	0	4	0	0	73	236	7	0	318	270	8	0	0	386	1	0	0	0	0	572	0	0	24	0	6	0	707	0	0	18
AB606249.1.1546	276	265	18	33	177	318	231	99	85	442	318	340	174	280	191	136	234	16	336	65	374	254	599	104	128	131	678	423	373	68	27	50	388	48	167	23	528	64	60	362	50	47	113	237	170	272	188	472	193	63	275	824	76	8
New.ReferenceOTU132	0	0	0	0	0	0	0	0	332	0	621	0	0	0	0	0	0	0	0	0	138	0	114	0	0	0	0	0	0	761	0	0	0	0	0	0	0	0	0	0	0	0	136	0	0	0	1	0	266	0	0	0	0	0
New.ReferenceOTU704	0	0	0	1	350	0	254	0	0	444	0	11	0	33	0	132	1	1	0	0	0	0	0	78	0	0	66	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	50	0	0	0	0	0	0	0
CDZJ01005172.1679.3161	0	11	34	0	0	0	42	0	8	0	0	1	0	0	0	0	5	136	0	0	2	0	38	11	10	33	0	9	8	230	262	0	16	24	34	0	0	384	329	0	389	149	0	17	469	172	0	321	258	1	780	486	0	80
AB606338.1.1472	5	0	197	0	5	0	12	0	81	0	24	0	1	17	0	38	22	4	0	0	1	0	1	21	3	134	5	51	0	0	71	47	0	0	313	0	0	70	267	0	730	0	2	1	166	0	72	162	12	82	0	43	22	104
EF406549.1.1479	1754	0	9	1	0	0	4	0	12	0	19	1	238	1	0	4	58	2	5	998	0	5	2	0	17	25	3	2	692	3	1	325	1	44	0	20	4	9	60	14	0	0	14	4	239	3	1	26	0	10	3	0	33	19
AJ308386.1.1462	8	3	0	23	7	3	0	0	0	0	10	0	82	1	0	1	8	4	0	14	0	0	1	0	16	12	0	65	1123	0	11	0	0	1	11	1	0	6	0	0	0	0	7	226	2	1	1	1	0	2	1	0	52	1
AB606393.1.1460	75	15	8	2	0	0	0	1	0	1	7	0	566	0	0	0	57	0	0	6	1	0	0	0	89	3	0	76	0	0	1	0	0	0	20	3	0	6	0	0	0	14	50	1	1	64	0	1	0	1	0	0	29	2
AB702909.1.1496	233	2	1	1096	58	1275	615	961	0	261	0	456	55	368	293	247	0	1	548	308	0	388	41	2331	0	0	112	0	1222	241	0	4	7	39	1	0	44	1	0	0	0	0	1	144	0	0	364	2	0	0	75	0	0	0
EU791043.1.1531	138	48	29	541	6	80	12	109	70	576	172	42	47	10	669	6	128	83	26	178	212	29	50	11	47	51	11	9	26	165	27	661	48	22	12	69	191	19	39	130	31	108	19	118	158	50	41	6	121	85	52	98	362	45
EU790984.1.1546	489	103	100	2029	19	287	28	417	245	1361	517	146	156	32	2505	19	436	167	78	357	679	115	181	25	121	132	26	49	101	450	87	2014	168	128	63	197	717	57	123	430	55	384	69	399	501	110	133	23	537	221	224	254	1372	111
AATF01001055.137.1680	187	46	27	667	6	110	17	133	76	575	208	46	61	13	814	11	156	77	21	209	311	49	60	13	49	76	16	23	18	221	45	778	48	49	18	64	244	23	45	120	53	112	32	142	198	49	43	6	185	103	87	122	529	57
LG042641.312212.313771	236	0	653	326	11	275	102	968	575	79	179	409	181	42	96	265	236	949	148	461	124	872	48	259	191	0	1	176	158	471	61	5	434	116	0	235	0	0	0	1	2	0	411	2	84	34	5	27	19	0	1	11	0	385
AB559728.1.1569	2	159	626	123	40	85	38	479	261	1	161	233	53	41	2	135	177	400	83	224	100	367	29	95	84	21	1	97	66	39	159	3	165	28	83	100	0	118	1	7	2	0	241	0	8	1	4	30	4	0	1	3	0	44
AB727348.1.1509	153	0	0	98	23	347	13	1422	1573	0	406	45	13	27	26	32	1076	87	59	32	120	93	12	25	17	1	2	0	372	295	1	15	29	108	1	0	408	3	0	2	0	0	451	3	15	0	7	0	430	0	230	0	2	1344
U07573.1.1470	22	0	0	385	12	78	48	0	8	80	90	170	122	8	11	29	24	0	36	41	79	150	38	70	0	0	481	1	69	8	0	7	10	33	1	12	12	0	0	90	0	0	30	45	53	0	53	0	23	120	12	0	87	17
CCPS01000022.154.1916	7	0	5	0	0	0	0	0	1	6	0	514	0	0	8	0	1	0	0	0	6	2	0	0	0	2	2	0	0	2	0	8	0	0	0	1	1123	5	0	1	34	0	5	6	142	0	248	713	31	37	25	1	5	2
AB606354.1.1519	96	321	161	115	692	260	98	295	62	26	544	26	492	317	71	79	196	160	100	107	495	956	341	191	603	161	726	880	269	78	91	22	172	164	683	50	129	34	22	286	83	103	681	20	139	469	37	746	105	674	106	475	59	1322
AYGZ01000001.327.1863	48	44	34	0	0	2	39	3	0	0	94	1	653	0	16	1	209	29	8	59	228	1293	37	159	177	59	0	1182	0	1	10	4	30	5	161	4	31	9	220	50	90	111	17	0	40	400	0	133	0	70	18	1741	81	759
DQ777900.1.1512	61	87	134	8	298	262	148	128	81	10	67	0	9	218	52	197	308	36	113	58	44	90	104	76	88	172	255	75	74	203	68	4	44	84	71	73	10	150	43	125	29	638	310	79	31	70	18	70	15	683	74	12	422	313
JX198592.1.1491	0	156	237	0	9	183	91	1496	0	27	0	56	0	40	92	257	1	8	10	0	39	210	22	440	7	9	15	1	0	125	116	0	0	0	26	0	49	22	10	0	11	0	0	0	7	0	51	18	0	0	0	0	0	27
EF097505.1.1302	78	1	13	854	7	54	16	280	237	16	0	18	277	55	43	14	0	109	48	594	1	114	0	48	55	1	6	2	33	1	8	124	1	10	1	4	10	2	2	9	0	61	1	4	19	13	41	0	314	14	19	2	2	119
AWFQ01000007.543.2090	62	1504	977	961	162	766	311	791	1905	399	0	683	3	649	640	758	2	6	672	39	164	332	303	299	2	539	73	0	42	1162	893	24	7	145	473	229	206	2	9	455	174	8	0	60	178	11	478	63	407	0	520	101	0	18
CP007457.1550473.1552017	2	32	647	112	167	392	225	305	49	290	0	315	1	486	200	910	0	0	444	0	24	132	6	295	0	533	38	1	0	121	653	1	0	64	405	86	57	2	6	197	319	5	0	0	0	8	276	29	3	1	278	14	0	1
AB627487.1.1525	30	114	210	59	0	41	13	0	0	0	160	10	11	10	1	6	156	127	18	39	3	2	3	37	6	110	7	1019	20	0	12	9	4	5	13	23	0	10	58	409	0	83	67	16	27	888	3	122	1	540	11	64	223	217
EU622689.1.1499	0	0	1	0	0	237	11	33	135	1	13	2	0	0	177	0	29	0	67	0	2	5	1	0	0	4	4	0	0	0	18	0	557	275	3	15	14	1	13	76	2	2	0	3	0	5	1	0	30	4	25	1	12	0
AB626908.1.1497	0	5	1	0	0	0	0	0	0	0	0	0	0	0	0	0	34	0	0	32	0	0	0	0	0	68	0	6	0	0	21	0	0	0	20	0	0	7	0	0	0	10	7	0	0	2	0	1	0	4	0	7	431	0
New.ReferenceOTU52	0	0	15	1	0	1	0	1	0	0	725	0	0	0	3	0	0	0	0	0	4	0	0	0	0	0	0	0	0	0	22	0	0	0	134	0	9	0	0	0	0	0	245	0	0	0	1	15	0	3	1	0	17	0
KC163095.1.1526	0	0	23	1	0	0	1	2	0	0	1181	0	0	0	6	0	1	0	0	0	8	0	0	0	0	0	0	0	0	0	39	0	0	0	272	0	14	0	0	0	0	0	334	0	0	0	5	16	0	2	0	0	24	0
JQ084996.1.1499	0	12	15	61	48	15	97	14	91	6	264	20	256	59	88	34	134	7	84	0	384	556	116	53	47	19	112	2	32	3	49	33	53	12	125	49	58	25	0	240	0	8	283	12	46	262	22	96	11	14	0	17	7	0
JX198656.1.1515	23	0	8	29	0	907	187	13	24	42	45	17	419	101	25	5	4	1	642	18	154	243	12	28	815	2	12	5	136	2	4	226	90	14	5	7	471	0	38	57	2	0	48	21	90	17	68	10	55	24	18	25	45	4
JQ085018.1.1515	94	24	19	2	11	78	36	219	14	66	35	99	144	63	47	168	14	76	87	42	931	231	289	128	57	43	18	329	159	11	15	23	122	156	10	20	99	3	40	56	8	5	34	130	16	8	42	40	46	11	128	310	162	23
AB627586.1.1525	0	1	2	0	0	0	0	0	0	8	0	10	0	0	0	0	0	1	0	0	0	0	0	0	197	0	0	0	0	0	0	0	0	0	8	7	0	0	0	2	0	80	2	0	1625	0	55	4	0	1	0	2	0	0
AB622818.1.1552	0	41	6	204	40	339	4	47	37	4	51	15	36	23	125	3	3	4	17	0	74	531	120	0	33	146	9	191	12	141	26	26	4	14	33	2	102	11	97	4	35	25	7	14	61	64	79	39	54	5	50	210	7	11
HM123880.1.1494	1	0	5	1296	0	115	225	1	0	0	1	0	0	0	0	0	45	0	255	0	0	0	0	28	68	0	0	0	7	0	0	0	0	10	15	58	0	3	1	51	0	0	6	1	0	0	0	0	0	1	12	0	12	0
HM124103.1.1494	0	0	3	439	0	43	77	0	1	0	2	0	0	1	0	0	19	0	101	0	0	0	0	16	22	0	0	0	0	0	0	0	0	0	0	12	0	0	0	14	0	0	2	1	0	0	0	0	0	0	6	0	0	0
CP015403.1369456.1371009	0	28	32	2	1	941	11	428	432	0	13	0	0	0	377	0	3	0	9	0	9	180	14	26	0	3	0	1	446	75	23	0	0	43	1	6	78	0	2	6	8	8	8	185	0	2	0	3	135	4	19	17	15	0
CDYJ01026408.72.1593	169	579	1446	1726	0	26	0	40	239	1	58	1	252	0	25	0	145	423	1	995	3	16	0	0	86	99	1	36	0	0	466	111	0	0	360	151	8	46	53	523	136	30	145	1	127	26	0	45	95	51	0	110	158	130
EU791028.1.1512	0	3	0	0	140	0	978	0	0	123	0	258	0	610	0	796	0	2	704	1	133	0	29	1305	0	0	261	0	0	143	0	0	47	357	0	0	0	0	0	0	0	0	0	0	0	0	160	0	0	0	200	1	0	0
EU791198.1.1518	0	0	0	0	40	0	312	0	0	41	0	79	0	182	0	300	0	0	200	0	50	0	9	509	0	0	81	0	0	53	0	0	5	95	0	0	0	0	0	0	0	0	1	0	0	0	55	0	0	0	47	0	0	0
JQ085082.1.1487	472	178	52	21	35	0	48	0	2	7	108	42	193	32	1	42	239	561	1	554	14	0	3	2	554	417	17	2	0	2	117	2449	0	0	75	217	0	944	1416	71	309	482	188	0	743	237	287	75	0	266	0	17	88	820
JQ083727.1.1489	1286	387	51	149	1570	473	627	376	859	34	1159	27	1765	1214	531	793	130	25	787	1214	837	39	85	302	347	141	536	518	650	519	95	434	742	2926	358	388	21	211	68	39	8	487	745	1312	193	30	36	69	82	212	14	35	1795	135
EF603126.1.1475	273	78	231	48	0	4	0	0	23	0	258	0	181	0	0	0	338	25	6	186	0	36	0	0	165	129	1	89	1	0	127	45	1	5	404	167	2	62	27	21	16	117	316	0	40	92	0	54	6	58	1	4	612	57
AB702766.1.1491	785	187	21	3	2	0	157	2	81	0	276	2	96	1	5	0	117	260	0	77	168	1	146	132	24	287	1	114	23	475	547	2150	0	0	389	132	17	415	784	428	502	361	69	274	355	421	0	212	19	315	0	409	236	552
EF603419.1.1478	0	45	615	1	62	58	1	44	55	70	232	132	0	212	95	164	192	1	67	0	217	43	81	1	0	510	208	230	51	51	292	1	13	114	353	429	127	86	131	120	92	25	358	205	0	79	114	148	67	152	91	97	900	0
DQ815760.1.1447	413	176	110	52	95	55	86	125	106	123	225	114	52	86	223	122	135	332	86	242	142	85	420	89	434	190	74	156	109	263	100	261	24	56	78	311	166	155	94	217	135	357	286	99	531	212	259	270	403	175	145	258	135	212
EU791171.1.1508	214	913	0	534	156	564	323	770	603	60	0	17	183	279	215	180	0	2	471	266	266	705	191	448	0	0	213	0	343	569	2	141	100	250	0	0	51	1	0	0	1	1291	1	199	0	299	41	0	59	1	89	1	1	1
JQ084893.1.1488	0	56	26	0	0	0	0	0	0	1	104	0	0	0	0	0	256	0	77	0	0	0	0	0	0	327	0	128	46	0	199	0	237	94	43	19	0	123	203	24	140	164	42	257	0	175	0	111	0	89	452	264	36	0
EU791086.1.1496	93	1870	1317	710	320	866	468	824	483	144	110	82	744	421	482	456	471	1735	952	617	341	699	75	435	1413	1738	533	1421	452	283	972	134	220	427	1042	104	59	400	422	101	441	264	19	22	304	20	30	28	1	6	16	8	27	12
EF603719.1.1492	202	418	301	897	877	1486	1050	1371	765	355	77	174	1637	1146	781	1357	780	156	1653	2500	892	1099	122	1268	101	439	1414	205	676	702	197	289	347	709	213	18	90	79	87	19	151	46	34	36	28	6	34	4	0	1	20	6	6	2
HM123899.1.1495	0	122	0	0	141	406	342	318	926	552	0	567	0	419	512	348	0	2	502	1	321	211	419	331	0	1	301	102	315	181	0	0	1070	131	1	0	469	0	0	1	0	0	0	389	0	0	488	0	682	0	472	228	0	0
EF406773.1.1501	109	58	98	11	4	0	1	0	15	16	97	39	20	6	0	4	132	150	0	58	0	1	0	0	164	141	13	101	0	0	104	50	1	2	84	239	13	120	110	163	113	170	217	0	458	142	25	198	123	312	11	42	143	113
EU791052.1.1511	320	259	266	753	228	275	242	148	310	201	111	167	243	646	231	456	507	49	451	231	151	252	190	214	376	436	595	286	156	307	133	241	999	201	220	1081	102	26	167	432	30	888	21	272	106	157	60	117	18	17	4	61	280	34
JN012682.1.1467	127	0	0	8	413	164	306	254	634	381	16	278	14	693	231	699	8	0	269	12	221	163	491	512	22	0	353	483	192	1017	0	36	317	83	0	48	265	0	0	31	0	1	17	182	254	0	384	0	369	1	452	98	0	4
AGZO01000001.452.1966	751	37	2	10	5	13	85	41	30	25	10	38	248	26	153	10	38	97	31	446	39	14	112	53	193	6	2	22	13	89	20	1982	13	8	3	58	339	0	14	68	75	47	1	3	2141	61	109	198	9	142	31	483	5	3517
JQ083904.1.1489	0	155	138	0	55	0	116	0	0	174	114	125	1	60	0	40	138	1	68	0	73	0	140	131	0	129	60	74	236	248	178	0	240	58	87	200	0	317	468	243	1163	266	300	638	0	706	231	331	3	294	149	134	232	0
AB606312.1.1491	1	812	733	2	181	0	715	0	1	872	687	813	0	273	1	206	913	3	419	3	308	0	862	540	1	547	242	435	1711	1215	1115	0	1636	364	481	1074	0	2066	3118	1467	835	2004	1469	3724	0	2872	1534	2231	3	1361	1016	673	1497	2
CP015402.2294917.2296466	145	439	332	328	142	393	296	293	844	301	376	246	170	293	461	358	175	310	409	304	232	319	898	365	206	109	415	148	799	770	331	354	24	348	322	667	605	46	8	586	109	463	335	259	130	98	425	444	1289	272	509	167	228	176
EU791056.1.1502	124	25	119	239	36	108	46	41	72	61	54	66	154	119	18	73	35	86	68	227	116	50	42	37	32	36	120	54	86	54	111	90	7	15	46	907	40	3	8	105	11	84	238	49	38	19	45	63	49	48	56	10	242	56
EF604562.1.1497	112	37	116	207	8	24	5	8	53	16	103	6	114	17	4	22	59	82	39	212	50	8	26	12	35	87	27	71	17	33	168	95	2	7	39	1451	12	3	8	138	27	95	430	10	30	25	12	80	32	96	28	20	526	37
AB606317.1.1500	686	134	150	129	24	86	88	119	126	48	219	67	63	117	67	70	193	884	71	202	75	130	168	109	810	63	44	115	184	91	295	611	20	52	244	882	274	78	163	782	278	475	248	132	1295	723	199	1414	270	2500	292	800	198	1927
EU451397.1.1414	0	0	106	0	0	0	0	8	0	0	0	0	1	0	20	0	0	620	0	0	2	9	1	2	209	204	0	0	0	0	445	1	1	0	199	0	17	965	2229	0	540	0	0	0	232	0	0	117	0	0	0	0	0	1065
AB606322.1.1499	0	0	0	0	90	92	281	659	506	54	172	7	0	367	146	317	11	0	204	0	455	780	30	127	0	20	210	0	509	58	0	1	608	306	335	496	0	122	218	32	0	1	524	337	0	0	68	0	170	75	149	0	790	0
New.ReferenceOTU559	0	0	0	0	53	46	126	315	629	30	225	3	0	149	81	113	4	0	69	1	440	387	40	66	0	20	82	0	228	71	0	0	222	119	127	6	1	70	89	1	0	0	597	175	0	0	42	0	160	0	59	1	21	0
EU791133.1.1509	131	241	148	305	7	298	106	147	98	5	88	2	406	41	26	10	128	265	124	600	142	165	5	113	243	193	26	87	166	42	207	145	35	32	68	759	67	14	0	255	19	193	207	58	50	28	67	51	91	331	19	56	428	46
EU791170.1.1521	400	897	555	1259	27	927	366	526	299	12	309	10	1333	141	98	23	330	759	430	1758	412	609	20	314	763	524	95	329	580	131	646	502	117	105	230	2839	220	23	4	797	37	756	639	245	170	107	209	142	335	924	49	171	1518	137
EF406604.1.1522	75	991	625	4	0	2	1	0	79	0	324	10	1	1	0	0	177	480	1	8	54	1	57	2	58	784	0	38	5	321	690	64	0	0	956	3715	0	281	63	1577	137	1075	46	16	25	1026	2	280	86	171	0	187	137	238
EF603122.1.1484	104	320	414	213	37	179	163	192	106	92	391	99	117	156	177	89	170	63	245	259	197	210	205	72	34	80	183	53	168	196	301	123	30	107	341	963	275	108	80	458	42	211	374	119	96	217	93	193	93	442	237	180	128	42
EU791037.1.1514	0	0	221	0	0	0	0	0	0	1	5	0	0	0	0	0	137	595	0	2	9	0	0	0	267	116	0	111	0	0	263	0	0	0	37	4	0	36	42	269	29	0	2	1	93	0	0	114	0	60	1	121	0	116
EF406671.1.1517	0	0	189	0	0	0	0	0	0	0	9	0	0	0	0	0	121	641	1	0	8	0	0	4	269	105	0	80	0	0	215	0	0	0	38	1	0	30	30	233	39	0	2	0	75	0	0	65	0	52	0	86	1	149
EU791186.1.1495	2	81	296	0	4	0	22	0	3	22	22	23	0	10	0	19	231	2	42	0	2	1	0	25	0	82	9	106	5	5	453	1	0	91	172	374	0	164	199	299	83	87	0	57	0	94	26	167	0	180	50	59	2	2
JQ084057.1.1501	223	239	937	20	18	108	68	114	47	67	58	64	36	32	0	50	696	677	121	82	310	78	268	128	336	288	39	346	33	783	1530	182	4	279	577	1170	117	599	660	1076	441	294	3	202	281	345	76	469	7	623	123	161	1	467
EF406453.1.1508	204	41	67	123	34	66	56	95	182	91	118	86	109	87	115	73	54	388	179	191	52	89	127	102	271	28	139	42	118	254	59	267	94	76	54	275	136	49	24	229	31	198	60	198	212	86	73	71	173	219	273	71	53	451
New.ReferenceOTU495	0	0	96	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	1	44	0	0	0	0	426	1	0	0	1193	0	0	650	183	0	177	0	0	0	0	0	0	179	0	0	0	0	0	0
EU791120.1.1513	56	206	121	173	12	271	91	138	84	5	101	6	220	35	20	3	83	126	101	308	131	189	4	89	142	111	25	58	130	37	142	92	37	27	58	489	45	6	1	155	9	132	199	52	21	20	52	32	88	182	10	29	255	36
EF406763.1.1510	1777	11	682	381	22	68	48	157	201	45	1604	54	1228	74	20	33	984	1829	278	1142	15	489	447	49	2456	645	74	1	464	100	157	982	6	448	728	3340	181	49	161	178	48	169	1710	295	373	121	36	240	57	1225	596	0	1064	1400
EU791051.1.1505	807	110	87	45	0	0	0	1	2	1	134	1	47	0	0	2	142	418	0	253	1	0	0	0	248	99	0	114	97	0	71	24	0	41	73	908	0	302	293	566	439	426	227	291	3115	330	0	543	1	933	37	446	364	1072
EU791168.1.1511	32	29	82	597	36	208	91	306	21	64	144	58	155	117	143	111	69	539	142	403	88	233	104	128	362	73	185	25	132	441	56	21	14	102	39	436	189	3	1	14	0	82	411	274	112	24	106	4	63	186	165	7	605	164
EF406783.1.1508	100	96	213	1631	94	582	223	851	34	149	439	127	510	268	388	289	170	1533	288	1419	206	743	244	299	1077	239	512	107	298	1142	153	68	28	252	109	1202	538	9	4	41	8	252	1565	720	330	60	197	9	123	468	417	11	1811	622
EU791091.1.1520	32	20	36	482	9	125	69	158	17	35	89	61	98	46	81	38	61	201	85	126	38	155	83	38	255	23	98	33	100	189	29	14	6	83	26	303	152	6	2	5	1	64	219	172	108	12	87	8	35	78	107	2	434	88
AB702764.1.1495	4	1286	3166	1	910	1	1840	4	5	251	1802	88	4	334	3	317	2564	3	1051	8	5	6	2	915	1	2853	3816	628	862	1	1861	1	35	3311	1511	361	0	529	3	4073	2311	2174	2008	988	2	427	1464	3365	1	2707	954	2911	2341	0
EU622719.1.1482	2	9	23	0	8	0	173	0	0	3	0	1	2	0	18	8	61	2	2	0	0	0	1	543	4	1550	0	75	107	0	411	63	7	1	167	0	0	1967	805	0	1154	876	0	2829	0	530	1	437	1	0	1	2509	0	0
New.ReferenceOTU235	10	0	0	1	74	1	17	4	4	55	0	2	6	7	108	31	0	0	21	0	0	1	0	8	0	0	4	0	0	0	0	517	70	0	1	0	0	0	0	0	0	0	0	0	0	0	37	0	14	0	0	0	0	0
New.ReferenceOTU172	114	1	0	1	434	3	167	36	23	280	0	58	36	23	827	188	0	0	266	0	0	0	0	33	0	0	39	0	0	0	0	3768	486	0	0	0	11	0	1	0	0	0	0	0	0	0	347	1	100	0	0	0	0	0
MOYD01000062.205.1704	138	0	0	0	417	2	217	52	31	336	0	63	55	54	1152	258	0	0	386	0	0	0	0	41	0	0	53	0	0	0	0	4555	689	0	1	1	20	0	0	0	0	0	0	0	0	0	449	0	172	0	0	0	0	0
AB606394.1.1473	169	41	0	0	18	0	398	0	185	163	218	645	0	31	0	76	110	67	29	1	44	0	112	641	0	32	0	30	0	471	600	0	37	69	8	8	0	590	967	163	348	645	22	0	794	148	561	93	364	732	369	0	32	73
EU451957.1.1393	24	1	17	1	16	0	27	15	3	3	5	1	2	27	18	55	8	27	6	144	5	6	23	31	129	43	27	14	8	20	96	24	0	14	29	19	2	6	47	21	20	88	7	18	48	8	11	116	1	487	12	13	8	185
KM244845.1.1500	135	0	1	144	213	24	2	1999	1385	4	329	16	57	395	2	341	32	1982	1	929	824	1	168	30	280	0	7	0	756	1497	1	1	0	0	1	0	14	0	0	0	0	0	291	4	4	1	34	0	71	0	2	0	4	726
ASTF01000015.433637.435168	3	142	7	2	1	1474	785	494	1	53	10	104	29	4	810	1	118	1	59	5	19	38	565	128	28	26	2	44	16	251	1	7	320	271	46	0	355	6	8	8	0	2	11	12	0	5	19	5	3	9	104	219	6	1
EF406503.1.1505	0	78	3	8	0	3	0	0	0	6	544	0	0	0	0	0	70	0	1	0	0	2	0	0	3	3	1	19	2	0	0	0	3	0	1	0	1	1	1	0	0	0	665	1	0	0	0	2	2	0	0	284	0	1
AB606242.1.1511	0	0	0	0	2	176	0	1	28	0	0	0	0	3	15	0	260	559	2	0	79	0	1	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	2	0	13	0	98	0	0	0	0	121
AB627560.1.1532	81	1	179	114	8	39	3	27	19	26	153	22	1	16	127	14	100	40	7	144	5	8	24	2	558	2	4	87	128	1	54	2	13	5	125	206	340	0	0	178	0	30	92	149	88	1	26	461	212	47	13	77	19	6
FJ879725.1.1573	2	12	10	4	89	10	28	2	0	49	15	31	0	66	20	96	6	3	0	5	0	0	40	14	1	2	17	389	6	37	3	1	122	74	3	2	0	2	105	3	0	6	2	3	0	13	34	16	3	2	33	433	2	5
AB969479.1.1526	0	1	1	0	0	6	30	12	25	0	77	1	0	2	10	0	0	77	34	0	39	0	386	21	0	0	0	0	0	436	18	0	3	0	84	0	0	0	222	247	93	0	0	0	27	0	41	1	18	340	3	0	11	82
HM124124.1.1484	0	1	0	0	308	0	1209	1	11	0	0	6	1	1630	0	15	1	0	0	0	3611	0	2837	394	1	372	0	1	2	937	0	0	0	1	0	1	0	1338	0	2	0	22	0	0	0	23	627	0	16	0	0	11	0	2
New.ReferenceOTU510	0	134	109	0	0	0	3	0	2	0	14	0	0	0	0	0	3	53	1	0	0	0	0	0	2	4	0	50	0	0	1	0	0	0	470	7	0	3	0	20	0	3	13	2	0	1	0	40	1	0	0	0	118	0
AB606270.1.1503	0	262	207	0	0	0	3	0	7	0	13	0	0	0	0	0	2	50	1	0	0	0	0	0	1	10	0	107	1	0	4	0	2	0	659	12	0	1	0	23	0	3	28	0	0	0	0	110	0	0	0	0	267	0
EF604541.1.1490	0	573	76	54	36	0	0	0	1	1	41	0	0	1	0	0	606	27	1	0	2	0	3	0	9	326	0	8	0	1	2	1	0	7	505	65	0	2	0	495	0	97	113	0	0	127	0	35	1	4	1	0	6	3
EU504917.1.1400	54	274	1	1	0	6	2	48	49	1	265	10	369	2	14	0	310	10	224	14	5	1	0	0	1	32	4	181	98	0	0	221	11	6	0	2	25	377	9	0	3	706	15	277	58	1239	1	4	156	0	37	79	0	40
FJ879561.1.1497	2	623	0	0	0	0	1	0	0	1	672	0	0	0	0	0	735	0	0	0	0	0	0	0	0	0	0	433	0	0	0	0	0	0	0	0	0	0	1	0	0	5	42	0	0	0	0	0	0	1	1	221	0	0
EF406575.1.1505	87	360	2	0	0	0	0	1	1	0	376	0	401	0	1	1	456	15	113	32	0	2	0	0	1	35	0	244	2	2	0	256	1	1	0	4	8	319	25	2	4	689	21	7	44	1239	0	2	7	0	2	156	0	49
AB702782.1.1499	0	4	5	0	5	0	29	14	13	0	5	0	0	2	17	0	13	99	0	0	0	0	0	0	525	0	0	366	1161	1	3	0	0	0	2	0	0	18	2	0	0	3	76	177	10	74	36	0	3	1	3	8	7	3
AB606327.1.1485	6	71	54	2	4805	1	4	25	193	0	506	7	0	35	4	22	543	2	3	1	185	7	662	5	2060	3	0	8	0	192	27	0	2	1478	374	44	4	0	0	116	0	950	1321	1	24	187	33	15	10	0	554	0	0	1
AB622844.1.1532	2	13	4	0	438	0	2	2	18	0	41	2	2	5	1	2	44	0	0	6	11	2	55	0	171	4	0	0	1	17	4	0	0	99	37	4	0	0	0	9	1	93	129	0	0	16	3	1	0	0	42	0	0	0
AB622841.1.1531	0	12	14	0	772	0	0	1	32	0	82	2	0	4	1	3	94	2	0	0	22	0	112	1	346	0	0	0	0	24	6	0	0	227	43	5	0	0	0	19	0	198	220	0	7	18	2	4	0	0	112	0	0	0

plot_heatmap(ps.ng.tax_abund, method = "NMDS", distance = "bray")

5 Taxonomic summary

5.1 Bar plots in phylum level

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))

# aes(color='Phylum', fill='Phylum') --> aes() ggplot(data=data, aes(x=Sample,
# y=Abundance, fill=Phylum))
plot_bar(ps.ng.tax_most_, 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 = 5,
    colour = "black")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 2))  #6 instead of theme.size

Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.

ps.ng.tax_most_pre_post_stroke <- merge_samples(ps.ng.tax_most_, "pre_post_stroke")
ps.ng.tax_most_pre_post_stroke_ = transform_sample_counts(ps.ng.tax_most_pre_post_stroke,
    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.ng.tax_most_pre_post_stroke_, 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"))

Use color according to phylum. Do separate panels Stroke and Sex_age.

ps.ng.tax_most_copied <- data.table::copy(ps.ng.tax_most_)
# FITTING7: regulate the bar height if it has replicates
otu_table(ps.ng.tax_most_)[, c("1")] <- otu_table(ps.ng.tax_most_)[, c("1")]/6
otu_table(ps.ng.tax_most_)[, c("2")] <- otu_table(ps.ng.tax_most_)[, c("2")]/6
otu_table(ps.ng.tax_most_)[, c("4")] <- otu_table(ps.ng.tax_most_)[, c("4")]/6
otu_table(ps.ng.tax_most_)[, c("5")] <- otu_table(ps.ng.tax_most_)[, c("5")]/6
otu_table(ps.ng.tax_most_)[, c("6")] <- otu_table(ps.ng.tax_most_)[, c("6")]/6
otu_table(ps.ng.tax_most_)[, c("7")] <- otu_table(ps.ng.tax_most_)[, c("7")]/6

otu_table(ps.ng.tax_most_)[, c("8")] <- otu_table(ps.ng.tax_most_)[, c("8")]/7
otu_table(ps.ng.tax_most_)[, c("9")] <- otu_table(ps.ng.tax_most_)[, c("9")]/7
otu_table(ps.ng.tax_most_)[, c("10")] <- otu_table(ps.ng.tax_most_)[, c("10")]/7
otu_table(ps.ng.tax_most_)[, c("11")] <- otu_table(ps.ng.tax_most_)[, c("11")]/7
otu_table(ps.ng.tax_most_)[, c("12")] <- otu_table(ps.ng.tax_most_)[, c("12")]/7
otu_table(ps.ng.tax_most_)[, c("13")] <- otu_table(ps.ng.tax_most_)[, c("13")]/7
otu_table(ps.ng.tax_most_)[, c("14")] <- otu_table(ps.ng.tax_most_)[, c("14")]/7

otu_table(ps.ng.tax_most_)[, c("15")] <- otu_table(ps.ng.tax_most_)[, c("15")]/6
otu_table(ps.ng.tax_most_)[, c("16")] <- otu_table(ps.ng.tax_most_)[, c("16")]/6
otu_table(ps.ng.tax_most_)[, c("17")] <- otu_table(ps.ng.tax_most_)[, c("17")]/6
otu_table(ps.ng.tax_most_)[, c("18")] <- otu_table(ps.ng.tax_most_)[, c("18")]/6
otu_table(ps.ng.tax_most_)[, c("19")] <- otu_table(ps.ng.tax_most_)[, c("19")]/6
otu_table(ps.ng.tax_most_)[, c("20")] <- otu_table(ps.ng.tax_most_)[, c("20")]/6

otu_table(ps.ng.tax_most_)[, c("21")] <- otu_table(ps.ng.tax_most_)[, c("21")]/8
otu_table(ps.ng.tax_most_)[, c("22")] <- otu_table(ps.ng.tax_most_)[, c("22")]/8
otu_table(ps.ng.tax_most_)[, c("23")] <- otu_table(ps.ng.tax_most_)[, c("23")]/8
otu_table(ps.ng.tax_most_)[, c("24")] <- otu_table(ps.ng.tax_most_)[, c("24")]/8
otu_table(ps.ng.tax_most_)[, c("25")] <- otu_table(ps.ng.tax_most_)[, c("25")]/8
otu_table(ps.ng.tax_most_)[, c("26")] <- otu_table(ps.ng.tax_most_)[, c("26")]/8
otu_table(ps.ng.tax_most_)[, c("27")] <- otu_table(ps.ng.tax_most_)[, c("27")]/8
otu_table(ps.ng.tax_most_)[, c("28")] <- otu_table(ps.ng.tax_most_)[, c("28")]/8

otu_table(ps.ng.tax_most_)[, c("29")] <- otu_table(ps.ng.tax_most_)[, c("29")]/4
otu_table(ps.ng.tax_most_)[, c("30")] <- otu_table(ps.ng.tax_most_)[, c("30")]/4
otu_table(ps.ng.tax_most_)[, c("31")] <- otu_table(ps.ng.tax_most_)[, c("31")]/4
otu_table(ps.ng.tax_most_)[, c("32")] <- otu_table(ps.ng.tax_most_)[, c("32")]/4

otu_table(ps.ng.tax_most_)[, c("33")] <- otu_table(ps.ng.tax_most_)[, c("33")]/7
otu_table(ps.ng.tax_most_)[, c("34")] <- otu_table(ps.ng.tax_most_)[, c("34")]/7
otu_table(ps.ng.tax_most_)[, c("35")] <- otu_table(ps.ng.tax_most_)[, c("35")]/7
otu_table(ps.ng.tax_most_)[, c("36")] <- otu_table(ps.ng.tax_most_)[, c("36")]/7
otu_table(ps.ng.tax_most_)[, c("37")] <- otu_table(ps.ng.tax_most_)[, c("37")]/7
otu_table(ps.ng.tax_most_)[, c("38")] <- otu_table(ps.ng.tax_most_)[, c("38")]/7
otu_table(ps.ng.tax_most_)[, c("39")] <- otu_table(ps.ng.tax_most_)[, c("39")]/7

otu_table(ps.ng.tax_most_)[, c("40")] <- otu_table(ps.ng.tax_most_)[, c("40")]/7
otu_table(ps.ng.tax_most_)[, c("41")] <- otu_table(ps.ng.tax_most_)[, c("41")]/7
otu_table(ps.ng.tax_most_)[, c("42")] <- otu_table(ps.ng.tax_most_)[, c("42")]/7
otu_table(ps.ng.tax_most_)[, c("43")] <- otu_table(ps.ng.tax_most_)[, c("43")]/7
otu_table(ps.ng.tax_most_)[, c("44")] <- otu_table(ps.ng.tax_most_)[, c("44")]/7
otu_table(ps.ng.tax_most_)[, c("45")] <- otu_table(ps.ng.tax_most_)[, c("45")]/7
otu_table(ps.ng.tax_most_)[, c("46")] <- otu_table(ps.ng.tax_most_)[, c("46")]/7

otu_table(ps.ng.tax_most_)[, c("47")] <- otu_table(ps.ng.tax_most_)[, c("47")]/9
otu_table(ps.ng.tax_most_)[, c("48")] <- otu_table(ps.ng.tax_most_)[, c("48")]/9
otu_table(ps.ng.tax_most_)[, c("49")] <- otu_table(ps.ng.tax_most_)[, c("49")]/9
otu_table(ps.ng.tax_most_)[, c("50")] <- otu_table(ps.ng.tax_most_)[, c("50")]/9
otu_table(ps.ng.tax_most_)[, c("51")] <- otu_table(ps.ng.tax_most_)[, c("51")]/9
otu_table(ps.ng.tax_most_)[, c("52")] <- otu_table(ps.ng.tax_most_)[, c("52")]/9
otu_table(ps.ng.tax_most_)[, c("53")] <- otu_table(ps.ng.tax_most_)[, c("53")]/9
otu_table(ps.ng.tax_most_)[, c("54")] <- otu_table(ps.ng.tax_most_)[, c("54")]/9
otu_table(ps.ng.tax_most_)[, c("55")] <- otu_table(ps.ng.tax_most_)[, c("55")]/9


# plot_bar(ps.ng.tax_most_swab_, x='Phylum', fill = 'Phylum', facet_grid =
# Patient~RoundDay) + geom_bar(aes(color=Phylum, fill=Phylum), stat='identity',
# position='stack') + theme(axis.text = element_text(size = theme.size,
# colour='black'))
plot_bar(ps.ng.tax_most_, x = "Phylum", fill = "Phylum", facet_grid = pre_post_stroke ~
    Sex_age) + 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 = 2))

5.2 Bar plots in class level

plot_bar(ps.ng.tax_most_copied, 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 = 5,
    colour = "black")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 3))

Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.

plot_bar(ps.ng.tax_most_pre_post_stroke_, 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"))

Use color according to class. Do separate panels Stroke and Sex_age.

#-- If existing replicates, to be processed as follows --
plot_bar(ps.ng.tax_most_, x = "Class", fill = "Class", facet_grid = pre_post_stroke ~
    Sex_age) + 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 order level

plot_bar(ps.ng.tax_most_copied, fill = "Order") + 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")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 4))

Regroup together pre vs post stroke and normalize number of reads in each group using median sequencing depth.

plot_bar(ps.ng.tax_most_pre_post_stroke_, fill = "Order") + 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")) + theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 4))

Use color according to order. Do separate panels Stroke and Sex_age.

# FITTING7: regulate the bar height if it has replicates
plot_bar(ps.ng.tax_most_, x = "Order", fill = "Order", facet_grid = pre_post_stroke ~
    Sex_age) + 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 = 4))

5.4 Bar plots in family level

plot_bar(ps.ng.tax_most_copied, fill = "Family") + geom_bar(aes(), stat = "identity",
    position = "stack") + scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF",
    "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080",
    "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000",
    "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF",
    "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9",
    "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017",
    "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117",
    "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117",
    "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = 5, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 8))

Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.

plot_bar(ps.ng.tax_most_pre_post_stroke_, fill = "Family") + geom_bar(aes(), stat = "identity",
    position = "stack") + scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF",
    "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080",
    "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000",
    "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF",
    "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9",
    "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017",
    "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117",
    "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117",
    "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = theme.size, colour = "black")) +
    theme(legend.position = "bottom") + guides(fill = guide_legend(nrow = 8))

Use color according to family. Do separate panels Stroke and Sex_age.

#-- If existing replicates, to be processed as follows --
plot_bar(ps.ng.tax_most_, x = "Family", fill = "Family", facet_grid = pre_post_stroke ~
    Sex_age) + geom_bar(aes(), stat = "identity", position = "stack") + scale_fill_manual(values = c("#FF0000",
    "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500",
    "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF",
    "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887",
    "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF",
    "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4",
    "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D",
    "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF",
    "#15317E", "#254117", "#FAAFBE", "#357EC7")) + 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 = 8))

6 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
#FITTING6: 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"))

Group	Chao-1	Shannon	OTU	Phylogenetic Diversity
Group1	6349.769	7.002665	2751	83.66721
Group1	7285.497	7.355100	2737	94.96087
Group1	5495.606	7.178491	2461	83.56788
Group1	6161.841	7.097847	2511	85.58772
Group1	5443.752	6.824459	2434	84.34122
Group1	4840.408	6.776871	2461	80.72348
Group2	7366.408	7.815119	3197	102.90263
Group2	6274.827	7.741948	2890	90.85852
Group2	5793.308	8.075933	2865	92.32601
Group2	5984.518	7.319433	2578	90.46149
Group2	5708.705	7.851736	2824	93.36881
Group2	6109.767	8.101662	2966	98.44462
Group2	6369.752	7.118622	2869	95.11685
Group3	6940.387	7.741135	3059	97.76165
Group3	6098.004	6.436854	2835	102.67316
Group3	6497.394	7.850194	2913	98.10649
Group3	5469.070	6.699081	2449	82.04387
Group3	6645.847	7.682096	2841	92.82798
Group3	6637.343	8.223364	3028	91.43127
Group4	6951.096	7.754592	3221	102.24768
Group4	6474.659	7.199445	2761	99.42394
Group4	6675.297	7.917291	3047	96.83298
Group4	6247.904	7.486238	2935	92.08430
Group4	6284.602	7.937227	3073	104.16131
Group4	5071.120	7.496080	2578	81.15671
Group4	7180.224	8.128899	3195	105.48484
Group4	6861.665	8.148061	3102	90.74225
Group5	6395.912	6.601610	2610	90.63905
Group5	6268.760	6.807457	2664	90.49587
Group5	5597.646	6.368809	2493	84.34328
Group5	6628.257	7.639552	2972	96.22455
Group6	6300.821	7.445990	2832	99.36449
Group6	6444.541	7.421449	2884	98.71807
Group6	6947.157	7.091146	2687	100.62399
Group6	5410.097	6.866969	2584	90.80504
Group6	5953.121	6.937662	2599	95.33328
Group6	6462.002	7.897205	3000	95.28576
Group6	6937.500	7.726382	2805	96.43301
Group7	6512.715	7.595379	2754	85.83723
Group7	5963.597	7.406097	2772	92.08358
Group7	6399.684	6.788279	2636	88.25698
Group7	4851.865	6.066666	2044	78.10996
Group7	5618.788	6.274673	2345	89.33655
Group7	5485.123	7.549293	2527	81.88175
Group7	6661.582	7.828921	2987	99.76522
Group8	5147.719	7.330667	2421	88.32128
Group8	7077.365	7.990243	3085	96.26960
Group8	7017.728	7.668941	3149	107.96902
Group8	6990.919	8.152600	3202	105.21691
Group8	6067.818	7.284291	2705	97.74334
Group8	7205.626	7.463663	3060	106.55682
Group8	7038.728	7.969380	3170	104.27469
Group8	5821.239	7.778258	2727	85.05796
Group8	4878.500	7.002098	2147	67.58582

#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	Group1	Group2	0.0022194	0.060	0.00222	**	T-test
Shannon	Group1	Group3	0.2335726	1.000	0.23357	ns	T-test
Shannon	Group1	Group4	0.0004373	0.012	0.00044	***	T-test
Shannon	Group1	Group5	0.5625942	1.000	0.56259	ns	T-test
Shannon	Group1	Group6	0.1132840	1.000	0.11328	ns	T-test
Shannon	Group1	Group7	0.9072871	1.000	0.90729	ns	T-test
Shannon	Group1	Group8	0.0023776	0.062	0.00238	**	T-test
Shannon	Group2	Group3	0.4111307	1.000	0.41113	ns	T-test
Shannon	Group2	Group4	0.8285950	1.000	0.82859	ns	T-test
Shannon	Group2	Group5	0.0425835	1.000	0.04258		T-test
Shannon	Group2	Group6	0.0895027	1.000	0.08950	ns	T-test
Shannon	Group2	Group7	0.0583335	1.000	0.05833	ns	T-test
Shannon	Group2	Group8	0.6382240	1.000	0.63822	ns	T-test
Shannon	Group3	Group4	0.3405858	1.000	0.34059	ns	T-test
Shannon	Group3	Group5	0.1831506	1.000	0.18315	ns	T-test
Shannon	Group3	Group6	0.7707906	1.000	0.77079	ns	T-test
Shannon	Group3	Group7	0.3689790	1.000	0.36898	ns	T-test
Shannon	Group3	Group8	0.5695443	1.000	0.56954	ns	T-test
Shannon	Group4	Group5	0.0379492	0.950	0.03795		T-test
Shannon	Group4	Group6	0.0490209	1.000	0.04902		T-test
Shannon	Group4	Group7	0.0436902	1.000	0.04369		T-test
Shannon	Group4	Group8	0.4628546	1.000	0.46285	ns	T-test
Shannon	Group5	Group6	0.1845088	1.000	0.18451	ns	T-test
Shannon	Group5	Group7	0.5837050	1.000	0.58370	ns	T-test
Shannon	Group5	Group8	0.0595454	1.000	0.05955	ns	T-test
Shannon	Group6	Group7	0.3969530	1.000	0.39695	ns	T-test
Shannon	Group6	Group8	0.1683518	1.000	0.16835	ns	T-test
Shannon	Group7	Group8	0.0917359	1.000	0.09174	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
#FITTING6: 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, ...))
    }
}

p2 <- p + 
stat_compare_means(
  method="t.test",
  #comparisons = L.pairs,    # L.pairs 
  comparisons = list(c("Group1", "Group2"), c("Group1", "Group3"), c("Group1", "Group4"), c("Group1", "Group6"), c("Group1", "Group8"), c("Group2", "Group5"),c("Group4", "Group5"),c("Group4", "Group6"),c("Group4", "Group7"),c("Group6", "Group7")), 
  label = "p.signif",
  symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
  #symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
)

[1] FALSE

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

p3 <- p + 
stat_compare_means(
  method="t.test",
  #comparisons = L.pairs,    # L.pairs 
  comparisons = list(c("Group1", "Group2"), c("Group1", "Group4"), c("Group3", "Group4"), c("Group5", "Group6"),c("Group7", "Group8")), 
  label = "p.signif",
  symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
  #symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
)

[1] FALSE

#stat_pvalue_manual
#print(p2)
#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)

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

Alpha diversity

Alpha diversity by pre_post_stroke: the alpha diversity of the post stroke samples is significantly different from that of the pre stroke samples.

div.df2 <- div.df[, c("pre_post_stroke", "chao1", "shannon", "observed_otus", "PD_whole_tree")]
colnames(div.df2) <- c("pre_post_stroke", "Chao-1", "Shannon", "OTU", "Phylogenetic Diversity")
#colnames(div.df2)
options(max.print=999999)
#27     H47 830.5000 5.008482 319               10.60177
#FITTING6: 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 ~ pre_post_stroke, 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	post	pre	0.0011532	0.0012	0.0012	**	T-test

div_df_melt <- reshape2::melt(div.df2)
p <- ggboxplot(div_df_melt, x = "pre_post_stroke", y = "value",
              facet.by = "variable", 
              scales = "free",
              width = 0.5,
              fill = "gray", legend= "right")
lev <- levels(factor(div_df_melt$pre_post_stroke))
L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i])
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, ...))
    }
}
p2 <- p + stat_compare_means(
  method="t.test",
  comparisons = L.pairs,    # L.pairs
  #CHANGE: 
  #comparisons = list(c("-", "21"), c("-", "28"),  c("00", "07"), c("00", "21"), c("00", "28")), 
  #comparisons = list(c("-", "021"), c("-", "028"),  c("000", "007"), c("000", "021"), c("000", "028")),
  label = "p.signif",
  symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
  hide.ns = TRUE
)

[1] FALSE

ggsave("./figures/alpha_diversity_pre_post_stroke.png", device="png", height = 10, width = 12)
ggsave("./figures/alpha_diversity_pre_post_stroke.svg", device="svg", height = 10, width = 12)
# send both png and svg files to Laura

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

Alpha diversity

Group1: f.aged and post
Group2: f.aged and pre
Group3: f.young and post
Group4: f.young and pre
Group5: m.aged and post
Group6: m.aged and pre
Group7: m.young and post
Group8: m.young and pre

Merge Group1 and Group5 as aged_post. Merge Group2 and Group6 as aged_pre. Merge Group3 and Group7 as young_post. Merge Group4 and Group8 as young_pre.

Then perform the statistical test between aged_post and aged_pre –> significant. the statistical test between young_post and young_pre –> not significant.

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")
div.df2[div.df2 == "Group1"] <- "aged.post"
div.df2[div.df2 == "Group5"] <- "aged.post"
div.df2[div.df2 == "Group2"] <- "aged.pre"
div.df2[div.df2 == "Group6"] <- "aged.pre"
div.df2[div.df2 == "Group3"] <- "young.post"
div.df2[div.df2 == "Group7"] <- "young.post"
div.df2[div.df2 == "Group4"] <- "young.pre"
div.df2[div.df2 == "Group8"] <- "young.pre"
#colnames(div.df2)
options(max.print=999999)
#27     H47 830.5000 5.008482 319               10.60177
#FITTING6: 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"))

Group	Chao-1	Shannon	OTU	Phylogenetic Diversity
aged.post	6349.769	7.002665	2751	83.66721
aged.post	7285.497	7.355100	2737	94.96087
aged.post	5495.606	7.178491	2461	83.56788
aged.post	6161.841	7.097847	2511	85.58772
aged.post	5443.752	6.824459	2434	84.34122
aged.post	4840.408	6.776871	2461	80.72348
aged.pre	7366.408	7.815119	3197	102.90263
aged.pre	6274.827	7.741948	2890	90.85852
aged.pre	5793.308	8.075933	2865	92.32601
aged.pre	5984.518	7.319433	2578	90.46149
aged.pre	5708.705	7.851736	2824	93.36881
aged.pre	6109.767	8.101662	2966	98.44462
aged.pre	6369.752	7.118622	2869	95.11685
young.post	6940.387	7.741135	3059	97.76165
young.post	6098.004	6.436854	2835	102.67316
young.post	6497.394	7.850194	2913	98.10649
young.post	5469.070	6.699081	2449	82.04387
young.post	6645.847	7.682096	2841	92.82798
young.post	6637.343	8.223364	3028	91.43127
young.pre	6951.096	7.754592	3221	102.24768
young.pre	6474.659	7.199445	2761	99.42394
young.pre	6675.297	7.917291	3047	96.83298
young.pre	6247.904	7.486238	2935	92.08430
young.pre	6284.602	7.937227	3073	104.16131
young.pre	5071.120	7.496080	2578	81.15671
young.pre	7180.224	8.128899	3195	105.48484
young.pre	6861.665	8.148061	3102	90.74225
aged.post	6395.912	6.601610	2610	90.63905
aged.post	6268.760	6.807457	2664	90.49587
aged.post	5597.646	6.368809	2493	84.34328
aged.post	6628.257	7.639552	2972	96.22455
aged.pre	6300.821	7.445990	2832	99.36449
aged.pre	6444.541	7.421449	2884	98.71807
aged.pre	6947.157	7.091146	2687	100.62399
aged.pre	5410.097	6.866969	2584	90.80504
aged.pre	5953.121	6.937662	2599	95.33328
aged.pre	6462.002	7.897205	3000	95.28576
aged.pre	6937.500	7.726382	2805	96.43301
young.post	6512.715	7.595379	2754	85.83723
young.post	5963.597	7.406097	2772	92.08358
young.post	6399.684	6.788279	2636	88.25698
young.post	4851.865	6.066666	2044	78.10996
young.post	5618.788	6.274673	2345	89.33655
young.post	5485.123	7.549293	2527	81.88175
young.post	6661.582	7.828921	2987	99.76522
young.pre	5147.719	7.330667	2421	88.32128
young.pre	7077.365	7.990243	3085	96.26960
young.pre	7017.728	7.668941	3149	107.96902
young.pre	6990.919	8.152600	3202	105.21691
young.pre	6067.818	7.284291	2705	97.74334
young.pre	7205.626	7.463663	3060	106.55682
young.pre	7038.728	7.969380	3170	104.27469
young.pre	5821.239	7.778258	2727	85.05796
young.pre	4878.500	7.002098	2147	67.58582

#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	aged.post	aged.pre	0.0021984	0.01100	0.0022	**	T-test
Shannon	aged.post	young.post	0.2366071	0.64000	0.2366	ns	T-test
Shannon	aged.post	young.pre	0.0000953	0.00057	9.5e-05	****	T-test
Shannon	aged.pre	young.post	0.2122282	0.64000	0.2122	ns	T-test
Shannon	aged.pre	young.pre	0.2685374	0.64000	0.2685	ns	T-test
Shannon	young.post	young.pre	0.0502394	0.20000	0.0502	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
#FITTING6: 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, ...))
    }
}
p2 <- p + 
stat_compare_means(
  method="t.test",
  #comparisons = L.pairs,    # L.pairs 
  comparisons = list(c("aged.pre", "aged.post"), c("young.pre", "young.post"), c("young.pre", "aged.post")), 
  label = "p.signif",
  symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
  #symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
)

[1] FALSE

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

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

Alpha diversity

7 Beta diversity

Do multivariate analysis based on unifrac distance and PCoA ordination.

# obtained beta diversity from distance matrices (DM) and PCoA plots (unifrac)
# https://forum.qiime2.org/t/tutorial-integrating-qiime2-and-r-for-data-visualization-and-analysis-using-qiime2r/4121?u=nicholas_bokulich
# https://www.rdocumentation.org/packages/qiimer/versions/0.9.4
# http://www.science.smith.edu/cmbs/wp-content/uploads/sites/36/2015/09/Tutorial-from-sample-to-analyzed-data-using-Qiime-for-analysis.pdf

# https://bioconductor.statistik.tu-dortmund.de/packages/3.8/bioc/manuals/phyloseq/man/phyloseq.pdf
# http://evomics.org/wp-content/uploads/2016/01/phyloseq-Lab-01-Answers.html#beta-diversity-distances

mp <- readRDS("./ps.ng.tax.rds")
mpra <- microbiome::transform(mp, "compositional")

# # Define taxa to keep.  keepTaxa = prevdt[(Prevalence >= 10 & TotalCounts >
# 3), TaxaID] # Define new object with relative abundance mpra =
# transform_sample_counts(mp, function(x) x / sum(x)) # Filter this new object
# mpraf = prune_taxa(keepTaxa, mpra) # Calculate distances #
# distance(esophagus, 'uunifrac') # Unweighted UniFrac # distance(esophagus,
# 'wunifrac') # weighted UniFrac DistUF = distance(mpra, method = 'wunifrac')


# IMPORT QIIME DISTANCE MATRIX TODO: !!!!!!!!!!!!!!!!!!!!!!!!!! HOW to
# construct a phyloseq::distance object?
# https://joey711.github.io/phyloseq/distance.html

# disttb = as.matrix(DistUF) disttb[1:10,1:10] DistUF2 <- as.dist(disttb)
# https://stackoverflow.com/questions/17875733/how-to-convert-a-symmetric-matrix-into-dist-object
# https://stackoverflow.com/questions/25845220/how-to-read-a-matrix-from-text-file-in-r
d <- read.table("./bdiv_even42434/weighted_unifrac_dm.txt", header = TRUE, sep = "\t",
    row.names = 1)
DistUF2 <- as.dist(d)

ordUF = ordinate(mpra, method = "PCoA", distance = DistUF2)
plot_scree(ordUF, "Scree Plot: Weighted UniFrac Multidimensional Scaling")

# plot_ordination(mpraf, ordUF, color = 'SampleType') + geom_point(mapping =
# aes(size = DaysSinceExperimentStart, shape = factor(Subject))) +
# ggtitle('PCoA: Unweigthed Unifrac')
unwt.unifrac <- plot_ordination(mpra, ordUF, color = "Group")
unwt.unifrac <- unwt.unifrac + ggtitle("Weighted UniFrac") + geom_point(size = 2)
# https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
# library(randomcoloR)
nb.cols <- 18
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

# install.packages('Polychrome')
library(Polychrome)
# https://stackoverflow.com/questions/9563711/r-color-palettes-for-many-data-classes
P36 = createPalette(36, c("#ff0000", "#00ff00", "#0000ff"))

# https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf
unwt.unifrac <- unwt.unifrac + theme_classic() + scale_color_brewer("Group", palette = "Paired")  #Set3
# print(unwt.unifrac)
ggsave("./figures/beta_diversity1.png", device = "png", height = 6, width = 12)

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

Beta diversity

# 31.4 and 15.6

d <- read.table("./bdiv_even42434/unweighted_unifrac_dm.txt", header = TRUE, sep = "\t",
    row.names = 1)
DistUF2 <- as.dist(d)

# method=c('DCA', 'CCA', 'RDA', 'CAP', 'DPCoA', 'NMDS', 'MDS', 'PCoA')
ordUF = ordinate(mpra, method = "PCoA", distance = DistUF2)
plot_scree(ordUF, "Scree Plot: Unweighted UniFrac Multidimensional Scaling")

# plot_ordination(mpraf, ordUF, color = 'SampleType') + geom_point(mapping =
# aes(size = DaysSinceExperimentStart, shape = factor(Subject))) +
# ggtitle('PCoA: Unweigthed Unifrac')
unwt.unifrac <- plot_ordination(mpra, ordUF, color = "Group")
unwt.unifrac <- unwt.unifrac + ggtitle("Unweighted UniFrac") + geom_point(size = 2)
unwt.unifrac <- unwt.unifrac + theme_classic() + scale_color_brewer("Group", palette = "Paired")  #    #mycolors
# print(unwt.unifrac)
ggsave("./figures/beta_diversity4.png", device = "png", height = 6, width = 12)
# Bug: the PCI component ratio is not identical to the results of QIIME1.

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

Beta diversity

8 Differential abundance analysis

Differential abundance analysis aims to find the differences in the abundance of each taxa between two groups of samples, assigning a significance value to each comparison.

8.1 Group1 vs Group2

library("DESeq2")
# ALTERNATIVE using ps.ng.tax_most_copied: ps.ng.tax (40594) vs.
# ps.ng.tax_most_copied (166)
ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("1", "2", "4", "5", "6", "7",
    "8", "9", "10", "11", "12", "13", "14")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "Group2")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)

[1] "Intercept"              "Group_Group1_vs_Group2"

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Domain	Phylum	Class	Order	Family	Genus	Species
EF603872.1.1478	12.26963	21.812585	2.7930949	7.809468	0.00e+00	0.0000000	D_0__Bacteria	D_1__Firmicutes	D_2__Clostridia	D_3__Clostridiales	D_4__Lachnospiraceae	D_5__Lachnospiraceae NK4A136 group	D_6__uncultured bacterium
New.ReferenceOTU132	19.71680	21.942035	2.8326453	7.746129	0.00e+00	0.0000000	D_0__Bacteria	D_1__Firmicutes	D_2__Clostridia	D_3__Clostridiales	D_4__Peptostreptococcaceae	D_5__Romboutsia	D_6__uncultured bacterium
MPKA01000044.182665.184200	77.50853	-3.630940	0.8749931	-4.149679	3.33e-05	0.0299107	D_0__Bacteria	D_1__Firmicutes	D_2__Erysipelotrichia	D_3__Erysipelotrichales	D_4__Erysipelotrichaceae	D_5__Dubosiella	D_6__Dubosiella newyorkensis
CCPS01000022.154.1916	65.95551	8.490080	1.7546207	4.838698	1.30e-06	0.0015655	D_0__Bacteria	D_1__Proteobacteria	D_2__Gammaproteobacteria	D_3__Enterobacteriales	D_4__Enterobacteriaceae	D_5__Escherichia-Shigella	NA
EU622719.1.1482	205.38847	9.137709	1.9088938	4.786913	1.70e-06	0.0017389	D_0__Bacteria	D_1__Cyanobacteria	D_2__Melainabacteria	D_3__Gastranaerophilales	D_4__uncultured bacterium	D_5__uncultured bacterium	D_6__uncultured bacterium

# rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)) subv %in% v
# returns a vector TRUE FALSE is.element(subv, v) returns a vector TRUE FALSE


library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))

# Error in checkForExperimentalReplicates(object, modelMatrix) : The design
# matrix has the same number of samples and coefficients to fit, so estimation
# of dispersion is not possible. Treating samples as replicates was deprecated
# in v1.20 and no longer supported since v1.22.

8.2 Group3 vs Group4

ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("15", "16", "17", "18", "19",
    "20", "21", "22", "23", "24", "25", "26", "27", "28")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "Group4")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)

[1] "Intercept"              "Group_Group3_vs_Group4"

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	baseMean	log2FoldChange	lfcSE	stat	padj	Domain	Phylum	Class	Order	Family	Genus	Species
HK240365.1.1492	1533.42683	9.540979	1.713161	5.569224	8.7e-06	D_0__Bacteria	D_1__Verrucomicrobia	D_2__Verrucomicrobiae	D_3__Verrucomicrobiales	D_4__Akkermansiaceae	D_5__Akkermansia	D_6__uncultured bacterium
EU510524.1.1382	28.46869	23.318939	3.004054	7.762489	0.0e+00	D_0__Bacteria	D_1__Firmicutes	D_2__Clostridia	D_3__Clostridiales	D_4__Ruminococcaceae	D_5__Ruminococcaceae UCG-014	D_6__uncultured bacterium
New.ReferenceOTU132	35.95017	23.643776	3.003872	7.871100	0.0e+00	D_0__Bacteria	D_1__Firmicutes	D_2__Clostridia	D_3__Clostridiales	D_4__Peptostreptococcaceae	D_5__Romboutsia	D_6__uncultured bacterium
New.ReferenceOTU495	23.53857	23.048897	2.518147	9.153117	0.0e+00	Unassigned	NA	NA	NA	NA	NA	NA

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))

8.3 Group5 vs Group6

ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("29", "30", "31", "32", "33",
    "34", "35", "36", "37", "38", "39")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "Group6")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)

[1] "Intercept"              "Group_Group5_vs_Group6"

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab),
    ], "matrix"))
sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)),
    ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Domain	Phylum	Class	Order	Family	Genus	Species
CCPS01000022.154.1916	72.4453	7.389872	1.6812744	4.395399	1.11e-05	0.0108856	D_0__Bacteria	D_1__Proteobacteria	D_2__Gammaproteobacteria	D_3__Enterobacteriales	D_4__Enterobacteriaceae	D_5__Escherichia-Shigella	NA
JQ083727.1.1489	530.9005	-4.174814	0.9088159	-4.593686	4.40e-06	0.0051448	D_0__Bacteria	D_1__Bacteroidetes	D_2__Bacteroidia	D_3__Bacteroidales	D_4__Rikenellaceae	D_5__Alistipes	D_6__uncultured bacterium

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels = names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels = names(x))
ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) + geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) + theme(axis.text.x = element_text(angle = -25,
    hjust = 0, vjust = 0.5))

8.4 Group7 vs Group8

ps.ng.tax_sel <- ps.ng.tax
otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[, c("40", "41", "42", "43", "44",
    "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55")]
diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "Group8")
diagdds = DESeq(diagdds, test = "Wald", fitType = "parametric")
resultsNames(diagdds)

[1] "Intercept"              "Group_Group7_vs_Group8"

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
# sigtab = cbind(as(sigtab, 'data.frame'),
# as(tax_table(ps.ng.tax_sel)[rownames(sigtab), ], 'matrix')) sigtab <-
# sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)), ]
kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

baseMean	log2FoldChange	lfcSE	stat	pvalue	padj

# library('ggplot2') theme_set(theme_bw()) scale_fill_discrete <-
# function(palname = 'Set1', ...) { scale_fill_brewer(palette = palname, ...) }
# x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x)) x =
# sort(x) sigtab$Order = factor(as.character(sigtab$Order), levels=names(x)) x
# = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x)) x =
# sort(x) sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))
# ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) +
# geom_point(aes(size=padj)) + scale_size_continuous(name='padj',range=c(8,4))+
# theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))