Tara Pacific metabarcoding 18S V9 SWARM table

Authors

Nicolas Henry

Cédric Berney

Frédéric Mahé

Colomban de Vargas

Published

March 18, 2024

Contact

  • Data analysis
    • Frédéric Mahé
    • Nicolas Henry
  • Reference database
    • Cédric Berney

OTUs or ASVs?

In the first version of the Tara Pacific metabarcoding 18S V9 SWARM table we used the term OTUs (Operational Taxonomic Units) to describe the genetic entities resulting from the bioinformatic analysis (Fred’s pipeline) and constituting each row of the SWARM table. The term OTU has a bad connotation because it mostly refers to deprecated approaches consisting of arbitrary clustering raw reads at an arbitrary global similarity threshold in order to eliminate errors. These approaches tend to cluster together sequences belonging to different species or genera when using genetic markers such as the 16S and 18S.

Fred’s pipeline relies on the use of SWARM, a local similarity clustering based approach, which is much more efficient at clustering sequencing errors only. For each cluster the sequence with the highest number of reads is considered as the true sequence or variant and its abundance is the sum of reads of all the sequences of the clusters (one true variant + errors). For more information about how SWARM works, see Mahé et al. (2014), Mahé et al. (2015) and Mahé et al. (2022).

We have decided to use the term ASV (Amplicon Sequence Variant) to avoid using an umpteenth term to describe the genetic entities resulting from Fred’s pipeline. This term is used to describe the output of another state of the art bioinformatic approach named DADA2 (Callahan et al. 2016).

Files

  1. tara_pacific_18s_v9_swarm_asvs.fasta.gz: ASV sequences
  2. tara_pacific_18s_v9_swarm_asvs.tsv.gz: ASV descriptions such as taxonomy
  3. tara_pacific_18s_v9_swarm_counts.tsv.gz: ASV read counts in Tara Pacific samples
  4. tara_pacific_18s_v9_swarm_counts_with_duplicates.tsv.gz: Same as tara_pacific_18s_v9_swarm_counts.tsv.gz but before summing the duplicated samples; see Sample dereplication section
  5. tara_pacific_18s_v9_swarm_run_corresp.tsv.gz: run origin and negative controls associated to each sample
  6. tara_pacific_18s_v9_swarm_taxo_pr2.tsv.gz: alternative taxonomic assignment using PR2

Content of files

  • Sequences (tara_pacific_18s_v9_swarm_asvs.fasta.gz): ASV sequences in FASTA format

  • ASVs (tara_pacific_18s_v9_swarm_asvs.tsv.gz): ASV descriptions, with 15 tab-separated fields:

    1. asv_id: ASV unique name / id.
    2. euk: Y if the ASV is considered to be a eukaryote, N if not; see Euk-Prok filtration section
    3. total: ASV total number of reads
    4. spread: number of samples the ASV appears in
    5. supergroup: supergroup from the EukRibo-specific annotation; see Taxogroup assignment section
    6. taxogroup1: taxogroup1 from the EukRibo-specific annotation; see Taxogroup assignment section
    7. taxogroup2: taxogroup2 from the EukRibo-specific annotation; see Taxogroup assignment section
    8. conf_supergroup: % confidence in the supergroup assignment; see Taxogroup assignment section
    9. conf_taxogroup1: % confidence in the taxogroup1 assignment; see Taxogroup assignment section
    10. conf_taxogroup2: % confidence in the taxogroup2 assignment ; see Taxogroup assignment section
    11. closest_hits_taxonomy
    12. closest_hits_similarity: percentage of similarity with closest reference sequence(s) (edit distance: (matching columns) / (alignment length)).
    13. closest_hits_references: GenBank accession number(s) of closest reference sequence(s).
    14. sequence: nucleic acid sequence.
    15. sequence_length: nucleic acid sequence base pairs length.

  • Counts: ASV counts (number of reads associated to each ASV in the samples), before (tara_pacific_18s_v9_swarm_counts_with_duplicates.tsv.gz) and after (tara_pacific_18s_v9_swarm_counts.tsv.gz) summing the duplicates, provided in long format, with 3 tab-separated fields:

    1. asv_id: ASV unique id, see tara_pacific_18s_v9_swarm_asvs.tsv.gz for ASV descriptions
    2. sample: sample unique id. Sample names starting with “CEB_” are negative controls. See tara_pacific_18s_v9_swarm_run_corresp.tsv.gz to know to which run they correspond to.
    3. nreads: ASV number of reads measured for the corresponding sample

  • Samples to run correspondence (tara_pacific_18s_v9_swarm_run_corresp.tsv.gz): details about run origin and negative controls associated to each sample, with 5 tab-separated fields:

    1. sample_code: sample id used in tara_pacific_18s_v9_swarm_counts.tsv.gz
    2. sample_code_duplicates: sample id used in tara_pacific_18s_v9_swarm_counts_with_duplicates.tsv.gz
    3. seq_run: run id
    4. neg_sample_pcr1: first PCR negative control
    5. neg_sample_pcr2: second PCR negative control

  • Alternative taxonomic assignment with PR2 (tara_pacific_18s_v9_swarm_taxo_pr2.tsv.gz) alternative ASV taxonomy using PR2 see PR2 taxonomic assignment section, with 6 fields:

    1. asv_id: ASV unique name / id.
    2. taxonomy: IDTAXA taxonomic assignment result
    3. confidence_scores: IDTAXA confidence score
    4. closest_hits_similarity: percentage of similarity with closest reference sequence(s) (edit distance: (matching columns) / (alignment length)).
    5. closest_hits_references: GenBank accession number(s) of closest reference sequence(s).
    6. closest_hits_taxonomy: taxonomic path to last common ancestor of closest reference sequence(s).

