import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

# Load the coverage data
positive_strand = np.load('positive_strand_coverage.npy')
negative_strand = np.load('negative_strand_coverage.npy') * -1  # Invert negative strand values

# Start-end coordinates
cus_start, cus_end = 0, 10000

# ---------- Prep x-axis and subset ----------
x_start, x_end = cus_start, cus_end
x = np.arange(x_start, x_end)

# Subset the data
pos_subset = positive_strand[x_start:x_end]
neg_subset = negative_strand[x_start:x_end]

# Calculate global maximum for rectangle height
global_max_coverage = max(pos_subset.max(), abs(neg_subset.min())) * 1.05

# ---------- Create the plot ----------
plt.figure(figsize=(12, 6))

# Plot positive strand (red) and negative strand (blue)
plt.plot(x, pos_subset, color='red', linewidth=1, label='Positive strand')
plt.plot(x, neg_subset, color='blue', linewidth=1, label='Negative strand')

# Fill between the lines and zero for visual effect
plt.fill_between(x, pos_subset, 0, color='red', alpha=0.3)
plt.fill_between(x, neg_subset, 0, color='blue', alpha=0.3)

regions = [
    #{'start': 1301182, 'end': 1302813, 'color': 'lightgreen',   'label': 'oppA'},
    #{'start': 1302899, 'end': 1303819, 'color': 'lightblue',    'label': 'oppB'},
    #{'start': 1301182, 'end': 1306768, 'color': "khaki",        'label': 'oppA-oppF'},
    #{'start': 1303834, 'end': 1304742, 'color': 'orange',         'label': 'oppC'},
    #{'start': 1300856, 'end': 1300915, 'color': 'blue',    'label': 'oppAp2'},
    #{'start': 1301025, 'end': 1301084, 'color': 'plum',    'label': 'oppAp'},
]

ax = plt.gca()
for r in regions:
    # Positive strand rectangle
    ax.add_patch(
        Rectangle(
            (r['start'], 0),
            r['end'] - r['start'],
            global_max_coverage,
            facecolor=r['color'],
            alpha=0.75,
            label=r['label']
        )
    )
    # Negative strand rectangle
    ax.add_patch(
        Rectangle(
            (r['start'], -global_max_coverage),
            r['end'] - r['start'],
            global_max_coverage,
            facecolor=r['color'],
            alpha=0.75
        )
    )

# ---------- Adjust plot cosmetics ----------
plt.ylim(bottom=-global_max_coverage, top=global_max_coverage)
plt.axhline(0, color='black', linewidth=0.5)  # Add zero line

plt.ticklabel_format(style='plain', axis='x', useOffset=False)
plt.xticks(np.arange(x_start, x_end+1, 500000), rotation=45)
plt.xlabel('Position')
plt.ylabel('Coverage (positive and negative strands)')
plt.title('Strand-specific coverage across genetic elements')

# Deduplicate legend entries
handles, labels = ax.get_legend_handles_labels()
unique_labels = dict(zip(labels, handles))
plt.legend(unique_labels.values(), unique_labels.keys(), loc='upper right')

plt.tight_layout()
plt.show()


Okay, I have positive_strand_coverage1-8 files with corresponding negative_strand_coverage1-8 files. Remove fill below lines. Plot for each type of line negative and positive in the same color. 1 and 2 plot in shades of blue, call them IVT, 3 green-blue, call IVT DNA 10x. 4+5+6 Plot shades of red, call "No polymerase", 6 "No polymerase 10x DNA". 7+8 plot shades of green, call them "DNA"