Exact steps to calculate the SNP index to identify locality of EMS induced mutation

Question

[n.anuragsharma]

We carried out EMS mutagenesis of wild type plants (WT) and identified a mutant (X1). Mutant X1 was crossed back to WT and allowed to self fertilize in the F1 generation. 150 segregating F2 individuals were pooled and sequenced. The WT was also sequenced. So I now had 2 x R1+R2 fastq files corresponding to WT and X1-F2 bulked samples respectively which were then subjected to the following workflow:

1. Trimming of adaptors and quality control using fastp
2. Alignment to reference genome using bwa mem
3. Bam files processing, variant calling using bcftools mpileup and bcftools call

I now need to plot the SNP index for WT and X1-F2 bulked sample and subsequently their delta-SNP index. For the life of me, I am unable to figure out which tool to use and how to do it. So far, this paper is the best lead I have: "The SNP index was calculated for a locus as the ratio of non-reference reads to the total reads mapped to a variant locus."

Other papers:

1. Paper 1: "To help visualize the patterns of SNP-index, we applied a sliding-window analysis by taking the average value of SNP-index over 4 Mb window and incrementing the window 10 Kb at a time."
2. Paper 2:"If we define the SNP index as the ratio between the number of reads of a mutant SNP and the total number of reads corresponding to the SNP, we expect that this index would equal 1 near the causal gene and 0.5 for the unlinked loci. SNP indices can be scanned across the genome to find the region with a SNP index of 1, harboring the gene responsible for the mutant phenotype."

Can someone help me figure out which tool I can use to analyse this? I am open to any solution using bash or R, but can use Python if needed.