Environmental samples overview

The samples correspond to DNA extracted from coral tissues (CORAL; n = 5044), seawater (PLANKTON; n = 1195), sediment (SEDIMENT; n = 219), aerosols (AEROSOL; n = 35) and fish (FISH; n = 6).

Negative controls

There are many ways of removing contaminants from an ASV table. So we decided to let the final user do it. The correspondence between samples and their respective PCR controls are listed in tara_pacific_18s_v9_swarm_run_corresp.tsv.gz

Method

Core pipeline (Fred’s pipeline)

Paired-end reads were merged with vsearch’s –fastq_mergepairs command (allowing for staggered reads). All merged fastq files were then trimmed with cutadapt, keeping only reads containing both forward and reverse primers. After trimming, the expected error per read was estimated with vsearch’s command –fastq_filter and the option –eeout. Each sample was then dereplicated, i.e. strictly identical reads were merged, using vsearch’s command –derep_fulllength, and converted into FASTA format.

Clustering was performed at the sample level with swarm using default parameters (Mahé et al. 2015), and results were set aside for later use. To prepare for global clustering, individual fasta files (one per sample) were pooled and further dereplicated with vsearch, resulting in 34,268,640 unique sequences, representing around 3.6 billion reads. Files containing per-read expected error values were also dereplicated to retain only the lowest expected error for each unique sequence. Clustering was performed with Swarm (with the fastidious option). Cluster representative sequences were then searched for chimeras with vsearch’s command –uchime_denovo (Edgar et al. 2011). In parallel, representative sequences were assigned to taxa using the stampa pipeline (https://github.com/frederic-mahe/stampa/) and the 18S database EukRibo (Berney et al. 2022; Berney 2022).

Clustering results, expected error values, taxonomic assignments, and chimera detection results were used to build a “raw” occurrence table, with all the reads except those without primers, shorter than 32 nucleotides, or reads with uncalled bases (“N”). In order to build a “filtered” occurrence table, additional filters were applied. Only non-chimeric sequences with an expected error per nucleotide below 0.0002 from a cluster containing at least 3 reads in 2 samples were retained. Primer trimming is not perfect and some sequences can still contain primer fragments or be excessively trimmed. These sub- or super-sequences were merged with their closest, most abundant, perfectly trimmed sequence. Finally, occurrence patterns throughout our sample collection were used to refine our results. Some clusters contain sub-clusters with only a single-nucleotide difference but with different ecological patterns (here defined as uncorrelated abundance values in at least 5% of the samples). These sub-clusters were turned into distinct clusters. On the other hand, some clusters with similar sequences present correlated abundance values in at least 95% of the samples. These correlated clusters were merged using mumu, a re-implementation of lulu’s method (Frøslev et al. 2017).

Fred’s pipeline scripts used for this table are here

Sample dereplication

Some samples have been sequenced twice, we call them duplicates. If the duplicates come from two different PCRs, the readset with the lowest number of reads is removed. If the duplicates are from the same PCR, they are summed together.

Filtration based on occurrence and abundance

ASVs present in only one sample or with less than 3 reads after dereplication were removed. Therefore the final ASV table contains 429,168 ASVs, 3,510,033,320 reads for 6,702 samples (counting 203 negative controls).

EukRibo

EukRibo (Berney et al. 2022) is a manually curated subset of sequences from PR2, with low sequence redundancy but high taxonomic accuracy and exhaustive coverage of the known diversity of eukaryotes. It is aimed specifically at taxonomic annotation of high-throughput metabarcoding datasets and to help discover new diversity at various taxonomic levels. A EukRibo-specific annotation into strictly monophyletic taxonomic groups at 3 levels that we called supergroup (e.g. Alveolata), taxogroup1 (e.g. Dinoflagellata), and taxogroup2 (e.g. MALV-II) was done.

Extra taxonomic information

Euk-Prok filtration

The primers used to amplify the V9 region of the 18S rRNA gene of eukaryotes also amplify the V9 region of the 16S rRNA from prokaryotes and organelles. In order to ease the selection of truly eukaryotic ASVs, we added a column (euk) in the table to indicate which ASVs are assigned to a eukaryote reference sequence with more than 80% similarity. This threshold has been defined by comparing V9 SSU rRNA gene reference sequences between eukaryotes, bacteria, archaea, mitochondria and chloroplasts. The figure below summarizes these comparisons. Each histogram represents the distribution of the best hits of each reference of the row group compared to the column group. For example, the top right histogram shows the distribution of the best hits of each eukaryotic reference sequence compared to all chloroplast reference sequences.

Taxogroup assignment

In complement to the taxonomic assignment made using stampa, ASVs were assigned to EukRibo taxonomic groups (i.e. supergroup, taxogroup1, taxogroup2) using IDTAXA (threshold = 50) (Murali, Bhargava, and Wright 2018). For each group, the associated percentage of confidence is reported in the columns conf_*group*.

PR2 taxonomic assignment

The taxonomic assignment was also performed with the reference database PR2 (version 5.0.0) in case you are not happy with EukRibo. Both classifier (IDTAXA with threshold = 50) and best hit (vsearch, stampa pipeline) approaches were used. More details here.

Références

Berney, Cédric. 2022. “EukRibo: a manually curated eukaryotic 18S rDNA reference database.” Zenodo. https://doi.org/10.5281/zenodo.6327891.
Berney, Cédric, Nicolas Henry, Frédéric Mahé, Daniel J. Richter, and Colomban de Vargas. 2022. “EukRibo: A Manually Curated Eukaryotic 18S rDNA Reference Database to Facilitate Identification of New Diversity.” bioRxiv. https://doi.org/10.1101/2022.11.03.515105.
Callahan, Benjamin J, Paul J McMurdie, Michael J Rosen, Andrew W Han, Amy Jo A Johnson, and Susan P Holmes. 2016. “DADA2: High-Resolution Sample Inference from Illumina Amplicon Data.” Nature Methods 13 (7): 581–83.
Edgar, Robert C, Brian J Haas, Jose C Clemente, Christopher Quince, and Rob Knight. 2011. “UCHIME Improves Sensitivity and Speed of Chimera Detection.” Bioinformatics 27 (16): 2194–2200.
Frøslev, Tobias Guldberg, Rasmus Kjøller, Hans Henrik Bruun, Rasmus Ejrnæs, Ane Kirstine Brunbjerg, Carlotta Pietroni, and Anders Johannes Hansen. 2017. “Algorithm for Post-Clustering Curation of DNA Amplicon Data Yields Reliable Biodiversity Estimates.” Nature Communications 8 (1): 1188.
Mahé, Frédéric, Lucas Czech, Alexandros Stamatakis, Christopher Quince, Colomban de Vargas, Micah Dunthorn, and Torbjørn Rognes. 2022. “Swarm V3: Towards Tera-Scale Amplicon Clustering.” Bioinformatics 38 (1): 267–69.
Mahé, Frédéric, Torbjørn Rognes, Christopher Quince, Colomban De Vargas, and Micah Dunthorn. 2015. “Swarm V2: Highly-Scalable and High-Resolution Amplicon Clustering.” PeerJ 3: e1420.
Mahé, Frédéric, Torbjørn Rognes, Christopher Quince, Colomban de Vargas, and Micah Dunthorn. 2014. “Swarm: Robust and Fast Clustering Method for Amplicon-Based Studies.” PeerJ 2: e593.
Murali, Adithya, Aniruddha Bhargava, and Erik S. Wright. 2018. “IDTAXA: A Novel Approach for Accurate Taxonomic Classification of Microbiome Sequences.” MICROBIOME 6 (140). https://doi.org/10.1186/s40168-018-0521-5.