SRA Data Download and Comparative Genome Analysis

Overview

In this exercise, we will download raw sequencing data from two published genome papers and analyse them using the training workflow.

The two papers are:

  1. Draft genome sequence of Leptospira yasudae strain Ss3a1f from Estero De Paco, Manila, Philippines
    DOI: https://journals.asm.org/doi/10.1128/mra.00547-26

  2. Draft genome sequence of multidrug-resistant Staphylococcus aureus strain MB02 isolated from the cesarean wound of a 32-year-old pregnant woman in Ogbomoso, Nigeria
    DOI: https://journals.asm.org/doi/10.1128/mra.01357-25

We will download the raw SRA data, perform read QC, trimming, genome assembly, assembly QC, taxonomy, MLST, annotation, AMR screening, and virulence factor screening.

Learning objectives

By the end of this exercise, you should be able to:

  • download raw sequencing reads from SRA
  • convert SRA files to FASTQ format
  • count raw and trimmed reads
  • perform read quality control
  • trim raw reads
  • assemble genomes
  • assess assembly quality
  • filter short contigs
  • estimate completeness and contamination
  • classify taxonomy
  • perform MLST
  • annotate genomes
  • screen for AMR genes and virulence factors
  • complete a comparative genome analysis results table

Input data

Paper Organism SRA accession Suggested sample name
Paper 1 Leptospira yasudae strain Ss3a1f SRR38327166 L_yasudae
Paper 2 Staphylococcus aureus strain MB02 SRR32454207 S_aureus
NoteImportant

Use the suggested sample names throughout this exercise:

L_yasudae
S_aureus

This will make the output files easier to compare.

Workflow by organism

Input information

Item Information
Organism Leptospira yasudae strain Ss3a1f
Paper https://journals.asm.org/doi/10.1128/mra.00547-26
SRA accession SRR38327166
Sample name used here L_yasudae

Step 1: Download SRA data

mkdir -p sra_comparative_analysis
cd sra_comparative_analysis

mkdir -p raw_sequence raw_reads qc_reports trimmed_reads assemblies quast checkm2 classify annotation abricate results

conda activate sra-tools

cd raw_sequence
prefetch SRR38327166
fasterq-dump SRR38327166/SRR38327166.sra
gzip *.fastq
cd ..

Step 2: Rename FASTQ files

cp raw_sequence/SRR38327166_1.fastq.gz raw_reads/L_yasudae_R1.fastq.gz
cp raw_sequence/SRR38327166_2.fastq.gz raw_reads/L_yasudae_R2.fastq.gz

Check files:

ls raw_reads

Step 3: Count total raw reads

echo "L_yasudae raw R1 reads:"
echo $(( $(zcat raw_reads/L_yasudae_R1.fastq.gz | wc -l) / 4 ))

echo "L_yasudae raw R2 reads:"
echo $(( $(zcat raw_reads/L_yasudae_R2.fastq.gz | wc -l) / 4 ))

echo "L_yasudae total raw reads:"
echo $(( $(zcat raw_reads/L_yasudae_R1.fastq.gz | wc -l) / 4 + $(zcat raw_reads/L_yasudae_R2.fastq.gz | wc -l) / 4 ))

Record this value in the results table.

Step 4: FastQC on raw reads

conda activate qc

fastqc \
  -t 4 \
  -o qc_reports \
  raw_reads/L_yasudae_R1.fastq.gz \
  raw_reads/L_yasudae_R2.fastq.gz

Step 5: Trim reads

trimmomatic PE \
  -threads 4 \
  raw_reads/L_yasudae_R1.fastq.gz \
  raw_reads/L_yasudae_R2.fastq.gz \
  trimmed_reads/L_yasudae_trim_R1.fastq.gz \
  trimmed_reads/L_yasudae_unpaired_R1.fastq.gz \
  trimmed_reads/L_yasudae_trim_R2.fastq.gz \
  trimmed_reads/L_yasudae_unpaired_R2.fastq.gz \
  SLIDINGWINDOW:4:20 \
  LEADING:3 \
  TRAILING:3 \
  HEADCROP:15 \
  MINLEN:40

Step 6: Count total reads after trimming

echo "L_yasudae trimmed R1 reads:"
echo $(( $(zcat trimmed_reads/L_yasudae_trim_R1.fastq.gz | wc -l) / 4 ))

echo "L_yasudae trimmed R2 reads:"
echo $(( $(zcat trimmed_reads/L_yasudae_trim_R2.fastq.gz | wc -l) / 4 ))

echo "L_yasudae total trimmed reads:"
echo $(( $(zcat trimmed_reads/L_yasudae_trim_R1.fastq.gz | wc -l) / 4 + $(zcat trimmed_reads/L_yasudae_trim_R2.fastq.gz | wc -l) / 4 ))

