BM3331 / BC5203

This dataset comes from the following paper;

Microbiota composition and evenness predict survival rate of oysters confronted to Pacific Oyster Mortality Syndrome

The full paper can be found at this link: Clerissi et al 2020

Pacific Oyster Mortality Syndrome (POMS) kills juvenile Pacific oysters (Crassostrea gigas) around the world. It is thought to be triggered by a herpesvirus, OsHV-1, which weakens the oyster’s defences so that opportunistic bacteria can take over. The authors bred three families of oysters in a hatchery. One family (RF21) came from parents selected over four generations for resistance to the disease. The other two (SF15 and SF32) came from wild parents and are susceptible.

Oysters from each family were either kept in the hatchery or placed on an oyster farm for 5 days during a disease outbreak (the field). Samples were taken at the end of this 5 day period, which the researchers assumed would be long enough to predispose animals to the disease but not enough for disease to develop. After this all of the oysters were returned to the hatchery and monitored for survival. By day 13 every SF15 and SF32 oyster had died and every RF21 oyster was alive.

Each sample is one whole oyster, and the V3–V4 region of the bacterial 16S rRNA gene was sequenced from it.

Family Phenotype Hatchery Field
RF21 resistant 14 14
SF15 susceptible 14 14
SF32 susceptible 13 14

Data Setup

This assignment starts from raw reads. You won’t need to demultiplex but you will need to perform all the steps after demultiplexing in the 16S microbiome workflow. This includes long-running steps such as removing primers and denoising the reads yourself. Start early because there will be times when you need to wait for jobs to finish on the server.

First create a directory raw_data for storing the original data for this project and another directory called cache for storing your own analysis outputs

mkdir raw_data
mkdir cache

The data for this assignment is located at the path /pvol/data/oyster on the class rstudio server. Run the following commands to create a symbolic link to this in the raw_data directory of your assignment

cd raw_data
ln -s /pvol/data/oyster .

The following files are present within the oyster data directory;

File Name Purpose
fastq/ Raw reads. There are two files per oyster, <sample-id>_R1.fastq.gz and <sample-id>_R2.fastq.gz. These are Illumina MiSeq 2 x 250 bp paired-end reads in compressed fastq format.
manifest.tsv A manifest listing the fastq files for each sample, for importing the reads with qiime tools import
classifier/silva-138-99-341f-805r-nb-classifier.qza Naive Bayesian classifier trained using SILVA 138. It was trained on the region amplified by the 341F and 805R primers used in this study and is appropriate for use with the --i-classifier argument required by qiime feature-classifier classify-sklearn.

The most important columns in the sample metadata file are;

Importing Reads using qiime

Unlike the metagenomics walkthrough exercise, this time we have pre-demultiplexed data (much more common these days) so the workflow for importing data into a qiime object is slightly different. In this case we will supply a file manifest.tsv with all the file paths for the raw fastq files. QIIME will read all these files and import them to a single object. This is the command you’ll need

qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path raw_data/oyster/manifest.tsv \
  --input-format PairedEndFastqManifestPhred33V2 \
  --output-path cache/demux.qza

Trimming reads

This data will need trimming to remove adapters. You should use the cutadapt plugin in qiime for this as follows;

qiime cutadapt trim-paired \
  --i-demultiplexed-sequences cache/demux.qza \
  --p-front-f CCTACGGGNGGCWGCAG \
  --p-front-r GACTACHVGGGTATCTAATCC \
  --p-discard-untrimmed \
  --p-cores 4 \
  --o-trimmed-sequences cache/trimmed.qza

After this finishes I recommend you create a visualisation for the trimmed reads as follows

qiime demux summarize \
  --i-data cache/trimmed.qza \
  --o-visualization cache/trimmed.qzv

Use trimmed.qzv to look at read quality before you choose truncation lengths for denoising.

Running qiime

Long-running commands should be run as slurm jobs, as you did in the metagenomics walkthrough. For these you’ll need to make sure you define the qiime alias as the top of your script. Here is an example script which performs the import command shown above.

#!/bin/sh
#SBATCH --time=60
#SBATCH --ntasks=1 --cpus-per-task=4 --mem=2gb

alias qiime='apptainer run -B /pvol/:/pvol --env "MPLCONFIGDIR=$HOME/matplotlib" /pvol/data/sif/qiime.sif qiime'

qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path raw_data/oyster/manifest.tsv \
  --input-format PairedEndFastqManifestPhred33V2 \
  --output-path cache/demux.qza

Although some qiime commands should run quickly there are four that take a long time and require a lot of memory. This table shows the approximate run times and memory. When you perform each of these steps you MUST do it using a slurm script.

Step Time Memory
Import, primer removal and summary 25 min 1 GB
qiime dada2 denoise-paired 2.5 h 5 GB
qiime feature-classifier classify-sklearn 2 h 9 GB
qiime phylogeny align-to-tree-mafft-fasttree 25 min 2 GB

Things to watch out for

This slurm script classifies your ASVs with 4 cores in about 2 hours;

#!/bin/sh
#SBATCH --time=4:00:00
#SBATCH --ntasks=1 --cpus-per-task=4 --mem=10gb

alias qiime='apptainer run -B /pvol/:/pvol --env "MPLCONFIGDIR=$HOME/matplotlib" /pvol/data/sif/qiime.sif qiime'

qiime feature-classifier classify-sklearn \
  --i-classifier raw_data/oyster/classifier/silva-138-99-341f-805r-nb-classifier.qza \
  --i-reads cache/rep-seqs.qza \
  --p-n-jobs 4 \
  --p-reads-per-batch 250 \
  --o-classification cache/taxonomy.qza

Suggested Analyses

As this is a 16S Metabarcoding dataset you should complete the following basic analyses in qiime2

Visualisation and statistical tests should focus on the question(s);