PipeRline

Bee Health (16S + ITS) example

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 workflow for the Bee Health project (16S + ITS primers)

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/
     

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/230421_M03633_0645_000000000-L29KJ #CHANGE TO YOUR SEQ RUN
outputdir=/group/pathogens/IAWS/Personal/Alexp/jessi_dataset/data/L29KJ #CHANGE TO YOUR DATA FOLDER RUN
samplesheet=/group/pathogens/IAWS/Personal/Alexp/jessi_dataset/SampleSheet_L29KJ.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 = "16SFwd515-16SRev806;gITS7ngs-ITS4ngsUni",
  for_primer_seq = "GTGCCAGCMGCCGCGGTAA;GTGARTCATCRARTYTTTG",
  rev_primer_seq = "GGACTACHVGGGTWTCTAAT;CCTSCSCTTANTDATATGC"
  )

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)

The genetic code parameter is used to check coding markers for stop codons and frameshifts that commonly indicate pseudogenes. For COI amplicons from insects, the invertebrate mitochondrial code is “SGC4”, while the standard genetic code is “SGC0”. If you are targetting other groups of organisms, all alternative genetic codes can be found by running Biostrings::GENETIC_CODE_TABLE in R.

Different or multiple primer sets per sample

# Params to add in step_'add_parameters
params <- tibble(
  # Run parameters
  pcr_primers = c("16SFwd515-16SRev806", "gITS7ngs-ITS4ngsUni"),
  target_gene=c("16S", "ITS2"),

  # Read filtering
  read_min_length = 50,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 225,
  read_trim_left = 0,
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = NA_integer_,
  asv_max_length = NA_integer_,
  genetic_code = NA_character_,
  coding = FALSE,
  phmm = NA_character_,
  
  # Taxonomic assignment
  idtaxa_db = c("reference/GTDB_bac120_arc122_ssu_r202_idtaxa.rds","reference/unite_general_fungi_20221129_idtaxa.rds"),
  ref_fasta = c("reference/GTDB_bac120_arc122_ssu_r202.fa.gz","reference/unite_general_fungi_20221129.fa.gz"),
  idtaxa_confidence = 60,
  run_blast=FALSE,
  blast_min_identity = 97,
  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 = 1000,
  min_taxa_reads= NA_integer_, # Minimum reads per ASV to retain
  min_taxa_ra = NA_real_, # 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

tar_make()

Results

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