Step 7: FastQC on trimmed reads

fastqc \
  -t 4 \
  -o qc_reports \
  trimmed_reads/L_yasudae_trim_R1.fastq.gz \
  trimmed_reads/L_yasudae_trim_R2.fastq.gz

conda deactivate

Step 8: Genome assembly

conda activate assembly

spades.py \
  -t 4 \
  -o assemblies/L_yasudae_assembly \
  -1 trimmed_reads/L_yasudae_trim_R1.fastq.gz \
  -2 trimmed_reads/L_yasudae_trim_R2.fastq.gz

cp assemblies/L_yasudae_assembly/contigs.fasta assemblies/L_yasudae_contigs.fasta

Step 9: Count contigs before filtering

grep -c "^>" assemblies/L_yasudae_contigs.fasta

Record this as:

Number of contigs >=0 bp

Step 10: QUAST on original assembly

quast \
  -t 4 \
  -o quast/L_yasudae_original \
  assemblies/L_yasudae_contigs.fasta

cat quast/L_yasudae_original/report.txt

Step 11: CheckM2 on original assembly

checkm2 predict \
  -t 4 \
  -x .fasta \
  -i assemblies/L_yasudae_contigs.fasta \
  -o checkm2/L_yasudae_original

cat checkm2/L_yasudae_original/quality_report.tsv

Step 12: Filter contigs smaller than 500 bp

seqkit seq \
  -m 500 \
  assemblies/L_yasudae_contigs.fasta \
  > assemblies/L_yasudae_contigs_filtered.fasta

Count filtered contigs:

grep -c "^>" assemblies/L_yasudae_contigs_filtered.fasta

Step 13: QUAST on filtered assembly

quast \
  -t 4 \
  -o quast/L_yasudae_filtered \
  assemblies/L_yasudae_contigs_filtered.fasta

cat quast/L_yasudae_filtered/report.txt

Record:

  • N50
  • L50

Step 14: CheckM2 on filtered assembly

checkm2 predict \
  -t 4 \
  -x .fasta \
  -i assemblies/L_yasudae_contigs_filtered.fasta \
  -o checkm2/L_yasudae_filtered

cat checkm2/L_yasudae_filtered/quality_report.tsv

conda deactivate

Record:

  • completeness
  • contamination

Step 15: Taxonomic classification

conda activate classify

kraken2 \
  --db ~/databases/minikrakendb \
  --threads 4 \
  --report classify/L_yasudae.kraken2.report \
  --output classify/L_yasudae.kraken2.output \
  assemblies/L_yasudae_contigs_filtered.fasta

head classify/L_yasudae.kraken2.report

conda deactivate

Record the main taxonomy result.

Step 16: MLST

conda activate mlst

mlst assemblies/L_yasudae_contigs_filtered.fasta > classify/L_yasudae_mlst.txt

cat classify/L_yasudae_mlst.txt

conda deactivate

If no MLST scheme is detected, write:

No MLST scheme detected

Step 17: Genome annotation

conda activate annotation

prokka \
  --outdir annotation/L_yasudae \
  --force \
  --cpus 4 \
  --prefix L_yasudae \
  --locustag LYAS \
  assemblies/L_yasudae_contigs_filtered.fasta

cat annotation/L_yasudae/L_yasudae.txt

conda deactivate

Record:

  • number of CDS
  • number of tRNA

Step 18: AMR gene screening

conda activate amr

abricate \
  --db ncbi \
  --threads 4 \
  assemblies/L_yasudae_contigs_filtered.fasta \
  > abricate/L_yasudae_amr.tsv

abricate --summary abricate/L_yasudae_amr.tsv > abricate/L_yasudae_amr_summary.tsv

cat abricate/L_yasudae_amr_summary.tsv

Record the AMR genes.

Step 19: Virulence factor screening

abricate \
  --db vfdb \
  --threads 4 \
  assemblies/L_yasudae_contigs_filtered.fasta \
  > abricate/L_yasudae_vf.tsv

abricate --summary abricate/L_yasudae_vf.tsv > abricate/L_yasudae_vf_summary.tsv

cat abricate/L_yasudae_vf_summary.tsv

conda deactivate

Record the virulence factors.

Input information

Item Information
Organism Staphylococcus aureus strain MB02
Paper https://journals.asm.org/doi/10.1128/mra.01357-25
SRA accession SRR32454207
Sample name used here S_aureus

Step 1: Download SRA data

mkdir -p sra_comparative_analysis
cd sra_comparative_analysis

mkdir -p raw_sequence raw_reads qc_reports trimmed_reads assemblies quast checkm2 classify annotation abricate results

