RosaBioinfo

Coverage visualization of alignment of Isolate reads on wPipPel and WO regions

January 16, 2026

Alice Brunner

Coverage visualization of alignment of Isolate reads on wPipPel and WO regions

Language: Python

#!pip3 install pandas openpyxl

import pandas as pd # 1.5.3 # numpy backend 1.26.0
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import FancyArrowPatch, Patch
from Bio import SeqIO
import os 
import re 

All wPipPel regions

The following input is required: Data/Coverage_Isolate_against_Pel.txt.

import numpy as np

def smooth_region(data, deb, end):
    """
    Smooths the region between `deb` and `end` in `data` by replacing it with the average of the values at `deb-10` and `end+10`.

    Args:
        data (np.array): Input data array.
        deb (int): Start index of the region to smooth.
        end (int): End index of the region to smooth.

    Returns:
        np.array: Data with the specified region smoothed.
    """
    # Calculate the average of the values at deb-10 and end+10
    left_val = data[max(0, deb - 10)]
    right_val = data[min(len(data) - 1, end + 10)]
    avg_val = (left_val + right_val) / 2

    # Replace the region between deb and end with the average value
    smoothed_data = data.copy()
    smoothed_data[deb:end + 1] = avg_val

    return smoothed_data

# Output path for the coverage figure
output_path = "Output/Cov_Total.pdf"

# Load per-base coverage file from read mapping to wPipPel
csv_path = "Data/Coverage_Isolate_against_Pel.txt"

df = pd.read_csv(
    csv_path,
    sep="\t",
    header=0,
    names=["contig", "position", "coverage"]
)

# Ensure numeric values and remove invalid entries
df["position"] = pd.to_numeric(df["position"], errors="coerce")
df["coverage"] = pd.to_numeric(df["coverage"], errors="coerce")
df = df.dropna()

# Smooth individual peaks of over-copied ribosomal-associated DNA for better visibility of coverage profiles
df["coverage"] = smooth_region(df["coverage"].values, 1135994, 1137497) # Smooth 16S
df["coverage"] = smooth_region(df["coverage"].values, 1235938, 1238876) # Smooth 23S
df["coverage"] = smooth_region(df["coverage"].values, 771502, 771929) # 50S Protein-Binding related sequence

# Downsample positions to reduce figure size and file weight
df = df.iloc[::10, :].copy()

# Generate coverage plot
fig, ax = plt.subplots(figsize=(12, 2))

ax.plot(df["position"], df["coverage"], linewidth=0.7)
ax.fill_between(df["position"], 0, df["coverage"], alpha=0.2)

# Highlight WOPip prophage regions along the genome
regions = {
    "WOPip\n1": (246932, 297466, "green"),
    "WOPip\n2": (297282, 346089, "orange"),
    "WOPip\n3": (346156, 375271, "blue"),
    "WOPip4": (429750, 502566, "yellow"),
    "WOPip5": (1369259, 1446799, "purple"),
}

for name, (start, end, color) in regions.items():
    ax.axvspan(start, end, color=color, alpha=0.2)
    ax.text(
        (start + end) / 2,
        max(df["coverage"]) * 0.95,
        name,
        ha="center",
        va="top",
        fontsize=6,
        fontweight="bold"
    )
    
# Axis formatting
ax.set_xlim(df["position"].min(), df["position"].max())
ax.set_xlabel("Position (bp)", fontsize=10)
ax.set_ylabel("Coverage", fontsize=10)
ax.set_title("Coverage of isolate reads mapped to wPipPel", fontsize=11)

# Aesthetic adjustments
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.grid(True, linestyle="--", alpha=0.2)

plt.tight_layout()

# Export figure in both raster-safe PDF and lightweight vector SVG formats
fig.savefig(output_path, bbox_inches="tight")
fig.savefig(output_path.replace(".pdf", ".svg"), bbox_inches="tight")

plt.show()
print("Figure saved to:", output_path)

png

Figure saved to: Output/Cov_Total.pdf

# Coverage WOPip

all_cov = pd.read_csv('Data/Coverage_Isolate_against_Pel.txt', sep='\t', header=None)

wopip1 = all_cov.iloc[246933-1:297466-1]
wopip2 = all_cov.iloc[297283-1:346089-1]
wopip3 = all_cov.iloc[346157-1:375271-1]
wopip4 = all_cov.iloc[429750-1:502566-1]
wopip5 = all_cov.iloc[1369259-1:1446799-1]

wopip1.to_csv('Data/WOPip1.csv', index=None, sep='\t')
wopip2.to_csv('Data/WOPip2.csv', index=None, sep='\t')
wopip3.to_csv('Data/WOPip3.csv', index=None, sep='\t')
wopip4.to_csv('Data/WOPip4.csv', index=None, sep='\t')
wopip5.to_csv('Data/WOPip5.csv', index=None, sep='\t')

WOPip1 coverage

The following inputs are required:

