Produce .vcf file of ALL mutations in .bam file

Question

I got tasked with a problem, that i thought would be quite simple to solve, but turned out to be quite tricky. Our lab is running targeted mutagenesis experiments in yeast using crispr base editors. What we have done is set up an experiment where the base editor selectively mutates a region inside of a gene of interest in a culture of yeast cells. We then extract the DNA of the whole culture and perform sequencing of our gene of interest (the size of which is around 1kb). What i would like to do is plot the mutation frequency of say, C -> G edits, across the length of the entire gene (around 300bp), with the hope that i see a spike in mutations at the site where the base editor binds. Or at least, higher mutation rates at this site compared to the background mutation rate.

These edits are quite rare (we think), and since it's in a culture of multiple cells, they don't occur at the exact same position on the gene.

I initially thought i could do this by quality filtering the reads, aligning them to the gene sequence in order to produce a .bam file, and then piping this file through variant calling algorithms such as vcftools. However, i noticed that i get very few variants out of this, even when using extremely relaxed settings (e.g. p-value < 1 for the variant call).

Im therefore wondering if anyone knows how to extract from a .bam file ALL the mutations across all reads that do not match the reference sequence?

I'm aware that this will probably lead to quite a huge .vcf file, but since the reference sequence is only 300bp, it should still be manageable? It could also be the the experiment does not work, but i would at least expect some mutations in the sequencing reads compared to the reference gene sequence.

Any help would be greatly appreciated!

Answer 1

Ensure that the variant caller stringency is relaxed all the way. For example in freebayes use --min-alternate-count 1 --min-alternate-fraction 0 .

Answer 2

You could perform an mpileup with samtools. I think this would be the most sensitive approach. Here's the man page: http://www.htslib.org/doc/samtools-mpileup.html

You might want to disable BAQ correction. If you specify a reference sequence, then any letter in the fifth column would represent a non reference base. Reference bases are coded as . or ,.

The reason why vcftools doesn't work is probably because it assumes a diploid individual, which is not what you have. You have lots of samples mixed together. Freebayes allows relaxation of the diploid assumption. Here is the relevant example from the documentation https://github.com/ekg/freebayes :

Generate frequency-based calls for all variants passing input thresholds. You'd do this in the case that you didn't know the number of samples in the pool.

freebayes -f ref.fa -F 0.01 -C 1 --pooled-continuous aln.bam >var.vcf

You could try a cancer/normal variant caller like mutect, but I don't think you'll have much luck.

The most straightforward approach is to use samtools mpileup and to process the output yourself. You could also use bcftools mpileup to get a VCF, and then extract the allele count from the INFO field. If you are a Python lover, then you could also use pysam to both perform the pileup and process the results.