PipeRline

Marine surveillance (COI only) on BASC

Introduction

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

This workflow is designed for the marine surveillance project using the mlCOIintF-jgHCO2198 COI 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

IF you already have the pipeline cloned and wish to update to the latest version, run the below code in the terminal.

git pull

# if an error occurs, run:
git stash
git pull

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
input_dir=/group/sequencing/240426_M01054_0007_000000000-LJRVM #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

Copy reference databases

Reference databases are stored in the referencedata directory on BASC

The Marine surveillance workflow first assigns ASVs against an internal database, then anything that remains unidentified to species is then assigned against a broader scope MARES database

# Change 'folder-name' to the directory you are running the analysis in
cp /group/referencedata/mspd-db/metabarcoding/marine/COI_deakin/* $working_dir/reference/.
cp /group/referencedata/mspd-db/metabarcoding/marine/COI_mares/* $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

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 mlCOIintF-jgHCO2198 \
  --for_primer_seq GGWACWGGWTGAACWGTWTAYCCYCC \
  --rev_primer_seq TAIACYTCIGGRTGICCRAARAAYCA \
  --target_gene COI \
  --max_primer_mismatch 0 \
  --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_min_length 303 \
  --asv_max_length 323 \
  --high_sensitivity TRUE \
  --concat_unmerged FALSE \
  --genetic_code SGC4 \
  --coding TRUE \
  --phmm NA \
  --idtaxa_db 'reference/COI_taxreturnmodel_idtaxa.rds;reference/mares_COI_taxreturnmodel_idtaxa.rds' \
  --ref_fasta 'reference/COI_taxreturnmodel_hierarchial.fa.gz;reference/mares_COI_taxreturnmodel_hierarchial.fa.gz' \
  --idtaxa_confidence 60 \
  --run_blast TRUE \
  --blast_min_identity 97 \
  --blast_min_coverage 90 \
  --target_kingdom NA \
  --target_phylum NA \
  --target_class NA \
  --target_order NA \
  --target_family NA \
  --target_genus NA \
  --target_species NA \
  --min_sample_reads 0 \
  --min_taxa_reads 0 \
  --min_taxa_ra 0 \
  --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