<jupyter_start><jupyter_text>Betabinomial mixture model for detecting mitochondria variants<jupyter_code>import bbmix
import numpy as np
from scipy import sparse
from scipy.io import mmread
import matplotlib.pyplot as plt
AD = mmread("../data/mitoDNA/cellSNP.tag.AD.mtx").tocsc().toarray()
DP = mmread("../data/mitoDNA/cellSNP.tag.DP.mtx").tocsc().toarray()
AD.shape, DP.shape<jupyter_output><empty_output><jupyter_text>Run on all variants<jupyter_code>from bbmix.models import MixtureBetaBinomial

model1 = MixtureBetaBinomial(n_components = 1, max_m_step_iter=250, tor=1e-10, n_init_searches=10)
model2 = MixtureBetaBinomial(n_components = 2, max_m_step_iter=250, tor=1e-10, n_init_searches=20)

for i in range(AD.shape[0]):
    a = AD[i, :]
    d = DP[i, :]
    idx = d >= 10
    a, d = a[idx], d[idx]
    
    params1 = model1.fit((a, d), max_iters=500, init_method="mixbin", early_stop=False)
    params2 = model2.fit((a, d), max_iters=1000, init_method="mixbin", early_stop=False)
    p_val = bbmix.models.LR_test(model1.losses[-1] - model2.losses[-1], df = 3)
    print(i, "mode1: %.2f\tmodel2:%.2f\tp: %.3e" %(model1.losses[-1], model2.losses[-1], p_val))<jupyter_output>0 mode1: 3762.45	model2:3760.84	p: 3.590e-01
1 mode1: 158787.75	model2:158755.41	p: 5.863e-14
2 mode1: 6631.89	model2:6620.92	p: 6.686e-05
3 mode1: 18800.42	model2:18800.36	p: 9.902e-01
4 mode1: 85371.67	model2:85333.99	p: 3.053e-16
5 mode1: 18068.62	model2:18022.85	p: 1.026e-19
6 mode1: 10041.48	model2:9984.59	p: 1.693e-24
7 mode1: 4618.31	model2:4617.40	p: 6.076e-01
8 mode1: 4436.28	model2:4435.06	p: 4.866e-01<jupyter_text>Test variant 1<jupyter_code>from bbmix.models import MixtureBetaBinomial

model1 = MixtureBetaBinomial(n_components = 1, max_m_step_iter=500, tor=1e-10, n_init_searches=10)
model2 = MixtureBetaBinomial(n_components = 2, max_m_step_iter=500, tor=1e-10, n_init_searches=20)
i = 1
idx = DP[i, :] >= 0
a = AD[i, idx]
d = DP[i, idx]
a, d
# Multiple initializations

for ite in range(20):
    params1 = model1.fit((a, d), max_iters=500, init_method="mixbin", early_stop=False)
    params2 = model2.fit((a, d), max_iters=1000, init_method="mixbin", early_stop=False)

    p_val = bbmix.models.LR_test(model1.losses[-1] - model2.losses[-1], df = 3)
    print(i, "mode1: %.2f\tmodel2:%.2f\tp: %.3e" %(model1.losses[-1], model2.losses[-1], p_val))<jupyter_output>1 mode1: 158787.77	model2:158764.79	p: 5.795e-10
