Skip to main content

Marine Genomics

Workshop 4: The Metagenome of Black Band Disease

Background

JCU PhD student Julia Hung studied metabolic processes within the black band disease (BBD) community. In 2023 she collected fresh samples of BBD from infected individuals of the coral, Acropora hyacinthus at Orpheus Island. She successfully extracted high molecular weight DNA from these samples and sequenced them using an oxford nanopore MinION sequencer.

This data formed part of a paper that she published in the journal, Environmental Microbiology

Microbial Metabolism and Disease Virulence Changes Across Day and Night in Coral Black Band Disease Lesions

You have been provided with a subset of Julia’s data which should be sufficient to generate a useful metagenome assembly. We will not be working with the full dataset because the computing resources required to do so are too large for a class setting.

In this workshop we will perform an analysis of the black band disease metagenome. At the end of this workshop you should have assembled reads from nanopore sequencing of black band disease, classified the assembled contigs into bins using metabat2, and then performed some functional analysis on the genomes to see what genes are prominent in the bacterial community of this disease.

Workshop Timeline

  • Do before the workshop: Steps 1-4. These take several hours to run. If you haven’t done them by the final workshop start directly at step 5
  • During the final workshop: Steps 5 onward

Setup

This final workshop is the one you should document on github for assessment. It has its own template repository, bbd-metagenome, which already contains an RMarkdown file and a folder for each step of the analysis. Just as you did with mb5370-workshops in workshop 1, first make your own repository from the template on github and then clone that into a new project using the git project type. The process is documented here

In these instructions we will assume that your repository is called bbd-metagenome and you cloned it into your home directory. After cloning the repository test this by using pwd as follows

pwd

It should print something like /home/jcXXXX/bbd-metagenome where jcXXXX is your student login id.

Step 1: Document your work

Your repository has a numbered RMarkdown document for each main logical step of the analysis, and a folder where you run the commands for that step. Every folder you need is already there, so you won’t need to make any.

Steps in this guide RMarkdown Folder
Step 2: Assembly 01_assembly.Rmd flye
Step 3: Binning 02_binning.Rmd metabat
Steps 4 and 5: Bin quality 03_bin_quality.Rmd checkm
Step 6: Taxonomy 04_taxonomy.Rmd gtdbtk
Step 7: Gene annotation 05_gene_annotation.Rmd prokka
Step 8: Metabolic pathways 06_metabolic_pathways.Rmd minpath

Take notes on each step in its RMarkdown document. Each one is set up to knit to a markdown (.md) file that github displays nicely, and your README.md already links to each of these. Edit the README.md so that it introduces your work. I do this for every paper I publish. If you want to see an example here is one from a paper on coral genetics

You will need to add these documentation files, and the files created when you knit them, to git and push them to github.

Your repository also has an index.Rmd file at the top level. You don’t write any notes in it. Instead it joins all six step files together, and when you knit it you get a single web page, index.html, with each step as a section and a table of contents down the side. Knit it whenever you want to read through your whole analysis in one place, for example before you write it up. The step files are knitted again as part of index.Rmd, so if it fails to knit, knit each step file on its own to find the one with the error.

Step 2: Assembly.

This step takes a long time to run. You will need to do it before the final workshop on Thursday. If you didn’t have time to do it that’s OK. You can use some pre-made results (see instructions below)

Take notes on this step in 01_assembly.Rmd.

Assuming you have time to perform this assembly use these instructions to start the job.

Change directory into flye and then make a symbolic link to the bbd read data

cd ~/bbd-metagenome/flye
ln -s /pvol/data/metagenome_assembly/bbd_0.1.fastq.gz .

Now make a new shell script called run_flye_bbd.sh and save it inside your flye directory. Edit it so it looks like this

#!/bin/bash
#SBATCH --time=720
#SBATCH --ntasks=12 --mem=20gb

echo "Starting flye bbd in $(pwd) at $(date)"

flye --nano-hq bbd_0.1.fastq.gz --threads 12 --out-dir bbd --meta

echo "Finished flye bbd in $(pwd) at $(date)"

And now submit your job using

sbatch run_flye_bbd.sh

What to add to git?

This is a good moment to think about what you want to add to git and then push to github. Obviously you will already be adding your RMarkdown docs. You might also want to add shell scripts such as the one you just created. There will soon be many large files within the flye directory. Your repository’s .gitignore file hides these from git, so the only file in flye that the git pane should offer you is your shell script. Even so, always check what you are adding before you commit.

Monitoring job progress

Since lots of people will be trying to run jobs on the server it is likely that yours will wait in the queue until it starts. You can check the queue with the squeue command

