nanopore SNV callers returns empty VCF

Question

[Meghan.T]

We have performed error prone PCR on a segment of a gene ( 800bp) and selected colonies with phage display. To identify the mutations that improve the function of the gene we decided to perform nanopore sequencing as we are interested in identifying mutations patterns on the same molecule. I use minimap to align to genome and then tried variant calling with clair3, However the VC returns no variants ( there is no error just the vcf file is empty) I was wondering how could I resolve this issue ( use diifferent tools, different option? ) I would appreciate your opinion on this.

I have attached the files here and would appreciate your opinion. bam files

Answer 1

[Istvan Albert]

well it looks like almost no reads aligned, out of 13K reads you have only a single primary alignment

hence the variant callers are unable to produce variants

samtools flagstat D2.bam

prints:

    13056 + 0 in total (QC-passed reads + QC-failed reads)
    13051 + 0 primary
    5 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    0 + 0 primary duplicates
    6 + 0 mapped (0.05% : N/A)
    1 + 0 primary mapped (0.01% : N/A)
    0 + 0 paired in sequencing
    0 + 0 read1
    0 + 0 read2
    0 + 0 properly paired (N/A : N/A)
    0 + 0 with itself and mate mapped
    0 + 0 singletons (N/A : N/A)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)

[Meghan.T]

Thank you so much for your response. I forgot to check the alignment statistics. I assume there were too many mutations during the process, I tried realiging with more lenient options (match =10, gap = -4, mismatch = -2) but still no read was mapped. Do you have any suggestions ?

[Istvan Albert]

it is very unlikely that the number of mutations caused the reads to fail to align. The number of mutations needed to cause that would be very high - say above 20-30% and when using long reads you would still have more regions that align.

You should not need to be more "lenient", beyond selecting the nanopore specific setting in minimap, which you have already done.

When reads do not align it typically means either that the sequenced DNA does not match the reference genome - for example there are other off targets or contaminants, or host related sequences that dominate.

or that the sequencing has failed.

[Meghan.T]

Thank you so much for your comprehensive response and tips. I reached out to the person who designed the experiment and turns out that the sequence was codon optimized for expression in yeast so I guess that's why alignment failed. After aligning to the modified sequence of gene, alignment performance was almost perfect.