app.py

import streamlit as st
import os
from analysis import load_data, compute_statistics, full_analysis
from plots import plot_summary_panels, plot_mena_delta_summary, plot_omim_medical_heatmaps, upset_plot_matplotlib, plot_pie_grid_3x5
from cyvcf2 import VCF
import tempfile
import subprocess
import globals
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
import glob

def save_and_index_vcf(uploaded_file, *, suffix=".vcf.gz"):


    with tempfile.NamedTemporaryFile(delete=False, suffix=suffix) as tmp:
        tmp.write(uploaded_file.getbuffer())
        tmp_path = tmp.name

    try:
        subprocess.run(
            ["bcftools", "index", "-f", tmp_path],
            check=True,
            capture_output=True,
            text=True,
        )
    except subprocess.CalledProcessError as e:
        for ext in ("", ".tbi", ".csi"):
            try:
                os.remove(tmp_path + ext)
            except FileNotFoundError:
                pass
        raise RuntimeError(
            f"bcftools index failed:\nSTDOUT:\n{e.stdout}\nSTDERR:\n{e.stderr}"
        )

    return tmp_path

st.set_page_config(page_title="A structural variant catalogue for MENA populations exposes critical gaps in global genomic databases", layout="centered")

st.title("SV diagnostic tool")
st.markdown(
"""
This is a web-based app that overlaps patient SVs with several databases to reduce the number of actionable SVs and enhance diagnostics. 

Patient SVs are overlapped with various databases based on the following parameters:  (Write the method of overlap). However, users have control over the similarity threshold when comparing with databases on the left-hand side. Generally, 50% reciprocal overlap is recommended to identify if an SV was reported previously. 

The databases used in the app to reduce the number of actionable SVs are;  

1- Human Genome Structural Variation Consortium, Phase 2 (HGSVC2). [https://www.internationalgenome.org/data-portal/data-collection/hgsvc2]

2- Sniffles-based calls from 1,019 individuals from the 1k-ONT project [https://www.nature.com/articles/s41586-025-09290-7#Sec2] merged with 61 individuals of MENA ancestry using  SURVIVOR [https://github.com/fritzsedlazeck/SURVIVOR] 

3- SVs from 61 MENA individuals, where an SV was called by at least three different SV callers, including  Svim, dellly, Sniffles and CuteSV. 

The app aims to identify previously reported SVs from publicly available databases and reduce the number of SVs in patients overlapping with OMIM genes exons (https://www.omim.org/) and GIAB medically relevant genes [https://www.nature.com/articles/s41467-024-53260-y#data-availability]. At the end of the anaysis users will be able to download patient SVs with various allele frequency metrics from the 1K-ONT and a list of patient SVs unique to that patient that were not reported in any of the databases, including SVs overlapping with OMIM exons and medically relevant genes
"""
)

if not os.path.exists("./output/"):
    os.makedirs("./output/")
    
    
if not os.path.exists("./figures/"):
    os.makedirs("./figures/")
    
    


st.sidebar.header("Input")
uploaded_files = st.sidebar.file_uploader(
    "Upload CVF file",
    type=["gz"],
    accept_multiple_files=True
)


if uploaded_files is None:
    st.info("Please upload a cvf file to get started.")
    st.stop()




with st.sidebar.form("parameters_form"):
    threshold = st.sidebar.number_input(
    "Similarity threshold",
    min_value=0.0,
    max_value=1.0,
    value=0.5,
    step=0.01,
    format="%.2f",
    )


    submitted = st.form_submit_button("Run analysis")

if not submitted:
    st.info("Adjust parameters and click *Run analysis*.")
    st.stop()


st.write(f"Running with threshold = {threshold}")


@st.cache_data
def load_cached(file):
    return load_data(file)

summary_tables = []
full_counts = np.array([[0, 0, 0, 0, 0, 0],[0, 0, 0, 0, 0, 0],[0, 0, 0, 0, 0, 0]])
sv_len = []