conda activate sra-tools

cd raw_sequence
prefetch SRR32454207
fasterq-dump SRR32454207/SRR32454207.sra
gzip *.fastq
cd ..

Step 2: Rename FASTQ files

cp raw_sequence/SRR32454207_1.fastq.gz raw_reads/S_aureus_R1.fastq.gz
cp raw_sequence/SRR32454207_2.fastq.gz raw_reads/S_aureus_R2.fastq.gz

Check files:

ls raw_reads

Step 3: Count total raw reads

echo "S_aureus raw R1 reads:"
echo $(( $(zcat raw_reads/S_aureus_R1.fastq.gz | wc -l) / 4 ))

echo "S_aureus raw R2 reads:"
echo $(( $(zcat raw_reads/S_aureus_R2.fastq.gz | wc -l) / 4 ))

echo "S_aureus total raw reads:"
echo $(( $(zcat raw_reads/S_aureus_R1.fastq.gz | wc -l) / 4 + $(zcat raw_reads/S_aureus_R2.fastq.gz | wc -l) / 4 ))

Step 4: FastQC on raw reads

conda activate qc

fastqc \
  -t 4 \
  -o qc_reports \
  raw_reads/S_aureus_R1.fastq.gz \
  raw_reads/S_aureus_R2.fastq.gz

Step 5: Trim reads

trimmomatic PE \
  -threads 4 \
  raw_reads/S_aureus_R1.fastq.gz \
  raw_reads/S_aureus_R2.fastq.gz \
  trimmed_reads/S_aureus_trim_R1.fastq.gz \
  trimmed_reads/S_aureus_unpaired_R1.fastq.gz \
  trimmed_reads/S_aureus_trim_R2.fastq.gz \
  trimmed_reads/S_aureus_unpaired_R2.fastq.gz \
  SLIDINGWINDOW:4:20 \
  LEADING:3 \
  TRAILING:3 \
  HEADCROP:15 \
  MINLEN:40

Step 6: Count total reads after trimming

echo "S_aureus trimmed R1 reads:"
echo $(( $(zcat trimmed_reads/S_aureus_trim_R1.fastq.gz | wc -l) / 4 ))

echo "S_aureus trimmed R2 reads:"
echo $(( $(zcat trimmed_reads/S_aureus_trim_R2.fastq.gz | wc -l) / 4 ))

echo "S_aureus total trimmed reads:"
echo $(( $(zcat trimmed_reads/S_aureus_trim_R1.fastq.gz | wc -l) / 4 + $(zcat trimmed_reads/S_aureus_trim_R2.fastq.gz | wc -l) / 4 ))

Step 7: FastQC on trimmed reads

fastqc \
  -t 4 \
  -o qc_reports \
  trimmed_reads/S_aureus_trim_R1.fastq.gz \
  trimmed_reads/S_aureus_trim_R2.fastq.gz

conda deactivate

Step 8: Genome assembly

conda activate assembly

spades.py \
  -t 4 \
  -o assemblies/S_aureus_assembly \
  -1 trimmed_reads/S_aureus_trim_R1.fastq.gz \
  -2 trimmed_reads/S_aureus_trim_R2.fastq.gz

cp assemblies/S_aureus_assembly/contigs.fasta assemblies/S_aureus_contigs.fasta

Step 9: Count contigs before filtering

grep -c "^>" assemblies/S_aureus_contigs.fasta

Step 10: QUAST on original assembly

quast \
  -t 4 \
  -o quast/S_aureus_original \
  assemblies/S_aureus_contigs.fasta

cat quast/S_aureus_original/report.txt

Step 11: CheckM2 on original assembly

checkm2 predict \
  -t 4 \
  -x .fasta \
  -i assemblies/S_aureus_contigs.fasta \
  -o checkm2/S_aureus_original

cat checkm2/S_aureus_original/quality_report.tsv

Step 12: Filter contigs smaller than 500 bp

seqkit seq \
  -m 500 \
  assemblies/S_aureus_contigs.fasta \
  > assemblies/S_aureus_contigs_filtered.fasta

Count filtered contigs:

grep -c "^>" assemblies/S_aureus_contigs_filtered.fasta

Step 13: QUAST on filtered assembly

quast \
  -t 4 \
  -o quast/S_aureus_filtered \
  assemblies/S_aureus_contigs_filtered.fasta

cat quast/S_aureus_filtered/report.txt

Record:

  • N50
  • L50

Step 14: CheckM2 on filtered assembly

checkm2 predict \
  -t 4 \
  -x .fasta \
  -i assemblies/S_aureus_contigs_filtered.fasta \
  -o checkm2/S_aureus_filtered

cat checkm2/S_aureus_filtered/quality_report.tsv

