<jupyter_start><jupyter_text>Performance assessment Determine binary classification statistics For the v3 site filter, vs:- phase 2 accessibility map- permissive (100% pass, 0 fail)based on:- Mendelian error in all crosses (autosomes)- male het calls (`n_male_het`, `all_males_called`) will be used instead of `mendel_error`<jupyter_code># Imports:
import gcsfs
import zarr
import dask.array as da
from dask.distributed import Client
import pandas as pd
pd.set_option('display.float_format', '{:.3E}'.format)
import numpy as np
import numba
from pathlib import Path
!pip install -q malariagen-data
import malariagen_data
from dask_kubernetes import KubeCluster
from dask.distributed import Client
import ag3, ag2
cluster = KubeCluster(n_workers=40)
cluster
client = Client(cluster)
client
# Settings:
cross_family_ids = ['18-5', '29-2', '36-9', '37-3', '42-4', '45-1', '46-9', '47-6', '73-2', '78-2', '80-2', 'B5', 'K2', 'K4', 'K6']
chrom_arms = ['2R', '2L', '3R', '3L']
eval_family_ids = ['B5', 'K6', '73-2', '45-1', '37-3']
cross_family_ids<jupyter_output><empty_output><jupyter_text>Determine which crosses stats (inc. Mendelian error) are predicted by each site filterPopulate the confusion matrix, such that:- If the crosses stats imply that a position is accessible ("good") and the site filter also says it passes, then mark the prediction as a `True Positive` (`TP`)- If the crosses stats imply that a position is accessible ("good") but the site filter says it fails, then mark the prediction as a `False Negative` (`FN`)- If the crosses stats imply that a position is inaccessible ("bad") and the site filter also says it fails, then mark the prediction as a `True Negative` (`TN`)- If the crosses stats imply that a position is inaccessible ("bad") but the site filter says it passes, then mark the prediction as a `False Positive` (`FP`)- If the crosses stats cannot determine whether a position is "good" or "bad" (e.g. there were no male het calls on the X chromosome, but not all samples were called), then mark the prediction result as undetermined and exclude from subsequent stats.<jupyter_code>@numba.njit
def cross_tab_values(eval_arr, test_arr):
    
    out = np.zeros((2, 2))
    
    # TN FP
    # FN TP
    
    # eval arr has 
    for i in range(eval_arr.shape[0]):
        if eval_arr[i] == -1:
            continue
        
        if test_arr[i]:
            out[eval_arr[i], 1] += 1
        else:
            out[eval_arr[i], 0] += 1
    
    return out.reshape((1, 2, 2))
ag3_release_data = ag3.release_data()

def load_phase3_accessibility(chrom, mask_id='gamb_colu'):
    return ag3_release_data.load_mask(chrom, mask_id=mask_id)
ag2_release_data = ag2.release_data()

def load_phase2_accessibility(chrom):
    
    phase2_is_accessible = ag2_release_data.load_is_accessible(chrom)
    filter_n = ag2_release_data.load_filter_n(chrom)
    phase2_is_accessible_nonN = da.compress(~filter_n, phase2_is_accessible, axis=0)

    # Compute chunk sizes avoids this error when concatenating arrays: ValueError: Arrays chunk sizes are unknown: (nan,)
    phase2_is_accessible_nonN.compute_chunk_sizes() 
    return phase2_is_accessible_nonN
def load_ground_truth(chrom, xid):
    
    if chrom == "X":
        raise ValueError("X Not supported as hemizygous")
    
    me = ag3_release_data.load_crosses(chrom, xid, 'mendel_error')
    called = ag3_release_data.load_crosses(chrom, xid, 'n_samples_called')

    training_negative = ag3_release_data.load_mask(chrom, mask_id='gamb_colu', field='training_negative')
    training_positive = ag3_release_data.load_mask(chrom, mask_id='gamb_colu', field='training_positive')
    training_mask = training_negative | training_positive

    # Mark training sites as undetermined (-1)
    return da.where(training_mask, -1, da.where(me > 0, 0, da.where(called == da.max(called), 1, -1)))

