Detection of pseudouridine, m6A, and m5C on individual molecules from direct RNA nanopore signals.
Detailed documentation is available at https://comprna.github.io/SWARM/
Demo dataset and commands available at https://github.com/comprna/SWARM/tree/main/demo
Benchmarking scripts available at https://github.com/comprna/SWARM/tree/main/SWARM_scripts/benchmark
SWARM was trained on signals basecalled with guppy 6.4.6 for RNA002 and with dorado 0.7.2 for RNA004.
Recommended parameters RNA002:
guppy_basecaller -i $INPUTDIR --recursive -s $output_path.fastq -c guppy/ont-guppy/data/rna_r9.4.1_70bps_hac.cfg --device cuda:all:100%
Recommended parameters RNA004:
MODEL=dorado-0.7.2-linux-x64/rna004_130bps_sup@v5.0.0
dorado basecaller $MODEL $INPUTDIR -r -x cuda:all --emit-fastq > $output_path.fastq
minimap 2.24 for alignment and samtools 1.12 for quality checks
Recommended parameters: -k 5 for sythetic IVTs and -k 14 for human transcriptomes
minimap2 -ax map-ont -k 5 ${fasta} ${input_path}/guppy_pass.fastq | samtools sort -o ${output_path}.bam
samtools index ${output_path}.bam
samtools view -b -F 2324 ${bam_file}.bam > ${bam_file}_pass_filtered.bam
samtools index ${bam_file}_pass_filtered.bam
This step is optional but highly recommended, especially for large datasets.
https://github.com/hasindu2008/slow5tools
Example conversion command:
#convert fast5 files to slow5 files using 8 I/O processes
slow5tools f2s $INPUT_DIR -d $TEMPDIR -p 8
#Merge all the slow5 files in to a single file using 8 threads
slow5tools merge $TEMPDIR -o $OUTDIR/${SAMPLE}.blow5 -t 8
#remove the temporary directory
rm -rf $TEMPDIR
Our workflow supports both f5c sam and nanopolish tsv formats. We highly recommend opting for f5c and sam files. This requires the slow5 conversion outlined in the previous step.
https://github.com/hasindu2008/f5c
Example event align command:
## f5c index
f5c index -t 48 $FASTQ_PATH --slow5 $SLOW5_PATH
## SQK-RNA002 event alignment
f5c eventalign -t 48 -r $FASTQ_PATH --rna -g $genome -b $BAM --slow5 $SLOW5_PATH --min-mapq 0 --secondary=yes --signal-index --scale-events --samples --print-read-names --sam > $OUT
## SQK-RNA004 event alignment
# first download the pore model used for SWARM training
wget https://raw.githubusercontent.com/hasindu2008/f5c/v1.3/test/rna004-models/rna004.nucleotide.5mer.model
# then run f5c giving the downloaded kmer-model path
f5c eventalign --kmer-model /path/to/rna004.nucleotide.5mer.model -t 48 -r $FASTQ_PATH --rna -g $genome -b $BAM --slow5 $SLOW5_PATH --min-mapq 0 --secondary=yes --signal-index --scale-events --samples --print-read-names --sam > $OUT
We used this format in earlier stages of the project, our workflow can still support it for SQK-RNA002. Note that our prediction workflow is optimised for f5c sam format.
https://github.com/jts/nanopolish
Example event align command:
nanopolish index -d ${fast5_path} -s ${guppy_files}/sequencing_summary.txt $fastq
nanopolish eventalign -t 48 --reads $fastq --bam $bam_file \
--genome $fasta --signal-index --scale-events --samples --print-read-names > $output_path
We provide a standalone containerised environment for running SWARM on Zenodo:https://zenodo.org/records/22123294
wget https://zenodo.org/records/22123294/files/SWARM.sif
singularity can be installed as per: https://docs.sylabs.io/guides/3.11/user-guide/quick_start.html#quick-installation-steps
You can run scripts from the SWARM repo using singularity and /opt/SWARM/path/to/script
# For example to run SWARM_read_level.py located at SWARM/SWARM_scripts/SWARM_read_level.py
singularity exec --nv SWARM.sif python3 /opt/SWARM/SWARM_scripts/SWARM_read_level.py --OPTIONS
*The image was built using singularity v3.11.0 and GO v1.18.2
*The image was tested on x86-64 Linux systems running CentOS and Ubuntu, and with NVIDIA Volta, Hopper, and Blackwell GPUs
For more advanced users and users with tensorflow already installed and GPU-configured. Skip if using singularity.
Clone SWARM from github (install lfs to download model h5 files)
git lfs install
git clone https://github.com/comprna/SWARM/ && cd SWARM
Compile SWARM_preprocess and prerequisites, should take ~2 minutes.
cd SWARM_scripts/preprocess/
#build and compile htslib, slow5tools, SWARM_preprocess
bash build.sh
If your HPC has a tensorflow module, simply load tensorflow and use the loaded python path for creating venv:
module load tensorflow/2.15.0
python3 -m venv swarm_env
source swarm_env/bin/activate
python3 -m pip install pysam==0.22.1 numpy==1.26.2 pandas==2.2.0 scikit-learn==1.4.0
# make sure to activate the venv before running SWARM read-level and site-level prediction
module load tensorflow/2.15.0
source /PATH/TO/swarm_env/bin/activate
If tensorflow is not pre-installed contact your HPC support or follow the official guide: https://www.tensorflow.org/install
Example bash code to run SWARM read-level prediction:
MOD=m6A # [<m6A> <m5C> <pU>]
FASTA=Homo_sapiens.GRCh38.cdna.fa
BLOW5=Hek293_mRNA.blow5
SAM=Hek293_mRNA_f5C.events.sam
OUT=Hek293_mRNA.$MOD.pred.tsv
# using singularity
SCRIPT=/opt/SWARM/SWARM_scripts/SWARM_read_level.py
singularity exec --nv SWARM.sif python3 $SCRIPT -m $MOD --sam $SAM --fasta $FASTA --raw $BLOW5 -o $OUT
# using pre-installed tensorflow
module load tensorflow/2.15.0
source /PATH/TO/swarm_env/bin/activate
SCRIPT=path/to/SWARM_read_level.py
python3 $SCRIPT -m $MOD --sam $SAM --fasta $FASTA --raw $BLOW5 -o $OUT
*Models for RNA002 or RNA004 chemistry are automatically selected based on the blow5 metadata.
First sort the read-level output, use cat if pooling multiple replicates.
cat Hek293_mRNA_rep1_pU.pred.tsv Hek293_mRNA_rep2_pU.pred.tsv > Hek293_mRNA_pooled_pU.pred.tsv
sort -k 1 Hek293_mRNA_pooled_pU.pred.tsv > Hek293_mRNA_pooled_pU.pred.tsv.sorted
Run site-level detection on sorted read-level data:
INPUT=Hek293_mRNA_pooled_pU.pred.tsv.sorted
OUT=Hek293_mRNA_pooled_pU.site.pred.tsv
# using singularity
SCRIPT=/opt/SWARM/SWARM_scripts/SWARM_site_level.py
singularity exec --nv SWARM.sif python3 $SCRIPT -i $INPUT -o $OUT
# using pre-installed tensorflow
module load tensorflow/2.15.0
source /PATH/TO/swarm_env/bin/activate
SCRIPT=/path/to/SWARM/SWARM_scripts/SWARM_site_level.py
python3 $SCRIPT -i $INPUT -o $OUT
*Site-level models are automatically selected based on the read-level tsv files.
Run differential modification test to find reference coordinates with stoichiometry changes across different conditions. Handles multiple replicates.
Set up a tab-separated config file with paths to inputs. Make sure that the header is present.
M2_file_path RepName Condition
WT_rep1.site.pred.tsv 1 WT
WT_rep2.site.pred.tsv 2 WT
KD_rep1.site.pred.tsv 1 KD
KD_rep2.site.pred.tsv 2 KD
Run SWARM_diff on site-level predictions (12 threads):
python3 SWARM_diff.py --data_file diff_config.tsv --output_file diff_out.tsv -n 12
Use --modsam tag with SWARM_read_level.py to get mod.sam output (pred.tsv is still produced too).
Note that this runs slower as multithreaded preprocessing is not implemented with modsam.
python3 SWARM_read_level.py -m $MOD --sam $SAM --fasta $FASTA --raw $BLOW5 -o $OUT --modsam
mod.sam can also be generated from sorted read-level pred.tsv files.
This should be faster on large datasets and also enables filtering of sites for cleaner visualisation.
prediction_sorted=Hek293_mRNA_pooled_pU.pred.tsv.sorted
BAM=Hek293_mRNA.bam
SITES=Hek293_mRNA_pooled_pU.site.pred.tsv
OUT=Hek293_mRNA_pooled_pU.site.pred.tsv
python3 convert_tsv_to_modsam.py -i $prediction_sorted -b $BAM --site_level $SITES -o $OUT
This optional step reduces the time to retrain models as preprocessing only a fraction of signals from a whole sample is usually enough for training. We trim for events comprising 500 signals per 9mer.
python3 train_models/read-level/trim_tsv_events.py -i <eventalign.tsv> -o <out_prefix> --limit-out 500
Preprocess trimmed files for read-level input features. Make sure to include --out_counter arg here!
python3 SWARM_read_level.py --preprocess -m <pU/m6A/m5C> --bam <BAM> //
--nanopolish <eventalign_trimmed.tsv> -o <out_prefix> --out_counter
Use this step for stratified sampling of the preprocessed signals. Splits equal number of signals per 9mer for positive/negative labels in each of train/validation/test set (60/20/20 split by default). Run on each sample, to make positive/negative data. Give same outpath if multiple samples are intended to be used under the same positive/negative label.
python3 train_models/read-level/split_training_by_9mers.py -i <preprocesed.pickle> //
--counts <preprocessed.counts> -o <outpath> --limit <signals_per_9mer> //
[--train_percent 0.6] [--validate_percent 0.2]
Use this step for finalising training/validation/testing data with matching positive/negative labels.
python3 train_models/read-level/assemble_data.py --input_positive <positive_prefix> //
--input_negative <negative_prefix> -o <outpath> [--positive_label 1] [--negative_label 0]
Use this step for training binary classifier of modification states with single-base single-molecule resolution.
python3 train_models/read-level/train_read-level.py -i <assembled_prefix> -o <outpath> //
[--vector_len 36] [--features 7] [--labels 2] [--epochs 100]
This optional step reduces the time to retrain models as preprocessing only a fraction of signals from a whole sample is usually enough for training. We trim for events comprising 1000 signals per 9mer (skipping first 500 used for training).
python3 train_models/site-level/trim_nanopolish_site-level.py -i <eventalign.tsv> -o <out_prefix> --limit_out 1500 --prev_limit 500 --out_counter
Preprocess trimmed files for read-level input features. Make sure to include --out_counter arg here!
python3 SWARM_read_level.py --preprocess -m <pU/m6A/m5C> --bam <BAM> //
--nanopolish <eventalign_trimmed.tsv> -o $out_prefix --out_counter
Predict read-level probabilities on preprocessed data
python3 SWARM_read_level.py --predict -m $MOD --pickle $out_prefix.pickle -o $out_prefix.pred.tsv
Create random mixtures of read-level probabilities with artificial coverages and stoichiometries for training site-level models. Need one pediction file for unmodified and one for target modification data. Provide cannonical nucleotide of the target modification to --center.
python3 train_models/site-level/make_site-level_features_parallel.py -NM unmodified.pred.tsv -M modified.pred.tsv --center <A/C/T> -o ${OutDir} -L 40 --kmer All --features Standard
Use this step for generating train/validation/test set (60/20/20 split by default). Run on each created mixture to make positive/negative data
# split positive data
python3 train_models/site-level/split_training_testing_site-level.py ${OutDir}/XWT_added_train.p WT_added ${OutDir}/
#split negative data
python3 train_models/site-level/split_training_testing_site-level.py ${OutDir}/XKO_train.p KO ${OutDir}/
Use this step for finalising training/validation/testing data with matching positive/negative labels.
python3 train_models/site-level/assemble_data_site-level.py --input_positive ${OutDir}/WT_added --input_negative ${OutDir}/KO -o $AssembledDataPrefix
Use this step for training binary classifier of site modification states.
python3 train_models/site-level/train_site-level.py -i $AssembledDataPrefix -o $OutDir --vector_len 100 --features 1 --labels 1 --arch Mini