squeue

Take a look for your username. When your job is running the ST column will show an R (running). You will also see a new slurm-XXX.out file appear in your flye directory.

This time flye will take quite a while to run. Depending on how many people are running it at the same time it could take many hours. The recommended workflow is to start your run, monitor it for a while to make sure it is working and then allow it to run overnight.

Check that flye worked

When your job has finished, check that it worked before you move on. Don’t rely on the last line of the slurm log. It will say Finished flye bbd even if flye failed, because the echo at the end of the script runs either way. Instead use sacct with your job number (the number in the name of your slurm-XXX.out file)

sacct -j XXX --format=JobID,State,MaxRSS,Elapsed

The State column should say COMPLETED. If it says something else, such as OUT_OF_MEMORY or TIMEOUT, ask an instructor for help or use the pre-made results (see below). You should also have a finished assembly

ls -l bbd/assembly.fasta

Step 3: Binning with metabat

After your flye assembly has run you should perform binning with metabat. Take notes on this step in 02_binning.Rmd.

First cd to the metabat folder in your project.

cd ~/bbd-metagenome/metabat

Extract depth information

cat ../flye/bbd/assembly_info.txt  | grep -v '#' | awk '{OFS="\t";print $1,$3}' > depths_bbd.txt

Reorder contigs based on depth

cat depths_bbd.txt | awk '{print $1}' | xargs samtools faidx ../flye/bbd/assembly.fasta > assembly_bbd_ordered.fasta

Run metabat

metabat2 -i assembly_bbd_ordered.fasta -a depths_bbd.txt --cvExt -o bins_bbd/metabat --seed 5

Step 4: Check bin quality with checkm

Now run checkm on the bins created by metabat. Take notes on this step and on step 5 in 03_bin_quality.Rmd.

As before cd to the folder for this step, checkm.

cd ~/bbd-metagenome/checkm

You will need to find the number of your largest contig in the bbd assembly. We need to manually add this contig because it usually isn’t included in the bins from metabat.

This command will display the lengths of all contigs sorted from shortest to longest.

cat ../flye/bbd/assembly.fasta | bioawk -c fastx '{print length($seq),$name}' | sort -n

Once you have found the longest contig substitute its name instead of contig_XX in the code below. The first command saves your longest contig as a bin of its own, in the same folder as the metabat bins, so that all the later steps treat it just like the other bins. The second command copies all of the bins, including your longest contig, into the bins_bbd folder inside checkm.

