Data Preparation & Quality Control¶
Experimental / Sampling Design¶
The sampling design is simple. There are 48 samples in total. The study compares apples from two different management groups (treatment). For each treatment, the authors took 6 tissue samples from four apples.
N(samples): 2 treatments x 4 replicates x 6 tissues = 48
The experimental design is well described in the paper.
Four apples, weighing 190 ± 5 g, were selected from each of the two management groups and each apple was divided into six tissues with the following weights: stem: 0.2 g, stem end: 2 g, peel: 9 g, fruit pulp: 12 g, seeds: 0.2 g, and calyx end: 3 g. Thus, each tissue was represented by four replicates, where each replicate consists of the respective tissue of one apple. source: Wassermann et al. 2019
The way I understand it, there are 4 data points for each treatment and tissue. This matters, because it does not correspond with the results of the diversity analysis shown in Figure 2 (see Alpha Diversity). The design also means the different tissue samples are not independent, since they derive from the same apple. This is relevant to the treatment comparison, where all samples are grouped together.
Heterogeneous Sample Preparation
There are three concerns with pooling samples from different tissues to study treatment effects. (i) Tissue sample weight differs by a factor of 60. This could be problematic, since sample collection (e.g., genomic DNA concentration) can influence the observed variability in microbial communities (Multinu et al. 2018). (ii) Data obtained using different DNA extraction protocols is problematic, since extraction method can influence observed bacterial diversity (e.g., Teng et al. 2018) and bias comparisons. (iii) A further reason for not combining exterior and interior tissue samples is the post-harvest treatment of non-organic apples (see below).
Management Groups¶
The study focused on Arlet (Swiss Gourmet) apples from two management groups: organic and conventional. The comparison between treatments is a central part of the study. The fact that the conventionally grown apples underwent further treatment after harvest, while the organic apples did not, is clearly described in the materials and methods but ignored in the interpretation.
In contrast to the organically produced apples, they underwent the following post-harvest treatments: directly after harvest, apples were short-term stored under controlled atmosphere (1–2°C, 1.5–2% CO2), washed and wrapped in polythene sheets for sale. Both apple management groups ("organic" and "conventional") were transported to the laboratory immediately and processed under sterile conditions. source: Wassermann et al. 2019
Treatment Bias
Based on Buchholz et al. (2018), it seems reasonable to assume that post-harvest handling of the non-organic apples could have influenced their microbiota composition, at least for exterior tissues such as peel, calyx, or stem. The authors do not address this possibility and conclude that the difference between treatments is caused by management practice. I do not think the available data justify attributing the observed differences solely to management practice.
Raw Data¶
The authors state that the raw sequence files supporting their study are available from the ENA.
The raw sequence files supporting the findings of this manuscript are available from the European Nucleotide Archive (ENA) at the study Accession Number: PRJEB32455. source: Wassermann et al. 2019
Because Illumina MiSeq paired-end 250nt was used for sequencing, I expected two files (R1 and R2) per sample. Instead, there was only one sequence file per sample. I contacted the authors, and received this response:
We receive our reads from the sequencing company in the form of one forward file and one reverse file containing all reads of the whole pool, which, in this case, contained other samples as well. Data are demultiplexed after joining reads and removing barcodes. Those sequences are then provided in ENA. Thus, I cannot send you the raw files from this pool. I'm sorry for that. source: personal communication with Birgit Wassermann
The raw MiSeq reads are therefore not available. The authors provide merged reads, with limited detail on the demultiplexing or merging process itself.
Raw sequence data preparation and data analysis was performed using QIIME 1.9.1 (Caporaso et al., 2010). After paired reads were joined and quality filtered (phred q20), chimeric sequences were identified using usearch7 (Edgar, 2010) and removed. source: Wassermann et al. 2019
The deposited data are processed rather than instrument-level sequencing reads, which makes the term "raw data" ambiguous in this context.
What are raw data?
Raw data from an Illumina MiSeq paired-end run would normally comprise a forward (R1) and reverse (R2) read, if standard Illumina adaptors were used and the system demultiplexed the samples. Customised barcoded primers may not allow automated demultiplexing, in which case uploading demultiplexed reads is justified, but this needs to be clearly stated and the process clearly described.
Why does this matter? Merging forward (R1) and reverse (R2) reads affects how much data is removed. Trimming the qualitatively weaker ends of reads commonly increases merging efficiency. Here is an example from a 16S amplicon dataset produced at the GDC:
# R1: 0nt; R2: 0nt => 50.68% (merging efficiency)
# R1: 5nt; R2: 10nt => 95.62%
# R1: 10nt; R2: 20nt => 98.55%
# R1: 15nt; R2: 30nt => 90.75%
Since the read-merging step is not clearly described, the adjusted quality scores of the merged data are not reproducible.
Data Preparation Principle
Throughout this re-analysis, I retained as much sequence data as possible while documenting every filtering step and its rationale.
Contamination Control and Experimental Controls¶
The study provides only a brief description of sample processing prior to DNA extraction. While the methods describe homogenization of the different apple tissues and physical disruption of seeds and stems, they provide little information on the precautions taken to prevent contamination during sample handling. For 16S amplicon sequencing, where thousands of bacterial taxa are detected and many downstream conclusions depend on low-abundance organisms, such details are essential. Information on instrument sterilization between samples, handling order, processing environment, and measures to avoid cross-contamination is largely absent.
This lack of detail is compounded by the absence of experimental controls. No negative controls are reported, making it impossible to assess contamination introduced during sample processing, DNA extraction, library preparation, or sequencing. This is particularly relevant because the samples appear to have been processed by multiple experimenters with different levels of laboratory experience, potentially increasing technical variability. Likewise, no positive controls or mock communities were included to verify that the sequencing and bioinformatics pipeline accurately recovered known bacterial communities.
The absence of these controls does not invalidate the dataset, but it limits the ability to distinguish genuine biological signals from technical artefacts. This is especially important in this study because several conclusions rely on differences in rare taxa, which are the most susceptible to contamination and sequencing artefacts.
Sequence Processing¶
I downloaded 48 fastq.gz files from the European Nucleotide Archive (ENA) under study number PRJEB32455. In total, there are 8,248,796 amplicon sequences (merged reads). I compared the amplicon length distribution with the expected distribution obtained from in-silico PCR results based on the SILVA 16S database (version 128) and NCBI bacterial genomes (version July 2019).
While the mean amplicon length is similar between expected (292nt) and observed (289nt), the distributions differ: there are noticeably more short fragments in the apple data than expected.
Primer Trimming¶
Because the downloaded sequences still contained primer sequences and the paper does not mention primer trimming, I processed the data under the assumption that primer removal had not been performed previously.
For culture-independent Illumina MiSeq v2 (250 bp paired end) amplicon sequencing, the primers 515f – 806r (Caporaso et al., 2010) were used to amplify the 16S rRNA gene using three technical replicates per sample. source: Wassermann et al. 2019
I assume the original 515F/806R primers targeting the V4 region of the 16S SSU rRNA (Caporaso et al., 2011) were used:
I used USEARCH to remove the primer region, allowing a maximum of 1 mismatch, not at the primer end.
usearch v11.0.667_i86linux64
Amplicon range: 100-2000
Number of mis-matches: 1 / not at the end
Coverage: full-length
Wildcards enabled: IUPAC codes
I allowed one mismatch to reduce data loss at this step:
# Mismatches: 0 => 7,713,726
# Mismatches: 1 => 335,973
# Mismatches: 2 => 16,135
# Mismatches: 3 => 7,082
Quality Filtering¶
I quality filtered the data using Prinseq:
PRINSEQ-lite 0.20.4
Size Range: 100-600
GC Range: 30-70
Min Q Mean: 20
Number of Ns: 0
Low Complexity: dust / 30
OTU Clustering / Amplicon Sequence Variants¶
I used two approaches to define operational taxonomic units (OTUs): UPARSE (Edgar 2016a) and UNOISE (Edgar 2016b), including additional clustering for zero-radius OTUs (ZOTUs) to account for possible early PCR errors.
After removing chimeric, mitochondrial and chloroplast sequences, the overall bacterial community of all apple samples, assessed by 16S rRNA gene amplicon sequencing, contained 6,711,159 sequences that were assigned to 92,365 operational taxonomic units (OTUs). source: Wassermann et al. 2019
Number of OTUs : 3,124 (308 Mitochondria / 52 Chloroplast)
Number of ZOTUs : 3,182 (210 Mitochondria / 70 Chloroplast)
Observation
The paper reports 92,365 OTUs. My pipeline, using USEARCH on the same deposited reads, produced 3,124 OTUs (or 3,182 ZOTUs) — roughly 30 times fewer.
Inference
Since both analyses used USEARCH, the discrepancy is unlikely to be explained by the choice of clustering software alone.
Unknown
The exact cause cannot be determined from the published methods. Candidate explanations include differences in clustering parameters, reference database handling, or upstream processing that isn't described in sufficient detail to reproduce.
Representative sequences were aligned, the open reference database SILVA (ver128_97_01.12.17) was used to pick operational taxonomic units (OTUs), and de novo clustering of OTUs was performed using usearch. source: Wassermann et al. 2019
Processing of the provided merged reads removed 2.5% of the data on average. Data loss and mean amplicon size were steady across samples. Sequencing depth ranged from 18k to 428k reads, with fruit pulp samples at the lower end of the range and calyx samples at the top.
Data-Prep Statistics¶
Per-sample read counts and quality stats (click to expand)
| Sample | Merged | Primer | Clean | MeanLength | D% |
|---|---|---|---|---|---|
| H-Stem-a | 205,945 | 201,311 | 201,311 | 251.7 | 97.7 |
| H-Stem-b | 231,409 | 223,661 | 223,661 | 251.7 | 96.7 |
| H-Stem-c | 217,307 | 211,969 | 211,969 | 250.1 | 97.5 |
| H-Stem-d | 218,388 | 214,743 | 214,743 | 251.5 | 98.3 |
| C-Stem-a | 181,279 | 174,667 | 174,667 | 252.1 | 96.4 |
| C-Stem-b | 144,427 | 141,113 | 141,113 | 252.4 | 97.7 |
| C-Stem-c | 372,408 | 363,099 | 363,099 | 249.9 | 97.5 |
| C-Stem-d | 385,316 | 375,813 | 375,813 | 252.5 | 97.5 |
| H-StemEnd-a | 221,403 | 206,834 | 206,834 | 253.0 | 93.4 |
| H-StemEnd-b | 320,614 | 310,926 | 310,926 | 251.9 | 97.0 |
| H-StemEnd-c | 86,899 | 84,075 | 84,075 | 252.3 | 96.8 |
| H-StemEnd-d | 197,981 | 194,176 | 194,176 | 251.9 | 98.1 |
| C-StemEnd-a | 197,662 | 193,093 | 193,093 | 251.6 | 97.7 |
| C-StemEnd-b | 311,501 | 303,738 | 303,738 | 252.3 | 97.5 |
| C-StemEnd-c | 251,718 | 245,561 | 245,561 | 251.0 | 97.6 |
| C-StemEnd-d | 180,369 | 176,591 | 176,591 | 251.4 | 97.9 |
| H-Seeds-a | 137,752 | 135,032 | 135,032 | 252.3 | 98.0 |
| H-Seeds-b | 91,782 | 89,809 | 89,808 | 252.3 | 97.8 |
| H-Seeds-c | 132,468 | 130,086 | 130,086 | 252.2 | 98.2 |
| H-Seeds-d | 95,599 | 93,835 | 93,834 | 252.2 | 98.2 |
| C-Seeds-a | 54,411 | 52,680 | 52,680 | 252.6 | 96.8 |
| C-Seeds-b | 163,095 | 159,090 | 159,090 | 252.1 | 97.5 |
| C-Seeds-c | 81,263 | 79,258 | 79,258 | 251.8 | 97.5 |
| C-Seeds-d | 95,157 | 92,992 | 92,992 | 252.0 | 97.7 |
| H-Peel-a | 157,919 | 155,094 | 155,094 | 238.6 | 98.2 |
| H-Peel-b | 242,554 | 238,273 | 238,273 | 239.4 | 98.2 |
| H-Peel-c | 111,213 | 108,852 | 108,852 | 245.5 | 97.9 |
| H-Peel-d | 76,360 | 74,889 | 74,889 | 238.7 | 98.1 |
| C-Peel-a | 18,693 | 18,266 | 18,266 | 252.3 | 97.7 |
| C-Peel-b | 39,117 | 38,176 | 38,176 | 252.4 | 97.6 |
| C-Peel-c | 68,716 | 67,016 | 67,016 | 252.3 | 97.5 |
| C-Peel-d | 75,548 | 73,492 | 73,492 | 249.1 | 97.3 |
| H-FruitPulp-a | 39,123 | 38,461 | 38,461 | 250.8 | 98.3 |
| H-FruitPulp-b | 33,287 | 31,395 | 31,395 | 247.9 | 94.3 |
| H-FruitPulp-c | 51,723 | 50,754 | 50,754 | 250.6 | 98.1 |
| H-FruitPulp-d | 30,431 | 29,913 | 29,913 | 248.7 | 98.3 |
| C-FruitPulp-a | 50,188 | 48,984 | 48,984 | 252.1 | 97.6 |
| C-FruitPulp-b | 33,734 | 32,922 | 32,922 | 251.7 | 97.6 |
| C-FruitPulp-c | 84,835 | 82,503 | 82,503 | 252.1 | 97.3 |
| C-FruitPulp-d | 44,167 | 42,822 | 42,822 | 252.3 | 97.0 |
| H-CalyxEnd-a | 225,921 | 222,124 | 222,124 | 242.7 | 98.3 |
| H-CalyxEnd-b | 273,104 | 268,632 | 268,632 | 239.9 | 98.4 |
| H-CalyxEnd-c | 349,714 | 342,698 | 342,698 | 246.0 | 98.0 |
| H-CalyxEnd-d | 158,404 | 155,816 | 155,816 | 250.5 | 98.4 |
| C-CalyxEnd-a | 289,295 | 283,139 | 283,139 | 251.7 | 97.9 |
| C-CalyxEnd-b | 427,153 | 414,207 | 414,207 | 252.6 | 97.0 |
| C-CalyxEnd-c | 350,489 | 340,059 | 340,059 | 248.3 | 97.0 |
| C-CalyxEnd-d | 440,955 | 428,514 | 428,514 | 251.7 | 97.2 |
| Total | 8,248,796 | 8,041,153 | 8,041,151 | ||
| Mean | 171,850 | 167,524 | 167,524 | 250.1 | 97.5 |
Literature¶
Multinu et al. (2018). Systematic Bias Introduced by Genomic DNA Template Dilution in 16S rRNA Gene-Targeted Microbiota Profiling in Human Stool Homogenates. mSphere, 3(2).
Buchholz et al. (2018) The potential of plant microbiota in reducing postharvest food loss. Microbial Biotechnology, 11(6), 971–975.
Caporaso et al. (2011). Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci USA 108, 4516–4522.
Krzywinski & Altman (2014) Visualizing samples with box plots. Nature Methods. Vol.11 No.2.
McDonald (2014) Handbook of Biological Statistics (3rd ed.). Sparky House Publishing, Baltimore, Maryland.
Edgar (2016a) UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing, https://doi.org/10.1101/081257
Edgar (2016b), UCHIME2: improved chimera prediction for amplicon sequencing, https://doi.org/10.1101/074252
Teng et al. (2018) Impact of DNA extraction method and targeted 16S-rRNA hypervariable region on oral microbiota profiling. Scientific Reports 8, 16321.