Marine Genomics Workshop 4

Workshop 4 - Functional analysis of the BBD metagenome

A flye assembly, binning and bin assessment were conducted on the BBD metagenome. The overall goal of the functional metagenome analysis is to find what taxonomic groups of microbes are present and their function.

Analyses

The flye assembly

Symbolic link to the bbd read data

ln -s /pvol/data/metagenome_assembly/bbd_0.1.fastq.gz .

Shell script to run flye

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

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)"

Binning with metabat

Extract depth information

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

Re-order 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

Bin assessment

Find the longest contigs and substitute it in ‘contig_xx’

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

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

Create a slurm script and run it.

#!/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 4 -x fa -f checkm_results_bbd.txt

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

Functional analysis of the BBD Metagenome

The analyses that were conducted were focused on the most complete bins, as these contain the majority of useful information within the BBD community, these will also be the most abundant and best assembled members of the community. We therefore will focus on bins with a completeness greater than 15% and perform taxonomic classification with gtdbtk on these.

Extract desired bins

The following code will extract the bins with a completeness greater than 15%, and send the output to a file text

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

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

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

Taxonomic classification

The following code will iterate over the list of topbins, for each of them we want to copy a file corresponding to each bin into the genomes folder.

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

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

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

Then we can run the classification workflow by creating the following shell script and running it.

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

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)"

Finding and annotating genes

Identify the genes present within the genomes of BBD taxa is essential for improving our understanding of metabolic processes. The following code will identify gene sequences, and identifies homologous genes with known functions, as genes with common ancestors will likely have similar functions.

Substitute with the first bin ID in your list

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

After this runs there will be a file with the name of the gene and its inferred name and functional information

Next we will run the following slurm script to run prokka on all the metagenomes.

#!/bin/bash
#SBATCH --time=60
#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 --outdir $bin --prefix $bin "../metabat/bins_bbd/${bin}.fa";
done < ../checkm/bbd_topbins.txt

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

Infer metabolic pathways

Inferring metabolic pathways is a complex and objective process, in this case we can refer to the general model of BBD developed by Sato et al. Which according to the model developed by Sato, we should expect cyanobacteria being the main organism capable of photosynthesis, and some other organisms capable of sulfate reduction.

The following code will take functional information for individual genes and use it to infer the presence of broader biological pathways.

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

Run minipath on the following file

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

Now run the following code to prepare minipath inputs for all the 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 minipath as follows

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

We will focus on two pathways, the first one is the photosynthesis light reactions and the metacyc code (PWY-101), it includes essential components of photosynthesis (Photosystem I and II)

grep 'PWY-101' *.ec.report

Modify the bin number to use the correct one

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

The second pathway is called superpathway of tetrathionate reduction and has metacyc code (PWY-5360), this pathways includes enzymes required to reduce sulfite to hydrogen sulfide.

Question - How many bins contain at least one gene from the PWY-101 pathway? and which bin has a complete pathway?

Four bins (contig_657, metabat.27, metabat.56, and metabat.66) contain at least one gene associated with photosynthesis light reactions (PWY-101). Of these, MinPath flags metabat.27 and metabat.66 as having the pathway confidently present.However, all four bins show only one expected gene family detected, therefore there is only partial evidence for a fully complete pathway.

Question - How many bins contain the pathway PWY-5360? and what family do these bins belong to?

Only one bin (metabat.67) contains at least one gene associated with superpathway of tetrathionate reduction (PWY-5360), however, the Minpath did not confidently confirm the pathway as complete. Just one gene family was detected, belonging to the Arcobacteraceae family.