analysis.py

from __future__ import annotations
import numpy as np
import pandas as pd
from scipy import stats
import globals
from cyvcf2 import VCF
from typing import Dict, Any, List, Optional, Tuple
from dataclasses import dataclass
from difflib import SequenceMatcher
import tqdm
from collections import defaultdict
import os
import math
import subprocess
import time
import subprocess
import os
from typing import Optional
import glob
from cyvcf2 import VCF
from typing import Dict, Any, List, Optional, Tuple


def _get_info_key(rec, keys: Tuple[str, ...]) -> Optional[str]:
    """Return first INFO value found among keys."""
    for k in keys:
        v = rec.INFO.get(k)
        if v is not None and v != "":
            return v
    return None


def attach_base_to_tp_comp(
    tp_base_vcf: str,
    tp_comp_vcf: str,
    *,
    base_fields: Tuple[str, ...] = ("ID", "CHROM", "POS", "END", "SVTYPE", "SVLEN", "REF", "ALT", "QUAL", "FILTER"),
    extra_base_info_keys: Tuple[str, ...] = ("AC", "AN", "AF", "SUPP", "SUPPORT", "RNAMES", "STRAND"),
) -> List[Dict[str, Any]]:

    matchid_keys = ("MatchId", "MATCHID", "TRUVARI_MatchId", "TruvariMatchId", "TRUVARI_MATCHID")


    base_by_match: Dict[str, Dict[str, Any]] = {}

    for b in VCF(tp_base_vcf):
        mid = _get_info_key(b, matchid_keys)
        if mid is None:
            continue

        base_rec: Dict[str, Any] = {"MatchId": mid}


        if "ID" in base_fields: base_rec["base_id"] = b.ID
        if "CHROM" in base_fields: base_rec["base_chrom"] = b.CHROM
        if "POS" in base_fields: base_rec["base_pos"] = int(b.POS)
        if "REF" in base_fields: base_rec["base_ref"] = b.REF
        if "ALT" in base_fields: base_rec["base_alt"] = ",".join(list(b.ALT or []))
        if "QUAL" in base_fields: base_rec["base_qual"] = b.QUAL
        if "FILTER" in base_fields: base_rec["base_filter"] = b.FILTER


        if "END" in base_fields: base_rec["base_end"] = b.INFO.get("END")
        if "SVTYPE" in base_fields: base_rec["base_svtype"] = b.INFO.get("SVTYPE")
        if "SVLEN" in base_fields: base_rec["base_svlen"] = b.INFO.get("SVLEN")


        for k in extra_base_info_keys:
            v = b.INFO.get(k)
            if v is not None:
                base_rec[f"base_INFO_{k}"] = v


        base_by_match.setdefault(mid, base_rec)

    if not base_by_match:
        raise RuntimeError(
            "No MatchId found in tp-base. "
            "Either Truvari didn't annotate MatchId (older version / settings), "
            "or MatchId is under a different INFO key."
        )


    merged: List[Dict[str, Any]] = []
    for c in VCF(tp_comp_vcf):
        mid = _get_info_key(c, matchid_keys)
        if mid is None:
            continue

        out: Dict[str, Any] = {"MatchId": mid}


        out["comp_id"] = c.ID
        out["comp_chrom"] = c.CHROM
        out["comp_pos"] = int(c.POS)
        out["comp_ref"] = c.REF
        out["comp_alt"] = ",".join(list(c.ALT or []))
        out["comp_qual"] = c.QUAL
        out["comp_filter"] = c.FILTER
        out["comp_end"] = c.INFO.get("END")
        out["comp_svtype"] = c.INFO.get("SVTYPE")
        out["comp_svlen"] = c.INFO.get("SVLEN")


        base_rec = base_by_match.get(mid)
        if base_rec is None:
            out["base_missing_for_matchid"] = True
        else:
            out.update(base_rec)

        merged.append(out)

    return merged


