Skip to content

Home

Jump to bottom
Zoe edited this page Dec 11, 2025 · 6 revisions

Baiting dataset workflow

This dataset consisted of 60 soil samples collected from around 60 kauri trees in the Waitakere Ranges, Auckland, New Zealand, in 2022.

Soil was collected from four cardinal points around each of 60 kauri trees to depth of 10 cm. Soil was pooled together to create one sample per tree before DNA was extracted. DNA was sequenced using amplicon sequencing (Illumina MiSeq) targeting the V3-V4 region of the 16S rRNA gene and the ITS1 gene region of the ITS gene. DNA was also sequenced using shotgun metagenome sequencing (Illumina NovaSeq 6000 S4). All 60 soil samples additionally underwent LAMP analysis to determine the presence and viability of Phytophthora agathidicida.

For this workflow we are looking at comparing amplicon and shotgun metagenome data to determine the taxonomy and functional potential of the soil microbial community surrounding kauri trees. Additionally, we want to assess the concordance of the two sequencing methods in their patterns for determining variation in the microbial communities between PA-detected and PA-not-detected samples (as determined by LAMP analysis).

Table of contents

Amplicon dataset

For a detailed description of the DADA2 pipeline used to generate ASVs, please refer to this GitHub repository.

After processing reads using the standard DADA2 pipeline to generate ASVs we will look at how we can identify the taxonomic composition and functional potential of this community.

Taxonomic data analysis

After generating the ASVs from the DADA2 pipeline, we can analyse the data using the scripts in the amplicon-data/01_taxonomy_scripts/ directory.

Functional inference

We used PICRUST2 to infer functional potential based on the taxonomic composition of the soil microbial community.

Installing PICRUSt2

Will be installing PICRUSt2 on z6. Initially I tried to install from bioconda, however, one of the default files was missing (as seen here). To fix this it is suggested to install PICRUSt2 from source, described here.

PICRUSt2 analysis

PICRUSt2 requires a few different input files, described here:

  • Sequence file (.fna) = sequence information of ASVs
  • Abundance file (BIOM or .tsv) = sequence abundance information

We can create a BIOM file or edit the output from DADA2 for the abundance information. As I already created a BIOM file for another analysis, I will just use that file. Instructions for creating a BIOM file from a phyloseq object can be found here. Instructions for formatting the DADA2 output to be compatible with PICRUSt2 can be found here.

Filtering low abundance taxa: More information on data pre-processing for PICRUSt2 here.

Running a 10% prevalence filter on the ASVs using 00_asv-prevalence-filtering.R script.

Running PICRUSt2

tmux new -s picrust
tmux attach -t picrust

conda activate picrust2

cd picrust2-2.6.1

scripts/picrust2_pipeline.py -s ../picrust-analysis/asv-16s-prev10.fna -i ../picrust-analysis/16S_prev10_biom.biom -o ../picrust-analysis/picrust-output-1404 --stratified -p 10 --verbose >> picrust_log.txt

KO information

Extract and processing KEGG KO information is discussed in the KEGG KO analysis section below.

Phytophthora identification

BLAST ITS ASVs against custom Phytophthora database

Extract ASV sequences (After completion of DADA2 workflow)

# Extract from phyloseq object
asv_seqs <- DNAStringSet(taxa_names(ps))  # Extracts DNA strings
names(asv_seqs) <- paste0("ASV", seq_along(asv_seqs))

# Save to FASTA
writeXStringSet(asv_seqs, filepath = "asvs.fasta")

Build reference database (on NeSI)

Download sequences from NCBI - ITS sequences from Phytophthra agathidicida species

makeblastdb -in phytophthora_refs.fasta -dbtype nucl -out phytophthora_db

Run BLASTn for higher sensitivity

blastn  -task blastn \
        -query its-asvs.fasta \
        -db phytophthora_refs.fasta \
        -strand both \
        -evalue 1e-10 \
        -word_size 24 \
        -perc_identity 97 \
        -qcov_hsp_perc 90 \
        -dust no \
        -outfmt "6 qseqid sseqid pident length mismatch gaps qstart qend sstart send evalue bitscore" \
        -max_target_seqs 25 > asv_vs_PA_blastn.tsv

Primer alignment analysis

BLAST analysis

Primer sequences:

