PipeRline

Marine Surveillance

Introduction

This metabarcoding pipeline is based around the targets package, which is a Make-like pipeline tool for R. The benefit of this is that all the code is automatically run, and the pipeline skips costly runtime for tasks that are already up to date.

This page lists the general workflow and all the options that can be adapted for your dataset.

Clone the pipeRline github repository

The first step is to clone this github repository, which contains the required code and directory structure to run the pipeline. To do this, you will need Git installed on your computer. If you are running the pipeline in Rstudio, it is best to create a new project from the github repository. This can be done by going file > new project > version control > git then adding https://github.com/alexpiper/piperline.git then change the name of the project to whatever you are working with.

Alternatively, the repository can be cloned using the command line, change the ‘folder-name’ to the desired name

# Change into the main directory you wish to make the project in
cd metabarcoding

# Clone the repository
git clone https://github.com/alexpiper/piperline.git folder-name

Prior to analysis

The PipeRline workflow assumes that your sequencing data meets certain criteria:

1. Samples have been demultiplexed, i.e. split into individual per-sample fastq files. If you want the pipeline to calculate the index-switching rate, the fastq files need to be re-demultiplexed as the miseq does not put indexes in fasta headers by default
2. For paired-end sequencing data, the forward and reverse reads are in separate files (ie. not interleaved) with reads arranged in matched order.

The output directory should be unique for each sequencing run, named as the flowcell id, within a directory called data

For example:

root/
  ├── data/
     ├── CJL7D/
     

You will also need BLAST+ installed on your computer, it can be downloaded from the NCBI website

Demultiplex sequencing reads

For this workflow to run, we need to first demultiplex the miseq run again as the miseq does not put indexes in fasta headers by default, and also obtain some necessary files from the sequencing folder. The below code is written for the Agriculture Victoria BASC server, and the locations will be different if you are using a different HPC cluster.

#load module
module load bcl2fastq2/2.20.0-foss-2018b

#raise amount of available file handles
ulimit -n 4000

###Run1

#Set up input and outputs
inputdir=/group/sequencing/230127_M01054_0817_000000000-KWD4N
outputdir=/group/pathogens/IAWS/Projects/Metabarcoding/SVBS_marine_metabarcoding/data/KWD4N #CHANGE TO YOUR DATA FOLDER RUN
samplesheet=/group/pathogens/IAWS/Projects/Metabarcoding/SVBS_marine_metabarcoding/SampleSheet_KWD4N.csv #CHANGE TO YOUR SAMPLESHEET

# convert samplesheet to unix format
dos2unix $samplesheet

#Demultiplex
bcl2fastq -p 12 --runfolder-dir $inputdir \
--output-dir $outputdir \
--sample-sheet $samplesheet \
--no-lane-splitting --barcode-mismatches 1