Data/WOPip1_gff.xlsx - Gene calling on WOPip1.fasta
Data/Orthogroups_mod_Pip1.tsv - Orthofinder launched against Wolbachia genes database launched on WOPip1 CDS
Data/WOPip1.csv - Coverage of Isolate reads on WOPip1 region of wPipPel

# WOPip1: gene synteny (functional categories) + isolate coverage track

# === Inputs / outputs ===
gff_path = "Data/WOPip1_gff.xlsx"
orthogroups_path = "Data/Orthogroups_mod_Pip1.tsv"
coverage_folder = "Data/WOPip1.csv"
output_path = "Output/WOPip1_synt.pdf"

# === Load gene coordinates and orthogroup-derived annotations ===
df = pd.read_excel(gff_path)
orthogroups_df = pd.read_csv(orthogroups_path, sep="\t")

df.columns = [col.strip().lower() for col in df.columns]
df["protein_id"] = df["proteines"].str.extract(r"___(\d+)$")[0]

# Map protein IDs to functional labels (from orthogroup annotation)
protein_function_map = {}
for _, row in orthogroups_df.iterrows():
    if pd.notna(row["PROTEINES_WOPip1"]):
        for pid in row["PROTEINES_WOPip1"].split(","):
            protein_function_map[pid.strip()] = str(row["WOqueries_041124"])

# === Functional category color mapping ===
def get_color(annotations):
    """Assign colors based on the primary functional label; prioritize cif annotations when present."""
    if annotations is None:
        return "grey"
    annotations = annotations.split(",")
    annotation = annotations[0]
    for a in annotations:
        if "cif" in a:
            annotation = a
    a = annotation

    return {
        "HEAD": "#660066",
        "ANK": "#ff0000",
        "CON_": "#0432ff",
        "TailFiber": "#ffc1c0",
        "Tail_": "#ff6666",
        "RECOMB": "#07f256",
        "LIGA": "#712ca3",
        "TOXIN": "#de8209",
        "WOPC1": "#4aeaed",
        "cifA": "#ffa4ff",
        "cifB": "#ffa4ff",
        "TR": "#7340f5",
        "PATATIN": "#593601",
        "TNP": "#fcfc2d",
        "MutL": "#272545",
        "METHYL": "#538201",
        "REP": "#88d4eb",
    }.get(
        next(
            (k for k in [
                "HEAD","ANK","CON_","TailFiber","Tail_","RECOMB","LIGA","TOXIN","WOPC1","cifA","cifB",
                "TR","PATATIN","TNP","MutL","METHYL","REP"
            ] if k in a),
            None
        ),
        "#c7c5c5"
    )

# === Figure layout: synteny track (top) + coverage track (bottom) ===
fig = plt.figure(figsize=(20, 3))
y_height_syn = 0.6
y_height_cov = 0.3
y_margin = 0.05

# This dataset corresponds to a single contig (WOPip1 region)
contig = df["contigs"].unique()[0]
subset = df[df["contigs"] == contig]
contig_len = subset["fin"].max()

# --- Synteny track ---
ax_syn = fig.add_axes([0.1, y_height_cov + y_margin, 0.8, y_height_syn])
ax_syn.set_xlim(0, contig_len)
ax_syn.set_ylim(-1, 1)
ax_syn.axis("off")
ax_syn.text(-1000, 0, "WOPip1", ha="right", va="center", fontsize=10, fontweight="bold")
ax_syn.hlines(0, 0, contig_len, color="black", linewidth=2, zorder=0)

# Draw gene arrows colored by functional category
for _, row in subset.iterrows():
    start, end = row["deb"], row["fin"]
    direction = row["sens"]
    prot_id = str(row["protein_id"])
    annotation = protein_function_map.get(prot_id)
    color = get_color(annotation)

    arrow = (
        FancyArrowPatch((start, 0), (end, 0), arrowstyle="simple", color=color,
                        mutation_scale=15, linewidth=1.5)
        if direction == "+"
        else FancyArrowPatch((end, 0), (start, 0), arrowstyle="simple", color=color,
                             mutation_scale=15, linewidth=1.5)
    )
    ax_syn.add_patch(arrow)

# --- Coverage track ---
ax_cov = fig.add_axes([0.1, 0, 0.8, y_height_cov])

df_cov = pd.read_csv(coverage_folder, header=0, sep="\t", names=["contig", "position", "coverage"])
df_cov = df_cov.dropna()
df_cov["position"] = pd.to_numeric(df_cov["position"], errors="coerce")
df_cov["coverage"] = pd.to_numeric(df_cov["coverage"], errors="coerce")

# Re-index positions to start at 1 (coverage was computed on wPipPel genomic coordinates)
offset = df_cov["position"].min() - 1
df_cov["position"] = df_cov["position"] - offset

