Add SER and PER #255

Author: XingerTang

src/utils.py

@@ -0,0 +1,68 @@
+import numpy as np
+
+
+def calc_switch_error_rate():
+    """
+    Run this function at root directory and input the path to the
+    assessed haplotype file.
+    The SER and PER_intra calculation satisfies the definition from:
+    https://www.cell.com/hgg-advances/fulltext/S2666-2477(25)00082-X
+    """
+    nLociAll = 2000
+    true_path = "tests/accuracy_tests/sim_for_alphapeel_accu_test/true-hap_0.5.txt"
+    file_name = input("Enter the path of the assessed haplotype file: ")
+    new_file = np.loadtxt(file_name, usecols=np.arange(1, nLociAll + 1))
+    true_file = np.loadtxt(true_path, usecols=np.arange(1, nLociAll + 1))
+
+    nInd = 1000
+    switch_error_count = 0
+    phase_error_count = 0
+    uncalled_count = 0
+    wrong_homo_count = 0
+    true_hetero_count = 0
+    homo_count = 0
+    hetero_count = 0
+    for ind in range(nInd):
+        hap_p_new = new_file[ind * 2]
+        hap_m_new = new_file[ind * 2 + 1]
+        hap_p_true = true_file[ind * 2]
+        hap_m_true = true_file[ind * 2 + 1]
+
+        switched = False
+        for loci in range(nLociAll):
+            if hap_p_true[loci] == hap_m_true[loci]:
+                homo_count += 1
+            else:
+                hetero_count += 1
+            if hap_p_new[loci] == 9 or hap_m_new[loci] == 9:
+                uncalled_count += 1
+                continue
+            if hap_p_true[loci] + hap_m_true[loci] == 1:
+                if hap_p_new[loci] == hap_m_new[loci]:
+                    wrong_homo_count += 1
+                    continue
+                true_hetero_count += 1
+                if hap_p_new[loci] != hap_p_true[loci]:
+                    phase_error_count += 1
+                if (hap_p_new[loci] != hap_p_true[loci] and switched is False) or (
+                    hap_p_new[loci] == hap_p_true[loci] and switched is True
+                ):
+                    switched = not switched
+                    switch_error_count += 1
+
+    print(f"Switch error rate: {switch_error_count / (nInd * (nLociAll - 1))}")
+    print(f"Phase error (intra) rate: {phase_error_count / (nInd * nLociAll)}")
+    print(f"Uncalled rate: {uncalled_count / (nInd * nLociAll)}")
+    print(f"Wrong homozygote rate: {wrong_homo_count / (nInd * nLociAll)}")
+    print(f"True heterozygote rate: {true_hetero_count / (nInd * nLociAll)}")
+    print(f"Homozygote count in true genotype: {homo_count}")
+    print(f"Heterozygote count in true genotype: {hetero_count}")
+    print(f"Homo to hetero ratio: {homo_count / hetero_count}")
+
+
+def main():
+    calc_switch_error_rate()
+
+
+if __name__ == "__main__":
+    main()