Computational protocol: A longitudinal study of the diabetic skin and wound microbiome

Similar protocols

Protocol publication

[…] The PCR amplicons from 264 samples (including positive and negative controls) were sequenced over two separate runs on an Illumina Miseq using 500 cycle V2 kits. Sequences were demultiplexed using phylosift () and read pairs merged using FLASH (). Sequences were quality filtered and processed into OTUs using USEARCH v 1.8.1 () (fastq_filter command with the fastq_maxee option set to ‘2’ to remove all sequences with two or more expected errors). Further quality filtering and operational taxonomic unit (OTU) clustering was carried out in QIIME () version 1.9.0. The split_libraries.py command was used with the –l and –L options set to 240 and 260 respectively, to remove sequences shorter than 240 and longer than 260 base pairs. Sequences were clustered into OTUs at 97% similarity using the pick_open_reference_otus.py script using default settings except that singleton OTUs were removed, and the usearch61 method was used for chimera filtering.Taxonomy was assigned to OTUs (assign_taxonomy.py) using the UCLUST method () against the Greengenes () database pre-clustered at 97% similarity, accessed from the QIIME website (ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz). Representative sequences from each OTU were aligned against the Greengenes alignment using Pynast () (align_seqs.py), OTUs which failed alignment were filtered from the final OTU table (filter_otus_from_otu_table.py). A phylogenetic tree was built from aligned representative OTU sequences (make_phylogeny.py script) using Fasttree2 (), with the –r option set to midpoint for tree rooting.Diabetic skin samples adjacent to wounds were found to be more similar to wound than contralateral skin samples (see ), and were removed so as not to confound comparisons between diabetic and non-diabetic skin. To ensure more even sample sizes between the diabetic and non-diabetic groups, only the right foot samples were included from the non-diabetic group for all downstream analyses. Alpha diversity was calculated using Phyloseq () for the observed number of OTUs, Chao 1 and Shannon diversity indices on data rarefied to 30,000 sequences per sample. Significance testing was carried out on alpha diversity estimates using the Wilcoxon rank sum test in R.Initial beta-diversity analysis was carried out in QIIME on a rarefied OTU table (30K sequences per sample) using the weighted unifrac metric, and the generate_boxplots.py script used to compare unifrac distances between groups of samples. Futher beta diversity analyses, were carried out in Phyloseq, using weighted unifrac distances calculated from an OTU table with raw counts subject to variance stabilising transformation implemented in DEseq2 () as described here (). Weighted unifrac distances matrices were also subject to principal coordinates analysis using the Phyloseq package, and significant differences in variance between groups (diabetic and control skin) were determined with PERMANOVA (adonis function) implemented in the Vegan package () in R, using a nested model formula (∼health/subject + subject) and the weighted unifrac distance matrix.The Wald test for differential abundance was used as implemented in the DESeq2 package in R. Multivariate correlation analysis was carried out against OTUs and wound duration and area using Pearson scores with Bonferroni correction, and p-values were determined via bootstrapping with 100 permutations (implemented in QIIME using the observation_metatdata_correlation.py command). OTU tables were filtered to remove OTUs present in less than 10% of samples for both differential abundance and correlation tests.A Random forest learning algorithm implemented in R () was used to determine if diabetic status could be predicted from the foot skin microbiome. Skin samples were randomly divided into two equal subsets (restricting samples from the same participant to the same subset) for training and testing of learning algorithms. The variance stabilizing transformed OTU table was filtered to include skin samples only, and to remove OTUs observed in less than 10% of samples, and used as the input matrix for the Random forest algorithm. The Random forest fitted on the training subset was created using bootstrapping of one third of the training samples with replacement. As a general practice the rest of the samples were used as a validation set in order to decrease the risk of over-fitting associated with classification algorithms. An optimisation to minimise the out of bag error (classification error on validation data) was used to obtain the optimal number of taxonomic units accessed at each iteration of decision tree creation. Two hundred decision trees consisting of 30 OTUs evaluated at each node of the tree were created. The Random forest model was then used to predict the health status of the subjects in the test subset.Analysis of the stability of skin microbial communities over time was carried out by comparing intrapersonal weighted unifrac distances between the diabetic and control skin samples, along with intrapersonal distances for all samples. Kruskal–Wallis tests were used to determine significant differences between groups.Pearson’s Product Moment Correlation was used to test for correlations between wound size or duration and OTU abundance in wound samples as implemented in QIIME (observation_metadata_correlation.py). P-values were calculated using bootstrapping with 100 permutations, and Bonferroni correction for multiple testing. Kruskal-Wallis tests for OTUs that were differentially abundant in healing vs non-healing wounds were implemented in QIIME (group_significance.py). Wounds were classified as healing or non-healing based on a reduction in wound area since the last sampling time (healing) or no change or greater wound size area since the last sampling (non-healing). OTU tables were filtered to remove OTUs present in less than 10% of samples prior to testing.Inter-visit weighted unifrac distances were compared to the overall degree of healing (1 − (final wound area/initial wound area)) using the lm function of the stats package in R.Quality filtered sequence data has been deposited in the European Nucleotide Archive under study accession number PRJEB17696. A script containing the code used to process the data in R is provided as supplementary data, along with all the necessary input files, including OTU table and phylogenetic tree. […]

Pipeline specifications

Software tools PhyloSift, USEARCH, QIIME, UCLUST, PyNAST, FastTree, phyloseq, DESeq2
Databases Greengenes
Applications Phylogenetics, 16S rRNA-seq analysis
Organisms Homo sapiens
Diseases Diabetes Mellitus, Diabetes Mellitus, Type 2, Skin Diseases