code: def consensus(args):
      
      p = OptionParser(consensus.__doc__)
      p.add_option("--fasta", default=False, action="store_true",
              help="Generate consensus FASTA sequences [default: %default]")
      p.add_option("--mask", default=0, type="int",
              help="Mask bases with quality lower than")
      opts, args = p.parse_args(args)
  
      if len(args) < 2:
          sys.exit(not p.print_help())
  
      fastafile, bamfile = args
      fasta = opts.fasta
      suffix = "fasta" if fasta else "fastq"
      pf = bamfile.rsplit(".", 1)[0]
      cnsfile = pf + ".cns.{0}".format(suffix)
      vcfgzfile = pf + ".vcf.gz"
      vcf([fastafile, bamfile, "-o", vcfgzfile])
      cmd += "zcat {0} | vcfutils.pl vcf2fq".format(vcfgzfile)
      if fasta:
          cmd += " | seqtk seq -q {0} -A -".format(opts.mask)
  
      sh(cmd, outfile=cnsfile)
docs: %prog consensus fastafile bamfile
  
      Convert bam alignments to consensus FASTQ/FASTA.