I ended up using samtools and bcftools like so:
samtools mpileup -uf human_g1k_v37.fasta sample.bam --region MT | bcftools call -m --ploidy 1 > sample_mt_only.vcf
This puts the entire MT in vcf. If I would supply "-v" parameter to bcftools, it would only store the variants.
samtools --version samtools 1.7 bcftools --version bcftools 1.7
We published a caller for mitochondrial ancient DNA but it is suitable for mitochondrial modern DNA as well. For mitochondrial modern DNA the problem is even simpler. The name of the program is endoCaller part of the (schmutzi package)
You can call it as such:
endoCaller -seq consensus.fa -log consensus.log -name [name sample here] chrMT.fa inputMT.bam
You might want to use the option -qual, this will filter bases according to a PHRED scale. e.g. -qual 50 means the confidence will be 1/100,000.