def bcftools_normalize_and_index(
    input_vcf: str,
    output_vcf: Optional[str] = None,
    *,
    force: bool = True,
    split_multiallelics: bool = True,
    threads: int = 1,
) -> str:

    if output_vcf is None:
        if input_vcf.endswith(".vcf.gz"):
            output_vcf = input_vcf.replace(".vcf.gz", ".norm.vcf.gz")
        else:
            output_vcf = input_vcf + ".norm.vcf.gz"

    norm_cmd = ["bcftools", "norm"]


    if split_multiallelics:
        norm_cmd.append("-m-any")  

    norm_cmd.extend(["-Oz", "-o", output_vcf])

    if threads and threads > 0:
        norm_cmd.extend(["--threads", str(threads)])

    if force:

        for ext in ("", ".tbi", ".csi"):
            p = output_vcf + ext
            if os.path.exists(p):
                os.remove(p)

    norm_cmd.append(input_vcf)

    try:
        subprocess.run(norm_cmd, check=True, capture_output=True, text=True)
    except subprocess.CalledProcessError as e:
        raise RuntimeError(
            f"bcftools norm failed for {input_vcf}\n"
            f"Command: {' '.join(norm_cmd)}\n\n"
            f"STDOUT:\n{e.stdout}\n\nSTDERR:\n{e.stderr}"
        )


    index_cmd = ["bcftools", "index"]
    if force:
        index_cmd.append("-f")
    if threads and threads > 0:
        index_cmd.extend(["--threads", str(threads)])
    index_cmd.append(output_vcf)

    try:
        subprocess.run(index_cmd, check=True, capture_output=True, text=True)
    except subprocess.CalledProcessError as e:
        raise RuntimeError(
            f"bcftools index failed for {output_vcf}\n"
            f"Command: {' '.join(index_cmd)}\n\n"
            f"STDOUT:\n{e.stdout}\n\nSTDERR:\n{e.stderr}"
        )

    return output_vcf
    

def run_truvari_unique(
    *,
    base_vcf: str,
    comp_vcf: str,
    outdir: str,
    reference_fasta: Optional[str] = None,
    pctovl: float = 0.5,
    pctsize: float = 0.5,
    refdist: int = 500,
    pctseq: float = 0.0,
    threads: int = 1,
    force: bool = True,
) -> Dict[str, str]:


    if pctseq > 0 and reference_fasta is None:
        raise ValueError("pctseq > 0 requires reference_fasta")

    os.makedirs(outdir, exist_ok=True)

    cmd = [
        "truvari", "bench",
        "-b", base_vcf,
        "-c", comp_vcf,
        "-o", outdir,
        "--pctovl", str(pctovl),
        "--pctsize", str(pctsize),
        "--refdist", str(refdist),
        "--pctseq", str(pctseq),
    ]

    if reference_fasta is not None:
        cmd.extend(["-f", reference_fasta])

    if force:
        
        for filename in os.listdir(outdir):
            file_path = os.path.join(outdir, filename)
            try:
                if os.path.isfile(file_path) or os.path.islink(file_path):
                    os.unlink(file_path)
                elif os.path.isdir(file_path):
                    shutil.rmtree(file_path)
            except Exception as e:
                print('Failed to delete %s. Reason: %s' % (file_path, e))
        os.rmdir(outdir)
        

    try:
        subprocess.run(
            cmd,
            check=True,
            capture_output=True,
            text=True,
        )
    except subprocess.CalledProcessError as e:
        raise RuntimeError(
            f"Truvari failed\n"
            f"STDOUT:\n{e.stdout}\n\n"
            f"STDERR:\n{e.stderr}"
        )

    results = {
        "tp_base": os.path.join(outdir, "tp-base.vcf.gz"),
        "tp_comp": os.path.join(outdir, "tp-comp.vcf.gz"),
        "fp": os.path.join(outdir, "fp.vcf.gz"),  # unique to comp
        "fn": os.path.join(outdir, "fn.vcf.gz"),  # unique to base
    }

    missing = [k for k, v in results.items() if not os.path.exists(v)]
    if missing:
        raise FileNotFoundError(
            "Expected Truvari output(s) missing:\n"
            + "\n".join(f"  {k}: {results[k]}" for k in missing)
        )

    return results
    

def _revcomp(seq: str) -> str:
    comp = str.maketrans("ACGTNacgtn", "TGCANtgcan")
    return seq.translate(comp)[::-1]
    






