I have a fastq file from Nanopore sequencing data that contains reads from the UK and South Africa variants of SARS-CoV-2. These variants are identified by three key mutations in the receptor-binding domain. I want to determine the number of reads then belongs to UK and South Africa variant. I will create a bam file from the fastq file using the SARS-CoV-2 reference genome. Then, how do I determine the reads belonging to each variants? Are there any tools available to classify the reads based on the pre-defined mutation criteria?
Based on a quick request from someone else, I wrote some code a couple of years ago for assigning read groups based on variants at particular locations; this was set up for HIV variant splitting, so I expect it should also work on SARS-CoV-2 nanopore reads as well:
It processes a SAM stream [i.e. in text mode] using Perl, so might be a bit slow, but should be fine for a few hundred thousand sequences. You'd use it as follows:
samtools view -h mapped_reads.bam | \ samVarSplitter.pl [-ref <refName>] [-pos <int>] | \ samtools sort > grouped_reads.bam
Then the reads can be split into different files based on their group using
Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>] [-f <format_string>] [-v] <merged.bam>