Data Prep Short-Reads (PE)

The full paired-end Illumina workflow, from quality control to taxonomic assignment.

The following describes our typical processing workflow for Illumina paired-end amplicon data. We assume a two-step PCR library preparation and paired-end sequencing with overlapping reads (amplicon length shorter than the total number of cycles).

The data preparation workflow is a constant construction site. We try to integrate promising new methods without eliminating established older steps, which results in similar-looking files (e.g. OTU and zOTU tables) and, occasionally, undocumented supporting files. Do not hesitate to ask if you have trouble understanding any of the files you receive.

Pipeline Philosophy: Modular, Not Fixed

Our pipeline is modular rather than fixed. Which modules run, and with what parameters, is assembled to fit the dataset in front of us, since read length and quality, insert size, primer set, sample number, and marker gene all shape the right choices. This also keeps HPC resource usage sensible: if a step fails or turns out to be poorly configured, we rerun that step, not the entire pipeline from scratch. One consequence is that a single comprehensive page describing exact, fixed parameters for every case is difficult, if not impossible, to write, and isn't really what you'd want anyway, since your data may call for different choices than someone else's. What you do get, regardless of how your specific run was assembled, is a full report of what was actually done at each step and why.

Data Loss Is Not the Problem, Understanding It Is

Every step below loses some data, and this is unavoidable. The percentages quoted for each step are typical, not targets to chase down to zero. What actually matters is not how much data you lose, but whether you understand why. If we lose reads at the filtering step because the GC or size range was set too tightly, that's worth knowing, because adjusting the parameter can make a real difference and recover data that should have been kept. If we lose reads because they turned out to be PhiX, there is little to adjust, that data was never going to be usable, however large the number looks. A large loss with a clear, benign explanation is nothing to worry about. A small loss with no clear explanation is the one that deserves your attention. Minimising data loss is not the goal; understanding it is.

The Steps, at a Glance

Raw FASTQ (R1 / R2)
      │
      ▼
  A. Quality control & parameter evaluation
      │
      ▼
  B. PhiX / low-complexity filtering & read merging
      │
      ▼
  C. Primer site trimming
      │
      ▼
  D. Quality, size & GC filtering
      │
      ▼
  E. Clustering (OTU / zOTU)
      │
      ▼
  F. Taxonomic assignment
      │
      ▼
  G. Extras (e.g. ITSx)
      │
      ▼
Annotated count table

Not every project runs every module, and parameters at each step are set to fit the data rather than fixed in advance, see the note on our modular approach above. The steps are explained in order in the next section, roughly in the order they happen: each step assumes the data has already passed through the ones before it, for example primer trimming (Step C) happens after merging (Step B) because the full primer sequence is much easier to locate reliably once forward and reverse reads have already been combined into one.

Typical Retention, by Step

Step	Typical retention	Watch for
B: Filtering & merging	~95%	Unusually poor merge rates, unexpected amplicon length
C: Primer trimming	>95%	Excessive primer mismatches, staggered-read issues
D: Filtering	>97%	Many reads failing on GC or size range specifically
E: Clustering	varies by community	Excessive singleton/rare sequences, unstable cluster counts across reruns
F: Taxonomic assignment	not a read-loss step	Many zOTUs with low-confidence or unresolved assignments

Retention naturally varies by project, marker gene, and sample type, so treat these as a rough sense of what "normal" looks like, not a pass/fail threshold. If a step falls well outside these ranges, that's the signal to ask why, per the note above, rather than to worry about the number itself.

Data Processing Steps - Explained

Step A1 - Quality Control

FastQC and MultiQC provide a first, convenient impression of run quality, and many sequencing centres already generate these reports for you. On top of that, we use USEARCH and our own scripts to better understand the data structure and set processing parameters accordingly.

QC Reports

Quality control applications are easy to use and intuitive, but their output is often context-sensitive, and the abbreviations they use aren't always self-explanatory. The best report is of no use if you can't understand it. Read the manual. See also our note on FastQC and amplicon data on the Data Download page.

• Andrews (2010). FastQC: a quality control tool for high throughput sequence data. bioinformatics.babraham.ac.uk/projects/fastqc
• Ewels et al. (2016) MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics

Step A2 - Data and Parameter Evaluation

For more detailed data analysis, we estimate the expected error profile of the run and use it to set appropriate processing parameters for the steps that follow. We currently do this with USEARCH's eestats and our own scripts. Only the relevant samples are included in this step; controls are ignored as far as possible. Numerous checks (e.g. unique sample identification, minimum number of reads) are an important part of this step.

Step B - Filtering & Read Merging

