biostars 141739

example.py

from pyfaidx import FastaVariant

consensus = FastaVariant('reference.fasta', 'sample1.vcf.gz', het=True, hom=True)
chrom = 'chr1'

seq = consensus[chrom][0:8]
print(seq)  # AGTGCG

# if you don't want to invariant sites masked, you're good to go. otherwise:

masked = []
for base in seq:
  if base.start in consensus.variant_sites:
   masked.append(base.seq)
  else:
   masked.append('N')
 
print(''.join(masked))  # ANTGCN