Plotting a ROC curve

[ACEnglish]

If your call data has quality scores, you can plot a ROC curve using truvari results. After running truvari bench and truvari vcf2df, run the following to make your ROC curve. I'll be using BioGraph results here.

import joblib
import seaborn as sb
import matplotlib.pyplot as plt
import sklearn.metrics as metrics
sb.set()

data = joblib.load("data.jl")
# Mark the subsets
tps = data["state"] == 'tp'
fps = data["state"] == 'fp'
dels = data["svtype"] == 'DEL'
ins = data["svtype"] == 'INS'
# Collect the subsets
all_data = data[tps | fps]
del_data = data[(tps | fps) & dels]
ins_data = data[(tps | fps) & ins]

# Plotting the total curve
fpr, tpr, threshold = metrics.roc_curve(all_data["state"] == 'tp', all_data["qual"])
auc = metrics.roc_auc_score(all_data["state"] == 'tp', all_data["qual"])
plt.plot(fpr, tpr, label="TOT auc=%.2f" % (auc))

# Plotting the DEL curve
fpr, tpr, threshold = metrics.roc_curve(del_data["state"] == 'tp', del_data["qual"])
auc = metrics.roc_auc_score(del_data["state"] == 'tp', del_data["qual"])
plt.plot(fpr, tpr, label="DEL auc=%.2f" % (auc))

# Plotting the INS curve
fpr, tpr, threshold = metrics.roc_curve(ins_data["state"] == 'tp', ins_data["qual"])
auc = metrics.roc_auc_score(ins_data["state"] == 'tp', ins_data["qual"])
plt.plot(fpr, tpr, label="INS auc=%.2f" % (auc))

plt.legend(loc=4)
plt.title("Truvari ROC curves for QUAL scores")
plt.show()

That flat tail on the curve indicates a quality score threshold could be applied to the calls to help boost the specificity. Let's look at the quality score distribution to see where that could happen.

p = sb.displot(data=data.reset_index(), x="qual", hue="state", hue_order=["tp", "fp"], multiple='stack')
p.set(title="QUAL distribution by state")
plt.show()

It seems that this result has SVs with QUAL==0 but FILTER==PASS. These calls should be removed. (Note that this turned out to be a bug in BioGraph and was fixed.)