docs: def plot_atacseq_insert_sizes(self, bam, plot, output_csv, max_insert=1500, smallest_insert=30):
          """
          Heavy inspiration from here:
          https://github.com/dbrg77/ATAC/blob/master/ATAC_seq_read_length_curve_fitting.ipynb
          """
          try:
              import pysam
              import numpy as np
              import matplotlib.mlab as mlab
              from scipy.optimize import curve_fit
              from scipy.integrate import simps
              import matplotlib
              matplotlib.use('Agg')
              import matplotlib.pyplot as plt
          except:
              print("Necessary Python modules couldn't be loaded.")
              return
  
          try:
              import seaborn as sns
              sns.set_style("whitegrid")
          except:
              pass
  
          def get_fragment_sizes(bam, max_insert=1500):
              frag_sizes = list()
  
              bam = pysam.Samfile(bam, 'rb')
  
              for i, read in enumerate(bam):
                  if read.tlen < max_insert:
                      frag_sizes.append(read.tlen)
              bam.close()
  
              return np.array(frag_sizes)
  
          def mixture_function(x, *p):
              """
              Mixture function to model four gaussian (nucleosomal)
              and one exponential (nucleosome-free) distributions.
              """
              m1, s1, w1, m2, s2, w2, m3, s3, w3, m4, s4, w4, q, r = p
              nfr = expo(x, 2.9e-02, 2.8e-02)
              nfr[:smallest_insert] = 0
  
              return (mlab.normpdf(x, m1, s1) * w1 +
                      mlab.normpdf(x, m2, s2) * w2 +
                      mlab.normpdf(x, m3, s3) * w3 +
                      mlab.normpdf(x, m4, s4) * w4 +
                      nfr)
  
          def expo(x, q, r):
              """
              Exponential function.
              """
              return q * np.exp(-r * x)
  
          # get fragment sizes
          frag_sizes = get_fragment_sizes(bam)
  
          # bin
          numBins = np.linspace(0, max_insert, max_insert + 1)
          y, scatter_x = np.histogram(frag_sizes, numBins, density=1)
          # get the mid-point of each bin
          x = (scatter_x[:-1] + scatter_x[1:]) / 2
  
          # Parameters are empirical, need to check
          paramGuess = [
              200, 50, 0.7,  # gaussians
              400, 50, 0.15,
              600, 50, 0.1,
              800, 55, 0.045,
              2.9e-02, 2.8e-02  # exponential
          ]
  
          try:
              popt3, pcov3 = curve_fit(
                  mixture_function, x[smallest_insert:], y[smallest_insert:],
                  p0=paramGuess, maxfev=100000)
          except:
              print("Nucleosomal fit could not be found.")
              return
  
          m1, s1, w1, m2, s2, w2, m3, s3, w3, m4, s4, w4, q, r = popt3
  
          # Plot
          plt.figure(figsize=(12, 12))
  
          # Plot distribution
          plt.hist(frag_sizes, numBins, histtype="step", ec="k", normed=1, alpha=0.5)
  
          # Plot nucleosomal fits
          plt.plot(x, mlab.normpdf(x, m1, s1) * w1, 'r-', lw=1.5, label="1st nucleosome")
          plt.plot(x, mlab.normpdf(x, m2, s2) * w2, 'g-', lw=1.5, label="2nd nucleosome")
          plt.plot(x, mlab.normpdf(x, m3, s3) * w3, 'b-', lw=1.5, label="3rd nucleosome")
          plt.plot(x, mlab.normpdf(x, m4, s4) * w4, 'c-', lw=1.5, label="4th nucleosome")
  
          # Plot nucleosome-free fit
          nfr = expo(x, 2.9e-02, 2.8e-02)
          nfr[:smallest_insert] = 0
          plt.plot(x, nfr, 'k-', lw=1.5, label="nucleosome-free")
  
          # Plot sum of fits
          ys = mixture_function(x, *popt3)
          plt.plot(x, ys, 'k--', lw=3.5, label="fit sum")
  
          plt.legend()
          plt.xlabel("Fragment size (bp)")
          plt.ylabel("Density")
          plt.savefig(plot, bbox_inches="tight")
  
          # Integrate curves and get areas under curve
          areas = [
              ["fraction", "area under curve", "max density"],
              ["Nucleosome-free fragments", simps(nfr), max(nfr)],
              ["1st nucleosome", simps(mlab.normpdf(x, m1, s1) * w1), max(mlab.normpdf(x, m1, s1) * w1)],
              ["2nd nucleosome", simps(mlab.normpdf(x, m2, s2) * w1), max(mlab.normpdf(x, m2, s2) * w2)],
              ["3rd nucleosome", simps(mlab.normpdf(x, m3, s3) * w1), max(mlab.normpdf(x, m3, s3) * w3)],
              ["4th nucleosome", simps(mlab.normpdf(x, m4, s4) * w1), max(mlab.normpdf(x, m4, s4) * w4)]
          ]
  
          try:
              import csv
  
              with open(output_csv, "w") as f:
                  writer = csv.writer(f)
                  writer.writerows(areas)
          except:
              pass