ax_cov.plot(df_cov["position"], df_cov["coverage"], color="blue", linewidth=0.7)
ax_cov.fill_between(df_cov["position"], 0, df_cov["coverage"], color="blue", alpha=0.2)
ax_cov.set_xlim(0, contig_len)
ax_cov.set_ylabel("Cov.", fontsize=6)
ax_cov.tick_params(axis="both", labelsize=6)
ax_cov.spines["top"].set_visible(False)
ax_cov.spines["right"].set_visible(False)

# === Legend ===
legend_elements = [
    mpatches.Patch(color="#660066", label="Head"),
    mpatches.Patch(color="#ff0000", label="Ankyrin"),
    mpatches.Patch(color="#0432ff", label="Connector/Baseplate"),
    mpatches.Patch(color="#ffc1c0", label="Tail Fiber"),
    mpatches.Patch(color="#ff6666", label="Tail"),
    mpatches.Patch(color="#ffa4ff", label="CifA/B"),
    mpatches.Patch(color="#07f256", label="Recombinase"),
    mpatches.Patch(color="#fcfc2d", label="Transposase"),
    mpatches.Patch(color="#4aeaed", label="WOPC1"),
    mpatches.Patch(color="#ff6600", label="Toxins"),
    mpatches.Patch(color="#712ca3", label="LigA"),
    mpatches.Patch(color="#7340f5", label="TR"),
    mpatches.Patch(color="#593601", label="Patatin"),
    mpatches.Patch(color="#272545", label="MutL"),
    mpatches.Patch(color="#88d4eb", label="Replication/Repair"),
    mpatches.Patch(color="#538201", label="Methylase"),
]
fig.legend(
    handles=legend_elements,
    loc="center left",
    bbox_to_anchor=(0.91, 0.5),
    title="Legend",
    fontsize=8,
    title_fontsize=9
)

# === Export ===
fig.savefig(output_path, bbox_inches="tight")
plt.show(fig)
print("Figure saved to:", output_path)

png

Figure saved to: Output/WOPip1_synt.pdf

WOPip2 coverage

The following inputs are required:

Data/WOPip2_gff.xlsx - Gene calling on WOPip2.fasta
Data/Orthogroups_mod_Pip2.tsv - Orthofinder launched against Wolbachia genes database launched on WOPip2 CDS
Data/WOPip2.csv - Coverage of Isolate reads on WOPip2 region of wPipPel

# WOPip2: gene synteny (functional categories) + isolate coverage track

# === Inputs / outputs ===
gff_path = "Data/WOPip2_gff.xlsx"
orthogroups_path = "Data/Orthogroups_mod_Pip2.tsv"
coverage_folder = "Data/WOPip2.csv"
output_path = "Output/WOPip2_synt.pdf"

# === Load gene coordinates and orthogroup-derived annotations ===
df = pd.read_excel(gff_path)
orthogroups_df = pd.read_csv(orthogroups_path, sep="\t")

df.columns = [col.strip().lower() for col in df.columns]
df["protein_id"] = df["proteines"].str.extract(r"___(\d+)$")[0]

# Map protein IDs to functional labels (from orthogroup annotation)
protein_function_map = {}
for _, row in orthogroups_df.iterrows():
    if pd.notna(row["PROTEINES_WOPip2"]):
        for pid in row["PROTEINES_WOPip2"].split(","):
            protein_function_map[pid.strip()] = str(row["WOqueries_041124"])

# === Functional category color mapping ===
def get_color(annotations):
    """Assign colors based on the primary functional label; prioritize cif annotations when present."""
    if annotations is None:
        return "grey"
    annotations = annotations.split(",")
    annotation = annotations[0]
    for a in annotations:
        if "cif" in a:
            annotation = a
    a = annotation

    return {
        "HEAD": "#660066",
        "ANK": "#ff0000",
        "CON_": "#0432ff",
        "TailFiber": "#ffc1c0",
        "Tail_": "#ff6666",
        "RECOMB": "#07f256",
        "LIGA": "#712ca3",
        "TOXIN": "#de8209",
        "WOPC1": "#4aeaed",
        "cifA": "#ffa4ff",
        "cifB_": "#ffa4ff",
        "TR": "#7340f5",
        "PATATIN": "#593601",
        "TNP": "#fcfc2d",
        "MutL": "#272545",
        "METHYL": "#538201",
        "REP": "#88d4eb",
    }.get(
        next(
            (k for k in [
                "HEAD","ANK","CON_","TailFiber","Tail_","RECOMB","LIGA","TOXIN","WOPC1","cifA","cifB_",
                "TR","PATATIN","TNP","MutL","METHYL","REP"
            ] if k in a),
            None
        ),
        "#c7c5c5"
    )

# === Figure layout: synteny track (top) + coverage track (bottom) ===
fig = plt.figure(figsize=(20, 3))
y_height_syn = 0.6
y_height_cov = 0.3
y_margin = 0.05

# This dataset corresponds to a single contig (WOPip2 region)
contig = df["contigs"].unique()[0]
subset = df[df["contigs"] == contig]
contig_len = subset["fin"].max()