def overlap_med_query_vcf(vcf, genes: Dict[str, Dict[str, str]], max_sv_size = 200):

    hits: List[Dict[str, Any]] = []
    ccounter = 0
    for gene_id, g in genes.items():
        coord = g["coords"]
        chrom = coord.split(":")[0]
        start = int(coord.split(":")[1].split("-")[0])
        stop = int(coord.split(":")[1].split("-")[1])
        # cyvcf2 region syntax is 1-based inclusive
        region = f"{chrom}:{start- max_sv_size}-{stop + max_sv_size}"
        try:
            counter = -1
            for rec in vcf(region):
                # Compute variant span
                counter = counter + 1
                v_start = rec.POS
                v_end = rec.POS + len(rec.REF) - 1
                end_info = rec.INFO.get("END")
                
                if end_info is not None:
                    try:
                        v_end = int(end_info)
                    except ValueError:
                        pass

                # Explicit overlap check (paranoia is good)
                if v_end < start or v_start > stop:
                    ccounter = ccounter +1 
                    continue
                    
                if v_start < start or v_start > stop:
                    continue
                    
                

                hits.append(
                    {
                        "chrom": chrom,
                        "pos": rec.POS,
                        "ref": rec.REF,
                        "alt": list(rec.ALT),
                        "gene": gene_id,
                        "coord": (g.get("coord") or "").strip(),
                        "in_mrg": (g.get("in_mrg") or "").strip(),
                        "counter": counter
                    }
                )
                

        except Exception as e:
            # Defensive: region queries can fail if contigs mismatch
            continue
    # 5027) and column (in_mrg)
    print("ccounter")
    print(ccounter)
    return hits
    
    
def overlap_genes_query_vcf(vcf, genes: Dict[str, Dict[str, str]], chrom_key = "#chr", 
start_key = "start", stop_key = "stop", name_key = "gene_id", add = True, add_af = False, max_sv_size = 200):

    hits: List[Dict[str, Any]] = []

    for gene_id, g in genes.items():
        chrom = g[chrom_key]
        start = int(g[start_key])
        stop = int(g[stop_key])

        region = f"{chrom}:{start- max_sv_size}-{stop + max_sv_size}"

        try:
            counter = -1
            for rec in vcf(region):
                # Compute variant span

                counter = counter + 1
                v_start = rec.POS
                v_end = rec.POS + len(rec.REF) - 1
                end_info = rec.INFO.get("END")
                if end_info is not None:
                    try:
                        v_end = int(end_info)
                    except ValueError:
                        pass

                # allow for any overlap
                if v_end < start or v_start > stop:
                    continue
                info = {
                        "chrom": chrom,
                        "pos": rec.POS,
                        "ref": rec.REF,
                        "alt": list(rec.ALT),
                        "info": list(rec.INFO),
                        name_key: gene_id,
                        "counter": counter,                        
                    }
                    
                if add == True:
                        info["gene_name"] = g.get("gene_name")
                        info["exon_number"] = (g.get("exon_number") or "").strip()
                        
                if add_af == True:
                        info["AMR_af"] = g.get("AMR_af")
                        info["AFR_af"] = g.get("AFR_af")
                        info["EAS_af"] = g.get("EAS_af")
                        info["EUR_af"] = g.get("EUR_af")
                        info["MENA_af"] = g.get("MENA_af")
                        info["SAS_af"] = g.get("SAS_af")

                hits.append(info)
        except Exception as e:
            # Defensive: region queries can fail if contigs mismatch
            continue

    return hits


def patient_vs_omim(patient_vcf, omim_exons):
    hits = overlap_genes_query_vcf(patient_vcf, omim_exons)
    return hits
    
def patient_vs_medical(patient_vcf, medical_data):
    hits = overlap_med_query_vcf(patient_vcf, medical_data)
    ### subselect medically difficult genes 
    hits_medically_dificult = [h for h in hits if h["in_mrg"] == 'TRUE']
    return hits, hits_medically_dificult
    
    

def load_background_statistics():
    
    background_vcf = VCF(globals.caller3)
    
    medical_data = {}
    with open(globals.medically_relevant_genes) as f:
        header = f.readline().rstrip("\n").split("\t")
        for line in f:
            fields = line.rstrip("\n").split("\t")
            row = dict(zip(header, fields))
            medical_data[row["gene"]] = row

    omin = {}
    with open(globals.coordinates_omin_exons) as f:
        skip = f.readline()
        header = f.readline().rstrip("\n").split("\t")
        for line in f:
            fields = line.rstrip("\n").split("\t")
            row = dict(zip(header, fields))
            omin[row["gene_id"]] = row
            
    frequencies = {}
    with open(globals.ONT_data) as f:
        header = f.readline().rstrip("\n").split("\t")
        for line in f:
            fields = line.rstrip("\n").split("\t")
            row = dict(zip(header, fields))
            frequencies[row["sv_id"]] = row

    return omin, frequencies, medical_data, background_vcf



