Make a pretty FSC plot

plotFSCwrp.py

import starfile
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np

# Read the data
warp_apo_df = starfile.read("warp_apo/apoferritin_fsc.star")
pyco_apo_df = starfile.read("pyco_apo/apoferritin_fsc.star")

# Create the plot
plt.figure(figsize=(12, 7))
sns.lineplot(y=warp_apo_df['wrpFSCMasked'], x=1/warp_apo_df['wrpResolution'], label='TM', color='#de62df')
sns.lineplot(y=pyco_apo_df['wrpFSCMasked'], x=1/pyco_apo_df['wrpResolution'], label='tomoCPT', color='#32d0e9')

# Add horizontal line at y=0.143
plt.axhline(y=0.143, color='#717172', linestyle='-', label='FSC=0.143')

# Find first intersection point where FSC drops below threshold
def find_first_intersection(x, y, threshold):
    for i in range(len(y) - 1):
        if y[i] >= threshold and y[i+1] < threshold:
            # Linear interpolation to find more precise intersection
            x_intersect = x[i] + (threshold - y[i]) * (x[i+1] - x[i]) / (y[i+1] - y[i])
            return x_intersect, 1/x_intersect
    return None, None  # Return None if no intersection found

warp_intersection, warp_resolution = find_first_intersection(1/warp_apo_df['wrpResolution'], warp_apo_df['wrpFSCMasked'], 0.143)
pyco_intersection, pyco_resolution = find_first_intersection(1/pyco_apo_df['wrpResolution'], pyco_apo_df['wrpFSCMasked'], 0.143)

# Add vertical lines at intersection points
if warp_intersection:
    plt.axvline(x=warp_intersection, color='#de62df', linestyle=':', label='TM intersection')
    plt.text(warp_intersection, 0.15, f'{warp_resolution:.2f}Å', 
             color='#de62df', ha='center', va='bottom', rotation=30)

if pyco_intersection:
    plt.axvline(x=pyco_intersection, color='#32d0e9', linestyle=':', label='tomoCPT intersection')
    plt.text(pyco_intersection, 0.15, f'{pyco_resolution:.2f}Å', 
             color='#32d0e9', ha='center', va='bottom', rotation=30)

# Customize the plot
plt.xlabel('Resolution (1/Å)')
plt.ylabel('FSC')
plt.title('FSC Curves for TM and tomoCPT')
plt.legend()

# Show the plot
plt.savefig("fsc.pdf", format="pdf")