for file_num, uploaded_file in enumerate(uploaded_files):

    tmp_path = save_and_index_vcf(uploaded_file, suffix=".vcf.gz")
    df, summary, counts, data_dicts = full_analysis(tmp_path, os.path.basename(uploaded_file.name) + "_summary.csv", threshold, file_num, threads = 10)
    full_counts = full_counts + np.array(counts)
    
    st.subheader(uploaded_file.name)
    st.write(f"Results for file {uploaded_file.name}")
    st.subheader("Figure")
    
    st.subheader("Preview of data")
    st.dataframe(df.head())
    st.subheader("Figure")
    
    titles = ["SVs overlapping with medical relevant genes", "Subset: has_medical_gene (medically difficult)", "SVs overlapping with OMIM exons"]
    figs = plot_summary_panels(df)
    for fig_number, fig in enumerate(figs):
        plt.close(fig)
        st.pyplot(fig)
        plt.savefig("./figures/" + "_".join(titles[fig_number].split(" ")) +"_" + os.path.basename(uploaded_file.name) + ".png")
    
    summary_tables.append(summary)
    sv_len.append(len(df))

    st.subheader("Reduction of actionable SVs after overlapping with databases")
    
    st.table(summary.applymap(
    lambda x: f"{x:.2f}" if isinstance(x, float) else x
    ))
    
    #fig7 = plot_mena_delta_summary(df, count_unique_ids=True)
    #st.pyplot(fig7, clear_figure=True)

    #fig8 = plot_omim_medical_heatmaps(df, count_unique_ids=True)
    #st.pyplot(fig8, clear_figure=True)
    
    data = [[summary["OMIM"][0] - counts[0][4], counts[0][4]],
        [summary["OMIM"][0] - counts[0][5], counts[0][5]],
        [summary["OMIM"][0] - counts[0][1], counts[0][1]],
        [summary["OMIM"][0] - counts[0][2], counts[0][2]],
        [summary["OMIM"][0] - counts[0][3], counts[0][3]],
        [summary["medical genes"][0] - counts[1][4], counts[1][4]],
        [summary["medical genes"][0] - counts[1][5], counts[1][5]],
        [summary["medical genes"][0] - counts[1][1], counts[1][1]], 
        [summary["medical genes"][0] - counts[1][2], counts[1][2]],
        [summary["medical genes"][0] - counts[1][3], counts[1][3]], 
        [summary["OMIM x medical gene"][0] - counts[2][4], counts[2][4]],
        [summary["OMIM x medical gene"][0] - counts[2][5], counts[2][5]],
        [summary["OMIM x medical gene"][0] - counts[2][1], counts[2][1]],
        [summary["OMIM x medical gene"][0] - counts[2][2], counts[2][2]],
        [summary["OMIM x medical gene"][0] - counts[2][3], counts[2][3]]]



    labels = [["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"], ["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"], ["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"]]


    titles = ["OMIM", "OMIM", "OMIM", "OMIM", "OMIM", "medical", "medical", "medical", "medical","medical", "medical-OMIM", "medical-OMIM", "medical-OMIM", "medical-OMIM","medical-OMIM"]


    fig9 = plot_pie_grid_3x5(data, labels, titles)
    st.pyplot(fig9, clear_figure=True)
    plt.savefig("./figures/" + "pie_summary_plot_" + os.path.basename(uploaded_file.name) + ".png")
        
    df2 = pd.DataFrame(data_dicts)
    
    fig10 = upset_plot_matplotlib(
    df2,
    set_columns=["medical_relevant_gene", "MENA","OMIM", "HGVSC2", "GLOBAL"],
    count_unique_ids=True,
    id_col="id",
    title="Overlap: medical vs OMIM vs MENA vs hgvsc2 vs GLOBAL",
    min_subset_size=2,
    )
    st.pyplot(fig10, clear_figure=True)
    plt.savefig("./figures/" + "UPSET_plot_" + os.path.basename(uploaded_file.name) + ".png")
    
    
    
    