@dataclass(frozen=True)
class Interval:
    chrom: str
    start0: int           
    end0: int           
    id: str
    ref: str
    alt: str
    svtype: str = "NA"   
    pos1: Optional[int] = None   
    ins_seq: Optional[str] = None 
    ac: int = None
    an:int = None
    af: float = None

    @property
    def length(self) -> int:
        return max(0, self.end0 - self.start0)



def vcf_to_intervals(
    vcf_path: str,
    *,
    assume_end_from_ref: bool = True,
    use_info_end: bool = True,
    skip_symbolic_no_end: bool = True,
    region: Optional[str] = None,
    anchor_ins_as_1bp: bool = True,
    skip_symbolic_ins_without_seq: bool = True,
    load_allele_stats: bool = False, 
) -> List["Interval"]:


    def _as_str(x) -> str:
        if x is None:
            return ""
        if isinstance(x, bytes):
            return x.decode("utf-8", errors="ignore")
        return str(x)

    def _as_int_scalar(x) -> Optional[int]:

        if x is None:
            return None

        if isinstance(x, (list, tuple)):
            if len(x) == 0:
                return None
            x = x[0]
        if isinstance(x, bytes):
            try:
                x = x.decode("utf-8", errors="ignore")
            except Exception:
                return None
        try:
            return int(x)
        except Exception:
            return None

    def _as_float_scalar(x) -> Optional[float]:
        if x is None:
            return None
        if isinstance(x, (list, tuple)):
            if len(x) == 0:
                return None
            x = x[0]
        if isinstance(x, bytes):
            try:
                x = x.decode("utf-8", errors="ignore")
            except Exception:
                return None
        try:
            return float(x)
        except Exception:
            return None

    def _infer_svtype(ref: str, alt: str, svtype_from_info: str) -> str:
        if svtype_from_info:
            return svtype_from_info
        if alt.startswith("<") and alt.endswith(">"):
            return alt.strip("<>").upper()
        if len(alt) > len(ref):
            return "INS"
        if len(alt) < len(ref):
            return "DEL"
        if len(ref) == 1 and len(alt) == 1:
            return "SNV"
        return "INDEL"

    vcf = VCF(vcf_path)
    it = vcf(region) if region else vcf

    out: List["Interval"] = []

    for rec in it:
        chrom = rec.CHROM
        pos1 = int(rec.POS)
        start0 = pos1 - 1

        rec_id = rec.ID if rec.ID is not None else "."
        ref = rec.REF if rec.REF is not None else "N"

        svtype_info = _as_str(rec.INFO.get("SVTYPE")).upper()


        ac: Optional[int] = None
        an: Optional[int] = None
        af: Optional[float] = None
        if load_allele_stats:
            ac = _as_int_scalar(rec.INFO.get("AC"))
            an = _as_int_scalar(rec.INFO.get("AN"))


            if an is not None and an != 0 and ac is not None:
                af = ac / an


        info_end: Optional[int] = None
        if use_info_end:
            info_end = _as_int_scalar(rec.INFO.get("END"))

        alts = list(rec.ALT) if rec.ALT is not None else ["."]
        if info_end is None:
            is_symbolic_any = any(str(a).startswith("<") and str(a).endswith(">") for a in alts)
            if is_symbolic_any and skip_symbolic_no_end:
                continue

        for alt_raw in alts:
            alt = str(alt_raw)
            svtype = _infer_svtype(ref, alt, svtype_info)


            if svtype != "INS":
                if info_end is not None and info_end > pos1:
                    end0 = info_end  # [POS-1, END)
                else:
                    end0 = start0 + (max(1, len(ref)) if assume_end_from_ref else 1)
            else:
                end0 = start0 + 1 if anchor_ins_as_1bp else start0


            ins_seq: Optional[str] = None
            if svtype == "INS":
                is_symbolic = alt.startswith("<") and alt.endswith(">")
                if is_symbolic:
                    if skip_symbolic_ins_without_seq:
                        continue
                    ins_seq = None
                else:
                    if alt.startswith(ref) and len(alt) > len(ref):
                        ins_seq = alt[len(ref):]
                    else:
                        ins_seq = alt

            out.append(
                Interval(
                    chrom=chrom,
                    start0=start0,
                    end0=end0,
                    id=rec_id,
                    ref=ref,
                    alt=alt,
                    svtype=svtype,
                    pos1=pos1,
                    ins_seq=ins_seq,

                    ac=ac if load_allele_stats else None,
                    an=an if load_allele_stats else None,
                    af=af if load_allele_stats else None,
                )
            )

    return out