# --- Synteny track ---
ax_syn = fig.add_axes([0.1, y_height_cov + y_margin, 0.8, y_height_syn])
ax_syn.set_xlim(0, contig_len)
ax_syn.set_ylim(-1, 1)
ax_syn.axis("off")
ax_syn.text(-1000, 0, "WOPip2", ha="right", va="center", fontsize=10, fontweight="bold")
ax_syn.hlines(0, 0, contig_len, color="black", linewidth=2, zorder=0)

# Draw gene arrows colored by functional category
for _, row in subset.iterrows():
    start, end = row["deb"], row["fin"]
    direction = row["sens"]
    prot_id = str(row["protein_id"])
    annotation = protein_function_map.get(prot_id)
    color = get_color(annotation)

    arrow = (
        FancyArrowPatch((start, 0), (end, 0), arrowstyle="simple", color=color,
                        mutation_scale=15, linewidth=1.5)
        if direction == "+"
        else FancyArrowPatch((end, 0), (start, 0), arrowstyle="simple", color=color,
                             mutation_scale=15, linewidth=1.5)
    )
    ax_syn.add_patch(arrow)

# --- Coverage track ---
ax_cov = fig.add_axes([0.1, 0, 0.8, y_height_cov])

df_cov = pd.read_csv(coverage_folder, header=0, sep="\t", names=["contig", "position", "coverage"])
df_cov = df_cov.dropna()
df_cov["position"] = pd.to_numeric(df_cov["position"], errors="coerce")
df_cov["coverage"] = pd.to_numeric(df_cov["coverage"], errors="coerce")

# Re-index positions to start at 1 (coverage was computed on wPipPel genomic coordinates)

offset = df_cov["position"].min() - 1
df_cov["position"] = (len(df_cov['position']) - (df_cov["position"] - offset) + 1) # reverse region


ax_cov.plot(df_cov["position"], df_cov["coverage"], color='blue', linewidth=0.7)
ax_cov.fill_between(df_cov["position"], 0, df_cov["coverage"], color='blue', alpha=0.2)
ax_cov.set_xlim(0, contig_len)
ax_cov.set_ylabel("Cov.", fontsize=6)
ax_cov.tick_params(axis='both', labelsize=6)
ax_cov.spines['top'].set_visible(False)
ax_cov.spines['right'].set_visible(False)


# === Legend ===
legend_elements = [
    mpatches.Patch(color="#660066", label="Head"),
    mpatches.Patch(color="#ff0000", label="Ankyrin"),
    mpatches.Patch(color="#0432ff", label="Connector/Baseplate"),
    mpatches.Patch(color="#ffc1c0", label="Tail Fiber"),
    mpatches.Patch(color="#ff6666", label="Tail"),
    mpatches.Patch(color="#ffa4ff", label="CifA/B"),
    mpatches.Patch(color="#07f256", label="Recombinase"),
    mpatches.Patch(color="#fcfc2d", label="Transposase"),
    mpatches.Patch(color="#4aeaed", label="WOPC1"),
    mpatches.Patch(color="#ff6600", label="Toxins"),
    mpatches.Patch(color="#712ca3", label="LigA"),
    mpatches.Patch(color="#7340f5", label="TR"),
    mpatches.Patch(color="#593601", label="Patatin"),
    mpatches.Patch(color="#272545", label="MutL"),
    mpatches.Patch(color="#88d4eb", label="Replication/Repair"),
    mpatches.Patch(color="#538201", label="Methylase"),
]
fig.legend(
    handles=legend_elements,
    loc="center left",
    bbox_to_anchor=(0.91, 0.5),
    title="Legend",
    fontsize=8,
    title_fontsize=9
)

# === Export ===
fig.savefig(output_path, bbox_inches="tight")
plt.show(fig)
print("Figure saved to:", output_path)

png

Figure saved to: Output/WOPip2_synt.pdf

WOPip3 coverage

The following inputs are required:

Data/WOPip3_gff.xlsx - Gene calling on WOPip3.fasta
Data/Orthogroups_mod_Pip3.tsv - Orthofinder launched against Wolbachia genes database launched on WOPip3 CDS
Data/WOPip3.csv - Coverage of Isolate reads on WOPip3 region of wPipPel

# WOPip3: gene synteny (functional categories) + isolate coverage track

# === Inputs / outputs ===
gff_path = "Data/WOPip3_gff.xlsx"
orthogroups_path = "Data/Orthogroups_mod_Pip3.tsv"
coverage_folder = "Data/WOPip3.csv"
output_path = "Output/WOPip3_synt.pdf"

# === Load gene coordinates and orthogroup-derived annotations ===
df = pd.read_excel(gff_path)
orthogroups_df = pd.read_csv(orthogroups_path, sep="\t")

df.columns = [col.strip().lower() for col in df.columns]
df["protein_id"] = df["proteines"].str.extract(r"___(\d+)$")[0]