conda deactivate

Record:

  • completeness
  • contamination

Step 15: Taxonomic classification

conda activate classify

kraken2 \
  --db ~/databases/minikrakendb \
  --threads 4 \
  --report classify/S_aureus.kraken2.report \
  --output classify/S_aureus.kraken2.output \
  assemblies/S_aureus_contigs_filtered.fasta

head classify/S_aureus.kraken2.report

conda deactivate

Record the main taxonomy result.

Step 16: MLST

conda activate mlst

mlst assemblies/S_aureus_contigs_filtered.fasta > classify/S_aureus_mlst.txt

cat classify/S_aureus_mlst.txt

conda deactivate

If no MLST scheme is detected, write:

No MLST scheme detected

Step 17: Genome annotation

conda activate annotation

prokka \
  --outdir annotation/S_aureus \
  --force \
  --cpus 4 \
  --prefix S_aureus \
  --locustag SAUR \
  assemblies/S_aureus_contigs_filtered.fasta

cat annotation/S_aureus/S_aureus.txt

conda deactivate

Record:

  • number of CDS
  • number of tRNA

Step 18: AMR gene screening

conda activate amr

abricate \
  --db ncbi \
  --threads 4 \
  assemblies/S_aureus_contigs_filtered.fasta \
  > abricate/S_aureus_amr.tsv

abricate --summary abricate/S_aureus_amr.tsv > abricate/S_aureus_amr_summary.tsv

cat abricate/S_aureus_amr_summary.tsv

Record the AMR genes.

Step 19: Virulence factor screening

abricate \
  --db vfdb \
  --threads 4 \
  assemblies/S_aureus_contigs_filtered.fasta \
  > abricate/S_aureus_vf.tsv

abricate --summary abricate/S_aureus_vf.tsv > abricate/S_aureus_vf_summary.tsv

cat abricate/S_aureus_vf_summary.tsv

conda deactivate

Record the virulence factors.

Final comparison table

Complete the table below using your analysis outputs.

Result Leptospira yasudae SRR38327166 Staphylococcus aureus SRR32454207
Number of total raw reads
Number of total reads after trimming
Number of contigs, >=0 bp
Number of contigs after filtering, >=500 bp
N50
L50
Completeness
Contamination
Taxonomy
MLST
Number of CDS
Number of tRNA
AMR genes
Virulence factors

Where to find each result

Result Command or file
Number of total raw reads Count reads from raw_reads/*fastq.gz
Number of total reads after trimming Count reads from trimmed_reads/*fastq.gz
Number of contigs, >=0 bp grep -c "^>" assemblies/*_contigs.fasta
Number of contigs after filtering grep -c "^>" assemblies/*_contigs_filtered.fasta
N50 quast/*_filtered/report.txt
L50 quast/*_filtered/report.txt
Completeness checkm2/*_filtered/quality_report.tsv
Contamination checkm2/*_filtered/quality_report.tsv
Taxonomy classify/*.kraken2.report
MLST classify/*_mlst.txt
Number of CDS annotation/*/*.txt
Number of tRNA annotation/*/*.txt
AMR genes abricate/*_amr_summary.tsv
Virulence factors abricate/*_vf_summary.tsv

Files generated

At the end of this exercise, you should have:

raw_reads/
trimmed_reads/
qc_reports/
assemblies/
quast/
checkm2/
classify/
annotation/
abricate/

Important output files:

assemblies/L_yasudae_contigs_filtered.fasta
assemblies/S_aureus_contigs_filtered.fasta

quast/filtered/report.txt
checkm2/filtered/quality_report.tsv

classify/L_yasudae.kraken2.report
classify/S_aureus.kraken2.report

classify/L_yasudae_mlst.txt
classify/S_aureus_mlst.txt

annotation/L_yasudae/L_yasudae.txt
annotation/S_aureus/S_aureus.txt

abricate/amr_summary.tsv
abricate/vf_summary.tsv

Discussion questions

NoteQuestion 1

Which genome produced more contigs after assembly?

NoteQuestion 2

Which genome had the higher N50?

NoteQuestion 3

Which genome had higher completeness?

NoteQuestion 4

Did both organisms have an MLST result?

NoteQuestion 5

Which organism had more AMR genes?

NoteQuestion 6

Which organism had more virulence factors?

Key points

NoteImportant
  • SRA data can be downloaded using prefetch.
  • SRA files can be converted to FASTQ using fasterq-dump.
  • Raw and trimmed reads can be counted from FASTQ files.
  • Assembly quality should be checked before and after filtering.
  • Taxonomy, MLST, annotation, AMR, and virulence screening help interpret genome content.
  • The final table should be completed using your own analysis results.
Back to top