This is the bioinformatics pipeline developed for the PacMAN (Pacific Islands Marine Bioinvasions Alert Network) project. This pipeline cleans and classifies sequences from eDNA samples. The steps in this pipeline are compiled from publicly available bioinformatic pipelines like ANACAPA, tourmaline, tagseq-qiime2-snakemake, pema, CASCABEL and MBARI-BOG. The pipeline is based on the snakemake workflow management system.
The pipeline has the following steps:
The following files are required to run the pipeline:
sample-id, file-path and direction (forward or reverse)occurrenceStatus as absent
--> The ASVs from these samples will be removed from all samples, before the occurrence table is madePROJECT name: Usually a specific sample setRUN name: the run with a specific combination of samples and/or parameters for the analysisSAMPLE_SET: manifest file pathsample-data-file: manifest file path:rocket: Quick start - If you want to try out the pipeline with example data, use the example config files and download the example reference databases. Also download the following data files: forward and reverse. Use the following folder structure:
└── data ├── config_files │ ├── config.yaml │ ├── manifest.csv │ └── sample_data.csv ├── raw_sequences │ ├── USP-24-01-172_S172_L001_R1_001.fastq.gz │ └── USP-24-01-172_S172_L001_R2_001.fastq.gz └── reference_databases ├── COI_ncbi_1_50000_pcr_pga_taxon_derep_clean_sintax.fasta └── COI_terrimporter ├── bergeyTrainingTree.xml ├── genus_wordConditionalProbList.txt ├── logWordPrior.txt ├── rRNAClassifier.properties └── wordConditionalProbIndexArr.txt
Install conda and snakemake. The pipeline is run with the --use-conda flag, where each rule has an isolated environment, that will be installed in the working directory using conda.
Download the repository and navigate to the root directory. Add your data files and reference databases to the data folder in the project root. Modify the example config file, manifest, and sample metadata file. Ensure that the reference databases are properly linked to from the config file, and that the file paths in the manifest match the location of your data files.
Perform a dry-run of the pipeline for testing with:
snakemake --use-conda --configfile data/config_files/config.yaml --rerun-incomplete --printshellcmds --cores 1 -npRemoving the -np flag will initiate the run.
The pipeline can also be run with Docker. In this case you will not need a copy of the pipeline code. After installing Docker, prepare a data folder in your workspace as explained in the data preparation section above. Then run the pipeline with the following command (change the paths to your workspace folder):
docker run --platform linux/amd64 -v $(pwd)/data:/pipeline/data -v $(pwd)/results:/pipeline/results -v $(pwd)/.snakemake:/pipeline/.snakemake pieterprovoost/pacman-pipelineOn Windows, you may have to use %cd% instead of $(pwd).
Download example output for a single sample here.
The snakemake workflow will run the following jobs:
The run will first initiate a folder structure in the results folder as follows
PROJECT
├── samples
| ├── sample_1
| │ ├── fastp
| │ ├── qc
| │ └── rawdata
| │ ├── forward (link to sample file)
| │ └── reverse (link to sample file)
| ├── sample_2
| ├── sample_3
| | .
| | .
| | .
| ├── sample_n
| └── multiqc_RUN.html
└── runs
└── RUN
├── 02-cutadapt
├── 03-dada2
├── 04-taxonomy
├── 05-dwca
└── 06-reportSamples will be linked to the file structure and their quality will be analysed with fastqc. All quality files of the raw sequence files will be summarized with multiqc, and can be found in /PROJECT/samples/multqc_RUN.html.
The sequences are trimmed and primers are removed utilising fastp and cutadapt. fastp automatically recognizes and removes illumina adapters. fastp is used as it also detects common artifacts of NovaSeq datasets. The primers must be added to the config file in both forward and reverse (reverse complement) directions.
ASVs are inferred with dada2, which is run in 2 steps. Initially filtering of samples is done based on user-defined parameters. The quality of sequences before and after this filtering is shown in aggregate in 2 plots (06-reports/dada2), and can be found separately for each samples in the 03-dada2/quality folder.
dada2 returns the ASV-table (03-dada2/seqtab-nochim.txt), as well as the sequences of each asv (03-dada2/rep-seqs.fna). In addition the number of reads filtered at each step and remaining after sample processing are returned in the 06-report/dada2_stats.txt file, and will be added to the report.
Taxonomy is annotated using two different methods. Annotation using the naive bayesian classifier rdp is the main method, and the results are filtered based on a user-defined probability threshold (i.e. 0.8). The output provides a confidence threshold for each taxonomic level. For rdp a trained classifier is required. Many trained classifiers can be found already online (see for example https://github.com/terrimporter/CO1Classifier), but a classifier can also be built from a custom reference database by following the steps in the tutorial of John Quensen.
To confirm the results with a local alignment to the reference database, also vsearch is used. The parameters of vsearch can be defined by the user. We use it to find matches with very high sequence similarity to the query (97-100%), to confirm the results of rdp as well as detect any disrepancies, or unclear results in the reference database. The top 10 hits are filtered with a built-in lca algorithm, to find the consensus taxonomy. The reference database format required for this process is the sintax format.
Note As two different taxonomic assignment methods are used to enable confirmation of the results, two different reference database formats are also necessary.
There is an option in the pipeline to further classify sequences that remained unclassified with BLASTn against the full ncbi nt database. However multiple resources are required for this to work. We recommend having a local copy of the full NCBI nt database available to run this step with the pipeline. If you have in total <10kbp of data (50 unknown sequences of 200 bp), you may also run the query in remote mode. We may include a loop to do this with more data at a later stage, but running the analysis remotely for more sequences will require a lot of time.
The user will need to also have access to NCBI-nt to TaxonID mapping files to get the scientific names of the sequences. The pipeline uses BASTA to filter and classify the Blast results based on an lca analysis. If a tax database is not provided in the configfile, the pipeline will prompt BASTA to download the tax_db (gb) to the resources folder. This will also take a long time.
In the final steps of the pipeline LSIDs are defined for the assigned taxonomic names, and the occurrence table and dna-derived data extension table are built for submitting into OBIS.
This step also returns a table 05-dwca/Taxa_not_in_worms.csv, containing the taxonomic names and linked asvs that were not given an lsid. This table will require manual inspection, and possibly contacting the WoRMS team.
In this step the unknown sequences are given the ID for 'Incertae sedis'. Non-marine species (most taxa with no lsid), and ASVs found in the control sample(s) (marked by occurrenceStatus: absent) are not added to the final dwca-tables. All of these can still be found in the table 05-dwca/Full_tax_table_with_lsids as well as the 05-dwca/phyloseq_object.rds, which can be read with the phyloseq R package for further analysis and visualization.
Note! With this strategy, sequences that are known but not marine, are not included in the occurrence tables, while sequences that are not known are always included (as 'Incertae sedis').
Fields that also still need to be added/modified based on the genetic data guidelines are:
identificationReferences: Website of this pipelineAn HTML report is made in the final steps with the statistics of the full run, to give an overview of what was done during the analysis and what the effect was on the results. Still more analysis will be added to this report.