<filename>src/cupcake/phasing/utils/evaluate_phase_switch.py
__author__ = "<EMAIL>"

"""
To be used to evaluate phase switching by assuming the genome is already fully "phased".
For example, in Falcon-Unzip, which is pseudo-phased, there could be switching errors.
Since IsoPhase only uses the primary contig (or one of the haplotypes), a switching error shows up
 by seeing isoforms that have SNPs on the same allele being switched between different phases.

ex: in SNP 1, isoform 1 has GT=0|1
    in SNP 2, isoform 1 has GT=1|0
    in SNP 3, isoform 1 has GT=0|1

this would be a possible switching error. This is two switching events. SNP1->2 one event, 2->3 one event.
For now we will only consider the diploid case. Fail as soon as we see more than 2 alleles.
We will use phased.nopartial.cleaned.vcf as the input VCF.

Output is:

   <directory>, <chrom>, <strand>, <start>, <end>, <num_isoforms>, <num_switch_event>

We would get the loci (chr:start-end strand) from the config file in each directory.
"""


import glob
from io import TextIOWrapper
from pathlib import Path
from typing import Tuple

# from collections import defaultdict
import vcfpy

from cupcake import cupcake_logger as logger


def read_config(config_filename: Path) -> Tuple[str, str, str, str]:
    """
    pbid=PB.9923
    ref_chr=000104F_0
    ref_strand=-
    ref_start=3484015
    ref_end=3519963

    :return: (chr, start, end, strand)
    """
    for line in open(config_filename, "r"):
        a, b = line.strip().split("=")
        if a == "ref_chr":
            _chr = b
        elif a == "ref_strand":
            _strand = b
        elif a == "ref_start":
            _start = b
        elif a == "ref_end":
            _end = b
    return _chr, _start, _end, _strand


def eval_isophase_phaseswitch(
    isophase_vcf: str, config_file: Path, out_f: TextIOWrapper, name: str = "NA"
) -> None:

    _chr, _start, _end, _strand = read_config(config_file)

    reader = vcfpy.Reader(isophase_vcf)
    # record the first SNP for each isoform
    prev = {}  # sample -> CallData.GT (ex: '0|1')
    r = next(reader)
    for c in r.samples:
        prev[c.sample] = c.data.GT

    num_switch = 0

    for r in reader:
        for c in r.samples:
            if c.data.GT.find("|") == -1:
                continue  # ignore those with just one allele
            a, b = c.data.GT.split("|")
            if a == b:
                continue  # for now, ignore IsoPhase results that only uses one allele
            if prev[c.sample] != c.data.GT:
                num_switch += 1
            prev[c.sample] = c.data.GT

    out_f.write(
        f"{name}\t{_chr}\t{_start}\t{_end}\t{_strand}\t{len(r.samples)}\t{num_switch}\n"
    )


def main_eval():

    with open("evaled.isophase_phase_switches.txt", "w") as out_f:

        out_f.write("dir\tchrom\tstart\tend\tstrand\tnum_iso\tnum_phase_switch\n")
        dirs = glob.glob("by_loci/*size*/")

        for d1 in dirs:
            configfile = Path(d1, "config")
            vcffile = Path(d1, "phased.nopartial.cleaned.vcf")
            nosnp = Path(d1, "phased.nopartial.NO_SNPS_FOUND")
            nohap = Path(d1, "phased.nopartial.NO_HAPS_FOUND")

            if not vcffile.exists():
                # no SNP, just skip
                assert nosnp.exists() or nohap.exists()
                logger.info(f"Skipping {d1} because no SNPs found.")
            else:
                logger.info(f"Evaluating {d1}.")
                name = d1.split("/")[1]
                eval_isophase_phaseswitch(vcffile, configfile, out_f, name)


main_eval()

<filename>src/cupcake/phasing/utils/evaluate_phase_switch.py
__author__ = "<EMAIL>"

"""
To be used to evaluate phase switching by assuming the genome is already fully "phased".
For example, in Falcon-Unzip, which is pseudo-phased, there could be switching errors.
Since IsoPhase only uses the primary contig (or one of the haplotypes), a switching error shows up
 by seeing isoforms that have SNPs on the same allele being switched between different phases.

ex: in SNP 1, isoform 1 has GT=0|1
    in SNP 2, isoform 1 has GT=1|0
    in SNP 3, isoform 1 has GT=0|1

this would be a possible switching error. This is two switching events. SNP1->2 one event, 2->3 one event.
For now we will only consider the diploid case. Fail as soon as we see more than 2 alleles.
We will use phased.nopartial.cleaned.vcf as the input VCF.

Output is:

   <directory>, <chrom>, <strand>, <start>, <end>, <num_isoforms>, <num_switch_event>

We would get the loci (chr:start-end strand) from the config file in each directory.
"""