>ITS1F
CTTGGTCATTTAGAGGAAGTAA
>ITS2
GCTGCGTTCTTCATCGATGC

The primer binding sites are in the 18S and 5.8S regions.

  • ITS1-F: 18S binding (SSU rRNA)
  • ITS2: 5.8S binding (LSU rRNA)

  1. Check forward primer binding to the SSU (18S rRNA) region

    • Download the PA SSU sequence from SILVA (accession number LGTR01001481).
    • Run an NCBI pairwise BLAST using the forward primer sequence.
    • The BLAST result shows a short 13 bp alignment between the forward primer and the PA SSU sequence, located near the end of the sequence.
    • The reverse primer does not bind to this SSU sequence.

  2. Check reverse primer binding to the LSU region

    • Download LSU rRNA sequence from UNITE database.
      • UNITE does not contain specific entries for PA but it does include entries from other Phytophthora species, e.g., P. cinnamomi.
      • One relevant entry links to a GenBank accession (L76535.1).

  3. Locate LSU in the PA genome

    • Use the P. cinnamomi LSU sequence as a query and BLAST it against the PA whole genome in NCBI.
    • Results of this search provide coordinates for the LSU region in the PA genome.
    • Extend this genomic region slightly to capture the adjacent SSU sequence as well.

  4. Final primer alignment test

    • Extract the extended LSU–SSU region from the PA genome.
    • Run pairwise BLASTs of both the forward and reverse primers against this sequence to determine their binding potential.
Run BLASTn-short on z6

To check the primer sequences that do align with the PA genome are in the correct place, I ran a pairwise BLAST with 4 sequences:

  1. ITS1-F primer
  2. ITS2 primer
  3. PA-18S sequence from SILVA (LGTR01001481)
  4. PA-ITS sequence from NCBI (JX122749.1)
blastn -task blastn-short \
       -query its-primers.fasta \
       -subject phyt-wg-CP106980.fasta \
       -strand both \
       -evalue 1000 \
       -word_size 7 \
       -outfmt "6 qseqid sseqid pident length mismatch gaps qstart qend sstart send evalue bitscore qseq sseq" \
       > primer_vs_PA_WG.tsv
  • blastn-short = optimised for short sequences like primers.
  • -strand both = searches both forward and reverse-complement orientations.
  • -evalue 1000 = relaxed for short sequences so nothing is filtered out.
  • -word_size 7 = smaller word size improves sensitivity for short sequences

Shotgun metagenome dataset

Firstly, we will run a series of QC steps to clean and filter the reads before assembly and gene prediction.

  1. Run FastQC and MultiQC on reads
    • 00_00a-meta-fastqc-slurm.slurm
    • Transfer outputs to central directory before running multiqc
  2. Trim adapters and PhiX reads
    • 00_01a-00_01a-trim-adapters-slurm.slurm and 00_01b-trim-adapters-workflow.sh
    • 00_02a-trim-phix-slurm.slurm and 00_02b-trim-phix-workflow.sh
  3. Filter out human associated DNA
    • 00_03a-human-filt-bbmap-slurm.slurm and 00_03b-human-filt-bbmap-workflow.sh

Taxonomic composition

To determine the taxonomic composition of the reads in this dataset we can use a few different tools/methods, described here. We will use Kraken2 and Bracken in this workflow.

In addition to running the Kraken2 analysis, we will also use Bracken for better abundance estimates.

The standard kraken2 bacterial and fungal databases have already been downloaded and built on the NeSI HPC. We now need to create the associated Bracken databases. If the database generation fails, you need to remove the files it created before resubmitting the script, if you don't the script will fail with a segmentation fault error - described here.

Build bracken script - 01_00-bracken-build.slurm.

Run kraken2 using scripts 01_01a-kraken-classify-slurm.slurm and 01_01b-kraken-classify-workflow.sh.

Running bracken takes seconds to run and no memory so we can run as a shell script - 01_02a-bracken-report.sh and 01_02b-bracken-workflow.sh.

Now that we have generate the bracken outputs for each sample, we can use the combine_bracken_outputs.py script to combine reports across samples.

# Download script
wget https://raw.githubusercontent.com/jenniferlu717/Bracken/refs/heads/master/analysis_scripts/combine_bracken_outputs.py

# Make script executable 
chmod +x combine_bracken_outputs.py