1 mode1: 158787.92	model2:158762.15	p: 3.775e-11
1 mode1: 158787.75	model2:158751.25	p: 9.699e-16
1 mode1: 158787.69	model2:158771.14	p: 3.064e-07
1 mode1: 158787.71	model2:158767.59	p: 9.424e-09
1 mode1: 158787.69	model2:158763.23	p: 1.354e-10
1 mode1: 158787.87	model2:158761.93	p: 3.155e-11
1 mode1: 158787.82	model2:158763.39	p: 1.398e-10
1 mode1: 158788.64	model2:158754.37	p: 8.819e-15
1 mode1: 158788.57	model2:158764.95	p: 3.101e-10
1 mode1: 158787.69	model2:158767.37	p: 7.747e-09
1 mode1: 158788.85	model2:158767.48	p: 2.797e-09
1 mode1: 158787.71	model2:158752.80	p: 4.671e-15
1 mode1: 158787.69	model2:158764.91	p: 7.015e-10
1 mode1: 158788.51	model2:158753.06	p: 2.731e-15
1 mode1: 158787.84	model2:158759.74	p: 3.816e-12
1 mode1: 158787.90	model2:158761.80	p: 2.724e-11
1 mode1: 158787.70	model2:158764.93	p: 7.126e-10
1 mode1: 158788.32	model2:158768.70	p: 1.552e-08
1 mode1: 158788.11	model2:158766.92	p: 3.345e-09<jupyter_text>Calibration of p value - mixture of betabinomial<jupyter_code>from bbmix.models import MixtureBetaBinomial

model1 = MixtureBetaBinomial(n_components = 1, max_m_step_iter=250, tor=1e-10, n_init_searches=10)
model2 = MixtureBetaBinomial(n_components = 2, max_m_step_iter=250, tor=1e-10, n_init_searches=20)

model0 = MixtureBetaBinomial(n_components = 1, max_m_step_iter=250, tor=1e-10, n_init_searches=10)
params0 = model0.fit((a, d), max_iters=500, init_method="mixbin", early_stop=True)
n_iter = 100
p_vals = np.ones(n_iter)

for ite in range(n_iter):
    a_sim = model0.sample(d)
    params1 = model1.fit((a_sim, d), max_iters=500, init_method="mixbin", early_stop=False)
    params2 = model2.fit((a_sim, d), max_iters=500, init_method="mixbin", early_stop=False)
    p_vals[ite] = bbmix.models.LR_test(model1.losses[-1] - model2.losses[-1], df = 3)
plt.plot(-np.log10(np.linspace(0, 1, len(p_vals)+1)[1:]), 
         -np.sort(np.log10(p_vals)), 'o')
plt.plot([0, np.log10(len(p_vals))], 
         [0, np.log10(len(p_vals))], '--', color='grey')
plt.xlabel("-log10(p value), expected")
plt.ylabel("-log10(p value), observed")
plt.title("QQ plot from Null model")
c = np.zeros((3,3))
c[c == 0]<jupyter_output><empty_output><jupyter_text>Calibration of p value - mixture of binomial<jupyter_code>from bbmix.models import MixtureBinomial

model1 = MixtureBinomial(n_components = 1, tor=1e-20)
model2 = MixtureBinomial(n_components = 2, tor=1e-20)

model0 = MixtureBinomial(n_components = 1, tor=1e-20)
params0 = model0.fit((a, d), max_iters=500, early_stop=True)
n_iter = 100
p_vals = np.ones(n_iter)

for ite in range(n_iter):
    a_sim = model0.sample(d)
    params1 = model1.fit((a_sim, d), max_iters=500, early_stop=True)
    params2 = model2.fit((a_sim, d), max_iters=500, early_stop=True)
    p_vals[ite] = bbmix.models.LR_test(model1.losses[-1] - model2.losses[-1], df = 3)
model1.losses[-1], model2.losses[-1]
plt.plot(-np.log10(np.linspace(0, 1, len(p_vals)+1)[1:]), 
         -np.sort(np.log10(p_vals)), 'o')
plt.plot([0, np.log10(len(p_vals))], 
         [0, np.log10(len(p_vals))], '--', color='grey')
plt.xlabel("-log10(p value), expected")
plt.ylabel("-log10(p value), observed")
plt.title("QQ plot from Null model")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1612041325258283, 'fraction_numerical': 0.16636978981118633, 'mean_word_length': 4.973404255319149, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '83042', 'license': 'Apache-2.0', 'n_tokens_mistral': 3096, 'n_tokens_neox': 2576, 'n_words': 478, 'path': 'examples/Mito_mutations.ipynb', 'repo_name': 'StatBiomed/BBMix', 'size': '5614'}