# Map protein IDs to functional labels (from orthogroup annotation)
protein_function_map = {}
for _, row in orthogroups_df.iterrows():
    if pd.notna(row["PROTEINES_WOPip3"]):
        for pid in row["PROTEINES_WOPip3"].split(","):
            protein_function_map[pid.strip()] = str(row["WOqueries_041124"])

# === Functional category color mapping ===
def get_color(annotations):
    """Assign colors based on the primary functional label; prioritize cif annotations when present."""
    if annotations is None:
        return "grey"
    annotations = annotations.split(",")
    annotation = annotations[0]
    for a in annotations:
        if "cif" in a:
            annotation = a
    a = annotation

    return {
        "HEAD": "#660066",
        "ANK": "#ff0000",
        "CON_": "#0432ff",
        "TailFiber": "#ffc1c0",
        "Tail_": "#ff6666",
        "RECOMB": "#07f256",
        "LIGA": "#712ca3",
        "TOXIN": "#de8209",
        "WOPC1": "#4aeaed",
        "cifA": "#ffa4ff",
        "cifB": "#ffa4ff",
        "TR": "#7340f5",
        "PATATIN": "#593601",
        "TNP": "#fcfc2d",
        "MutL": "#272545",
        "METHYL": "#538201",
        "REP": "#88d4eb",
    }.get(
        next(
            (k for k in [
                "HEAD","ANK","CON_","TailFiber","Tail_","RECOMB","LIGA","TOXIN","WOPC1","cifA","cifB",
                "TR","PATATIN","TNP","MutL","METHYL","REP"
            ] if k in a),
            None
        ),
        "#c7c5c5"
    )

# === Figure layout: synteny track (top) + coverage track (bottom) ===
fig = plt.figure(figsize=(20, 3))
y_height_syn = 0.6
y_height_cov = 0.3
y_margin = 0.05

# This dataset corresponds to a single contig (WOPip3 region)
contig = df["contigs"].unique()[0]
subset = df[df["contigs"] == contig]
contig_len = subset["fin"].max()

# --- Synteny track ---
ax_syn = fig.add_axes([0.1, y_height_cov + y_margin, 0.8, y_height_syn])
ax_syn.set_xlim(0, contig_len)
ax_syn.set_ylim(-1, 1)
ax_syn.axis("off")
ax_syn.text(-1000, 0, "WOPip3", ha="right", va="center", fontsize=10, fontweight="bold")
ax_syn.hlines(0, 0, contig_len, color="black", linewidth=2, zorder=0)

# Draw gene arrows colored by functional category
for _, row in subset.iterrows():
    start, end = row["deb"], row["fin"]
    direction = row["sens"]
    prot_id = str(row["protein_id"])
    annotation = protein_function_map.get(prot_id)
    color = get_color(annotation)

    arrow = (
        FancyArrowPatch((start, 0), (end, 0), arrowstyle="simple", color=color,
                        mutation_scale=15, linewidth=1.5)
        if direction == "+"
        else FancyArrowPatch((end, 0), (start, 0), arrowstyle="simple", color=color,
                             mutation_scale=15, linewidth=1.5)
    )
    ax_syn.add_patch(arrow)

# --- Coverage track ---
ax_cov = fig.add_axes([0.1, 0, 0.8, y_height_cov])

df_cov = pd.read_csv(coverage_folder, header=0, sep="\t", names=["contig", "position", "coverage"])
df_cov = df_cov.dropna()
df_cov["position"] = pd.to_numeric(df_cov["position"], errors="coerce")
df_cov["coverage"] = pd.to_numeric(df_cov["coverage"], errors="coerce")

# Re-index positions to start at 1 (coverage was computed on wPipPel genomic coordinates)
offset = df_cov["position"].min() - 1
df_cov["position"] = df_cov["position"] - offset

ax_cov.plot(df_cov["position"], df_cov["coverage"], color="blue", linewidth=0.7)
ax_cov.fill_between(df_cov["position"], 0, df_cov["coverage"], color="blue", alpha=0.2)
ax_cov.set_xlim(0, contig_len)
ax_cov.set_ylabel("Cov.", fontsize=6)
ax_cov.tick_params(axis="both", labelsize=6)
ax_cov.spines["top"].set_visible(False)
ax_cov.spines["right"].set_visible(False)

# === Legend ===
legend_elements = [
    mpatches.Patch(color="#660066", label="Head"),
    mpatches.Patch(color="#ff0000", label="Ankyrin"),
    mpatches.Patch(color="#0432ff", label="Connector/Baseplate"),
    mpatches.Patch(color="#ffc1c0", label="Tail Fiber"),
    mpatches.Patch(color="#ff6666", label="Tail"),
    mpatches.Patch(color="#ffa4ff", label="CifA/B"),
    mpatches.Patch(color="#07f256", label="Recombinase"),
    mpatches.Patch(color="#fcfc2d", label="Transposase"),
    mpatches.Patch(color="#4aeaed", label="WOPC1"),
    mpatches.Patch(color="#ff6600", label="Toxins"),
    mpatches.Patch(color="#712ca3", label="LigA"),
    mpatches.Patch(color="#7340f5", label="TR"),
    mpatches.Patch(color="#593601", label="Patatin"),
    mpatches.Patch(color="#272545", label="MutL"),
    mpatches.Patch(color="#88d4eb", label="Replication/Repair"),
    mpatches.Patch(color="#538201", label="Methylase"),
]
fig.legend(
    handles=legend_elements,
    loc="center left",
    bbox_to_anchor=(0.91, 0.5),
    title="Legend",
    fontsize=8,
    title_fontsize=9
)