def _parasail_identity(a: str, b: str) -> float:

    import parasail

    if not a or not b:
        return 0.0


    match = 2
    mismatch = -2
    gap_open = 3
    gap_extend = 1
    matrix = parasail.matrix_create("ACGTN", match, mismatch)

    # Global alignment (Needleman-Wunsch)
    res = parasail.nw_stats_striped_16(a, b, gap_open, gap_extend, matrix)
    matches = res.matches
    length = res.length
    if length <= 0:
        return 0.0
    return matches / length

def _best_ins_identity(a: str, b: str) -> float:
    return max(_parasail_identity(a, b), _parasail_identity(a, _revcomp(b)))


def reciprocal_overlap_pairs(
    a_intervals: List[Interval],
    b_intervals: List[Interval],
    *,
    ro_threshold: float = 0.5,
    max_ins_pos_delta: int = 25,
    min_ins_identity: float = 0.80,
    min_ins_len: int = 1,
) -> List[Tuple[Interval, Interval, int, Optional[float], Optional[float], Optional[float]]]:

    if ro_threshold <= 0 or ro_threshold > 1:
        raise ValueError("ro_threshold must be in (0, 1].")


    B_by_chr: Dict[str, List[Interval]] = {}
    for b in b_intervals:
        B_by_chr.setdefault(b.chrom, []).append(b)
    for chrom in B_by_chr:
        B_by_chr[chrom].sort(key=lambda x: (x.start0, x.end0))

    matches: List[Tuple[Interval, Interval, int, Optional[float], Optional[float], Optional[float]]] = []
    ins_cont_counter = 0

    for a in tqdm.tqdm(a_intervals):

        b_list = B_by_chr.get(a.chrom)
        
        if not b_list:
            continue

        
        if a.svtype == "INS":

            if a.pos1 is None or not a.ins_seq or len(a.ins_seq) < min_ins_len:
                continue

            for b in b_list:
                if b.svtype != "INS":
                    continue
                if b.pos1 is None or not b.ins_seq or len(b.ins_seq) < min_ins_len:
                    continue
                if abs(a.pos1 - b.pos1) > max_ins_pos_delta:
                    continue
                if b.end0 <= a.start0:
                    continue
                if b.start0 >= a.end0:
                    continue
                
                
                ident = _best_ins_identity(a.ins_seq, b.ins_seq)
                if ident >= min_ins_identity:
                    matches.append((a, b, 1, 1.0, 1.0, ident))
                    
            
            continue

        for b in b_list:
            if a.svtype == "INS":
                continue
            if b.end0 <= a.start0:
                continue
            if b.start0 >= a.end0:
                continue

            ov = min(a.end0, b.end0) - max(a.start0, b.start0)
            if ov <= 0:
                continue

            lenA = a.length
            lenB = b.length
            if lenA <= 0 or lenB <= 0:
            
                continue

            roA = ov / lenA
            roB = ov / lenB
            
            if roA < ro_threshold or roB< ro_threshold:
                ins_cont_counter = ins_cont_counter +1
            if roA >= ro_threshold and roB >= ro_threshold:
                matches.append((a, b, ov, roA, roB, None))

    return matches


def write_bed(intervals: List[Interval], bed_path: str) -> None:

    with open(bed_path, "w") as f:
        for iv in intervals:
            name = iv.id if iv.id != "." else f"{iv.chrom}:{iv.start0+1}"
            f.write(f"{iv.chrom}\t{iv.start0}\t{iv.end0}\t{name}\t{iv.ref}\t{iv.alt}\n")

    
