parseFreyjaBootstraps.py

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import sys
import pandas as pd
import numpy as np
from getDisplayName import *


if len(sys.argv) != 4:
    raise Exception('Incorrect call to the script.')


# Data regarding the current sample passed by UNIX
freyjaOutputFile = sys.argv[1]
freyjaBootFile = sys.argv[2]
outfilename = sys.argv[3]


# Process the bootstrap result table generated by Freyja
freyja_boot = pd.read_table(freyjaBootFile, index_col=0, header=0, sep=',')
boot_names = [ getDisplayName(x) for x in list(freyja_boot) ]

percentiles = list(freyja_boot.index)
index25 = percentiles.index(0.25)
index50 = percentiles.index(0.50)
index75 = percentiles.index(0.75)
percentile25 = list(freyja_boot.iloc[index25])
percentile50 = list(freyja_boot.iloc[index50])
percentile75 = list(freyja_boot.iloc[index75])


# Import the freyja file
(lineages, abundances, freyja_names) = import_freyja_demix(freyjaOutputFile)


freyja_hits = []
for dname in set(freyja_names):
    pct = 100 * sum([ y for (x,y) in zip(freyja_names, abundances) if x==dname ])
    median = 100 * sum([ y for (x,y) in zip(boot_names, percentile50) if x==dname ])
    spread = 100 * np.sqrt(sum([ (z-y)**2 for (x,y,z) in zip(boot_names, percentile25, percentile75) if x==dname ]))   
    freyja_hits.append( (dname, pct, median, spread) )


# Make a scatter plot, value vs. estimated spread
plt.rcParams.update({'font.size': 14})

for (dname, pct, median, spread) in freyja_hits:
    color2plot = getColor(dname)
    plt.plot(pct, median, 'o', color=color2plot, markersize=20)
    plt.errorbar(pct, median, yerr=3*spread, color='k')
    r = np.random.uniform(5,8)
    phi = np.random.uniform(0,np.pi)
    plt.text(pct+r*np.sin(phi), median+r*np.cos(phi), dname, color=color2plot)


plt.plot([0,100], [0,100], 'k--')
plt.axis('tight')
plt.xlim(0, 100)
plt.ylim(0, 100)
plt.xlabel('Freyja demix prediction')
plt.ylabel('Freyja boot median +/- 3 FWHM')
plt.title('Freyja bootstrapping (Experimental)')

plt.savefig(outfilename, dpi=300)
plt.close()