# Move bracken outputs (not reports) to a central directory
mkdir bracken-output

# Create list of directories
find -name "*-bacteria-bracken" > dir-list.txt
# Sanity check for correct number of directories
wc -l dir-list.txt
# For loop to move outputs
for dirnames in $(cat dir-list.txt)
do
  cd $dirnames || continue
  cp *.bracken ../../../filter-output/bracken-output/bacteria/
  cd "$OLDPWD"
done

# Load python module
module purge
module load Python/3.11.6-foss-2023a

# Run combine script
python combine_bracken_outputs.py --files *.bracken --output all-bracken.tsv

With the combined bracken reports, we can add lineage information to the taxIDs using a script from bit.

Firstly, we need to create a list of the taxIDs identified by the kraken2 analysis, we can then run the bit-get-lineage-from-taxids script to obtain the linage information.

# Extract column of taxIDs from bracken analysis
cut -f2 all-fungi-bracken.tsv > fungi-taxid.txt

# Remove the header row
sed '1d' -i fungi-taxid.txt

# Sanity check that number of rows are the same
wc -l all-fungi-bracken.tsv 
#93 all-fungi-bracken.tsv (incl header row)

wc -l fungi-taxid.txt 
#92 fungi-taxid.txt

# Get lineage information
bit-get-lineage-from-taxids -i fungi-taxid.txt -o fungi-lineages.tsv

# Sanity check that number of rows are the same
wc -l fungi-lineages.tsv 
#93 fungi-lineages.tsv (incl header row)

# Combine tables together
paste all-fungi-bracken.tsv fungi-lineages.tsv > fungi-bracken.tsv

This output can now be analysed with the scripts within the folder metagenome-data/01_taxonomy-scripts/R-scripts/.

Functional annotation

For functional annotation, we will first filter out non-target DNA such as plant and animal DNA that may be present in our samples.

Using kraken2 we can filter out the DNA we are not interested in against a set of both standard and custom kraken2 databases. These databases have all be built on the NeSI HPC and encompass the following taxonomic groups:

  • Plant (standard database)
  • Annelida
  • Arthropoda
  • Chordata
  • Mollusca
  • Nematoda
  • Platyhelminthes
  • Tardigrada
  • Oomycota

Filtering of DNA can be done using the scripts: 02_01a-kraken-exclude-reads-slurm.slurm and 02_01b-kraken-exclude-reads-workflow.sh.

We can then start functional gene prediction and annotation.

  1. Assembly of reads in contigs using MEGAHIT

    • 02_02a-megahit-slurm.slurm and 02_02b-megahit-workflow.sh

    Generate some stats on the assemblies

    module purge
module load BBMap/39.01-GCC-11.3.0