def create_output_table2(patients_svs,
    omim_hits: List[Dict[str, Any]],
    medical_hits: List[Dict[str, Any]],
    mena_labels: List[Dict[str, Any]],
    global_hits: List[Dict[str, Any]],
    hvg_labels: List[Tuple[Any, Any]],
    *,
    out_path: Optional[str] = None,
) -> pd.DataFrame:

    def key(chrom: str, pos: int) -> Tuple[str, int]:
        return (chrom, int(pos))
        

    omim_by = defaultdict(list)
    for o in omim_hits:
        omim_by[key(o["chrom"], o["pos"])].append(o)
        
    medical_by = defaultdict(list)
    for m in medical_hits:
        medical_by[key(m["chrom"], m["pos"])].append(m)
        
        
    global_by = defaultdict(list)
    for g in global_hits:
        global_by[key(g["chrom"], g["pos"])].append(g)
        
        
    hvg_by = defaultdict(list)
    for g in hvg_labels:
        hvg_by[key(g.chrom, g.pos1)].append(g)


    mena_by = defaultdict(list)
    for g in mena_labels:
        mena_by[key(g["comp_chrom"], g["comp_pos"])].append(g)
        
        
    final_rows: List[Dict[str, Any]] = []
    enum = 0
    

    for patient_sv in tqdm.tqdm(patients_svs):
    

        k = key(patient_sv.chrom, patient_sv.pos1)
        relevant_global_hits = global_by.get(k, [])
        relevant_omim = omim_by.get(k, [])
        relevant_medical = medical_by.get(k, [])
        relevant_hvg = hvg_by.get(k, [])
        relevant_mena = mena_by.get(k, [])
        
        
        info: Dict[str, Any] = {
            "id": patient_sv.id,
            "chrom": patient_sv.chrom,
            "pos": int(patient_sv.pos1),
            "ref": patient_sv.ref,
            "alt": patient_sv.alt,
            "sv_type": patient_sv.svtype,
        }
        

        if relevant_omim:
            info["OMIM"] = True
            info["OMIM_genes"] = ",".join([o.get("gene_name", "") for o in relevant_omim if o.get("gene_name")])
            info["OMIM_gene_ids"] = ",".join([o.get("gene_id", "") for o in relevant_omim if o.get("gene_id")])

        else:
            info["OMIM"] = False

        if relevant_medical:
            genes = []
            in_mrg_vals = []

            for m in relevant_medical:
                genes.append(m.get("gene", ""))
                in_mrg_vals.append(str(m.get("in_mrg", "")))

            info["medical_relevant_gene"] = ",".join([g for g in genes if g])
            info["in_mrg"] = ",".join([v for v in in_mrg_vals if v])
            
            
        if relevant_global_hits:

            sv_ids = []
            strands = []
            afs = []
            
            afr, eas, eu, sas, amr,menaaf = [],[],[],[],[],[]
            info["global_population"] = True
            for gh in relevant_global_hits:
                sv_ids.append(str(gh.get("sv_id", "")))
                sv_ids.append(str(gh.get("AMR_af", "")))
                afr.append(str(gh.get("AFR_af", "")))
                eas.append(str(gh.get("EAS_af", "")))
                eu.append(str(gh.get("EU_af", "")))
                sas.append(str(gh.get("SAS_af", "")))
                menaaf.append(str(gh.get("MENA_af", "")))
                try:
                    strands.append(str(gh["info"][7][1]))
                except Exception:
                    strands.append("")
                try:
                    afs.append(str(gh["info"][8][1]))
                except Exception:
                    afs.append("")

            info["sv_id"] = ",".join([s for s in sv_ids if s])
            info["strand"] = ",".join([s for s in strands if s])
            info["af"] = ",".join([a for a in afs if a])
            
            info["AMR_af"] = ",".join([a for a in amr if a])
            info["AFR_af"] = ",".join([a for a in afr if a])
            info["EAS_af"] = ",".join([a for a in eas if a])
            info["EU_af"] = ",".join([a for a in eu if a])
            info["SAS_af"] = ",".join([a for a in sas if a])
            info["MENA_af"] = ",".join([a for a in menaaf if a])
            
            
        if relevant_hvg: 
            info["hgvsc"] = False
        else:
            info["hgvsc"] = True
            
        if relevant_mena: 
            info["mena"] = True
            info["mena_af_3caller"] = []
            info["mena_pos"] = []

            for relevant in relevant_mena:

                info["mena_af_3caller"].append(relevant['base_INFO_AC']/max(relevant['base_INFO_AN'],1))
                info["mena_pos"].append(relevant['base_pos'])
                
            info["mena_af_3caller"] = ",".join([str(a) for a in info["mena_af_3caller"] if a])
            info["mena_pos"] = ",".join([str(a) for a in info["mena_pos"] if a])
        else:
            info["mena"] = False
        
        final_rows.append(info)
            
    df = pd.DataFrame(final_rows)
    if out_path.endswith(".csv"):
        df.to_csv(out_path, index=False)
    else:
        df.to_parquet(out_path, index=False)

    return df
    