def calculate_metrics_of_interest(true_neg, false_pos, false_neg, true_pos):
    
    r = {}
    r["tpr"] = true_pos / (true_pos + false_neg)
    r["fpr"] = false_pos / (false_pos + true_neg)
    r["fdr"] = false_pos / (true_pos + false_pos)
    r["tnr"] = 1 - (false_pos / (false_pos + true_neg))
    r["youden"] = r["tpr"] + r["tnr"] - 1
    r["frac_pass"] = (true_pos + false_pos) / (true_pos + false_pos + true_neg + false_neg)
    
    return pd.Series(r)
chrom_arms
cross_family_ids
desired_chunks = (500_000,)
%%time
data_holder = {}

for chrom_arm in chrom_arms:
        
    # varies by chrom only, but need to be chunked depending on ground_truth
    alternative_calls = {
        "phase2": load_phase2_accessibility(chrom_arm).rechunk(desired_chunks),
        "phase3": load_phase3_accessibility(chrom_arm).rechunk(desired_chunks)
    }
    
    alternative_calls["r100"] = da.ones(alternative_calls["phase2"].shape[0], dtype=np.bool, chunks=desired_chunks)
    
    chunks = ((1, ) * alternative_calls["phase2"].numblocks[0], 2, 2)
#     x = da.map_blocks(
#         cross_tab_values, 
#         alternative_calls["phase2"].astype("int"), 
#         alternative_calls["phase3"], 
#         chunks=chunks,
#         dtype=np.int32,
#         new_axis=[1, 2]).sum(axis=0).compute()

    for cross_id in cross_family_ids:
        
        print("processing", chrom_arm, cross_id)

        # varies by cross and chrom
        ground_truth = load_ground_truth(chrom_arm, cross_id).rechunk(desired_chunks)


        for key, eval_data in alternative_calls.items():
            
            if (key, chrom_arm, cross_id) in data_holder:
                print("skipping", (key, chrom_arm, cross_id))
                continue

            chunks = ((1, ) * ground_truth.numblocks[0], 2, 2)
            q = da.map_blocks(
                cross_tab_values, 
                ground_truth, 
                eval_data, 
                chunks=chunks, 
                dtype=np.int32, 
                new_axis=[1, 2]).sum(axis=0)

            (tn, fp), (fn, tp) = q.compute()

            # frac accessible is simply TPs + FPs. / all
            #res["frac_accessible"] = eval_data.mean().compute()

            data_holder[key, chrom_arm, cross_id] = pd.Series(
                [tn, fp, fn, tp], dtype=np.int32, index=["TN", "FP", "FN", "TP"])

df = pd.concat(data_holder, axis=0, sort=False, names=["description", "chrom", "cross_id", "metric"])
df.name = "value"<jupyter_output><empty_output><jupyter_text>List of tables1. Autosomes: Phase 3 vs unfiltered. FDR, TNR, ie 4 cols x (5 x 4) rows2. X: Phase 3 vs unfiltered3. Autosomes: Phase 2 vs phase 3. FDR / % accessible / Youden.4. X: Phase 2 va phase 35. Accessibility summaryAccessibility is shifted to separate section, where we can pull arab + gamb_colu_arab.<jupyter_code>output = pd.pivot_table(
    pd.DataFrame(df).reset_index(level="metric"), 
    index=["description", "chrom", "cross_id"], 
    columns="metric", 
    values="value")
output
output_all_chroms = output.groupby(level=[0, 2]).agg(np.sum)
output_summary = output.apply(lambda y: calculate_metrics_of_interest(*y), axis=1)
output_all_chroms_summary = output_all_chroms.apply(lambda y: calculate_metrics_of_interest(*y), axis=1)<jupyter_output><empty_output><jupyter_text>generate table 1.<jupyter_code>f = output_summary.reorder_levels([2, 0, 1]).loc[eval_family_ids]
f = f.drop("phase2", level=1)

f.rename({
    "fdr": "FalseDiscoveryRate", 
    "tpr": "TruePositiveRate"}, axis=1, inplace=True)
