Tephritid metabarcoding workshop

Bioinformatic analysis

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 an example pipeline for insect COI sequencing.

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

Updating the pipeline

To update to the latest version of the pipeline, run the below code in the terminal.

git pull

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 MiSeq run

For this workflow to run, we will need some sequencing runs to work with. If you are working with MiSeq data, it is recommended that the data is demultiplexed again using bcl2fastq, as the miseq does not put indexes in fasta headers by default which is required for the index swtiching calculation.

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/210219_M03633_0489_000000000-JDYG3 #CHANGE TO YOUR SEQ RUN
outputdir=/group/pathogens/IAWS/Projects/Metabarcoding/dros_surveillance/data/JDYG3 #CHANGE TO YOUR DATA FOLDER RUN
samplesheet=/group/pathogens/IAWS/Projects/Metabarcoding/dros_surveillance/SampleSheet_JDYG3.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

This pipeline depends on a number of other R packages which this step will install. The versions of these packages are managed using renv to ensure they match the versions the pipeline was developed on. The first time installing packages may take some time, but they should be quick to load for any future runs.

# 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")

Reference databases

The pipeRline workflow uses the IDTAXA and BLAST software to assign taxonomy to the sequence reads. These require a pre-trained IDTAXA model, and a fasta file with heirarchial taxonomy. For insects, these

Download reference database

Reference databases for insects and arachnids have been hosted on Zenodo, either download them manually, or with R below:

# Download files from zenodo
download_zenodo(
  doi = "10.5281/zenodo.7655352",
  path = "reference"
)

Create sample tracking sheet

In order to track samples and relevant QC statistics throughout the metabarcoding pipeline, we will first create a new sample tracking sheet from our input illumina 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.

For this workflow a single primer set was used across all samples. For multiple primers see the options in the General introduction to the pipeRline workflow

# Add PCR primers to sample sheet
samdf <- samdf %>%
  mutate(
    pcr_primers = "fwhF2-fwhR2nDac;EIF3LminiF4-EIF3lminiR4",
    for_primer_seq = "GGDACWGGWTGAACWGTWTAYCCHCC;GATGCGYCGTTATGCYGATGC",
    rev_primer_seq = "GTRATWGCHCCIGCTAADACHGG;TTRAAYACTTCYARATCRCC"
    )

#Write out sample CSV for use in pipeline
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.

For this workflow a single primer set was used across all samples. For multiple primers see the options in the General introduction to the pipeRline workflow

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 filtering:
  • 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
  • blast_min_coverage - Minimum query coverage 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

Create parameters file

The pipeline also requires a locus parameters file which lists important information about the target barcodes, as well as the PHMM model used to clean the loci, as well as the different reference databases that are to be used for taxonomic assignment.

If multiple reference databases are to be used iteratively for assignment, they should be split with a semicolon, and in the order they are to be used for assignment.

The below code generates the loci_params.csv file, however this can also be done in excel.

# Params to add in step_add_parameters
params <- tibble(
  # Primer parameters
  pcr_primers = c("fwhF2-fwhR2nDac", "EIF3LminiF4-EIF3lminiR4"),
  target_gene=c("COI", "EIF3L"),
  max_primer_mismatch=1,

  # Read filtering
  read_min_length = 20,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 150,
  read_trim_left = 0,
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = c(195, 207),
  asv_max_length = c(215, 227),
  genetic_code = c("SGC4", "SGC0"),
  coding = c(TRUE, TRUE),
  phmm = c("reference/phmm/Bactrocera_COI.rds", "reference/phmm/Bactrocera_EIF3L.rds"),
  
  # Taxonomic assignment
  idtaxa_db = c("reference/COI_idtaxa.rds","reference/EIF3L_idtaxa.rds"),
  ref_fasta = c("reference/COI_hierarchial.fa.gz", "reference/EIF3L_hierarchial.fa.gz"),
  idtaxa_confidence = 60,
  run_blast=TRUE,
  blast_min_identity = 97,
  blast_min_coverage = 90,
  target_kingdom = "Metazoa",
  target_phylum = "Arthropoda",
  target_class = "Insecta",
  target_order = "Diptera",
  target_family = NA, 
  target_genus = NA,  
  target_species = NA,  
  
  # Sample & Taxon filtering
  min_sample_reads = c(1000, 1000),
  min_taxa_reads= NA, 
  min_taxa_ra = c(1e-4, 1e-4),
  
  # General pipeline parameters
  threads = 1
)

# Write out the parameters file
write_csv(params, "sample_data/loci_params.csv")

Run the pipeline

Now that the sample data and parameters file is ready, its time to run the pipeline!

Run the targets pipeline

tar_make()

Results

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