It is often useful to start by cleaning up the raw data and removing problematic sequences before merging, rather than after, so that contaminating or low-complexity reads don't distort the merge statistics or get carried forward into every later step. First, we remove all PhiX-related reads. PhiX is added to the library to increase base complexity (see the Data Download page for more on why). The MiSeq Control Software (MCS) should already have filtered these out, but it is better to be safe than sorry. We also apply a low-complexity filter, which cleans up the data further, quickly and easily.

Once the worst of the data has been removed, we merge the read pairs. It is often worth trimming the 3' end of the reads first, to reduce mismatches and improve merging efficiency. How much to trim depends on the quality profile of the reads, the expected amplicon length, and the amount of overlap. Data loss of around 5% at this step is normal. Be careful with negative controls, which often have a much higher loss rate here and can skew the mean if included in summary statistics.

• Illumina PhiX Support Bulletin
• Seqtk
• Edgar (2010) Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26(19):2460-2461.
• Magoc and Salzberg (2011) FLASH: Fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27(21):2957-2963.

Step C - Primer Site Trimming

Primer regions are not part of the biological amplicon and must be removed. We do this after merging rather than before, since the full primer sequence is much easier to identify reliably once the forward and reverse reads have already been combined into a single, longer sequence. This step also lets us ensure that only amplicons with a recognisable primer site are carried forward. We require the complete primer site to be found with as few mismatches as possible, normally 1-2, and never at the 3' end. The exact parameters depend on your primer sequences, including their length and number of wobble bases. Expect data loss below 5% here; anything above that warrants a closer look. Pay particular attention when the expected amplicon is shorter than the read length (staggered reads).

Sample Carryover on the MiSeq

The MiSeq uses a fixed template-loading fluidic pathway to load sequencing template onto the flow cell before a run. A small number of template molecules can remain in this pathway after loading and get washed onto the flow cell in a later run, appearing there as clusters. On instruments maintained according to Illumina's recommendations, carryover is typically below 0.1% (one read in a thousand). For most applications these rare reads are simply discarded during trimming and filtering and do not pose a problem.

• Martin (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17(1):10-12.
• Edgar (2010) Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26(19):2460-2461.

Step D - Filtering

Before clustering the amplicons, we filter the merged reads further: an average quality filter (e.g. Q20), a size range, and a GC-content range selection. Data loss at this step is typically below 3%.

• Schmieder and Edwards (2011) Quality control and preprocessing of metagenomic datasets. Bioinformatics 27:863-864. (prinseq.sourceforge.net)
• Edgar (2010) Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26:2460-2461.
• BBTools, a suite of fast, multithreaded bioinformatics tools for DNA and RNA sequence data.

Coding Loci Look Different Here Too (e.g. COI)

Most amplicon markers, 16S, ITS, 18S, are non-coding, and their length can vary fairly freely between taxa without any functional consequence, which is part of why their length distributions tend to be comparatively broad. Protein-coding loci like COI behave differently: because an insertion or deletion that isn't a multiple of 3 bases shifts the reading frame and typically breaks the protein, indels away from a multiple of 3 are strongly selected against in real biological sequences. In practice this means coding loci usually show a much sharper, more tightly constrained peak in their length distribution than non-coding markers do (see the Data Download page's read-length check), and a length outside that tight peak is itself a mild warning sign at this stage. This constraint also enables an additional filtering step once clustering is done, see Step G below, rather than at this raw-read stage.

Step E - Clustering

At this stage, each sequence should represent a high-confidence biological amplicon, cleaned of adapters, primers, and low-quality bases. The remaining challenge is figuring out which sequences represent the same underlying biological variant and which differ only because of leftover sequencing error, which is what clustering does.

We currently use two clustering approaches, UPARSE and UNOISE, both implemented in USEARCH (see the Overview Workflows page for the conceptual difference between them). In UPARSE, dereplicated amplicons are clustered at 97% identity; singletons (below an abundance threshold of 2) and chimeras are removed. In UNOISE, dereplicated amplicons are first error-corrected, and each remaining sequence becomes a potential zero-radius OTU (zOTU) if it clears its own minimum abundance threshold. Both methods, in other words, use an abundance threshold to discard low-confidence sequences, UPARSE before clustering, UNOISE after denoising. This approach often yields a high number of rare zOTUs and can therefore overestimate diversity if taken at face value. To get a better sense of the underlying data structure, we additionally cluster the zOTUs at 99%, 98%, and 97% identity.

To obtain count/OTU tables, amplicons are mapped back to their clusters. Sequence identity and coverage thresholds during this mapping step balance efficiency against accuracy.

For Experienced Users: Mapping Algorithms

We map amplicons to (z)OTUs using a heuristic approach that finds a high-scoring alignment quickly. This can be replaced with a full dynamic-programming algorithm (Needleman-Wunsch, global alignment), which is guaranteed to find the highest-scoring alignment but is much slower, days rather than hours, than the heuristics USEARCH uses by default.