for sample in LIC034-*/*megahit-assembly/final.contigs.fa
do
statswrapper.sh in=$sample, format=6 >> assembly-stats.txt
done

# Tidy up output file to remove extra header rows
sed -n '3~2!p' assembly-stats-adj-kmer-settings.txt > assembly-stats-adj-kmer-settings-tidy.txt

  2. Gene prediction using prodigal

    • 02_03a-prodigal-assembly-slurm.slurm and 02_03b-prodigal-assembly-workflow.sh.
    • After running gene prediction we need to prepare the files for creating a non-redundant gene catalogue.

      • Firstly, concatenate all predicted sequences into one master .fna file.

        # Create list of directories
find . -type d -name "LIC*-prodigal" > dir-list.txt

# Concatenate all .fna files
for dirname in $(cat dir-list.txt)
do
cd $dirname || continue
find *.fna | xargs cat >> ../../all-cds.fna
cd "$OLDPWD"
done
        • Next, strip the metadata from the prodigal output for easier handling.
        cut -f1 -d ' ' all-cds.fna > all-cds-no-meta.fna
        • Adjust ID names using seqkit as there are likely entries with the same name but different sequences (due to the same contig names being used across the different assemblies).
        module purge
module load SeqKit/2.4.0

seqkit rename all-cds-no-meta.fna -o all-cds-no-meta-rename.fna

  3. Create a non-redundant gene catalogue

    • 02_04-cd-hit-est.slurm

  4. Functional annotation using eggNOG database

    • 02_05-eggnog-annotation.slurm

  5. Generate coverage information

    • 02_06a-cds-map-index.slurm
    • 02_06b-cds-map-reads-slurm.slurm and 02_06c-cds-map-reads-workflow.sh
    • 02_06d-coverm-cds.slurm

  6. Normalise coverage information

    • 02_07a-norm-gene-length.py and 02_07b-norm-gene-length.slurm
    • Before running the scripts for this section you need to create the seq-lengths.txt file
    module purge
module load SeqKit/2.4.0

seqkit fx2tab --length --name non-redundant-cds.fna > seq-lengths.txt

# Optional to convert these base pair counts to kilbase counts
# If you do this you have to carry through the rest of your analysis using kilobases
# Mathematically using either is fine, its a matter of scaling that will be adjusted
#awk '{printf "%s\t%.3f\n", $1, $2 / 1000}' seq-lengths.txt > gen-lengths-kb.txt

# Sanity check that all genes were processed
wc -l seq-lengths.txt
# 8519250 

#wc -l non-redundant-cds.fna
grep ">" -c non-redundant-cds.fna
# 8519250

  7. Merge coverage and annotation tables

    • 02_08a-merge-tables.py and 02_08b-merge-tables.slurm

  8. Filter out genes that don't have a KO number

    • 02_09a-filter-kegg.py and 02_09b-filter-kegg.slurm

  9. Sum genes by KO number

    • 02_10a-sum-abund-by-ko.py and 02_10b-sum-abund-by-ko.slurm

KEGG KO analysis

To compare functional annotations between PICRUSt2 outputs and eggNOG annotations we can compare the identified KOs by each method. We want to group KOs based on their BRITE hierarchy information. To do this we will pull down information from the list of KOs using kegg_pull, parse the .txt flat files to extract the BRITE text block, then process this parsed output in R to generate a dataframe that has the KOs and their BRITE hierarchy levels split. We can then combine this dataframe with abundance information, normalise using TSS and plot some barcharts.

Firstly, extract a list of KOs identified by eggNOG and PICRUSt2 (two separate lists).

# KO abundance info per sample
shot_dat <- read_tsv("../metagenome/summed_by_ko.tsv") %>% 
  mutate(KEGG_ko = str_remove(KEGG_ko, "^ko:")) %>% 
  data.frame() %>% 
  column_to_rownames(var = "KEGG_ko")

# List of KOs
kos <- rownames(shot_dat)
# Write to file with no quotes or row numbers
writeLines(kos, "kegg-pull/ko_ids.txt")

Run kegg_pull on z6

conda activate kegg_env

cat ko_ids.txt | kegg_pull pull entry-ids - --multi-process --output=ko_entries --sleep-time=8.0
# 1:00:12 to complete

Unfortunately, this does not pull down all our KOs in one go and a large portion tend to fail each time. We can extract the names of the failed searches using the 02_11a-extract-failed-kos.py script. The output of this script then becomes the input for kegg_pull. To keep track of the failed and retried KO searches I would rename the pull-results.json and failed_ids.txt files each time I resubmitted. I also move the pull KOs in the ko_entries/ directory into a new central directory so they don't get overwritten.

Some KOs will always fail and this is likely due to depreciated KOs that are no longer present on the KEGG website but are still in the eggNOG database that uses an older version of KEGG. Once the failed list is small enough you can manually check the KOs on the KEGG website to see if they are depreciated.

# Some quick code to find missing KOs (likely depreciated)
ls ko_entries/*.txt | xargs -n1 basename > found_files.txt
sed 's/\.txt$//' found_files.txt > found_files_stripped.txt

comm -23 <(sort ko_ids.txt) <(sort found_files_stripped.txt) > missing_kos.txt

Once all the KO info files have been pulled from KEGG, we can combine the BRITE info into one file using the 02_11b-extract-brite-batch.py script. Copy the output of this script back to a local computer for further analysis with R.

We can extract the mapping information using the kegg-parse-info.R script (in the misc-scripts directory). The output of this script is essentially a mapping file that you can use to link annotated KOs to their function.

This mapping file can now be used with the coverage/abundance information to analyse this data using the scripts in metagenome-data/02_function-scripts/R-scripts/.

Comparative analysis

After generating all of the outputs from the code described above we can run the scripts within the R-combined-scripts directory to combine the plots from the amplicon and shotgun datasets to compare outputs.

Clone this wiki locally