I am trying to use variant calling using the metrics of DP>=10, %QUAL>=30 and mapping MQ>30 and prepare a whole genome alignment. Everything was working well. However, I had an issue replacing "variants that did not pass filtering" with "N". Could anyone here help me please?

Here is the code:

Filtering variants:

bcftools mpileup -Ou -f 0140J.fasta  $sorted_bam_files | bcftools call -O b -cv --threads 24 --ploidy 1 -V indels | bcftools filter --threads 24 -i 'MIN(INFO/DP)>=10 & %QUAL>=30||MQ>30'> Suberis.vcf

Mask the bases that failed filtering parameter with “N”

bgzip -c Suberis.vcf > Suberis.vcf.gz
bcftools query -f'%CHROM\t%POS0\t%END\n'  Suberis.vcf.gz  > variants.bed
bedtools genomecov -bga -ibam $sorted_bam_files | awk '$4 < 10' > low_coverage_sites.bed

bedtools subtract -a low_coverage_sites.bed -b variants.bed > mask.bed
bcftools consensus -f 0140J.fasta -m mask.bed -i 'FORMAT/VAF > 0.90 & INFO/MQ >= 20 & FORMAT/DP >= 10' Suberis.vcf.gz > bcftools.consensus.fa

Note: I need to replace all the bases that failed the filtering metrics of DP>=10, QUAL>=30 and MQ>30, with the character “N”.