f.rename({"phase3": "site_filters", "r100": "unfiltered"}, level=1, inplace=True)

tab1 = f.reset_index().pivot(
    index=["chrom", "cross_id"], 
    columns=["description"], 
    values=["FalseDiscoveryRate", "TruePositiveRate"])
tab1.columns = [".".join(s) for s in tab1.columns.to_flat_index()]<jupyter_output><empty_output><jupyter_text>generate table 3<jupyter_code>f = output_summary.reorder_levels([2, 0, 1]).loc[eval_family_ids]
f = f.drop("r100", level=1)

f.rename({
    "fdr": "FalseDiscoveryRate", 
    "tpr": "TruePositiveRate",
    "youden": "Youden_J"}, axis=1, inplace=True)
f.rename({"phase3": "Phase3_filters", "phase2": "Phase2_filters"}, level=1, inplace=True)

tab3 = f.reset_index().pivot(
    index=["chrom", "cross_id"], 
    columns=["description"], 
    values=["FalseDiscoveryRate", "TruePositiveRate", "Youden_J"])


tab3.columns = [".".join(s) for s in tab3.columns.to_flat_index()]
tab1.head()
tab1.to_csv(
    "../content/tables/site-filters/table1_performance.csv", float_format='%.4E')
tab3.to_csv(
    "../content/tables/site-filters/table2_phase2vsphase3.csv", float_format='%.4E')<jupyter_output><empty_output><jupyter_text>Numbers for paper<jupyter_code>pd.set_option('display.float_format', '{:.5f}'.format)
output_all_chroms_summary.loc["phase3"].loc[eval_family_ids].apply([np.max, np.min]) * 100
output_all_chroms_summary.loc["r100"].loc[eval_family_ids].apply([np.max, np.min]) * 100
output_all_chroms_summary.loc["r100"].loc[eval_family_ids].fdr / output_all_chroms_summary.loc["phase3"].loc[eval_family_ids].fdr
output_all_chroms_summary.loc["phase2"].loc[eval_family_ids]<jupyter_output><empty_output><jupyter_text>Check fraction accessible matches up<jupyter_code>chromosomes = chrom_arms + ["X",]
contig_shape = pd.DataFrame(
    0,
    index=chromosomes, 
    columns=["phase2", "phase3:gamb_colu", "phase3:arab", "phase3:gamb_colu_arab"],
    dtype=int)

contig_sum_accessible = pd.DataFrame(
    0,
    index=chromosomes, 
    columns=["phase2", "phase3:gamb_colu", "phase3:arab", "phase3:gamb_colu_arab"],
    dtype=int)
for c in chromosomes:

    for mid in ["arab", "gamb_colu", "gamb_colu_arab"]:
        a = load_phase3_accessibility(c, mid)
        contig_sum_accessible.at[c, f"phase3:{mid}"] = a.sum().compute()
        contig_shape.at[c, f"phase3:{mid}"] = a.shape[0]
        
    a = load_phase2_accessibility(c)
    contig_sum_accessible.at[c, "phase2"] = a.sum().compute()
    contig_shape.at[c, "phase2"] = a.shape[0]
cluster.adapt()
accessibility_sumz = contig_sum_accessible / contig_shape
accessibility_sumz
accessibility_sumz.to_csv("../content/tables/accessibility/table3_summary.csv")<jupyter_output><empty_output><jupyter_text>Phase 2 vs phase 3<jupyter_code>p2stats = output_all_chroms_summary.loc["phase2"].loc[eval_family_ids]
p2stats
p3stats = output_all_chroms_summary.loc["phase3"].loc[eval_family_ids]
p3stats
p2stats.fdr - p3stats.fdr<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13467687239608564, 'fraction_numerical': 0.023641120046507123, 'mean_word_length': 4.277096114519427, 'pattern_counts': {'":': 12, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '37877', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3793, 'n_tokens_neox': 3574, 'n_words': 802, 'path': 'notebooks/site-filters_autosome-evaluation.ipynb', 'repo_name': 'malariagen/ag1000g-phase3-data-paper', 'size': '10321'}