$ samtools view input.bam | cut -f 5 | sort | uniq  -c | sort -n | awk '{printf("MAPQ:%s\t%d\n",$2,$1);}' | gnuplot -e " set terminal dumb ;set nokey; plot '-' using 2:xtic(1) with boxes"



  400 ++---+---------+---------+---------+--------+---------+----***********
      |    +         +         +         +        +         +    *    +    *
  350 ++                                                         *        +*
      |                                                          *         *
      |                                                          *         *
  300 ++                                                         *        +*
      |                                                          *         *
  250 ++                                                         *        +*
      |                                                          *         *
  200 ++                                                         *        +*
      |                                                          *         *
      |                                                          *         *
  150 ++                                                         *        +*
      |                                                          *         *
  100 ++                                                         *        +*
      |                                                          *         *
      |                                                          *         *
   50 ++                                               ***********        +*
      |    +         +         +         +   ***********    +    *    +    *
    0 **********************************************************************
        MAPQ:7    MAPQ:25   MAPQ:0    MAPQ:70  MAPQ:29   MAPQ:37   MAPQ:60

This may be too simple for what you need, but I have a script here that uses pysam to extract MapQ scores from a bam file and then generates a kernel density plot of the scores.

Here is the actual code:

#!/usr/bin/env python

import pysam
import seaborn as sns
import argparse
from matplotlib import pyplot as plt

parser = argparse.ArgumentParser(description="Plot mapping quality score from bam file")

parser.add_argument('--bam', required=True, help='path to input bam file', action='store')
parser.add_argument('--png', required=True, help='path to output png file', action='store')

args=parser.parse_args()

output = args.png
input = args.bam 

bam = pysam.AlignmentFile(input, "rb")

quals = [read.mapping_quality for read in bam.fetch()]

sns.set_style("dark")
sns.set_context("paper")

plt.figure(figsize=(10,6))

ax = sns.kdeplot(quals, shade = True)
ax.set_ylabel('Density')
ax.set_xlabel('Mapping Quality Score')

plt.savefig(output)

I'm pretty sure that Qualimap BamQC (not to be confused with the BamQC tool from Babraham mentioned above) calculates and plots this. Qualimap is supported by MultiQC, however off the top of my head, I don't think that MultiQC replicates this plot in reports. I doubt it would be particularly difficult to add if you fancy it (or opening a GitHub issue with some example data), or you could probably sidestep that by using MultiQC's "custom content" feature.