if len(summary_tables)>1:
    st.subheader(f"Combined results for all files")
    st.subheader("Reduction of actionable SVs after overlapping with databases (MEAN)")

    sv_num_mean = np.mean(sv_len)
    sv_num_sum = np.sum(sv_len)
    
    summary =  pd.DataFrame({
        f"total # of patient SVs: {sv_num_mean}": ["patient", "patient x MENA", "patient x global", "hgvsc2 x patient", "hgvsc2 x global x patient","hgvsc2 x global x MENA x patient", "percentage reduction"],
        "OMIM": [0,0, 0,0,0,0,0],
        "medical genes": [0,0, 0,0,0,0,0],
        "OMIM x medical gene": [0,0, 0,0,0,0,0],
        "rate of reduction from OMIM": [0,0, 0,0,0,0,""],

    })
    
    summary["OMIM"] = [np.mean([s["OMIM"][0] for s in summary_tables]), np.mean([s["OMIM"][1] for s in summary_tables]), np.mean([s["OMIM"][2] for s in summary_tables]), np.mean([s["OMIM"][3] for s in summary_tables]), np.mean([s["OMIM"][4] for s in summary_tables]), np.mean([s["OMIM"][5] for s in summary_tables]), (np.mean([s["OMIM"][0] for s in summary_tables]) - np.mean([s["OMIM"][5] for s in summary_tables]))/np.mean([s["OMIM"][0] for s in summary_tables])]

    summary["medical genes"] = [np.mean([s["medical genes"][0] for s in summary_tables]), np.mean([s["medical genes"][1] for s in summary_tables]), np.mean([s["medical genes"][2] for s in summary_tables]), np.mean([s["medical genes"][3] for s in summary_tables]), np.mean([s["medical genes"][4] for s in summary_tables]), np.mean([s["medical genes"][5] for s in summary_tables]), (np.mean([s["medical genes"][0] for s in summary_tables]) - np.mean([s["medical genes"][5] for s in summary_tables]))/np.mean([s["medical genes"][0] for s in summary_tables])]
    
    summary["OMIM x medical gene"] = [np.mean([s["OMIM x medical gene"][0] for s in summary_tables]), np.mean([s["OMIM x medical gene"][1] for s in summary_tables]), np.mean([s["OMIM x medical gene"][2] for s in summary_tables]), np.mean([s["OMIM x medical gene"][3] for s in summary_tables]), np.mean([s["OMIM x medical gene"][4] for s in summary_tables]), np.mean([s["OMIM x medical gene"][5] for s in summary_tables]), (np.mean([s["OMIM x medical gene"][0] for s in summary_tables]) - np.mean([s["OMIM x medical gene"][5] for s in summary_tables]))/np.mean([s["OMIM x medical gene"][0] for s in summary_tables])]

    summary["rate of reduction from OMIM"] = [1-(np.mean([s["OMIM x medical gene"][0] for s in summary_tables])/np.mean([s["OMIM"][0] for s in summary_tables])), 1-(np.mean([s["OMIM x medical gene"][1] for s in summary_tables])/np.sum([s["OMIM"][1] for s in summary_tables])), 1-(np.mean([s["OMIM x medical gene"][2] for s in summary_tables])/np.mean([s["OMIM"][2] for s in summary_tables])), 1-(np.mean([s["OMIM x medical gene"][3] for s in summary_tables])/np.mean([s["OMIM"][3] for s in summary_tables])), 1-(np.mean([s["OMIM x medical gene"][4] for s in summary_tables])/np.mean([s["OMIM"][4] for s in summary_tables])), 1-(np.mean([s["OMIM x medical gene"][5] for s in summary_tables])/np.mean([s["OMIM"][5] for s in summary_tables])), ""]
     
     
    st.table(summary.applymap(
    lambda x: f"{x:.2f}" if isinstance(x, float) else x
    ))

    st.subheader("Reduction of actionable SVs after overlapping with databases (SUM)")
    
    summary =  pd.DataFrame({
        f"total # of patient SVs: {sv_num_sum}": ["patient", "patient x MENA", "patient x global", "hgvsc2 x patient", "hgvsc2 x global x patient","hgvsc2 x global x MENA x patient", "percentage reduction"],
        "OMIM": [0,0, 0,0,0,0,0],
        "medical genes": [0,0, 0,0,0,0,0],
        "OMIM x medical gene": [0,0, 0,0,0,0,0],
        "rate of reduction from OMIM": [0,0, 0,0,0,0,""],
    })
    
    summary["OMIM"] = [np.sum([s["OMIM"][0] for s in summary_tables]), np.sum([s["OMIM"][1] for s in summary_tables]), np.sum([s["OMIM"][2] for s in summary_tables]), np.sum([s["OMIM"][3] for s in summary_tables]), np.sum([s["OMIM"][4] for s in summary_tables]), np.sum([s["OMIM"][5] for s in summary_tables]), (np.sum([s["OMIM"][0] for s in summary_tables]) - np.sum([s["OMIM"][5] for s in summary_tables]))/np.sum([s["OMIM"][0] for s in summary_tables])]

    summary["medical genes"] = [np.sum([s["medical genes"][0] for s in summary_tables]), np.sum([s["medical genes"][1] for s in summary_tables]), np.sum([s["medical genes"][2] for s in summary_tables]), np.sum([s["medical genes"][3] for s in summary_tables]), np.sum([s["medical genes"][4] for s in summary_tables]), np.sum([s["medical genes"][5] for s in summary_tables]), (np.sum([s["medical genes"][0] for s in summary_tables]) - np.sum([s["medical genes"][5] for s in summary_tables]))/np.sum([s["medical genes"][0] for s in summary_tables])]
    
    summary["OMIM x medical gene"] = [np.sum([s["OMIM x medical gene"][0] for s in summary_tables]), np.sum([s["OMIM x medical gene"][1] for s in summary_tables]), np.sum([s["OMIM x medical gene"][2] for s in summary_tables]), np.sum([s["OMIM x medical gene"][3] for s in summary_tables]), np.sum([s["OMIM x medical gene"][4] for s in summary_tables]), np.sum([s["OMIM x medical gene"][5] for s in summary_tables]), (np.sum([s["OMIM x medical gene"][0] for s in summary_tables]) - np.sum([s["OMIM x medical gene"][5] for s in summary_tables]))/np.sum([s["OMIM x medical gene"][0] for s in summary_tables])]
    
    summary["rate of reduction from OMIM"] = [1-(np.sum([s["OMIM x medical gene"][0] for s in summary_tables])/np.sum([s["OMIM"][0] for s in summary_tables])), 1-(np.sum([s["OMIM x medical gene"][1] for s in summary_tables])/np.sum([s["OMIM"][1] for s in summary_tables])), 1-(np.sum([s["OMIM x medical gene"][2] for s in summary_tables])/np.sum([s["OMIM"][2] for s in summary_tables])), 1-(np.sum([s["OMIM x medical gene"][3] for s in summary_tables])/np.sum([s["OMIM"][3] for s in summary_tables])), 1-(np.sum([s["OMIM x medical gene"][4] for s in summary_tables])/np.sum([s["OMIM"][4] for s in summary_tables])), 1-(np.sum([s["OMIM x medical gene"][5] for s in summary_tables])/np.sum([s["OMIM"][5] for s in summary_tables])), ""]
    st.table(summary.applymap(
    lambda x: f"{x:.2f}" if isinstance(x, float) else x
    ))



    
    data = [[summary["OMIM"][0] - counts[0][4], counts[0][4]],
        [summary["OMIM"][0] - counts[0][5], counts[0][5]],
        [summary["OMIM"][0] - counts[0][1], counts[0][1]],
        [summary["OMIM"][0] - counts[0][2], counts[0][2]],
        [summary["OMIM"][0] - counts[0][3], counts[0][3]],
        [summary["medical genes"][0] - counts[1][4], counts[1][4]],
        [summary["medical genes"][0] - counts[1][5], counts[1][5]],
        [summary["medical genes"][0] - counts[1][1], counts[1][1]], 
        [summary["medical genes"][0] - counts[1][2], counts[1][2]],
        [summary["medical genes"][0] - counts[1][3], counts[1][3]], 
        [summary["OMIM x medical gene"][0] - counts[2][4], counts[2][4]],
        [summary["OMIM x medical gene"][0] - counts[2][5], counts[2][5]],
        [summary["OMIM x medical gene"][0] - counts[2][1], counts[2][1]],
        [summary["OMIM x medical gene"][0] - counts[2][2], counts[2][2]],
        [summary["OMIM x medical gene"][0] - counts[2][3], counts[2][3]]]



    labels = [["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"], ["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"], ["patients w.o. MENA", "MENA"], ["patients w.o. global", "global"], ["patients w.o. hgvsc2", "hgvsc2"], ["patients w.o. hgvsc2 or global", "hgvsc2 x global"], ["patients w.o. hgvsc2 or global or MENA", "hgvsc2 x global x MENA"]]


    titles = ["OMIM", "OMIM", "OMIM", "OMIM", "OMIM", "medical", "medical", "medical", "medical","medical", "medical-OMIM", "medical-OMIM", "medical-OMIM", "medical-OMIM","medical-OMIM"]

        
    df2 = pd.DataFrame(data_dicts)
    
    fig10 = upset_plot_matplotlib(
    df2,
    set_columns=["medical_relevant_gene", "MENA","OMIM", "HGVSC2", "GLOBAL"],
    count_unique_ids=True,
    id_col="id",
    title="Overlap: medical vs OMIM vs MENA vs hgvsc2 vs GLOBAL",
    min_subset_size=2,
    )
    st.pyplot(fig10, clear_figure=True)
    plt.savefig("./figures/" + "UPSET_plot_all"+ ".png")