def full_analysis(input_file, output_name, threshold, file_num,threads = 10):


    A_norm = bcftools_normalize_and_index(globals.caller3)
    B_norm = bcftools_normalize_and_index(input_file)
    
    ###########################################################################################################
    ################################# compare to hgvsc cvfs ##################################################
    ###########################################################################################################

    B_hgvsc_selection = B_norm
    start = time.time()
    for enum, hvgsc_file in enumerate(globals.hgvsc_selection):
        hvgsc_file_norm = bcftools_normalize_and_index(hvgsc_file)
    
        results = run_truvari_unique(
            base_vcf=hvgsc_file_norm,
            comp_vcf=B_hgvsc_selection,
            outdir=f"patient_vs_hvg{enum}_{file_num}",
            threads = threads,
            pctovl = threshold,
        )

        B_hgvsc_selection = results["fp"]
        B_hgvsc_selection = bcftools_normalize_and_index(B_hgvsc_selection)

        os.remove(hvgsc_file_norm)


    ###########################################################################################################
    ################################# compare to global  ##################################################
    ###########################################################################################################

    
    #### frequencies stuff
    omim_exons, frequencies, medical_data, background_vcf = load_background_statistics()
    patient_vcf = VCF(B_norm)
    global_hits = overlap_genes_query_vcf(patient_vcf, frequencies, chrom_key = "chr", start_key = "start", stop_key = "end", name_key = "sv_id", add = False, add_af = True)

    ###########################################################################################################
    ################################# compare to MENA  ##################################################
    ###########################################################################################################

    
    results = run_truvari_unique(
        base_vcf=A_norm,
        comp_vcf=B_norm,
        outdir=f"patient_vs_mena_{file_num}",
        pctovl = threshold,
    )

    shared_base = results["tp_base"] 
    shared_comp = results["tp_comp"]

    pairs = attach_base_to_tp_comp(
    tp_base_vcf=shared_base,
    tp_comp_vcf=shared_comp,
    )


    ###########################################################################################################
    #################### match full files to medical and OMIM #############################
    ###########################################################################################################
    

    omim_hits = patient_vs_omim(patient_vcf, omim_exons)
    medical_hits, hits_medically_dificult = patient_vs_medical(patient_vcf, medical_data)
    patient = vcf_to_intervals(input_file, use_info_end=True)
    hgvcsc = vcf_to_intervals(B_hgvsc_selection, use_info_end=True)
    df = create_output_table2(patient, omim_hits, medical_hits, pairs, global_hits, hgvcsc, out_path = "./output/" + output_name)

    
    
    omim_raw = 0
    omim_and_hgvcs = 0
    omim_and_hgvcs_and_global = 0
    omim_and_hgvcs_or_global_or_mena = 0
    omim_and_mena = 0
    omim_and_global = 0

    med_raw = 0
    med_and_hgvcs = 0
    med_and_hgvcs_and_global = 0
    med_and_hgvcs_or_global_or_mena = 0
    med_and_mena = 0
    med_and_global = 0
    
    
    omim_and_med = 0
    omim_and_med_and_hgvcs = 0
    omim_and_med_and_hgvcs_and_global = 0
    omim_and_med_and_hgvcs_or_global_or_mena = 0
    omim_and_med_and_mena = 0
    omim_and_med_and_global = 0
    
    
    data_dicts = {"id":[],"medical_relevant_gene": [], "MENA": [],"OMIM": [], "HGVSC2": [], "GLOBAL": []}
    
        
    for enum, gene_data in enumerate(df["medical_relevant_gene"]):
        
        omim_available = df["OMIM"][enum]
        mena_hit = df["mena"][enum]
        hgvsc = df["hgvsc"][enum]
        overlap_global_population = df["global_population"][enum]
        id_ = df["id"][enum] + str(enum)
            
            
            
        if pd.isnull(gene_data) == False: data_dicts["medical_relevant_gene"].append(True)
        else: data_dicts["medical_relevant_gene"].append(False)
        data_dicts["OMIM"].append(omim_available)
        data_dicts["MENA"].append(mena_hit)
        data_dicts["HGVSC2"].append(hgvsc)
        if overlap_global_population == True: data_dicts["GLOBAL"].append(overlap_global_population)
        else: data_dicts["GLOBAL"].append(False)
        data_dicts["id"].append(id_)


        if omim_available == True: omim_raw = omim_raw + 1
        if omim_available == True and hgvsc == True: omim_and_hgvcs = omim_and_hgvcs +1
        if omim_available == True and (hgvsc == True or overlap_global_population == True): omim_and_hgvcs_and_global = omim_and_hgvcs_and_global +1
        if omim_available == True and (hgvsc == True or overlap_global_population == True or mena_hit == True): omim_and_hgvcs_or_global_or_mena = omim_and_hgvcs_or_global_or_mena +1
            
        if omim_available == True and mena_hit == True: omim_and_mena = omim_and_mena + 1
        if omim_available == True and overlap_global_population == True: omim_and_global = omim_and_global + 1
            
            
        
        if pd.isnull(gene_data) == False: med_raw = med_raw + 1
        if pd.isnull(gene_data) == False and hgvsc == True:  med_and_hgvcs = med_and_hgvcs + 1
        if pd.isnull(gene_data) == False and (hgvsc == True or overlap_global_population == True) == True:  med_and_hgvcs_and_global = med_and_hgvcs_and_global + 1
        if pd.isnull(gene_data) == False and (hgvsc == True or overlap_global_population == True or mena_hit == True): med_and_hgvcs_or_global_or_mena = med_and_hgvcs_or_global_or_mena +1
            
        if pd.isnull(gene_data) == False and mena_hit == True: med_and_mena = med_and_mena + 1
        if pd.isnull(gene_data) == False and overlap_global_population == True: med_and_global = med_and_global + 1


        if pd.isnull(gene_data) == False and  omim_available == True: omim_and_med = omim_and_med + 1
        if pd.isnull(gene_data) == False and  omim_available == True and hgvsc == True:  omim_and_med_and_hgvcs = omim_and_med_and_hgvcs + 1
        if pd.isnull(gene_data) == False and  omim_available == True and (hgvsc == True or overlap_global_population == True) == True:  omim_and_med_and_hgvcs_and_global = omim_and_med_and_hgvcs_and_global + 1
        if pd.isnull(gene_data) == False and  omim_available == True and (hgvsc == True or overlap_global_population == True or mena_hit == True): omim_and_med_and_hgvcs_or_global_or_mena = omim_and_med_and_hgvcs_or_global_or_mena +1
            
        if omim_available == True and pd.isnull(gene_data) == False and mena_hit == True: omim_and_med_and_mena = omim_and_med_and_mena + 1
        if omim_available == True and pd.isnull(gene_data) == False and overlap_global_population == True: omim_and_med_and_global = omim_and_med_and_global + 1
            
            
        
    counts = [[omim_raw, omim_and_hgvcs, omim_and_hgvcs_and_global, omim_and_hgvcs_or_global_or_mena, omim_and_mena, omim_and_global],[med_raw, med_and_hgvcs, med_and_hgvcs_and_global, med_and_hgvcs_or_global_or_mena, med_and_mena, med_and_global],[omim_and_med, omim_and_med_and_hgvcs, omim_and_med_and_hgvcs_and_global, omim_and_med_and_hgvcs_or_global_or_mena, omim_and_med_and_mena, omim_and_med_and_global]]
        
    summary =  pd.DataFrame({
        f"total # of patient SVs: {len(df)}": ["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": [counts[0][0], counts[0][0] - counts[0][4], counts[0][0] - counts[0][5], counts[0][0] - counts[0][1],  counts[0][0] - counts[0][2], counts[0][0] - counts[0][3], (counts[2][3])/counts[0][0]],
        "medical genes": [counts[1][0], counts[1][0] - counts[1][4], counts[1][0] - counts[1][5], counts[1][0] - counts[1][1], counts[1][0] - counts[1][2], counts[1][0] - counts[1][3], (counts[1][3])/counts[1][0]],
        "OMIM x medical gene" : [counts[2][0], counts[2][0] - counts[2][4], counts[2][0] - counts[2][5], counts[2][0] - counts[2][1], counts[2][0] - counts[2][2], counts[2][0] - counts[2][3], (counts[2][3])/counts[2][0]],
        "rate of reduction from OMIM": [1-(counts[2][0]/counts[0][0]), 1-((counts[2][0] - counts[2][4])/(counts[0][0] - counts[0][4])), 1-((counts[2][0] - counts[2][5])/(counts[0][0] - counts[0][5])), 1-((counts[2][0] - counts[2][1])/(counts[0][0] - counts[0][1])), 1-((counts[2][0] - counts[2][2])/(counts[0][0] - counts[0][2])),  1-((counts[2][0] - counts[2][3])/(counts[0][0] - counts[0][3])), ""]

    })
    
    return df, summary, counts, data_dicts