PEMA investigating metabarcoding

Looking for PEMA v2.x? See Basics for running PEMA (v2.x).

Basics for running PEMA

Running PEMA takes two things: a folder with your sequence files and a parameters file. Then, it is a single command.

Step 1: Your sequence files

Put the paired-end fastq files of your samples in a folder; .fastq, .fq, .fastq.gz and .fq.gz files are accepted, and any other file in the folder is ignored.

PEMA recognises the two files of each sample by a 1/2 (or R1/R2, in any case) between separators (_, . or -) in their names; the part of the name before it becomes the sample name. For example:

Read 1 Read 2 Sample
S1_1.fastq.gz S1_2.fastq.gz S1
S1_R1.fq.gz S1_R2.fq.gz S1
S1.R1.fastq S1.R2.fastq S1
S1_S1_L001_R1_001.fastq.gz S1_S1_L001_R2_001.fastq.gz S1_S1_L001
site_1_1.fastq.gz site_1_2.fastq.gz site_1

So files downloaded from ENA (ERR654526_1.fastq.gz) and files straight from an Illumina sequencer (..._R1_001.fastq.gz) both work as they are.

Before running anything, PEMA checks the files and lists all the problems it finds, if:

  • the folder does not exist or has no fastq files;
  • a file has no 1/2 in its name, or no sample name before it;
  • a file has no mate, or two files are the same read of a sample;
  • two pairs end up with the same sample name (e.g. S1_R1_001 and S1_R1_002).

Sample names may not contain spaces. Samples sequenced on several lanes (_L001, _L002, …) are treated as separate samples.

Step 2: The parameters file

All the parameters of an analysis are in a single YAML file. Copy pema.yaml from the PEMA repository and fill it in; every parameter is described in it. A parameter left empty is not passed to the tool, so the tool’s own default applies.

The main ones:

Parameter Description Default
raw_reads The folder with your fastq files required
marker_gene 16S, 18S, 12S, COI or ITS required
outdir The output folder results
threads Threads for fastp, cutadapt, Swarm, VSEARCH and crest4 3
clustering_algo swarm or dada2 (ASVs), or vsearch (OTUs) swarm
classifier crest or rdp (RDPClassifier) crest
min_oligotons ASVs (Swarm) or unique sequences (VSEARCH) with a lower total abundance are removed; 0 keeps all 5

and then one section per tool:

Section Parameters
cutadapt ITS only: forward_primer, reverse_primer, and any other cutadapt option
fastp Any fastp option, without the leading --
swarm Any Swarm option, e.g. differences, fastidious
vsearch threshold: the identity threshold of the OTUs
dada2 pool, min_overlap, max_mismatch, chimera_method, nbases
crest database, mount_crest_dbs, search_algo, min_score, … and train for your own database
rdp database, databases_dir, confidence, memory, and train for your own database
ncbi_taxonomy get_ncbi_taxonomy, ncbi_api_key, max_parallel: map the taxonomies to NCBI Taxonomy (see Mapping to NCBI Taxonomy)

PEMA keeps the tools consistent: for example, with Swarm and DADA2, fastp’s n_base_limit is set to 0, as they do not accept Ns, and the Swarm options that only apply to some values of differences are left out otherwise.

There is no single best set of parameters; see Tuning tuning tuning!. For each marker gene, see 16/18S, COI, 12S and ITS.

Step 3: Running PEMA

From the folder where you cloned PEMA:

nextflow run pema.nf --params-file my_analysis.yaml

The parameters file is given with --params-file (two dashes); Nextflow’s own -params-file (one dash) is a different option.

By default the containers run with Docker. To use another container engine, add a profile:

nextflow run pema.nf --params-file my_analysis.yaml -profile singularity

The profiles are docker (default), podman, singularity and apptainer, and slurm to submit the steps as SLURM jobs (see Running on HPC). Adding reports (e.g. -profile docker,reports) writes a Nextflow report of the run, with the time and memory of each step, in pipeline_info/.

While it runs, Nextflow shows the progress of each step. The results are written in outdir; see PEMA output. Nextflow keeps its intermediate files in a work/ folder, in the folder where you ran PEMA; you can delete it once you are happy with the results.

Resuming a run

Add -resume to run PEMA again, re-running only the steps whose input or parameters changed:

nextflow run pema.nf --params-file my_analysis.yaml -resume

This replaces the checkpoints of PEMA v2.x. It is how you tune an analysis: change, for example, the taxonomy assignment parameters and resume; quality control, merging and clustering are not run again. It is also useful when a run was interrupted, for example by “too many requests” errors of the NCBI Taxonomy step.