# === Export ===
fig.savefig(output_path, bbox_inches="tight")
plt.show(fig)
print("Figure saved to:", output_path)

png

Figure saved to: Output/WOPip3_synt.pdf

WOPip4 coverage

The following inputs are required:

Data/WOPip4_gff.xlsx - Gene calling on WOPip4.fasta
Data/Orthogroups_Pip4.tsv - Orthofinder launched against Wolbachia genes database launched on WOPip4 CDS
Data/WOPip4.csv - Coverage of Isolate reads on WOPip4 region of wPipPel

# WOPip4: gene synteny (functional categories) + isolate coverage track

# === Inputs / outputs ===
gff_path = "Data/WOPip4_gff.xlsx"
orthogroups_path = "Data/Orthogroups_Pip4.tsv"
coverage_folder = "Data/WOPip4.csv"
output_path = "Output/WOPip4_synt.pdf"

# === Load gene coordinates and orthogroup-derived annotations ===
df = pd.read_excel(gff_path)
orthogroups_df = pd.read_csv(orthogroups_path, sep="\t")

df.columns = [col.strip().lower() for col in df.columns]
df["protein_id"] = df["proteines"].str.extract(r"___(\d+)$")[0]

# Map protein IDs to functional labels (from orthogroup annotation)
protein_function_map = {}
for _, row in orthogroups_df.iterrows():
    if pd.notna(row["PROTEINES_WOPip4"]):
        for pid in row["PROTEINES_WOPip4"].split(","):
            protein_function_map[pid.strip()] = str(row["WOqueries_041124"])

# === Functional category color mapping ===
def get_color(annotations):
    """Assign colors based on the primary functional label; prioritize cif annotations when present."""
    if annotations is None:
        return "grey"
    annotations = annotations.split(",")
    annotation = annotations[0]
    for a in annotations:
        if "cif" in a:
            annotation = a
    a = annotation

    return {
        "HEAD": "#660066",
        "ANK": "#ff0000",
        "CON_": "#0432ff",
        "TailFiber": "#ffc1c0",
        "Tail_": "#ff6666",
        "RECOMB": "#07f256",
        "LIGA": "#712ca3",
        "TOXIN": "#de8209",
        "WOPC1": "#4aeaed",
        "cifA": "#ffa4ff",
        "cifB": "#ffa4ff",
        "TR": "#7340f5",
        "PATATIN": "#593601",
        "TNP": "#fcfc2d",
        "MutL": "#272545",
        "METHYL": "#538201",
        "REP": "#88d4eb",
    }.get(
        next(
            (k for k in [
                "HEAD","ANK","CON_","TailFiber","Tail_","RECOMB","LIGA","TOXIN","WOPC1","cifA","cifB",
                "TR","PATATIN","TNP","MutL","METHYL","REP"
            ] if k in a),
            None
        ),
        "#c7c5c5"
    )

# === Figure layout: synteny track (top) + coverage track (bottom) ===
fig = plt.figure(figsize=(20, 3))
y_height_syn = 0.6
y_height_cov = 0.3
y_margin = 0.05

# This dataset corresponds to a single contig (WOPip4 region)
contig = df["contigs"].unique()[0]
subset = df[df["contigs"] == contig]
contig_len = subset["fin"].max()

# --- Synteny track ---
ax_syn = fig.add_axes([0.1, y_height_cov + y_margin, 0.8, y_height_syn])
ax_syn.set_xlim(0, contig_len)
ax_syn.set_ylim(-1, 1)
ax_syn.axis("off")
ax_syn.text(-1000, 0, "WOPip4", ha="right", va="center", fontsize=10, fontweight="bold")
ax_syn.hlines(0, 0, contig_len, color="black", linewidth=2, zorder=0)

# Draw gene arrows colored by functional category
for _, row in subset.iterrows():
    start, end = row["deb"], row["fin"]
    direction = row["sens"]
    prot_id = str(row["protein_id"])
    annotation = protein_function_map.get(prot_id)
    color = get_color(annotation)

    arrow = (
        FancyArrowPatch((start, 0), (end, 0), arrowstyle="simple", color=color,
                        mutation_scale=15, linewidth=1.5)
        if direction == "+"
        else FancyArrowPatch((end, 0), (start, 0), arrowstyle="simple", color=color,
                             mutation_scale=15, linewidth=1.5)
    )
    ax_syn.add_patch(arrow)

