Metagenomic Assembly and Read Mapping

This walkthrough follows paired-end shotgun reads from host screening (when needed) through metaSPAdes assembly and read-back mapping. It focuses on file contracts and checks that let you judge whether the outputs are usable; mapping support alone does not validate an assembly.

What you need and what you will produce

Start with quality-controlled paired-end FASTQ files, sample_R1.fastq.gz and sample_R2.fastq.gz, with matching read names and mates in the same order. Adapter contamination matters to metaSPAdes, so inspect the upstream QC report before assembly. For host-associated samples, obtain an appropriate host reference FASTA and record its source and checksum. Host screening is optional for environmental samples when host contamination is not expected.

Install SPAdes (including metaSPAdes), BWA, and SAMtools using a package manager or your site's supported software environment; installation is not covered here. From the project working directory, confirm that each executable is available and record the versions reported by spades.py --version, bwa, and samtools --version. Before analysis, test compressed inputs with gzip -t sample_R1.fastq.gz sample_R2.fastq.gz and inspect the QC report and mate-name/order agreement. Stop if a check fails or inputs are not the same validated pair. The commands assume these files are in the current directory and adequate memory, storage, and CPU capacity; choose thread and memory values for your system. This is a command template, not a completed biological run.

The workflow creates spades_output/contigs.fasta, a coordinate-sorted and indexed BAM of reads aligned to those contigs, summary files, and an optional IGV view. Record tool versions, input checksums, reference provenance, and commands in a project log; proceed only when each command exits successfully and its expected output exists and is non-empty.

Screen host reads when appropriate

Map the paired reads to the host reference, then retain templates for which both mates are unmapped. The SAM flags 0x4 and 0x8 represent an unmapped read and an unmapped mate; -f 12 requires both. Excluding secondary and supplementary records avoids writing alternate alignment records as if they were additional read templates.

# Download the human GRCh38.p14 reference only when appropriate for the sample.
# Verify the checksum against a trusted source before using the reference.
set -euo pipefail
curl -L -o host_ref.fa.gz "https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.fna.gz"
sha256sum host_ref.fa.gz

gunzip host_ref.fa.gz
bwa index host_ref.fa

bwa mem -t 8 host_ref.fa sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools sort -@ 4 -o host_mapped.bam
samtools index host_mapped.bam
test -s host_mapped.bam && test -s host_mapped.bam.bai

samtools view -b -f 12 -F 2304 host_mapped.bam \
  | samtools sort -n -@ 4 -o nonhost.name_sorted.bam
samtools fastq -n \
  -1 microbial_R1.fastq \
  -2 microbial_R2.fastq \
  -0 /dev/null -s /dev/null \
  nonhost.name_sorted.bam
test -s microbial_R1.fastq && test -s microbial_R2.fastq

# Count FASTQ records as lines / 4 (counts reads, not pairs).
gzip -cd sample_R1.fastq.gz | awk 'END {print NR/4}'
awk 'END {print NR/4}' microbial_R1.fastq

samtools fastq requires reads grouped by name; the name-sort step supplies that ordering. The -1 and -2 outputs contain paired reads, while -0 and -s send singletons and other reads to /dev/null. Since the filter asks for both mates unmapped, such discarded reads should not ordinarily be produced, but explicit destinations avoid mixing them into the paired outputs. Confirm the two output files have equal FASTQ record counts and that the starting pair names correspond before proceeding.

Assemble the paired reads

metaSPAdes is SPAdes's metagenomic mode. Its documented short-read input contract supports one paired-end short-read library; this example therefore assumes one sample/library represented by the R1/R2 pair. The command does not combine independent libraries. The official manual warns that adapter read-throughs can affect metaSPAdes, so upstream QC remains important.

spades.py --meta \
  -1 microbial_R1.fastq \
  -2 microbial_R2.fastq \
  -o spades_output \
  -t 16 -m 64

On successful completion, inspect the run log and the output directory. Use spades_output/contigs.fasta as the contig reference for the mapping step; retain params.txt and the logs with the analysis record. Do not assume that scaffolds were produced or that the assembly succeeded simply because a command was submitted. If SPAdes reports an input, memory, or output-directory error, resolve that before continuing.

Map reads back to contigs

Build the BWA index from the same contig FASTA that you will later load as the reference. Align the same paired reads used for assembly (screened reads here), sort the resulting alignments by genomic coordinate, then index the BAM. Coordinate order is required for efficient regional access and IGV viewing.