• Edgar (2013) UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nature Methods 10:996-998.
• Edgar (2016) UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing. bioRxiv 081257; doi.org/10.1101/081257

Step F - Predict Taxonomic Associations

For taxonomic assignment, we use SINTAX together with an appropriate reference database (see Choosing a Reference Database for the full walkthrough). SINTAX reports a confidence value at each taxonomic level; we typically use cutoffs between 70 and 90% to control sensitivity. Choosing the right reference isn't always easy: it doesn't need to be the largest one available, but it does need to fit your data. We commonly use SILVA for bacterial 16S data and UNITE for fungal ITS data.

For Experienced Users: Custom Reference Databases

Beyond these standard, general-purpose databases, we sometimes build custom reference databases for individual projects, typically where no suitable standard reference exists or where a project needs finer resolution than a general-purpose database can offer. These are usually project-specific rather than reusable across projects. We've developed a couple of approaches for building them in a SINTAX-compatible format, including curating candidate sequences from BLAST searches and applying a lowest common ancestor (LCA) approach. If you think your project might need a custom database, talk to us.

• Edgar (2016) SINTAX, a simple non-Bayesian taxonomy classifier for 16S and ITS sequences. bioRxiv 074161; doi.org/10.1101/074161
• Blackman, Walser, Rüber, Brantschen, Villalba, Brodersen, Seehausen, and Altermatt (2023) General principles for assignments of communities from eDNA: Open versus closed taxonomic databases. Environmental DNA 5(2):326-342.

Annotation Error Rates

Taxonomic annotation, regardless of method or database, carries a non-trivial error rate, for a number of reasons covered in detail on the Choosing a Reference Database page. Automatic annotation with SINTAX is simply great, and it's the only realistic way to add taxonomic labels to thousands of zOTUs, but treat the results as a strong starting hypothesis rather than a final answer. A species name assigned to a zOTU is a taxonomic hypothesis supported by the available reference data, not direct proof that the organism was present in your sample. zOTUs that turn out to matter for your conclusions should always be checked by hand.

Step G - Extra(s)

Some datasets need an additional step. For ITS amplicon data, for example, we use ITSx to identify and remove non-fungal clusters.

• Bengtsson-Palme et al. (2013) ITSx: Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for use in environmental sequencing. Methods in Ecology and Evolution 4:914-919.

ORF Checking for Protein-Coding Loci (e.g. COI)

For protein-coding markers like COI, we can apply an analogous filter to ITSx, but working on the zOTU representative sequences rather than raw reads: translate each zOTU in the appropriate reading frame, using the correct genetic code for the target group (for COI in most invertebrates, this is the NCBI invertebrate mitochondrial code, translation table 5, not the standard code), and check for an intact open reading frame, no premature stop codons, no frameshift.

Doing this after clustering rather than before has real practical advantages. There are far fewer zOTUs than raw reads, so the check is cheap to run, and since clustering and mapping are already complete by this point, dropping a failing zOTU is just removing a row from the finished count table. It doesn't change how the remaining zOTUs were clustered or how reads were mapped to them, so mapping efficiency for everything else is unaffected.

zOTUs that fail this check are often NUMTs (nuclear copies of mitochondrial genes), chimeras, or sequences carrying an uncorrected indel error, none of which represent a genuine, functional mitochondrial sequence. This isn't foolproof: some NUMTs are old enough, or young enough, to still translate cleanly, so passing the ORF check is reassuring but not definitive proof a zOTU is genuine. Talk to us if you'd like this applied to your COI project.

---

GDC Report Files, by Step

This section is specific to how we organise processing at the GDC and won't be relevant if you're running your own pipeline elsewhere. If you're working with data we've processed for you, this is where to find the report for each step above.

We divide processing from raw data to annotated count tables into discrete steps, both to make better use of computer resources and to speed up reprocessing when only part of the pipeline needs to be rerun.

Step	What happens	Report
A	Quality control, data and parameter evaluation	y_help/A_*.report
B	Filtering and read merging	y_help/B_*_CTM.report
C	Primer site trimming	y_help/C_*_PCR.report
D	Filtering	y_help/D_*_SQF.report
E	Clustering (OTU/zOTU)	y_help/E*.report
F	Taxonomic prediction	y_help/F*.report
G	Extras (e.g. ITSx for ITS amplicons)	y_help/G_*.report

The * in each filename is a placeholder for your project and sample identifiers.

Reproducibility

Every processing run is parameterised and logged. The reports listed above document the exact software versions, parameters, and modules used for your project, so the analysis can be reproduced or repeated later if needed, by us or by you.

---

Individual projects differ in read quality, marker gene, and the exact parameters used, and that's by design, not an inconsistency. The goal stays the same across all of them: turning raw sequencing reads into a reliable, biologically interpretable count table, while keeping genuine biological variation intact and technical artefacts out.