# --- Coverage track ---
ax_cov = fig.add_axes([0.1, 0, 0.8, y_height_cov])

df_cov = pd.read_csv(coverage_folder, header=0, sep="\t", names=["contig", "position", "coverage"])
df_cov = df_cov.dropna()
df_cov["position"] = pd.to_numeric(df_cov["position"], errors="coerce")
df_cov["coverage"] = pd.to_numeric(df_cov["coverage"], errors="coerce")

# Re-index positions to start at 1 (coverage was computed on wPipPel genomic coordinates)
offset = df_cov["position"].min() - 1
df_cov["position"] = df_cov["position"] - offset

ax_cov.plot(df_cov["position"], df_cov["coverage"], color="blue", linewidth=0.7)
ax_cov.fill_between(df_cov["position"], 0, df_cov["coverage"], color="blue", alpha=0.2)
ax_cov.set_xlim(0, contig_len)
ax_cov.set_ylabel("Cov.", fontsize=6)
ax_cov.tick_params(axis="both", labelsize=6)
ax_cov.spines["top"].set_visible(False)
ax_cov.spines["right"].set_visible(False)

# === Legend ===
legend_elements = [
    mpatches.Patch(color="#660066", label="Head"),
    mpatches.Patch(color="#ff0000", label="Ankyrin"),
    mpatches.Patch(color="#0432ff", label="Connector/Baseplate"),
    mpatches.Patch(color="#ffc1c0", label="Tail Fiber"),
    mpatches.Patch(color="#ff6666", label="Tail"),
    mpatches.Patch(color="#ffa4ff", label="CifA/B"),
    mpatches.Patch(color="#07f256", label="Recombinase"),
    mpatches.Patch(color="#fcfc2d", label="Transposase"),
    mpatches.Patch(color="#4aeaed", label="WOPC1"),
    mpatches.Patch(color="#ff6600", label="Toxins"),
    mpatches.Patch(color="#712ca3", label="LigA"),
    mpatches.Patch(color="#7340f5", label="TR"),
    mpatches.Patch(color="#593601", label="Patatin"),
    mpatches.Patch(color="#272545", label="MutL"),
    mpatches.Patch(color="#88d4eb", label="Replication/Repair"),
    mpatches.Patch(color="#538201", label="Methylase"),
]
fig.legend(
    handles=legend_elements,
    loc="center left",
    bbox_to_anchor=(0.91, 0.5),
    title="Legend",
    fontsize=8,
    title_fontsize=9
)

# === Export ===
fig.savefig(output_path, bbox_inches="tight")
plt.show(fig)
print("Figure saved to:", output_path)

png

Figure saved to: Output/WOPip4_synt.pdf

WOPip5 coverage

The following inputs are required:

Data/WOPip5_gff.xlsx - Gene calling on WOPip5.fasta
Data/Orthogroups_Pip5.tsv - Orthofinder launched against Wolbachia genes database launched on WOPip5 CDS
Data/WOPip5.csv - Coverage of Isolate reads on WOPip5 region of wPipPel

# WOPip5: gene synteny (functional categories) + isolate coverage track

# === Inputs / outputs ===
gff_path = "Data/WOPip5_gff.xlsx"
orthogroups_path = "Data/Orthogroups_Pip5.tsv"
coverage_folder = "Data/WOPip5.csv"
output_path = "Output/WOPip5_synt.pdf"

# === Load gene coordinates and orthogroup-derived annotations ===
df = pd.read_excel(gff_path)
orthogroups_df = pd.read_csv(orthogroups_path, sep="\t")

df.columns = [col.strip().lower() for col in df.columns]
df["protein_id"] = df["proteines"].str.extract(r"___(\d+)$")[0]

# Map protein IDs to functional labels (from orthogroup annotation)
protein_function_map = {}
for _, row in orthogroups_df.iterrows():
    if pd.notna(row["PROTEINES_WOPip5"]):
        for pid in row["PROTEINES_WOPip5"].split(","):
            protein_function_map[pid.strip()] = str(row["WOqueries_041124"])

# === Functional category color mapping ===
def get_color(annotations):
    """Assign colors based on the primary functional label; prioritize cif annotations when present."""
    if annotations is None:
        return "grey"
    annotations = annotations.split(",")
    annotation = annotations[0]
    for a in annotations:
        if "cif" in a:
            annotation = a
    a = annotation

    return {
        "HEAD": "#660066",
        "ANK": "#ff0000",
        "CON_": "#0432ff",
        "TailFiber": "#ffc1c0",
        "Tail_": "#ff6666",
        "RECOMB": "#07f256",
        "LIGA": "#712ca3",
        "TOXIN": "#de8209",
        "WOPC1": "#4aeaed",
        "cifA": "#ffa4ff",
        "cifB_": "#ffa4ff",
        "TR": "#7340f5",
        "PATATIN": "#593601",
        "TNP": "#fcfc2d",
        "MutL": "#272545",
        "METHYL": "#538201",
        "REP": "#88d4eb",
    }.get(
        next(
            (k for k in [
                "HEAD","ANK","CON_","TailFiber","Tail_","RECOMB","LIGA","TOXIN","WOPC1","cifA","cifB_",
                "TR","PATATIN","TNP","MutL","METHYL","REP"
            ] if k in a),
            None
        ),
        "#c7c5c5"
    )