import glob
from io import TextIOWrapper
from pathlib import Path
from typing import Tuple

# from collections import defaultdict
import vcfpy

from cupcake import cupcake_logger as logger


def read_config(config_filename: Path) -> Tuple[str, str, str, str]:
    """
    pbid=PB.9923
    ref_chr=000104F_0
    ref_strand=-
    ref_start=3484015
    ref_end=3519963

    :return: (chr, start, end, strand)
    """
    for line in open(config_filename, "r"):
        a, b = line.strip().split("=")
        if a == "ref_chr":
            _chr = b
        elif a == "ref_strand":
            _strand = b
        elif a == "ref_start":
            _start = b
        elif a == "ref_end":
            _end = b
    return _chr, _start, _end, _strand


def eval_isophase_phaseswitch(
    isophase_vcf: str, config_file: Path, out_f: TextIOWrapper, name: str = "NA"
) -> None:

    _chr, _start, _end, _strand = read_config(config_file)

    reader = vcfpy.Reader(isophase_vcf)
    # record the first SNP for each isoform
    prev = {}  # sample -> CallData.GT (ex: '0|1')
    r = next(reader)
    for c in r.samples:
        prev[c.sample] = c.data.GT

    num_switch = 0

    for r in reader:
        for c in r.samples:
            if c.data.GT.find("|") == -1:
                continue  # ignore those with just one allele
            a, b = c.data.GT.split("|")
            if a == b:
                continue  # for now, ignore IsoPhase results that only uses one allele
            if prev[c.sample] != c.data.GT:
                num_switch += 1
            prev[c.sample] = c.data.GT

    out_f.write(
        f"{name}\t{_chr}\t{_start}\t{_end}\t{_strand}\t{len(r.samples)}\t{num_switch}\n"
    )


def main_eval():

    with open("evaled.isophase_phase_switches.txt", "w") as out_f:

        out_f.write("dir\tchrom\tstart\tend\tstrand\tnum_iso\tnum_phase_switch\n")
        dirs = glob.glob("by_loci/*size*/")

        for d1 in dirs:
            configfile = Path(d1, "config")
            vcffile = Path(d1, "phased.nopartial.cleaned.vcf")
            nosnp = Path(d1, "phased.nopartial.NO_SNPS_FOUND")
            nohap = Path(d1, "phased.nopartial.NO_HAPS_FOUND")

            if not vcffile.exists():
                # no SNP, just skip
                assert nosnp.exists() or nohap.exists()
                logger.info(f"Skipping {d1} because no SNPs found.")
            else:
                logger.info(f"Evaluating {d1}.")
                name = d1.split("/")[1]
                eval_isophase_phaseswitch(vcffile, configfile, out_f, name)


main_eval()

To be used to evaluate phase switching by assuming the genome is already fully "phased". For example, in Falcon-Unzip, which is pseudo-phased, there could be switching errors. Since IsoPhase only uses the primary contig (or one of the haplotypes), a switching error shows up  by seeing isoforms that have SNPs on the same allele being switched between different phases.  ex: in SNP 1, isoform 1 has GT=0|1     in SNP 2, isoform 1 has GT=1|0     in SNP 3, isoform 1 has GT=0|1  this would be a possible switching error. This is two switching events. SNP1->2 one event, 2->3 one event. For now we will only consider the diploid case. Fail as soon as we see more than 2 alleles. We will use phased.nopartial.cleaned.vcf as the input VCF.  Output is:     <directory>, <chrom>, <strand>, <start>, <end>, <num_isoforms>, <num_switch_event>  We would get the loci (chr:start-end strand) from the config file in each directory. # from collections import defaultdict pbid=PB.9923     ref_chr=000104F_0     ref_strand=-     ref_start=3484015     ref_end=3519963      :return: (chr, start, end, strand) # record the first SNP for each isoform # sample -> CallData.GT (ex: '0|1') # ignore those with just one allele # for now, ignore IsoPhase results that only uses one allele # no SNP, just skip