<jupyter_start><jupyter_text>Find principal component pairs with a good enough separation line between `maternal` and `paternal` haplotypes:<jupyter_code>N_to_best_pair_and_info_dict = {}

for N in ['All', 'XY'] + list(range(1, 23)):
    print(N)
    HG002_df = pd.read_csv(
        os.path.join(dir_base, f'CMS.LongAlignments.NoUtgs.NoAltExceptPeregrine.chr{N}.PCA.tsv'),
        sep = '\t',
        index_col=0
    )

    #HG002_df = HG002_df.replace('maternal', '#FF0000')
    #HG002_df = HG002_df.replace('paternal', '#00FF00')

    pair_and_info_list = []

    for i, j in itertools.combinations(range(0, HG002_df.shape[0]), 2):
        X = HG002_df.iloc[:, [i, j]].values
        y = HG002_df.Haplotype.values
        
        clf = SVC(kernel='linear', C=1.0).fit(X, y)
        
        score = clf.score(X, y)
        print('\t', i, j, score)
        
        w = clf.coef_[0]
        slope = -w[0] / w[1]
        intercept = -clf.intercept_[0] / w[1]
        
        pair_and_info_list.append([(i, j), score, slope, intercept])
        
        if score == 1:
            break

    pair_and_info_list = sorted(pair_and_info_list, key=lambda x: x[1], reverse=True)
    
    N_to_best_pair_and_info_dict[N] = pair_and_info_list[0]
N_to_best_pair_and_info_dict<jupyter_output><empty_output><jupyter_text>Write the results:<jupyter_code>import os
base_dir = '/home/guarracino/Downloads/Pangenomics/HG002_bakeoff'
with open(os.path.join(base_dir, 'BestPairsMatPatSeparationInPCA.tsv'), 'w') as fw:
    fw.write('\t'.join(['N', 'PCx', 'PCy', 'slope', 'intercept']) + '\n')
    
    for N, ((PCx, PCy), score, slope, intercept) in N_to_best_pair_and_info_dict.items():
        # PC + 1 because in R array indexes start from 1
        fw.write('\t'.join([str(N), str(PCx + 1), str(PCy + 1), str(slope), str(intercept)]) + '\n')<jupyter_output><empty_output><jupyter_text>Debugging: plot the last (or one of the) solutions:<jupyter_code>def plot_svc_decision_function(model, ax=None, plot_support=True):
    """Plot the decision function for a 2D SVC"""
    if ax is None:
        ax = plt.gca()
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    
    # create grid to evaluate model
    x = np.linspace(xlim[0], xlim[1], 30)
    y = np.linspace(ylim[0], ylim[1], 30)
    Y, X = np.meshgrid(y, x)
    xy = np.vstack([X.ravel(), Y.ravel()]).T
    P = model.decision_function(xy).reshape(X.shape)
    
    # plot decision boundary and margins
    ax.contour(X, Y, P, colors='k',
               levels=[-1, 0, 1], alpha=0.5,
               linestyles=['--', '-', '--'])
    
    # plot support vectors
    if plot_support:
        ax.scatter(model.support_vectors_[:, 0],
                   model.support_vectors_[:, 1],
                   s=100, linewidth=1, facecolors='none',edgecolors='black');
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)


N = '22'
i, j = 0, 8
    
HG002_df = pd.read_csv(
    f'/home/guarracino/Downloads/Pangenomics/HG002_bakeoff/CMS.LongAlignments.NoUtgs.NoAltExceptPeregrine.chr{N}.PCA.tsv',
    sep = '\t',
    index_col=0
)
HG002_df = HG002_df.replace('maternal', '#FF0000')
HG002_df = HG002_df.replace('paternal', '#00FF00')

X = HG002_df.iloc[:, [i, j]].values
y = HG002_df.Haplotype.values

clf = SVC(kernel='linear', C=5).fit(X, y)
score = clf.score(X, y)
print('\t', i, j, score)

plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plot_svc_decision_function(clf);
plt.xlabel(f"PC {i+1}");
plt.ylabel(f"PC {j+1}");
plt.title(f'No Utgs and No alt except Peregrine - chr {N}')<jupyter_output>0 8 1.0
{'fraction_non_alphanumeric': 0.1662425784563189, 'fraction_numerical': 0.03505795872208086, 'mean_word_length': 3.191943127962085, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '379080', 'license': 'MIT', 'n_tokens_mistral': 1537, 'n_tokens_neox': 1421, 'n_words': 268, 'path': 'scripts/HG002_bakeoff_SVM_on_PCA.ipynb', 'repo_name': 'AndreaGuarracino/HG002_assemblies_assessment', 'size': '3537'}