mSWEEP User Manual
This document details how to use the rust port / rewrite of mSWEEP, available at https://codeberg.org/themaklin/msweep-rs.
msweep-rs incorporates the functionality of both mSWEEP and mGEMS. The tool also provides a method for quality checking the binned reads, inspired by demix_check.
Tools comprising the pipeline
See installation for recommended ways to obtain these.
Abundance estimation and binning
- msweep-rs
This manual does not cover the C++ mSWEEP/mGEMS packages.
Pseudoalignment
One of
Themisto is recommended for indexes that are sparse, eg. those containing multiple species.
Fulgor is recommended for indexes that are dense, eg. containing many assemblies from a single species.
Both tools can be used for either case. The results may vary slightly due to differences in the algorithms they use.
Metagraph and Bifrost are supported but not recommended. Alignments in the SAM format are also supported but must be processed with ahda.
Using the pseudoalignment outputs
msweep-rs can natively read input from either, but it is recommended to use ahda to compress the results for better I/O performance and to save disk space.
Defining target groups/clusters
PopPUNK is recommended for clustering the reference assemblies.
Installation
Requirements
Exact requirements may vary depending on what channels are used to install the tools.
- Rust compiler, [https://rust-lang.org/learn/get-started/].
- C++ compiler
- git
- conda
How-to
msweep-rs
From source:
git clone https://codeberg.org/themaklin/msweep-rs
cd msweep-rs
git clone https://github.com/tmaklin/sketchlib.rust vendor/sketchlib.rust -b no-panic-from-rollhash-set_k
git clone https://github.com/tmaklin/sparrowhawk vendor/sparrowhawk -b assembly-api
cd vendor/sketchlib.rust
git apply ../../patches/sketchlib.patch
cd ../../
cd vendor/sparrowhawk
git apply ../../patches/sparrowhawk.patch
cd ../../
cargo build --release
ahda
From crates.io:
cargo install ahda
From source:
git clone https://codeberg.org/themaklin/ahda -b v0.1.1
cd ahda
cargo build --release
Pseudoaligners
Themisto
Download a precompiled binary from GitHub or see https://github.com/algbio/themisto for other options.
Fulgor
Install from bioconda or see https://github.com/jermp/fulgor for other options.
PopPUNK
Install from bioconda or see https://github.com/bacpop/poppunk for other options.
Usage example
This will demonstrate how to use msweep-rs to analyse a (simulated) sequencing read set with multiple E. coli strains in it.
Reference data
Setup
Download ~1900 hybrid Nanopore+Illumina E. coli assemblies described in Arredondo-Alonso et al. (2025) from: https://www.ncbi.nlm.nih.gov/datasets/genome/?bioproject=PRJEB57633. In general, it’s OK to use short-read only assemblies as the reference, the hybrid dataset is just used as an example here.
Unzip the assemblies and create the input list needed for poppunk:
unzip ncbi_dataset.zip
paste <(ls ncbi_dataset/data/*/*.fna | xargs basename -a | cut -d. -f1) <(ls -d ncbi_dataset/data/*/*.fna) > poppunk_input_list.tsv
The input list popppunk_input_list.tsv will look like the following:
GCA_947669175 ncbi_dataset/data/GCA_947669175.1/GCA_947669175.1_30345_1_256_genomic.fna
GCA_964035765 ncbi_dataset/data/GCA_964035765.1/GCA_964035765.1_30134_6_107_genomic.fna
GCA_964036675 ncbi_dataset/data/GCA_964036675.1/GCA_964036675.1_30348_2_324_genomic.fna
GCA_964036685 ncbi_dataset/data/GCA_964036685.1/GCA_964036685.1_30345_1_313_genomic.fna
GCA_964036695 ncbi_dataset/data/GCA_964036695.1/GCA_964036695.1_30602_1_150_genomic.fna
GCA_964036715 ncbi_dataset/data/GCA_964036715.1/GCA_964036715.1_30224_1_362_genomic.fna
GCA_964036725 ncbi_dataset/data/GCA_964036725.1/GCA_964036725.1_30345_1_127_genomic.fna
GCA_964036735 ncbi_dataset/data/GCA_964036735.1/GCA_964036735.1_30348_2_197_genomic.fna
GCA_964036745 ncbi_dataset/data/GCA_964036745.1/GCA_964036745.1_30345_1_45_genomic.fna
GCA_964036755 ncbi_dataset/data/GCA_964036755.1/GCA_964036755.1_30348_1_161_genomic.fna
Creating the poppunk clustering
Create the poppunk database:
poppunk --create-db --r-files poppunk_input_list.tsv --output NORM_pp_db
Create the poppunk clusters:
poppunk --fit-model dbscan --ref-db NORM_pp_db --output NORM_pp_dbscan
For more detailed information on using poppunk, see the documentation at: https://poppunk-docs.bacpop.org.
The clustering is located at NORM_pp_dbscan/NORM_pp_dbscan_clusters.csv and looks something like this:
Taxon,Cluster
GCA_964040825,1
GCA_964044985,1
GCA_964055635,1
GCA_964042775,1
GCA_964055745,1
GCA_964039205,1
GCA_964051785,1
GCA_964045395,1
GCA_964048425,1
To inspect the quality of the poppunk clustering, we can download the ST labels from another study and merge them with the poppunk results:
wget https://zenodo.org/records/17608175/files/NORM_1991_colibactin_coevolution_metadata.tsv
sed '1d' NORM_1991_colibactin_coevolution_metadata.tsv | tr -d '["]' | cut -f1,7 >> NORM_ST_assignments.tsv
join -1 1 -2 1 <(sed '1d' NORM_pp_dbscan/NORM_pp_dbscan_clusters.csv | tr ',' '\t' | sort -k 1,1) <(sort -k1,1 NORM_ST_assignments.tsv) | sort -gk2,2 | tr ' ' '\t' > NORM_clusters_with_st_labels.tsv
We can see that the poppunk clusters (2nd column) roughly correspond to clonal complexes of the different STs (third column):
> shuf NORM_clusters_with_st_labels.tsv | head | sort -gk2,2
GCA_964045755 1 73
GCA_964049165 1 73
GCA_964051775 1 73
GCA_964042555 2 95
GCA_964042625 2 95
GCA_964042865 2 95
GCA_964043495 2 95
GCA_964046195 3 69
GCA_964051275 5 404
GCA_964042815 7 12
This is generally a good clustering resolution for msweep if the species is clonal.
Indexing the reference data with Fulgor
Extract the paths to the assemblies from the poppunk input list we created earlier:
cut -f2 poppunk_input_list.tsv > NORM_input_paths.txt
Build the fulgor index with
fulgor build -l NORM_input_paths.txt -o fulgor_NORM_index -k 31 -m 15 -d tmp -g 8 -t 10 --meta --diff
This will produce three indexes
fulgor_NORM_index.fur
fulgor_NORM_index.mdfur
fulgor_NORM_index.mfur
The index we want to use is fulgor_NORM_index.mdfur which is the smallest in size (185M compared to 403M and 277M).
Note we use the --meta and --diff flags, which instruct Fulgor to build a
compressed index. These two flags provide good compression for indexes that are
dense, eg. sequences from the same species. The size of the index determines the
amount of memory required to align reads against it, so using these flags is
recommended when applicable although they increase the amount of time required
to construct the index.
Matching the poppunk clusters with the fulgor index
Fulgor reorders the input sequences, so we need to match the poppunk clusters with the order they are stored in the index file fulgor_NORM_index.mdfur:
fulgor print-filenames -i fulgor_NORM_index.mdfur | sed '1,3d' | cut -f2 | xargs basename -a | cut -d. -f1 > fulgor_filenames.txt
while read -r accession; do grep "^$accession\t" NORM_clusters_with_st_labels.tsv; done < fulgor_filenames.txt | cut -f1,2 > fulgor_clusters.tsv
Setting up the test query
To set up a simulated strain mixture, we mix Illumina data from three isolate sequencing runs from the same cohort that the reference data was obtained from. These mixed samples were not included in our hybrid assembly data set.
Download the read files:
wget -O ERR4034704-ST144_1.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/004/ERR4034704/ERR4034704_1.fastq.gz
wget -O ERR4034704-ST144_2.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/004/ERR4034704/ERR4034704_2.fastq.gz
wget -O ERR4036743-ST95_1.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/003/ERR4036743/ERR4036743_1.fastq.gz
wget -O ERR4036743-ST95_2.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/003/ERR4036743/ERR4036743_2.fastq.gz
wget -O ERR4037355-ST127_1.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/005/ERR4037355/ERR4037355_1.fastq.gz
wget -O ERR4037355-ST127_2.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR403/005/ERR4037355/ERR4037355_2.fastq.gz
Concatenate the reads from these files:
gunzip -c ERR4034704-ST144_1.fastq.gz ERR4036743-ST95_1.fastq.gz ERR4037355-ST127_1.fastq.gz | gzip --fast > mixed_1.fastq.gz
gunzip -c ERR4034704-ST144_2.fastq.gz ERR4036743-ST95_2.fastq.gz ERR4037355-ST127_2.fastq.gz | gzip --fast > mixed_2.fastq.gz
Aligning the reads
We’re now ready to run the msweep pipeline.
Running fulgor
To start, align the reads in mixed_1.fastq.gz and mixed_2.fastq.gz against the fulgor index fulgor_NORM_index.mdfur that we generated earlier
fulgor pseudoalign -i fulgor_NORM_index.mdfur -q mixed_1.fastq.gz -r 1.0 -o mixed_1.txt
fulgor pseudoalign -i fulgor_NORM_index.mdfur -q mixed_2.fastq.gz -r 1.0 -o mixed_2.txt
The -r 1.0 controls what percentage of the k-mers in the reads should be present in the index to report an alignment. Generally, it’s recommended to use -r 1.0 or a value very close to 1.0 as this also functions as a filter for reads that originate from eg. host DNA. The equivalent argument for Themisto is --threshold 1
Compressing the output alignment
The file mixed_1.txt output in the previous step is rather large (11G). We can make it much smaller by piping the output from fulgor to ahda for compression:
fulgor pseudoalign -i fulgor_NORM_index.mdfur -q mixed_1.fastq.gz -r 1.0 -o /dev/stdout | ahda encode -q mixed_1.fastq.gz -t fulgor_filenames.txt > mixed_1.txt.ahda
fulgor pseudoalign -i fulgor_NORM_index.mdfur -q mixed_2.fastq.gz -r 1.0 -o /dev/stdout | ahda encode -q mixed_2.fastq.gz -t fulgor_filenames.txt > mixed_2.txt.ahda
This will compress the output alignment and encode the target and fragment names in the output file mixed_1.txt.ahda so they can be used later to extract binned reads from the original file. The binary format produced by ahda is also faster to parse for msweep than the original text file.
Compared to gzip --fast, the ahda file is 5% of the size (83M vs. 1.6G). Compared to gzip --best, the ahda file is 12% of the size (83M vs. 689M).
(Optional) intersect the alignment files
We can further reduce the size of the alignment files mixed_1.txt.ahda and mixed_2.txt.ahda by merging them:
ahda set mixed_1.txt.ahda mixed_2.txt.ahda > mixed.txt.ahda
rm mixed_1.txt.ahda
rm mixed_2.txt.ahda
This intersection is performed anyway when running msweep so it does not change the results.
Estimating abundances and binning
Run msweep-rs with:
mSWEEP2 estimate -i fulgor_clusters.tsv --bin-reads -o mixed mixed.txt.ahda
This will generate the mixed.abundances.tsv and mixed.bins.ahda files which contain the relative abundances and the read assignments to the reference clusters.
We can view the relative abundances file mixed.abundances.tsv with sort -gk2 mixed.abundances.tsv | tail:
5 0.0000026 0
50 0.0000026 0
7 0.0000033 0
82 0.0000110 0
118 0.0000195 0
1 0.0002956 1714
111 0.0003842 89
14 0.3202538 1159651
6 0.3254141 1081402
2 0.3535427 1199918
The columns contain the following information:
- First column is the name of the cluster supplied via
fulgor_clusters.tsv. - Second column is the relative abundance of this cluster.
- Third column is the number of reads assigned to this cluster.
We can check what sequence types the clusters correspond to by looking at the `` file:
grep "\t14\t" NORM_clusters_with_st_labels.tsv | head -n1
> GCA_964038655 14 144
grep "\t6\t" NORM_clusters_with_st_labels.tsv | head -n1
> GCA_964037375 6 127
grep "\t2\t" NORM_clusters_with_st_labels.tsv | head -n1
> GCA_947669175 2 95
ie. we found ST144, ST127, and ST95 in roughly equal proportions. This matches the number of reads in the inputs that we originally mixed, so the results appear correct.
Verifying the relative abundances and read assignments
Sometimes, msweep can falsely report the presence of a cluster in a sample because of ambiguity in the input alignments. msweep-rs includes a tool to verify the results based on comparing distances from the assigned reads to the reference assemblies.
To perform this verification, we need to perform some additional setup on the reference assemblies. The following command creates a list that matches each input assemble to its poppunk cluster:
join -1 1 -2 1 <(sort -k1,1 poppunk_input_list.tsv) <(sed '1d' NORM_pp_dbscan/NORM_pp_dbscan_clusters.csv | tr ',' '\t' | sort -k1,1) | tr ' ' '\t' | cut -f2,3 > msweep_setup_input_list.tsv
This file looks something like the following:
ncbi_dataset/data/GCA_947669175.1/GCA_947669175.1_30345_1_256_genomic.fna 2
ncbi_dataset/data/GCA_964035765.1/GCA_964035765.1_30134_6_107_genomic.fna 4
ncbi_dataset/data/GCA_964036675.1/GCA_964036675.1_30348_2_324_genomic.fna 5
ncbi_dataset/data/GCA_964036685.1/GCA_964036685.1_30345_1_313_genomic.fna 18
ncbi_dataset/data/GCA_964036695.1/GCA_964036695.1_30602_1_150_genomic.fna 21
ncbi_dataset/data/GCA_964036715.1/GCA_964036715.1_30224_1_362_genomic.fna 1
ncbi_dataset/data/GCA_964036725.1/GCA_964036725.1_30345_1_127_genomic.fna 2
ncbi_dataset/data/GCA_964036735.1/GCA_964036735.1_30348_2_197_genomic.fna 1
ncbi_dataset/data/GCA_964036745.1/GCA_964036745.1_30345_1_45_genomic.fna 3
ncbi_dataset/data/GCA_964036755.1/GCA_964036755.1_30348_1_161_genomic.fna 7
Next, we pass this file to mSWEEP2 setup to sketch the sequences belonging to each cluster:
mkdir sketches
mSWEEP2 setup -i msweep_setup_input_list.tsv -o sketches
This will create a number of reference sketches in the sketches/ directory that can be used to compare the reads assigned to each cluster.
The comparison can be performed by running:
mSWEEP2 check -s sketches/sketch_list.tsv -b mixed.bins.ahda mixed_1.fastq.gz mixed_2.fastq.gz > mixed.abundances.check.tsv
We can view the quality of the assignments in mixed.abundances.check.tsv with sort -gk2 mixed.abundances.check.tsv | tail:
5 0.0000026 0 0
50 0.0000026 0 0
7 0.0000033 0 0
82 0.0000110 0 0
118 0.0000195 0 0
1 0.0002956 1714 0.05
111 0.0003842 89 0
14 0.3202538 1159651 0.93
6 0.3254141 1081402 0.93
2 0.3535427 1199918 0.93
Compared to mixed.abundances.tsv, this new file has a fourth column that
contains a likelihood ratio of the distribution of k-mer distances between the
reads compared to the distribution of k-mer distances between the reference
assemblies. Values close or above 1.0 indicate a good match between the two
and imply that the cluster was identified correctly.
Internally, mSWEEP2 check uses sketchlib.rust to compute the distances.
Verifying assignments with a low read count
The standard workflow for mSWEEP2 check uses the sequencing reads to compute
the k-mer distances. This approach may fail in some cases, for example if the
sequencing depth is low, as there are not enough reads to correctly compute the
distribution of the distances.
In this case, we can perform a draft assembly of the reads and compare the distances from this draft against the references by adding the --assemble flag:
mSWEEP2 check -s sketches/sketch_list.tsv -b mixed.bins.ahda --assemble mixed_1.fastq.gz mixed_2.fastq.gz > mixed.abundances.check.tsv
Looking at the quality in mixed.abundances.check.tsv now shows some additional information:
5 0.0000026 0 0 0 0
50 0.0000026 0 0 0 0
7 0.0000033 0 0 0 0
82 0.0000110 0 0 0 0
118 0.0000195 0 0 0 0
1 0.0002956 1714 0.05 0.02 5079 3
111 0.0003842 89 0 0 0
14 0.3202538 1159651 0.93 0.77 4838502 570
6 0.3254141 1081402 0.93 1.02 4835524 500
2 0.3535427 1199918 0.93 0.80 4816063 623
The first four columns are the same as before. The new columns contain the following information:
- Fifth column contains the likelihood ratio of distances computed from the draft assembly against the reference assemblies.
- Sixth column contains the length of the draft assembly, computed from contigs of length 1000bp or more that were used to comptue the value in the fifth column.
- Seventh column contains the number of contigs that are 1000bp or longer.
Similarly to the reads only checking, the likelihood ratio in the fifth column is close or above 1 when the assignment is good. However, the distance from the assembly may falsely indicate a good fit if the draft assembly is substantially shorter than a real genome from this target would be. In this case, it’s useful to look at the draft assembly length and number of contigs to verify that the assembly used is approximately the length expected, 5M bases in this case for an E. coli genome.
Internally, the assembly checking uses sparrowhawk to perform the draft assembly. This is a relatively low-resource assembly method that tends to produce more fragmented assemblies compared to something more computationally intense. If you’d like to access these, add the --write-assemblies flag to mSWEEP2 check.
Thresholding the likelihood ratio values
For general use, it’s recommended to use something similar to the following to consider an assignment good:
- Likelihood ratio from the reads above 0.60.
- Likelihood ratio from the assembly above 0.60.
- Assembly length within ~75%-125% of expected.
If you run further QC on the extracted data, these thresholds may be relaxed to eg. to 0.20 for the likelihood ratios and the length may be ignored.
If you plan to assemble the reads assigned to each target group, it’s recommended to relax the thresholds and then use the assembly quality as a further indicator for which detections to keep.
Extracting the assigned reads
After running estimation and the optional (but recommended) verification, you can use mSWEEP2 extract to get the reads assigned to a target:
mSWEEP2 extract -q mixed_1.fastq.gz -t 14 mixed.bins.ahda | gzip --fast > mixed_cluster14_1.fastq.gz
mSWEEP2 extract -q mixed_2.fastq.gz -t 14 mixed.bins.ahda | gzip --fast > mixed_cluster14_2.fastq.gz
This will write the reads assigned to the target cluster 14 to the files mixed_cluster14_1.fastq.gz and mixed_cluster14_2.fastq.gz
You can change the cluster via the -t 14 argument.
These reads are good to go for downstream analyses, for example we can assemble them using shovill:
shovill --R1 mixed_cluster14_1.fastq.gz --R2 mixed_cluster14_2.fastq.gz --gsize 5M --outdir cluster14 --minlen 1000
This generates an assembly of length 5 126 071 in 106 contigs which is comparable to the isolate assembly from the unmixed reads (length 5 130 696 in 95 contigs).
Contributing to mWEEP
Manual
Any type of contributions to the manual are welcome. You can submit contributions to the git repository at https://codeberg.org/themaklin/msweep-user-manual.
Code
Source code for msweep-rs is available at https://codeberg.org/themaklin/msweep-rs. Contributions are welcome.
msweep-rs relies on the following code for certain key parts in the pipeline:
- Pseudoalignment processing: ahda.
- Statistical inference: mixt.
- Computing k-mer distances: sketchlib.rust.
- Sequence assembly: sparrowhawk-asm.
If your contribution concerns one of these parts, prefer making it to the particular project’s repository.
LLM policy
Contributions that are mostly/entirely produced by LLMs (“generative AI”) are not accepted in msweep-rs nor in the manual.
You must disclose any LLM usage when making a pull request.
If you believe that a feature should be added, prefer opening an issue instead of submitting an LLM generated pull request.