get SNPs from draft genomes

Question

[anubhav9]

Hi, I wanted to compare a draft genome of a bacterial strain to its reference genome. The idea is to see what SNPs or changes have occurred in the new strain compared to its parent. The strain was sequenced using Oxford Nanopore long read sequencing and de novo assembled using Flye and annotated using bakta.

In the past, I have used breseq for this kind of work as it gives a list of all changes in a nice format. However, it requires short reads to align to the reference genome. Hence, I was looking for any tool which will take either the fasta or genbank files and output SNPs and other changes in a tabular format. It would be nice to also know which SNPs lie inside coding regions.

Any help will be appreciated!

Regards, Anubhav

Answer 1

[Charles Yin]

You may perform sequence alignments of the reference genome and new genome sequences, then parse the SNPs from the alignment files.

Reference paper (my paper): https://www.sciencedirect.com/science/article/pii/S0888754320303189

[anubhav9]

Hi, thanks a lot for your answer! I think I can use minimap2 to align the draft genome to the reference genome to get a sam file. But how do I parse the SNPs from that file? I am sure there are a million ways to do this but I have never done it. Could you please suggest a tools and which commands to use for this?

Thanks in advance! :)

[Charles Yin]

Sorry I do not have a tool to parse the SNPs from alignments.

Answer 2

[colindaven]

For reference based - so reads vs genome, which is not really what you are trying to do

Samtools sort, index -> BAMfile -> longshot or clair3

For your use case I would use PGGB to create a pangenome of the two, and this will give you a VCF of the differences between the two.

Other options are Freebayes or Samtools on a BAM file giving out all differences (coverage will only be 1, which is not typical).

Further options - mummer