# Copy other necessary files and move fastqs
cd $outputdir
cp -r $inputdir/InterOp $outputdir
cp $inputdir/RunInfo.xml $outputdir
cp $inputdir/[Rr]unParameters.xml $outputdir
cp $samplesheet $outputdir
mv **/*.fastq.gz $outputdir

# Append fcid to start of sample names if missing
fcid=$(echo $inputdir | sed 's/^.*-//')
for i in *.fastq.gz; do
  if ! [[ $i == $fcid* ]]; then
  new=$(echo ${fcid} ${i}) #append together
  new=$(echo ${new// /_}) #remove any white space
  mv -v "$i" "$new"
  fi
done

Install and load R packages and setup directories

# Load the targets and renv packages
library(targets)
library(tarchetypes)

# Load all packages using renv::restore() and source package list
renv::restore()
source("_targets_packages.R")

# Source ancillary functions
source("R/functions.R")
source("R/themes.R")

Create sample tracking sheet

In order to track samples and relevant QC statistics throughout the metabarcoding pipeline, we will first create a new samplesheet from our input samplesheets. This function requires both the SampleSheet.csv used for the sequencing run, and the runParameters.xml, both of which should have been automatically obtained from the demultiplexed sequencing run folder in the bash step above

runs <- dir("data/") #Find all directories within data
SampleSheet <- list.files(paste0("data/", runs), pattern= "SampleSheet", full.names = TRUE)
runParameters <- list.files(paste0("data/", runs), pattern= "[Rr]unParameters.xml", full.names = TRUE)

# Create samplesheet containing samples and run parameters for all runs
samdf <- create_samplesheet(SampleSheet = SampleSheet, runParameters = runParameters, template = "V4") %>%
  distinct()

# Check that sample_ids contain fcid, if not; attatch
samdf <- samdf %>%
  mutate(sample_id = case_when(
    !str_detect(sample_id, fcid) ~ paste0(fcid,"_",sample_id),
    TRUE ~ sample_id
  ))

# Check that samples match samplesheet
fastqFs <- purrr::map(list.dirs("data", recursive=FALSE),
                      list.files, pattern="_R1_", full.names = TRUE) %>%
  unlist() %>%
  str_remove(pattern = "^(.*)\\/") %>%
  str_remove(pattern = "(?:.(?!_S))+$")

# Filter undetermined reads from sample sheet
fastqFs <- fastqFs[!str_detect(fastqFs, "Undetermined")]

# Check for fastq files that are missing from samplesheet
if (length(setdiff(fastqFs, samdf$sample_id)) > 0) {warning("The fastq file/s: ", setdiff(fastqFs, samdf$sample_id), " are not in the sample sheet") }

# Check for sample_ids that dont have a corresponding fastq file
if (length(setdiff(samdf$sample_id, fastqFs)) > 0) {
  warning(paste0("The fastq file: ",
                 setdiff(samdf$sample_id, fastqFs),
                 " is missing, dropping from samplesheet \n")) 
  samdf <- samdf %>%
    filter(!sample_id %in% setdiff(samdf$sample_id, fastqFs))
}

# Write out sample tracking sheet
write_csv(samdf, "sample_data/Sample_info.csv")

Add PCR primers to sample sheet

This can either be done manually by editing the sample_data/Sample_info.csv file, or it can be done in R as below.

If a single primer set was used across all samples, these can simply be added using a mutate call (tab 1).

If different primers were used for different samples, these can be set using pattern matching on the sample names with case_when (tab 2).

If multiple primer sets are used per sample, before indexing, the pipeline will conduct an extra round of demultiplexing to split each sample by primer. This option can be set by splitting each primer set with a semicolon (tab 3)

Multiple primer sets per sample

# Add PCR primers to sample sheet
samdf <- samdf %>%
  mutate(
    pcr_primers = "mlCOIintF-jgHCO2198;Uni18S-Uni18SR",
    for_primer_seq = "GGWACWGGWTGAACWGTWTAYCCYCC;AGGGCAAKYCTGGTGCCAGC",
    rev_primer_seq = "TAIACYTCIGGRTGICCRAARAAYCA;GRCGGTATCTRATCGYCTT"
    )

write_csv(samdf, "sample_data/Sample_info.csv")

Create parameters file

The parameters file table the respective target gene, reference databases, and filtering parameters for each primer set used to amplify the samples. As the pcr_primers column is used in the pipeline to match the parameters to the respective sample, it is critical that the primer names match those set in the previous step.

If a single primer set was used across all samples, these can simply be added to the table as below (tab 1).

If different primers were used for different samples, these can be set as different rows in the table using c() when creatign the tibble (tab 2).

If multiple primer reference databases are to be used for each sample, these can be set by splitting each a semicolon. The taxonomic assignment will be conducted sequentially through the databases from left to right, with the second reference database only being used for those ASVs that couldnt be assigned to species level using the first databases (tab 3)

Parameter options

• Primer parameters:
  • pcr_primers - Name of PCR primers - must match samdf file
  • target_gene - Name of target gene
  • max_primer_mismatch - How much mismatch to allow when detecting primer sequences
• Read filterign:
  • read_min_length - Minimum length of primer trimmed reads
  • read_max_length - Maximum length of primer trimmed reads
  • read_max_ee - Maximum expected errors of primer trimmed reads
  • read_trunc_length - Length to cut all longer reads to
  • read_trim_left - Remove this many bp from left side of primer trimmed reads
  • read_trim_right - Remove this many bp from right side of primer trimmed reads
• ASV filtering
  • asv_min_length - Minimum length of amplicon
  • asv_max_length - Maximum length of amplicon
  • high_sensitivity - Option to turn on pseudo pooling, setting this to FALSE greatly speeds up ASV inference at the expense of detecting fewer rare (<5 reads) taxa.
  • concat_unmerged - Retain any unmerged read pairs by concatenating them together separated by a string of 10 N bases
  • genetic_code - Genetic code for amplicon - see Biostrings::GENETIC_CODE_TABLE
  • coding - Is the amplicon from a protein coding gene
  • phmm - Path to profile hidden markov model (Optional)
• Taxonomic assignment:
  • idtaxa_db - Path to trained IDTAXA model
  • ref_fasta - Path to fasta file of reference database
  • idtaxa_confidence - Minimum bootstrap confidence for IDTAXA
  • run_blast - Whether a blast top hit search should be conducted in addition to IDTAXA
  • blast_min_identity - Minimum nucleotide identity for BLAST
  • target_kingdom - Subset to target kingdom
  • target_phylum - Subset to target phylum
  • target_class - Subset to target class
  • target_order - Subset to target order
  • target_family - Subset to target family
  • target_genus - Subset to target genus
  • target_species - Subset to target species
• Sample & Taxon filtering:
  • min_sample_reads - Minimum reads per sample after filtering
  • min_taxa_reads - Minimum reads per ASV to retain
  • min_taxa_ra - Minimum relative abundance per ASV to retain. 1e-4 is 0.01%
• General pipeline parameters
  • threads - Number of CPU threads to run the pipeline across

Different or multiple primer sets per sample

# Params to add in step_add_parameters
params <- tibble(
  # Amplicon parameters
  pcr_primers = c("mlCOIintF-jgHCO2198", "Uni18S-Uni18SR"),
  target_gene=c("COI", "18S"),
  max_primer_mismatch=0,

  # Read filtering
  read_min_length = 20,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 224,
  read_trim_left = 0,
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = c(303, 50), #MtCOIintF-jgHCO2198 = 313bp
  asv_max_length = c(323, 500),
  high_sensitivity = TRUE,
  concat_unmerged = TRUE,
  genetic_code = c("SGC4", "SGC0"),
  coding = c(TRUE, FALSE),
  phmm = NA_character_,
  
  # Taxonomic assignment
  idtaxa_db = c("reference/COI_idtaxa.rds","reference/18S_idtaxa.rds"),
  ref_fasta = c("reference/COI_hierarchial.fa.gz", "reference/18S_hierarchial.fa.gz"),
  idtaxa_confidence = 60,
  run_blast=TRUE,
  blast_min_identity = 97,
  blast_min_coverage = 90,
  target_kingdom = NA_character_,
  target_phylum = NA_character_,
  target_class = NA_character_,
  target_order = NA_character_,
  target_family = NA_character_, 
  target_genus = NA_character_,  
  target_species= NA_character_,  
  
  # Sample & Taxon filtering
  min_sample_reads = 0,
  min_taxa_reads= 0, # Minimum reads per ASV to retain
  min_taxa_ra = 0, # Minimum relative abundance per ASV to retain. 1e-4 is 0.01%
  
  # General pipeline parameters
  threads = 1
)

write_csv(params, "sample_data/loci_params.csv")

Run pipeline

Now that the sample data sheet and parameters are defined, the pipeline steps can now be visualised, and then run automatically.

tar_make()

Results

See the general introduction vignette for an in-depth explanation of the results and quality control outputs