This link provides the tool to create strain specific reference genome for effective mapping of the reads.

http://alleleseq.gersteinlab.org/tools.html

Download personnel genome constructor. This script will take vcf file and reference genome as input and will create enhanced reference genome as output. As this script is meant for human, two versions including maternal and paternal for each chromosome will be created. But as reference genome of mouse is an inbred organism, it has the same allele on both the chromosomes. so you can consider any copy either paternal or maternal. But make sure your vcf file has the homozygous variant tagged as 1/1 for GT tag and you should remove any heterozygous variants if any in the vcf file to skip the unnecessary complications.

I also have this script that I wrote a long ago. If you know very basic python you should be able to modify it to work for your case. It exactly does what you need. Here is the link: https://github.com/ashutoshkpandey/SimplePrograms/blob/master/psuedogenome_SNPs.py

Loading a full chromosome as an array is never a good idea, you can load the full chromosome as a "string" and modify the positions with "substr", something like this:

#!/usr/bin/perl

my %pos = loadPositions($file); # returns a hash like this $pos{'chr1'}{121726} = 'G'
my %seq = loadSequences($fasta); # returns a hash like this $seq{'chr1'} = 'ATCGATCATCTACTA....';

while (my ($chr, $seq) = each %seq) {
    foreach my $pos (keys %{ $pos{$chr} }) {
        substr($seq, $pos - 1, 1) = $pos{$chr}{$pos}; # be careful with 0-based and 1-based coordinates
    }
    printFasta($seq);
}

In this case, I'm loading the full genome, each chromosome as a "string", but you can optimize the code to do it chromosome by chromosome.

Anyone still coming to this page, if you have a fasta reference genome and a bam file that you want to turn into the reference file by changing SNP's and N's, you may try this one-liner using samtools, bcftools and vcfutils.pl (ps for beginners: both samtools and bcftools can be compiled in a computing cluster or in Linux, if so just add the locations of each before the software name; vcfutils is already a perl script from bcftools)

samtools mpileup -d8000 -q 20 -Q 10 -uf  REFERENCE.fasta Your_File.bam | bcftools call -c | vcfutils.pl vcf2fq  > OUTPUT.fastq

-d, --max-depth == -q, -min-MQ Minimum mapping quality for an alignment to be used == -Q, --min-BQ Minimum base quality for a base to be considered == (You can use different values of course, see http://www.htslib.org/doc/samtools.html)

Which generates a weird format that looks like fastq but isn't, so you can't convert it using a converter, but you can use the following sed command, which I wrote specific for this output:

sed -i -e '/^+/,/^\@/{/^+/!{/^\@/!d}}; /^+/ d; s/@/>/g' OUTPUT.fastq

In the end, make sure to compare your new fasta files to your reference to be sure that everything is fine.

EDIT BE CAREFUL WITH THE SED COMMAND IT MAY DELETE SOME OF YOUR READS IN DIFFERENT CASES OF QUALITY SCORING THAN I HAD