set -euo pipefail
bwa index spades_output/contigs.fasta
bwa mem -t 8 spades_output/contigs.fasta microbial_R1.fastq microbial_R2.fastq \
  | samtools sort -@ 4 -o reads_to_contigs.bam
samtools index reads_to_contigs.bam
test -s reads_to_contigs.bam && test -s reads_to_contigs.bam.bai
samtools flagstat reads_to_contigs.bam > reads_to_contigs.flagstat.txt
samtools idxstats reads_to_contigs.bam > reads_to_contigs.idxstats.txt
samtools depth -a reads_to_contigs.bam > contig_depth.tsv

flagstat summarizes alignment-record categories, while idxstats reports reference sequence lengths and mapped/unmapped read counts per reference sequence. depth -a reports per-position depth, including positions with zero depth; its output can be large for large assemblies. Both summary commands should be interpreted with their counting units in mind: the BAM contains alignment records, and paired reads can contribute two records.

Calculate a simple mean of the emitted per-position depth for each contig:

awk '{sum[$1] += $3; n[$1]++} END {for (c in sum) print c, sum[c]/n[c]}' \
  contig_depth.tsv > contig_mean_coverage.tsv

Because depth -a includes zero-depth positions, this mean includes uncovered positions in the reported reference contigs. It is a descriptive average, not a genome abundance estimate or a measure of assembly correctness. Keep the per-position file if you need to examine variation along a contig.

Inspect the evidence before interpreting contigs

1. Check the run and read accounting. Confirm SPAdes completed, compare starting and post-screening FASTQ record counts, and check that the two retained mate files remain paired. Unexpected count changes can indicate filtering or file-handling problems.
2. Review mapping summaries. Read flagstat and idxstats alongside the read counts. The fraction mapping back depends on sequencing depth, community complexity, filtering, and sequence divergence. It is not a universal assembly pass threshold.
3. Look at depth and contig lengths together. Very low support, zero-coverage stretches, or abrupt depth changes are reasons to investigate. Repeats, conserved genes, collapsed sequences, strain mixtures, and uneven organism abundance can also create nonuniform coverage, so depth alone cannot diagnose a misassembly.
4. Inspect candidate joins in IGV. Create the FASTA index with samtools faidx spades_output/contigs.fasta, then load that exact FASTA and the indexed coordinate-sorted BAM. Examine paired-read placement around suspicious regions. A visual pattern is a prompt for further validation, not proof of a chimera.

There is no general coverage cutoff such as 5× that separates valid from artifactual metagenomic contigs across samples. Likewise, contig length or a high mapping percentage alone does not establish organism identity, genome completeness, correct binning, or biological function. Interpret those questions with appropriate downstream evidence and sample-specific context.

Common problems and next actions

• Few reads map back: verify that the correct R1/R2 files were assembled and mapped, check counts before and after host filtering, and inspect SPAdes logs for completion or input warnings. Low depth or high community complexity can also leave reads unassembled; the mapping percentage by itself does not identify the cause.
• Many short contigs: consider read quality, sequencing depth, community complexity, and strain variation before changing assembler settings. Do not alter k-mer options without checking the version-specific SPAdes documentation and the study design.
• Unexpectedly few retained host-screened pairs: verify the host reference, sample identity, and filtering flags. Host screening can remove microbial reads that resemble host sequence; compare counts and examine the reference choice rather than treating a retained-read total as a target.
• IGV cannot display the tracks: confirm BAM and BAI names match, BAM is coordinate sorted, and the loaded reference is the exact contig FASTA used for mapping. Index that FASTA with samtools faidx.

Limits and responsible handling

Short-read community assembly can collapse related strains or repeats, and host screening depends on the reference used. Read-back mapping documents sequence support under the aligner's settings; it cannot establish that a contig is complete, correctly assembled, correctly binned, or free of host sequence. Retain commands, tool versions, checksums, logs, and reference provenance so the analysis can be reviewed.

For human-derived material, follow institutional approvals, consent, and data-sharing requirements. Do not treat host-read removal alone as proof that a dataset is safe to release.

References and software documentation

• SPAdes command-line manual and SPAdes output files
• BWA source and usage examples
• SAMtools manual
• IGV desktop documentation
• Quince et al., shotgun metagenomics and assembly review