Metabarcoding workflow for COI amplicon sequencing

The COI region is commonly used for metabarcoding practices and consequently there are many primer options to choose from. The Fisheries team at GMGI has optimized the Leray Geller set (outlined in red box below). Citation: Leray et al 2013.

We have primarily use this set for invertebrate targets and 12S for vertebrate communities.

Analysis workflow:
1. Assess quality of fastq files using FASTQC and MultiQC
2. Identify sequencing outliers and sub-sample if needed
3. Run nf-core/ampliseq to remove adapters, predict ASVs, assign taxonomy, and generate counts tables
4. Run blastn and least common ancestor analysis to supplement assignment with DADA2

Taxonomic identification is completed using 1) NCBI at 90% percent identity + Least Common Ancestor (LCA) and 2) RDP Classifier via DADA2 with MetaZooGene. After filtering and annotation confirmation, the remaining unassigned sequences are *3) assigned to family-level at best from NCBI at 80% identity.

*This step is option for users. E.g., Instead users can use 80% identity for Class level instead of Family.

Project set-up

Typical project folder structure (once analysis is run) for Fisheries team. Create scripts, results (for ampliseq output), fastqc, blast, and metadata. E.g., mkdir results.

blast  fastqc  metadata  multiqc_raw_data  multiqc_raw.html  raw_data  results  scripts

Create/edit a slurm script: nano name.sh
View a file without editing: less name.sh
Start an interactive node: srun --pty bash

Confirm conda environment is available

The conda environment is started within each slurm script, but to activate conda environment outside of the slurm script to update packages or check what is installed:

# Activate conda
source /projects/gmgi/miniconda3/bin/activate

# Activate fisheries eDNA conda environment 
conda activate fisheries_eDNA

# List all available environments 
conda env list 

# List all packages installed in fisheries_eDNA
conda list 

# Update a package
conda update [package name]

# Update nextflow ampliseq workflow 
nextflow pull nf-core/ampliseq

Assess quality of raw data

Background information on FASTQC. This should take ~10 seconds per file for a typical full MiSeq run. If <2 seconds, then check output/error files.

Create slurm script: nano 01-fastqc.sh. Copy below script into file and save (Ctrl+X; Y; Enter).

#!/bin/bash
#SBATCH --error=output/fastqc_output/"%x_error.%j" #if your job fails, the error report will be put in this file
#SBATCH --output=output/fastqc_output/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=20:00:00
#SBATCH --job-name=fastqc
#SBATCH --mem=3GB
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=2

### USER TO-DO ### 
## 1. Set paths for your project

# Activate conda environment
source /projects/gmgi/miniconda3/bin/activate fisheries_eDNA

## SET PATHS (USER EDITS)
raw_path="/projects/gmgi/example/COI/raw_data"
out_dir="/projects/gmgi/example/COI/fastqc"

