PipeRline

BASC example using SLURM

Introduction

The below code is written for the Agriculture Victoria BASC computing cluster.

This workflow currently only works with insect COI sequencing using a single primer set

Clone the pipeRline github repository

# Define the directory you will be running the analysis in
working_dir=/group/pathogens/IAWS/Personal/Alexp/metabarcoding/marine_surveillance #CHANGE TO YOUR DIRECTORY

git clone https://github.com/alexpiper/piperline.git $working_dir
cd $working_dir

Updating the pipeline

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

cd $working_dir
git pull

# if an error occurs, run:
git stash
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/
     
     

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 computing cluster, 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

fcid=$(echo $input_dir | sed 's/^.*-//')
outputdir=$working_dir/data/$fcid
samplesheet=$input_dir/SampleSheet.csv

# convert samplesheet to unix format
dos2unix $samplesheet

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

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

# Append fcid to start of sample names if missing
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

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

Copy reference databases

Reference databases are stored in the referencedata directory on BASC

# Change 'folder-name' to the directory you are running the analysis in
cp /group/referencedata/mspd-db/metabarcoding/arthropod/imappests_coi_18_08_2020/* $working_dir/reference/.

Submit job with parameters

The parameters for the run are parsed along with the slurm script.

if you wish to use multiple values (i.e. multiple reference databases, multiple primers per index) encapsulate them in a ’’ and seperate them with a ;

You will need to change your email address and account code below, see BASC account code lists

# Change 'folder-name' to the directory you are running the analysis in
cd $working_dir

# Submit slurm job - CHANGE EMAIL TO YOUR OWN EMAIL
sbatch --mail-user=your.name@email.com --account=your.account --cpus-per-task=10 \
supplementary_scripts/basc_shifter.slurm \
  --pcr_primers fwhF2-fwhR2n \
  --for_primer_seq GGDACWGGWTGAACWGTWTAYCCHCC \
  --rev_primer_seq GTRATWGCHCCDGCTARWACWGG \
  --target_gene COI \
  --max_primer_mismatch 0 \
  --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_min_length 195 \
  --asv_max_length 215 \
  --high_sensitivity TRUE \
  --concat_unmerged FALSE \
  --genetic_code SGC4 \
  --coding TRUE \
  --phmm reference/folmer_fullength_model.rds \
  --idtaxa_db reference/idtaxa_bftrimmed.rds \
  --ref_fasta reference/insecta_hierarchial_bftrimmed.fa.gz \
  --idtaxa_confidence 60 \
  --run_blast TRUE \
  --blast_min_identity 97 \
  --blast_min_coverage 90 \
  --target_kingdom Metazoa \
  --target_phylum Arthropoda \
  --target_class NA \
  --target_order NA \
  --target_family NA \
  --target_genus NA \
  --target_species NA \
  --min_sample_reads 1000 \
  --min_taxa_reads NA \
  --min_taxa_ra 1e-4 \
  --start_fresh FALSE

Parameter descriptions

• Primer parameters:
  • pcr_primers - Name of PCR primers
  • for_primer_seq - Forward primer sequence
  • rev_primer_seq - Reverse primer sequence
  • 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 - NOT USED FOR BASC
  • start_fresh - Whether to remove all intermediate files and start the pipeline fresh