# === Figure layout: synteny track (top) + coverage track (bottom) ===
fig = plt.figure(figsize=(20, 3))
y_height_syn = 0.6
y_height_cov = 0.3
y_margin = 0.05

# This dataset corresponds to a single contig (WOPip5 region)
contig = df["contigs"].unique()[0]
subset = df[df["contigs"] == contig]
contig_len = subset["fin"].max()

# --- Synteny track ---
ax_syn = fig.add_axes([0.1, y_height_cov + y_margin, 0.8, y_height_syn])
ax_syn.set_xlim(0, contig_len)
ax_syn.set_ylim(-1, 1)
ax_syn.axis("off")
ax_syn.text(-1000, 0, "WOPip5", ha="right", va="center", fontsize=10, fontweight="bold")
ax_syn.hlines(0, 0, contig_len, color="black", linewidth=2, zorder=0)

# Draw gene arrows colored by functional category
for _, row in subset.iterrows():
    start, end = row["deb"], row["fin"]
    direction = row["sens"]
    prot_id = str(row["protein_id"])
    annotation = protein_function_map.get(prot_id)
    color = get_color(annotation)

    arrow = (
        FancyArrowPatch((start, 0), (end, 0), arrowstyle="simple", color=color,
                        mutation_scale=15, linewidth=1.5)
        if direction == "+"
        else FancyArrowPatch((end, 0), (start, 0), arrowstyle="simple", color=color,
                             mutation_scale=15, linewidth=1.5)
    )
    ax_syn.add_patch(arrow)

# --- Coverage track ---
ax_cov = fig.add_axes([0.1, 0, 0.8, y_height_cov])

df_cov = pd.read_csv(coverage_folder, header=0, sep="\t", names=["contig", "position", "coverage"])
df_cov = df_cov.dropna()
df_cov["position"] = pd.to_numeric(df_cov["position"], errors="coerce")
df_cov["coverage"] = pd.to_numeric(df_cov["coverage"], errors="coerce")

# Re-index positions to start at 1 (coverage was computed on wPipPel genomic coordinates)

offset = df_cov["position"].min() - 1
df_cov["position"] = (len(df_cov['position']) - (df_cov["position"] - offset) + 1) # reverse region


ax_cov.plot(df_cov["position"], df_cov["coverage"], color='blue', linewidth=0.7)
ax_cov.fill_between(df_cov["position"], 0, df_cov["coverage"], color='blue', alpha=0.2)
ax_cov.set_xlim(0, contig_len)
ax_cov.set_ylabel("Cov.", fontsize=6)
ax_cov.tick_params(axis='both', labelsize=6)
ax_cov.spines['top'].set_visible(False)
ax_cov.spines['right'].set_visible(False)


# === Legend ===
legend_elements = [
    mpatches.Patch(color="#660066", label="Head"),
    mpatches.Patch(color="#ff0000", label="Ankyrin"),
    mpatches.Patch(color="#0432ff", label="Connector/Baseplate"),
    mpatches.Patch(color="#ffc1c0", label="Tail Fiber"),
    mpatches.Patch(color="#ff6666", label="Tail"),
    mpatches.Patch(color="#ffa4ff", label="CifA/B"),
    mpatches.Patch(color="#07f256", label="Recombinase"),
    mpatches.Patch(color="#fcfc2d", label="Transposase"),
    mpatches.Patch(color="#4aeaed", label="WOPC1"),
    mpatches.Patch(color="#ff6600", label="Toxins"),
    mpatches.Patch(color="#712ca3", label="LigA"),
    mpatches.Patch(color="#7340f5", label="TR"),
    mpatches.Patch(color="#593601", label="Patatin"),
    mpatches.Patch(color="#272545", label="MutL"),
    mpatches.Patch(color="#88d4eb", label="Replication/Repair"),
    mpatches.Patch(color="#538201", label="Methylase"),
]
fig.legend(
    handles=legend_elements,
    loc="center left",
    bbox_to_anchor=(0.91, 0.5),
    title="Legend",
    fontsize=8,
    title_fontsize=9
)

# === Export ===
fig.savefig(output_path, bbox_inches="tight")
plt.show(fig)
print("Figure saved to:", output_path)

png

Figure saved to: Output/WOPip5_synt.pdf

RosaBioinfo

Coverage visualization of alignment of Isolate reads on wPipPel and WO regions

Table of Contents

All wPipPel regions

WOPip1 coverage

WOPip2 coverage

WOPip3 coverage

WOPip4 coverage

WOPip5 coverage