1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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:

https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl

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 samtools split:

Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>]
                      [-f <format_string>] [-v] <merged.bam>
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.