samtools faidx ../flye/bbd/assembly.fasta contig_XX > ../metabat/bins_bbd/contig_XX.fa
cp ../metabat/bins_bbd/*.fa bins_bbd/

Next create a slurm script for your checkm run. Your script should look like this and you should save it within your checkm directory as run_checkm_bbd.sh

#!/bin/bash
#SBATCH --time=60
#SBATCH --ntasks=6 --mem=40gb

echo "Starting checkm in $(pwd) at $(date)"

shopt -s expand_aliases
alias checkm='apptainer run -B /pvol/:/pvol /pvol/data/sif/checkm.sif checkm'

checkm lineage_wf bins_bbd out -t 6 -x fa -f checkm_results_bbd.txt

echo "Finished checkm in $(pwd) at $(date)"

Submit your checkm job like this

sbatch run_checkm_bbd.sh

Functional Analysis of the BBD Metagenome

In theory each of the bins from our metagenomic assembly corresponds to a species or strain forming part of the black band disease community. These genomes are valuable resources that can be used to enhance analysis of other data such as RNA sequencing data. In addition (and relevant here), they can be used to make deductions about the metabolic processes occurring within a black band disease lesion. The full process of making these deductions involves quite a bit of biochemical knowledge, literature searching and cross referencing of outputs from many statistical tools. Here we will focus on some preliminary analyses that reveal basic processes within the mat and that provide information on the genes and pathways present.

The overall goal of functional analysis of metagenomes is to answer the following questions

  1. Who is there, or what taxonomic groups of microbes are present?
  2. What can they do? (what genes do they have and what biochemical pathways do these support)

Our starting point for functional analysis is a bbd assembly on which we have already performed binning with MetaBat and a quality check with checkm (steps 1-4 above). Files that you will need from these analysis and their expected locations are;

Analysis Step File or folder
Flye flye/bbd
MetaBAT metabat/bins_bbd
CheckM checkm/checkm_results_bbd.txt

Note: Expected locations are given as relative paths within your bbd-metagenome project.

If any of these steps did not work, or you didn’t get time to complete them, don’t worry. You can retrieve pre-analysed results from a shared folder

Step 5: Select bins based on checkm results

Change directory to the checkm folder. Keep taking notes in 03_bin_quality.Rmd.

cd ~/bbd-metagenome/checkm

Now take a look at the results from our checkm analysis. These are summarised in the file checkm_results_bbd.txt

You can view the results in the terminal using cat, head, more or less commands. These all represent different ways to view a text file.

For example

head checkm_results_bbd.txt

If your screen is wide enough this should display the top 10 lines in the file in a nicely formatted way. If your screen is too narrow the lines will wrap and it will be hard to read. Take a look at the Marker lineage column for the top few bins. This column isn’t a full taxonomic classification but it does seem to give family, class or phylum level classification for the top most complete bins. As you go further down the list you will eventually see the marker lineage shift to a very broad level such as k_Bacteria. Try viewing a few more results by telling the head command to show more lines

head -n 20 checkm_results_bbd.txt

In the analyses that follow we would like to focus on the most complete bins because these are likely to contain the majority of useful information on function within the bbd community. These bins will also usually correspond to the most abundant (and therefore best assembled) members of the community. Here we are just using the checkm completeness to assess these factors but a more complete analysis might also look at the size of contigs within the bins and their read coverage.

Let’s focus on bins with a completeness greater than 15%. To extract the bin IDs for these we can use awk. This command should print only the lines where the 13th column has a value greater than 15.

cat checkm_results_bbd.txt | awk '$13>15'

This seems to be the bins we want but we just want the bin ID not the entire line. To get this we can use the print command

cat checkm_results_bbd.txt | awk '$13>15{print $1}'

Now let’s send this output to a file so that we can use it in the next step

cat checkm_results_bbd.txt | awk '$13>15{print $1}' > bbd_topbins.txt

Step 6: Perform taxonomic classification with gtdbtk

Take notes on this step in 04_taxonomy.Rmd. First change directory to the gtdbtk folder.

cd ~/bbd-metagenome/gtdbtk

Inside the gtdbtk folder is a subfolder called genomes. This is where we will put a copy of each of the top bins we identified in step 5.

This next bit of code is a bit more advanced. Don’t worry if you don’t understand it all. The idea here is that we will iterate over our list of topbins. For each of the topbins we want to copy a file corresponding to that bin into our genomes folder.

We’ll build up the command in steps. First just get the while loop working

while read bin;do echo $bin;done < ../checkm/bbd_topbins.txt

We’re not doing any copying of files yet. Instead, we have used echo $bin to simply print the name of the bin. We need to turn this bin name into a path specifying the fasta file for the bin. So for every bin we need to prepend the path to its folder ../metabat/bins_bbd/ and then append the file extension .fa

while read bin;do echo "../metabat/bins_bbd/${bin}.fa";done < ../checkm/bbd_topbins.txt

Now we are ready to perform the copy

while read bin;do cp "../metabat/bins_bbd/${bin}.fa" genomes;done < ../checkm/bbd_topbins.txt

Check that it worked by listing files within the genomes directory. You should see one file for each of the bins

ls genomes

These bins form input for gtdbtk. We can now run the classification workflow. This is another very resource intensive step. We will need to run it by creating a slurm job script. Hopefully the workflow for this is familiar by now.

First create a blank shell script file and save it inside your gtdbtk directory. Name your file run_gtdbtk.sh

Copy the following code into your run_gtdbtk.sh file and save it.

#!/bin/bash
#SBATCH --time=120
#SBATCH --ntasks=2 --mem=56gb

echo "Starting gtdbtk in $(pwd) at $(date)"

shopt -s expand_aliases
alias gtdbtk='apptainer run -B /pvol:/pvol,/pvol/data/mg/gtdbtkdata/release207_v2/:/refdata /pvol/data/sif/gtdbtk.sif'

gtdbtk classify_wf -x fa --cpus 2 --genome_dir genomes --out_dir out --skip_ani_screen

echo "Finished gtdbtk in $(pwd) at $(date)"

Then submit your job using

sbatch run_gtdbtk.sh

This will take quite a while. While you are waiting you can continue with step 7.

When your job has finished, check that it produced a summary of its classifications

ls -l out/gtdbtk.bac120.summary.tsv

If this file doesn’t exist the job probably ran out of memory. Use sacct to check, as you did for flye in step 2.

Now copy the summary up into the gtdbtk folder. Git ignores everything inside out, but it can see a .tsv file saved directly in gtdbtk, so this lets you add the summary to your repository.

cp out/gtdbtk.bac120.summary.tsv .

If your gtdbtk job didn’t work and you are using the pre-made results from steps 2-4, copy the pre-made summary instead. This only works with the pre-made results, because the bins from your own assembly will have different numbers.

cp /pvol/data/mg/bbd/gtdbtk/out/gtdbtk.bac120.summary.tsv .

Step 7. Find and annotate genes with prokka

Our ultimate goal is to use the BBD metagenome to improve our understanding of metabolic processes within the mat. A crucial step toward this is identifying the genes present within the genomes of bbd taxa. We will be using a program called prokka for this purpose. Prokka will not only identify the gene sequences, it also uses a range of methods to identify homologous genes with known functions. The idea behind this is that if two genes share a common ancestor (ie are homologous) their functions are likely to be similar. This in turn allows Prokka to assign a putative name and function to some of the genes within our newly assembled metagenome bins.

Take notes on this step in 05_gene_annotation.Rmd. As usual we will change to the folder for our prokka work.

cd ~/bbd-metagenome/prokka

First we try running prokka for a single bin. This produces quite a bit of output. Watch the output as it goes by. This output gives a good idea of the various steps prokka goes through.

Next pick a bin to test with. Look at the list of bins you made in step 5

cat ../checkm/bbd_topbins.txt

These instructions use metabat.5. If you are using the pre-made results it will be in your list and you can use it too. If you are using your own assembly your bins will have different numbers, so pick any bin from your list and substitute its id where appropriate in the commands below.

bin=metabat.5
prokka --cpus 2 --outdir $bin --prefix $bin  ../metabat/bins_bbd/${bin}.fa

The --cpus 2 option stops prokka from using more than two CPUs. This command runs on the RStudio server itself, which everyone in the class shares.

This should run fairly quickly (about 3 minutes), and if it works you should end up with a folder called metabat.5. List the contents of this folder. It should look like this

ls metabat.5
metabat.5.err  metabat.5.faa  metabat.5.ffn  metabat.5.fna  metabat.5.fsa  metabat.5.gbk  metabat.5.gff  metabat.5.log  metabat.5.sqn  metabat.5.tbl  metabat.5.tsv  metabat.5.txt

There are many output files here. For our purposes the most important one is metabat.5.tsv which is a tab separated file with the name of each gene and its inferred name and functional info (if available)

Now we are ready to run prokka on all the metagenomes. This is another computationally intensive step, so it’s time for another slurm script.

Create a new blank text file and save it in your prokka directory. Call it run_prokka.sh and copy the following code into it.

#!/bin/bash
#SBATCH --time=90
#SBATCH --ntasks=2 --mem=2gb

echo "Starting prokka in $(pwd) at $(date)"

shopt -s expand_aliases
alias prokka='apptainer run -B /pvol/:/pvol /pvol/data/sif/prokka.sif prokka'

while read bin;do 
  prokka --cpus 2 --outdir $bin --prefix $bin "../metabat/bins_bbd/${bin}.fa";
done < ../checkm/bbd_topbins.txt

echo "Finished prokka in $(pwd) at $(date)"

Submit your job with

sbatch run_prokka.sh

Your slurm log will include a message like Folder 'metabat.5' already exists! Please change --outdir or use --force. This is expected. It just means prokka skipped the bin you already annotated in your test run.

While this is running you can keep working on step 8.

Step 8: Infer metabolic pathways

Infering metabolic pathways can be quite complex and is difficult to do in a completely objective or automated fashion. In our case we are lucky because we can refer to the general model of BBD developed by Sato et al.

img
img

According to this figure we expect to find that the cyanobacteria should be the main organism capable of photosynthesis. We also expect to find some organisms capable of sulfate reduction. Nominally this should be species of Desulfovibrio, however, there may be other organisms with this capability and since our assembled bins for Desulfovibrio are quite incomplete we might not be able to find sulfate-reducing enzymes in this group.

We will use a tool called minpath to take functional info for individual genes and use this to infer the presence of broader biological pathways (involving multiple genes/steps). We will run this on each of the genomes (bins) in our analysis. Minpath will look at all of the annotated genes in our data and infer the minimal set of biological pathways that can explain the presence of all these genes. It is important to note that this is just one approach (parsimony) to inferring pathways and other approaches also exist.

Take notes on this step in 06_metabolic_pathways.Rmd. First change into the minpath folder, where we will store minpath results.

cd ~/bbd-metagenome/minpath

Inputs for minpath will be prepared for each bin based on the prokka .tsv file. Start by looking at the .tsv file for a single bin

head ../prokka/metabat.5/metabat.5.tsv

You should see 9 genes displayed. Many will be missing annotations (called hypothetical protein) but some will have a gene name (eg manA), an EC (Enzyme Commission) number (eg 5.3.1.8), a COG (Clusters of Orthologous Genes) id (eg COG1482) and a protein name (eg Mannose-6-phosphate isomerase ManA).

For input to minpath we only want to keep a unique ID for each gene (column 1) and its EC number (column 5). We can transform the tsv into this format using awk like this

head ../prokka/metabat.5/metabat.5.tsv | awk -F '\t' '{print $1,$5}'

Note that here we are using -F '\t' to tell awk that our input is tab separated. We also want to filter the output to only include lines with an EC number

head ../prokka/metabat.5/metabat.5.tsv | awk -F '\t' '$5{print $1,$5}'

Next we need to ensure that the output is tab separated and we also need to remove the header line

head ../prokka/metabat.5/metabat.5.tsv | awk -F '\t' '$5{OFS="\t";print $1,$5}' | grep -v '^locus'

This of course only processes the first 10 lines … we need to change head to cat to include all genes

cat ../prokka/metabat.5/metabat.5.tsv | awk -F '\t' '$5{OFS="\t";print $1,$5}' | grep -v '^locus'

Finally we send the output to a file so that it can be used as input to minpath

cat ../prokka/metabat.5/metabat.5.tsv | awk -F '\t' '$5{OFS="\t";print $1,$5}' | grep -v '^locus' > metabat.5.ec

Now run minpath on this file as follows;

minpath -ec metabat.5.ec -report metabat.5.ec.report -details metabat.5.ec.details

Take a look at the top few lines of output

head metabat.5.ec.report

As an example consider the first line

path PWY-7321 any n/a  naive 1  minpath 1  fam0  4  fam-found  2  name  ecdysteroid metabolism (arthropods)

Important parts of this line are

  • PWY-7321 which is a pathway identifier in the metacyc database. You can use this to lookup the full pathway. Eg PWY-7321
  • naive 1 tells us that at least one gene from this pathway was found. Simply calling a pathway present whenever one of its genes is found is the “naive” approach
  • minpath 1 tells us that minpath kept this pathway in its minimal set of pathways. If this says minpath 0 then minpath decided that the genes found are better explained by other pathways, so the evidence for this one is weak
  • fam0 4 tells us that four gene families are involved in this pathway
  • fam-found 2 tells us that we found 2 gene families from the pathway in this genome
  • name ecdysteroid metabolism (arthropods) is a name for the pathway

Now go back to check on your prokka job. You can do this by using the command squeue. Look for a job named run_prok and your username. When your job is done it will not appear in this list. You can also look at the output log file

tail ../prokka/slurm-XXX.out

When your job is done it should show Finished prokka on the last line of output.

When prokka is done you can run the following code to prepare minpath inputs for all bins.

while read bin;do 
  cat ../prokka/$bin/${bin}.tsv | awk -F '\t' '$5{printf("%s\t%s\n", $1,$5)}' | grep -v '^locus' > ${bin}.ec 
done < ../checkm/bbd_topbins.txt

Now run minpath on every bin. This should run pretty fast.

while read bin;do 
  minpath -ec ${bin}.ec -report ${bin}.ec.report -details ${bin}.ec.details
done < ../checkm/bbd_topbins.txt

There is a huge amount of information in the outputs from minpath. To keep things tractable we are just going to look at two pathways.

The first pathway is named photosynthesis light reactions and has the metacyc code PWY-101. This pathway includes the essential components of photosynthesis that occur in light, Photosystem I and Photosystem II.

Use the following command to search for this pathway across all the bins

grep 'PWY-101' *.ec.report
Question 4.1
  • How many bins contain at least one gene from this pathway?
  • Which bin has all of the genes (ie a complete pathway)?
  • For which bins does minpath keep the pathway (minpath 1)?

Cross reference your results with the output from gtdbtk to see the taxonomy of bins containing this pathway. Note that you will need to modify the commands below to use the correct bin numbers (these will be different for everyone).

grep 'metabat.26' ../gtdbtk/gtdbtk.bac120.summary.tsv

The second pathway is a little less obviously named. It is called superpathway of tetrathionate reduction and has metacyc code PWY-5360. This pathway is important because it includes enzymes required to reduce sulfite to hydrogen sulfide. Note that sulfite acts as the terminal electron acceptor in this pathway, however, in a broader sense the term “sulfate reducing” refers to micro-organisms that can produce sulfide from a range of oxidised sulfur compounds.

Question 4.2
  • How many bins contain the pathway PWY-5360?
  • What order (in a taxonomic sense) do these bins belong to? Are they in the same family?
  • How strong is the evidence that these bins really have this pathway? Look at the fam-found and minpath columns.