## CREATE SAMPLE LIST FOR SLURM ARRAY
### 1. Create list of all .gz files in raw data path
ls -d ${raw_path}/*.gz > ${raw_path}/rawdata

### 2. Create a list of filenames based on that list created in step 1
mapfile -t FILENAMES < ${raw_path}/rawdata

### 3. Create variable i that will assign each row of FILENAMES to a task ID
i=${FILENAMES[$SLURM_ARRAY_TASK_ID]}

## RUN FASTQC PROGRAM 
fastqc ${i} --outdir ${out_dir}

To run:
- Start slurm array (e.g., with 138 files) = sbatch --array=0-137 01-fastqc.sh. Note that this starts with file 0.

Count the number of files by navigating to raw data folder and using ls *.gz | wc. The first column is the number of rows in the generated list (files).

Notes:

• This is going to output many error and output files. After job completes, navigate to cd scripts/output/fastqc_output and use cat *output.* > ../fastqc_output.txt to create one file with all the output and cat *error.* > ../fastqc_error.txt to create one file with all of the error message outputs.
• Within the out_dir output folder, use ls *html | wc to count the number of html output files (1st/2nd column values). This should be equal to the --array range used and the number of raw data files. If not, the script missed some input files so address this before moving on.

Visualize quality of raw data

Background information on MULTIQC.

Run the following code in an interactive node. This program will be quick for eDNA data (typically <2 minutes).

# Use srun to claim a node
srun --pty bash 

# Activate conda environment
source /projects/gmgi/miniconda3/bin/activate fisheries_eDNA

## SET PATHS 
## fastqc_output = output from 00-fastqc.sh; fastqc program
fastqc_output="" 
multiqc_dir="" 

## RUN MULTIQC 
### Navigate to your project
multiqc --interactive ${fastqc_output} -o ${multiqc_dir} --filename multiqc_raw.html

Path examples:
- fastqc_output="/projects/gmgi/Fisheries/eDNA/tidal_cycle/fastqc_output"
- multiqc_dir="/projects/gmgi/Fisheries/eDNA/tidal_cycle"

Identify sequencing outliers and sub-sample if needed

This step can be done two ways: 1) download read counts from Illumina basespace run information or 2) Download read counts from the multiqc report. Either way, you need an excel file with the sample ID and total read count.

Read counts can be downloaded from the General Statistics section on Multiqc by selecting Export as CSV > Data > Format (csv) > Download Plot Data. Save this file in the project repository.

We use a median absolute deviation (MAD) approach that works by first finding the median of the data, then measuring how far each value deviates from that median and using the median of those deviations as the typical spread. Values whose deviation is much larger than this typical spread are flagged as outliers.

Visit this page for the R script or download here.

Download seqtk program

This program has already been installed on NU with the below code:

# move to gmgi/packages folder 
cd /projects/gmgi/packages/

# download seqtk program from their github page; move into that folder; and enable (`make`) the program
# ; indicates a pipe much like %>% in R
git clone https://github.com/lh3/seqtk.git;
cd seqtk; make

Path to the program: /projects/gmgi/packages/seqtk/seqtk

Running seqtk with a single (or handful) of samples

Calculate the median number of reads prior to this code. This example will use 70,377.

# start an interactive session on 1 node
srun --pty bash 

# cd to your raw data folder
# make a new folder called outliers (`mkdir` = make directory)
mkdir outliers

# define a path to reference later using ${seqtk}
seqtk="/projects/gmgi/packages/seqtk/seqtk"

# reference the program path above using ${} and subsample the original fastq file 
# to 70,377 reads using a set seed of -s100 following by gzipping that file using a pipe |
# that is similar to %>% in R. I have two commands then: 1) sub-sample 2) gzip 
# > indicates a new output file 

${seqtk} sample -s100 INS2-FS-0011_S284_L001_R1_001.fastq.gz 70377 | gzip > INS2-FS-0011_S284_L001_R1_001_subset.fastq.gz
${seqtk} sample -s100 INS2-FS-0011_S284_L001_R2_001.fastq.gz 70377 | gzip > INS2-FS-0011_S284_L001_R2_001_subset.fastq.gz

# move original files into the outliers folders so they are out of the way of our analysis 
mv INS2-FS-0011_S284_L001_R1_001.fastq.gz outliers/
mv INS2-FS-0011_S284_L001_R2_001.fastq.gz outliers/

Running seqtk with a list of files

Create a samplesheet using 03a-metadata.R below as an example. If you have a csv with 'sampleID', 'forwardReads', and 'reverseReads' columns then turn that csv into a text file with:

tail -n +2 outlier_list_for_subsampling.csv | cut -d',' -f2,3 | tr ',' '\n' > outlier_paths.txt

The input to the below .sh is a txt file with a single column with no header that lists the full path of each fastq file to subset.

02-outlier.sh

#!/bin/bash
#SBATCH --error=output/"%x_error.%j" #if your job fails, the error report will be put in this file
#SBATCH --output=output/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=4:00:00
#SBATCH --job-name=02-outlier
#SBATCH --mem=1GB
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1

### Running seqtk for a list of files

# path to seqtk
seqtk="/projects/gmgi/packages/seqtk/seqtk"

# input list
input_list="/projects/gmgi/example/COI/outlier_paths.txt"

# directory for subset files 
output_directory="/projects/gmgi/example/COI/raw_data"

# directory to archive files
archive_directory="/projects/gmgi/oceanX/example/COI/raw_data/outliers"

# number of reads
nreads=70377
seed=100

# Run seqtk on the input_list files
while read -r fq; do
    [[ -z "$fq" ]] && continue

    base=$(basename "$fq" .fastq.gz)
    out="${output_directory}/${base}_subset.fastq.gz"

    echo "Subsampling: $fq -> $out"
    "${seqtk}" sample -s"${seed}" "$fq" "${nreads}" | gzip > "$out"

done < "$input_list"


# Move input_list files to the archive folder
mkdir -p "$archive_directory"

while read -r fq; do
    [[ -z "$fq" ]] && continue

    echo "Archiving: $fq -> $archive_directory"
    mv "$fq" "$archive_directory/"

done < "$input_list"

Downloading and updating reference databases

Download and/or update NBCI blast nt database

NCBI is updated daily and therefore needs to be updated each time a project is analyzed. This is the not the most ideal method but we were struggling to get the -remote flag to work within slurm because I don't think NU slurm is connected to the internet? NU help desk was helping for awhile but we didn't get anywhere.

Within /projects/gmgi/databases/ncbi, there is a update_nt.sh script with the following code. To run sbatch update_nt.sh. This won't take long as it will check for updates rather than re-downloading every time.

/projects/gmgi/databases/ncbi/update_nt.sh (this script is already in shared gmgi folder):

#!/bin/bash
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=24:00:00
#SBATCH --job-name=update_ncbi_nt
#SBATCH --mem=50G
#SBATCH --output=%x_%j.out
#SBATCH --error=%x_%j.err

# Activate conda environment
source /projects/gmgi/miniconda3/bin/activate fisheries_eDNA

# Create output directory if it doesn't exist
cd /projects/gmgi/databases/ncbi/nt

# Update BLAST nt database
update_blastdb.pl --decompress nt

# Print completion message
echo "BLAST nt database update completed"

This can take hours so plan ahead - start updating NCBI nt database when you know you are a couple days out from analyzing data.

View the update_ncbi_nt.out file to confirm the echo printed at the end.

Download Taxonkit

Download taxonkit to your home directory. You only need to do this once, otherwise the program is set-up and you can skip this step. Taxonkit program is in /projects/gmgi/databases/taxonkit but the program needs required files to be in your home directory.

# Move to home directory and download taxonomy information
cd ~
wget -c ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz

# Extract the downloaded file
tar -zxvf taxdump.tar.gz

# Create the TaxonKit data directory 
mkdir -p $HOME/.taxonkit

# Copy the required files to the TaxonKit data directory
cp names.dmp nodes.dmp delnodes.dmp merged.dmp $HOME/.taxonkit

Download MetaZooGene

Check the following webpages for the most recent release's version number. If that version number is the same as the database we have on our server, then you don't need to update the database. If not, we need to download the most recent version.

MetaZooGene North Altantic
MetaZooGene Global

Instructions to update:

Open RStudio in NU Discovery Cluster OOD and run the following R script to update both the Global and Atlantic versions of MetaZooGene. User needs to update the file name in base::write() to reflect the correct version number. This script is very quick, plan for ~3-5 minutes total to run the script, change the version number, and confirm output.

Path to R script: /projects/gmgi/databases/COI/MetaZooGene/MZG_to_DADA2.R (this script is already in gmgi shared folder).

## GLOBAL DB 

library(magrittr)
input_file = "https://www.st.nmfs.noaa.gov/copepod/collaboration/metazoogene/atlas/data-src/MZGdata-coi__MZGdbALL__o00__A.csv.gz"


Global_MZG <- 
  # Download
  readr::read_csv(input_file, 
                              col_names = FALSE, 
                              trim_ws = TRUE, 
                              na = c("", "NA", "N/A"),
                              col_types = cols(.default = "c"),
                              show_col_types = FALSE) %>%
  # Select ID, Sequence and taxa columns
  dplyr::select(X34, X31, X33, X2) 

Edited_Global_MZG <- Global_MZG %>% 
  # Separate taxa column into separate columns
  tidyr::separate(X34, into = as.character(1:21), sep = ";") %>%

  # Merge wanted taxa columns and add ">"
  tidyr::unite(col = "Taxa", 1,4,5,7,8,9,10,11,12,16,18,X2, sep = ";" ) %>%
  dplyr::mutate(Taxa = paste(">", Taxa, sep = "")) %>%

  # Restructure into one single column vector
  dplyr::select(X33, Taxa, X31) %>% 
  tidyr::pivot_longer(cols = c("Taxa", "X31")) %>%
  dplyr::pull(value)

# Write vector to fasta 
### USER TO CHANGE VERSION NUMBER
Edited_Global_MZG %>% base::write("/projects/gmgi/databases/COI/MetaZooGene/DADA2_MZG_v2023-m07-15_Global_modeA.fasta")

## ATLANTIC DB

input_file_ATL = "https://www.st.nmfs.noaa.gov/copepod/collaboration/metazoogene/atlas/data-src/MZGdata-coi__T4000000__o02__A.csv.gz"

ATL_MZG <- 
  # Download
  readr::read_csv(input_file_ATL, 
                  col_names = FALSE, 
                  trim_ws = TRUE, 
                  na = c("", "NA", "N/A"),
                  col_types = cols(.default = "c"),
                  show_col_types = FALSE) %>%
  # Select ID, Sequence and taxa columns
  dplyr::select(X34, X31, X33, X2) 

Edited_ATL_MZG <- ATL_MZG %>% 
  # Separate taxa column into separate columns
  tidyr::separate(X34, into = as.character(1:21), sep = ";") %>%

  # Merge wanted taxa columns and add ">"
  tidyr::unite(col = "Taxa", 1,4,5,7,8,9,10,11,12,16,18,X2, sep = ";" ) %>%
  dplyr::mutate(Taxa = paste(">", Taxa, sep = "")) %>%

  # Restructure into one single column vector
  dplyr::select(X33, Taxa, X31) %>% 
  tidyr::pivot_longer(cols = c("Taxa", "X31")) %>%
  dplyr::pull(value)

# Write vector to fasta
### USER TO CHANGE VERSION NUMBER
Edited_ATL_MZG %>% base::write("/projects/gmgi/databases/COI/MetaZooGene/DADA2_MZG_v2023-m07-15_NorthAtlantic_modeA.fasta")

Step 4: nf-core/ampliseq

Nf-core: A community effort to collect a curated set of analysis pipelines built using Nextflow.
Nextflow: scalable and reproducible scientific workflows using software containers, used to build wrapper programs like the one we use here.

[https://nf-co.re/ampliseq/2.11.0]: nfcore/ampliseq is a bioinformatics analysis pipeline used for amplicon sequencing, supporting denoising of any amplicon and supports a variety of taxonomic databases for taxonomic assignment including 16S, ITS, CO1 and 18S.

We use ampliseq for the following programs:

• Cutadapt is trimming primer sequences from sequencing reads. Primer sequences are non-biological sequences that often introduce point mutations that do not reflect sample sequences. This is especially true for degenerated PCR primer. If primer trimming would be omitted, artifactual amplicon sequence variants might be computed by the denoising tool or sequences might be lost due to become labelled as PCR chimera.
• DADA2 performs fast and accurate sample inference from amplicon data with single-nucleotide resolution. It infers exact amplicon sequence variants (ASVs) from amplicon data with fewer false positives than many other methods while maintaining high sensitivity. This program also assigns taxonomy using a RDP classifier.

COI primer sequences (required)

Below is what we used for COI amplicon sequencing. This results in ~313 bp expected ASV.

LG COI amplicon F: GGWACWGGWTGAACWGTWTAYCCYCC
LG COI amplicon R: TAIACYTCIGGRTGICCRAARAAYCA

Ampliseq will automatically calculate the reverse compliment and include this for us.

Metadata sheet (optional)

The metadata file has to follow the QIIME2 specifications. Below is a preview of the sample sheet used for this test. Keep the column headers the same for future use. The first column needs to be "ID" and can only contain numbers, letters, or "-". This is different than the sample sheet. NAs should be empty cells rather than "NA".

Create samplesheet sheet for ampliseq

This file indicates the sample ID and the path to R1 and R2 files. Below is a preview of the sample sheet used in this test. File created on RStudio Interactive on Discovery Cluster using (create_metadatasheets.R).

• sampleID (required): Unique sample IDs, must start with a letter, and can only contain letters, numbers or underscores (no hyphons!).
• forwardReads (required): Paths to (forward) reads zipped FastQ files
• reverseReads (optional): Paths to reverse reads zipped FastQ files, required if the data is paired-end
• run (optional): If the data was produced by multiple sequencing runs, any string

sampleID	forwardReads	reverseReads	run
sample1	./data/S1_R1_001.fastq.gz	./data/S1_R2_001.fastq.gz	A
sample2	./data/S2_fw.fastq.gz	./data/S2_rv.fastq.gz	A
sample3	./S4x.fastq.gz	./S4y.fastq.gz	B
sample4	./a.fastq.gz	./b.fastq.gz	B

This is an R script, not slurm script. Open RStudio interactive on Discovery Cluster to run this script.

Prior to running R script, use the rawdata file created for the fastqc slurm array from within the raw data folder to create a list of files. Below is an example from our Offshore Wind project but the specifics of the sampleID will be project dependent. This project had four sequencing runs with different file names.

On OOD, open RStudio. Create a new R script file:

03-metadata.R

## Creating ampliseq metadata sheet for GOBLER
## Load libraries 
library(dplyr)
library(stringr)
library(strex) 

### Read in sample sheet 
sample_list <- read.delim2("/projects/gmgi/ecosystem-diversity/Gobler/COI/raw_data/rawdata", header=F) %>% 
  dplyr::rename(forwardReads = V1) %>%
  mutate(sampleID = str_after_nth(forwardReads, "data/", 1),

         ### Confirm that _S is correct for your file names, sometimes this is _R
         sampleID = str_before_nth(sampleID, "_S", 1))

         ## add any other mutate functions as needed for sample names
         ## sampleID needs to match a column in your metadata file for later analyses 

# creating sample ID 
sample_list$sampleID <- gsub("-", "_", sample_list$sampleID)

# keeping only rows with R1
sample_list <- filter(sample_list, grepl("R1", forwardReads, ignore.case = TRUE))

# duplicating column 
sample_list$reverseReads <- sample_list$forwardReads

# replacing R1 with R2 in only one column 
sample_list$reverseReads <- gsub("R1", "R2", sample_list$reverseReads)

# rearranging columns 
sample_list <- sample_list[,c(2,1,3)]

sample_list %>% write.csv("/projects/gmgi/ecosystem-diversity/Gobler/COI/metadata/samplesheet.csv", 
                          row.names=FALSE, quote = FALSE)

Run nf-core/ampliseq (Cutadapt & DADA2)

Update ampliseq workflow if needed: nextflow pull nf-core/ampliseq.

Create slurm script: nano 04-ampliseq.sh. Copy below script into file and save (Ctrl+X; Y; Enter).

Set trunclenf and trunclenr based on data quality.

#!/bin/bash
#SBATCH --error=output/"%x_error.%j" #if your job fails, the error report will be put in this file
#SBATCH --output=output/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=48:00:00
#SBATCH --job-name=ampliseq_project_marker
#SBATCH --mem=30GB
#SBATCH --ntasks=12
#SBATCH --cpus-per-task=2

### USER TO-DO ### 
## 1. Set paths for project 
## 2. Adjust SBATCH options above (time, mem, ntasks, etc.) as desired  
## 3. Fill in F primer information based on primer type (no reverse compliment needed)
## 4. Adjust parameters as needed (below is Fisheries team default for COI)

# Load nextflow conda environment
source /projects/gmgi/miniconda3/bin/activate nextflow_env

# SET PATHS (PROJECT SPECIFIC)
metadata="" 
output_dir=""

# SET PATHS (EDIT VERSION # ON MZG IF NEEDED)
assignTaxonomy="/projects/gmgi/databases/COI/MetaZooGene/DADA2_MZG_v2023-m07-15_Global_modeA.fasta"
taxlevels="Kingdom,Phylum,Subphylum,Superclass,Class,Subclass,Infraclass,Superorder,Order,Family,Genus,Species"

nextflow run nf-core/ampliseq -resume \
   -c /projects/gmgi/Fisheries/scripts/ampliseq_minBoot.config \
   -profile singularity \
   --input ${metadata}/samplesheet.csv \
   --FW_primer "GGWACWGGWTGAACWGTWTAYCCYCC" \
   --RV_primer "TAIACYTCIGGRTGICCRAARAAYCA" \
   --outdir ${output_dir} \
   --trunclenf 220 \
   --trunclenr 220 \
   --trunc_qmin 25 \
   --max_ee 2 \
   --min_len_asv 300 \
   --max_len_asv 330 \
   --dada_ref_tax_custom ${assignTaxonomy} \
   --dada_assign_taxlevels ${taxlevels} \
   --sample_inference pseudo \
   --ignore_failed_trimming

To run:
- sbatch 04-ampliseq.sh

Path examples:
- metadata="/projects/gmgi/ecosystem-diversity/Gobler/COI/metadata"
- output_dir="/prokects/gmgi/ecosystem-diversity/Gobler/COI/results_Global"

Files generated by ampliseq

Pipeline summary reports:

• summary_report/
• summary_report.html: pipeline summary report as standalone HTML file that can be viewed in your web browser.
• *.svg*: plots that were produced for (and are included in) the report.
• versions.yml: software versions used to produce this report.

Preprocessing:

• FastQC: fastqc/ and *_fastqc.html: FastQC report containing quality metrics for your untrimmed raw fastq files.
• Cutadapt: cutadapt/ and cutadapt_summary.tsv: summary of read numbers that pass cutadapt
• MultiQC: multiqc, multiqc_data/, multiqc_plots/ with multiqc_report.html: a standalone HTML file that can be viewed in your web browser;

ASV inferrence with DADA2:

• dada2/, dada2/args/, data2/log/
• ASV_seqs.fasta: Fasta file with ASV sequences.
• ASV_table.tsv: Counts for each ASV sequence.
• DADA2_stats.tsv: Tracking read numbers through DADA2 processing steps, for each sample.
• DADA2_table.rds: DADA2 ASV table as R object.
• DADA2_table.tsv: DADA2 ASV table.
• dada2/QC/
• *.err.convergence.txt: Convergence values for DADA2's dada command, should reduce over several magnitudes and approaching 0.
• *.err.pdf: Estimated error rates for each possible transition. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. The estimated error rates (black line) should be a good fit to the observed rates (points), and the error rates should drop with increased quality.
• *_qual_stats.pdf: Overall read quality profiles: heat map of the frequency of each quality score at each base position. The mean quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position.
• *_preprocessed_qual_stats.pdf: Same as above, but after preprocessing.

We add an ASV length filter that will output asv_length_filter/ with:

• ASV_seqs.len.fasta: Fasta file with filtered ASV sequences.
• ASV_table.len.tsv: Counts for each filtered ASV sequence.
• ASV_len_orig.tsv: ASV length distribution before filtering.
• ASV_len_filt.tsv: ASV length distribution after filtering.
• stats.len.tsv: Tracking read numbers through filtering, for each sample.

Step 5: Running additional NCBI taxonomic assignment ID script

Ampliseq uses MetaZooGene Global to assign taxonomy but we also use NCBI Eukaryotic database sequences and blastn to find 90+% hits then use Least Common Ancestor approach. After filtering and confirming annotations, we then take NCBI hits at 80+% at family-level only.

Create new folder blast in project folder.

Blast NCBI at 80% identity

Create slurm script: nano 05-blast_NCBI.sh. Copy below script into file and save (Ctrl+X; Y; Enter).

#!/bin/bash
#SBATCH --error=output/"%x_error.%j" #if your job fails, the error report will be put in this file
#SBATCH --output=output/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=20:00:00
#SBATCH --job-name=tax_ID
#SBATCH --mem=30GB
#SBATCH --ntasks=12
#SBATCH --cpus-per-task=2

### USER TO-DO ### 
## 1. Set paths for project; change db path if not 12S

# Activate conda environment
source /projects/gmgi/miniconda3/bin/activate fisheries_eDNA

# SET PATHS (PROJECT SPECIFIC)
## ASV fasta file path, excluding the file name (already included below)
ASV_fasta=""
out=""

# SET PATHS (NO EDITS)
ncbi="/projects/gmgi/databases/ncbi/nt"
taxonkit="/projects/gmgi/databases/taxonkit"

#### DATABASE QUERY ####
### NCBI database 
blastn -db ${ncbi}/"nt" \
   -query ${ASV_fasta}/ASV_seqs.len.fasta -num_threads 16 \
   -taxids 2759 -culling_limit 50 -max_target_seqs 50 -qcov_hsp_perc 95 -evalue 1e-30 \
   -perc_identity 80 \
   -out ${out}/BLASTResults_NCBI.txt \
   -outfmt '6  qseqid   sseqid   sscinames   staxid pident   length   mismatch gapopen  qstart   qend  sstart   send  evalue   bitscore'

############################

#### TAXONOMIC CLASSIFICATION #### 
## creating list of staxids from all three files 
awk -F $'\t' '{ print $4}' ${out}/BLASTResults_NCBI.txt | sort -u > ${out}/NCBI_sp.txt

## annotating taxid with full taxonomic classification
cat ${out}/NCBI_sp.txt | ${taxonkit}/taxonkit reformat -I 1 -r "Unassigned" > ${out}/NCBI_taxassigned.txt

To run:
- sbatch 05-blast_NCBI.sh

Reformat taxonkit output for LCA

Create 06-blast_sort.R and use RStudio interface on NU OOD to run this script to create a taxid list to use as input for LCA.

## Sorting BLAST output for LCA 
library(tidyverse)

### Load data: user edits project path
NCBI_input="/projects/gmgi/example/COI/blast/BLASTResults_NCBI.txt"
TAXID_output="/projects/gmgi/example/COI/blast/taxids.txt"

blast <- read.table(NCBI_input, header=F,
                    col.names = c("ASV_ID", "sseqid", "sscinames", "staxid", "pident", "length", "mismatch",
                                  "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore"),
                    colClasses = c(rep("character", 3), "integer", rep("numeric", 9)))

### Print number of ASVs that had a BLAST output
length(unique(blast$ASV_ID))

### Filter to top pident; this will include ties so any entries with the top pident
blast_filtered <- blast %>% group_by(ASV_ID) %>% slice_max(pident, n=1) 

## nrow should = the number of unique ASVs 
staxid_list <- blast_filtered %>% 
  dplyr::select(ASV_ID, staxid) %>%
  group_by(ASV_ID) %>%

  ## combine tax ids into one row and separate with a space 
  summarise(tax_list = paste(staxid, collapse = " ")) 

### USER EDITS PATH
staxid_list %>% write_tsv(TAXID_output)

Run LCA on BLAST output

Running LCA on Blast output. At this point, the blast hits have been filtered down to the top hit(s). LCA is performed on those top hit(s) rather than everything above 80%. This way, we can take 90% in our first step then 80% later on.

Create slurm script: nano 07-LCA.sh. Copy below script into file and save (Ctrl+X; Y; Enter).

#!/bin/bash
#SBATCH --error=output/"%x_error.%j" #if your job fails, the error report will be put in this file
#SBATCH --output=output/"%x_output.%j" #once your job is completed, any final job report comments will be put in this file
#SBATCH --partition=short
#SBATCH --nodes=1
#SBATCH --time=20:00:00
#SBATCH --job-name=LCA
#SBATCH --mem=10GB
#SBATCH --ntasks=2
#SBATCH --cpus-per-task=2

### USER TO-DO ### 
## 1. Set paths for project 

# SET PATHS (PROJECT SPECIFIC)
taxid_list="/projects/gmgi/example/COI/blast/taxids.txt"
out="/projects/gmgi/example/COI/blast"

# SET PATHS (NO EDITS)
taxonkit="/projects/gmgi/databases/taxonkit"

###################

# Output file for LCA results
lca_results_file="${out}/lca_results.txt"
lca_reformat_file="${out}/lca_results_reformatted.txt"

##### LCA by TaxonKit #####
## Running LCA 
${taxonkit}/taxonkit lca ${taxid_list} -i 2 -o ${lca_results_file}

## Reformatting output 
${taxonkit}/taxonkit reformat ${lca_results_file} -I 3 -o ${lca_reformat_file}

Download files for R analysis

In working folder (for Fisheries team, this is your folder on Box), create the following folders: data, metadata, and scripts. Within data, make a folder for Taxonomic_assignments.

Download and place in data:
- results/summary_report/summary_report.html
- results/overall_summary.tsv
- results/pipeline_info/execution_report_*date*.html
- results/asv_length_filter/ASV_seqs.len.fasta
- results/ASV_table.len.tsv

Download and place in data/Taxonomic_assignments:
- results/phyloseq/dada2_phyloseq.rds
- blast/NCBI_taxassigned.txt
- blast/lca_results_reformatted.txt