<jupyter_start><jupyter_text>Compare similarity search strategies (OLD)In this notebook I am comparing the following similarity search strategies for finding genomic regions with similar peaks:1. kNN search in the autoencoded latent space(1b. kNN search in the UMAP 2D embedding of the autoencoded latent space)2. kNN search in the Euclidean space of the original data3. kNN search in the DTW space of the original data4. kNN search in the 2D UMAP embedding of the original data5. Template-matching with normalized cross correlation<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline
import numpy as np
import os
import sys
import warnings

# Ignore warnings as they just pollute the output
warnings.filterwarnings('ignore')

# Enable importing modules from the parent directory
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
module_path = os.path.abspath(os.path.join('../experiments'))
if module_path not in sys.path:
    sys.path.append(module_path)

# DNase-seq 2011, hg19
bw = 'data/ENCFF158GBQ.bigWig'<jupyter_output><empty_output><jupyter_text>Load and normalize the data Download the data<jupyter_code>from download import download_encode_file
from pathlib import Path

Path('data').mkdir(parents=True, exist_ok=True) 

download_encode_file('ENCFF158GBQ.bigWig')<jupyter_output>File already exist. To overwrite pass `overwrite=True`<jupyter_text>Load models, extract genomic windows, and normalize data<jupyter_code># Load model and data and normalize data

from keras.models import load_model
from server import bigwig

encoder_2_4kb = load_model('../../examples/encode-gm12878-dnase-seq-2011-encoder-2_4kb.h5')
decoder_2_4kb = load_model('../../examples/encode-gm12878-dnase-seq-2011-decoder-2_4kb.h5')

data_2_4kb = bigwig.chunk(
    bw,
    2400,
    100,
    2400 / 2,
    ['chr22'],
    verbose=True,
)
data_2_4kb_b = bigwig.chunk(
    bw,
    2400,
    100,
    2400 / 2,
    ['chr1', 'chr2', 'chr3'],
    verbose=True,
)<jupyter_output>Extracted 207708 windows from chr1 with a max value of 1.0.
Extracted 202666 windows from chr2 with a max value of 1.0.
Extracted 165018 windows from chr3 with a max value of 1.0.<jupyter_text>Encode data<jupyter_code>from ae.utils import predict

predicted_2_4kb, _, latent_2_4kb = predict(
    encoder_2_4kb,
    decoder_2_4kb,
    data_2_4kb.reshape(data_2_4kb.shape[0], data_2_4kb.shape[1], 1)
)<jupyter_output><empty_output><jupyter_text>Helper methods Simple kNN search<jupyter_code>from scipy.spatial.distance import cdist

def knn(data, target_idx, k, metric='euclidean', sax = None, ignore: int = 0):
    """K nearest neighbors
    
    Find the `k` nearest neighbors of a 
    """
    
    target = data[target_idx]
    
    if sax is None:
        dist = cdist(data, target.reshape((1, target.size)), metric='euclidean').flatten()

    else:
        N = data.shape[0]
        dist = np.zeros(N)
        for i in range(N):
            dist[i] = sax.distance_sax(target, data[i])

    # Ensure that the target is always first
    dist[target_idx] = -1
    for i in range(1, ignore + 1):
        dist[min(target_idx + i, data.shape[0] - 1)] = -1
        dist[max(target_idx - i, 0)] = -1
    
    return np.argsort(dist)[1 + (2 * ignore):k + 1 + (2 * ignore)]<jupyter_output><empty_output><jupyter_text>DTW distance<jupyter_code>from scipy.spatial.distance import euclidean
from fastdtw import fastdtw
from multiprocessing import Pool

def dtw(data, target_idx: int, n: int, target = None, return_all = False, print_progress: bool = False):
    N = data.shape[0]
    dist = np.zeros(N)
    
    if target is None:
        target = data[target_idx]
        
    p = ((np.arange(4) + 1) * (N // 4)).astype(int)
    
    for i in np.arange(N):
        if i in p and print_progress:
            print('.', end='', flush=True)
        d, _ = fastdtw(data[i], target, dist=euclidean)
        dist[i] = d
    
    if return_all:
        return dist
    
    return np.argsort(dist)[:n]

def pooled_dtw(data, target_idx: int, n: int, print_progress: bool = False, ignore: int = 0):
    target = data[target_idx]
    is_not_target = np.arange(data.shape[0]) != target_idx
    for i in range(0, ignore + 1):
        is_not_target[min(target_idx + i, data.shape[0] - 1)] = False
        is_not_target[max(target_idx - i, 0)] = False
    
    with Pool() as pool:
        args = [[d, 0, n, target, True, print_progress] for d in np.array_split(data[is_not_target], pool._processes)]
        dist = np.concatenate(pool.starmap(dtw, args))
        
        return np.argsort(dist)[:n]<jupyter_output><empty_output><jupyter_text>Normalized cross correlation search<jupyter_code>from scipy.signal import correlate

def norm(data, zero_norm: bool = False):
    mean = np.mean(data) if zero_norm else 0
    
    return (data - mean) / np.std(data)

def norm2d(data, zero_norm: bool = False):
    mean = np.mean(data, axis=1).reshape(-1, 1) if zero_norm else np.zeros((data.shape[0], 1))
    std = np.std(data, axis=1).reshape(-1, 1)
    
    return (data - mean) / std

def xcorrelation(data, template_idx, n, normalize=False, zero_normalize=False, ignore: int = 0):
    unknown = data
    template = data[template_idx]
    
    if norm:
        unknown = norm2d(unknown, zero_norm=zero_normalize)
        template = norm(template, zero_norm=zero_normalize)
        
    xcorr = np.apply_along_axis(lambda m: correlate(m, template, mode='full'), axis=1, arr=unknown)
    xcorr[np.where(np.isnan(xcorr))] = 0

    max_xcorr = np.nanmax(xcorr, axis=1)
    
    # Ensure that the target is always last
    max_xcorr[template_idx] = -1
    for i in range(1, ignore + 1):
        max_xcorr[min(template_idx + i, data.shape[0] - 1)] = -1
        max_xcorr[max(template_idx - i, 0)] = -1
    
    return np.argsort(max_xcorr)[::-1][:n]<jupyter_output><empty_output><jupyter_text>Generate UMAP embeddings of the latent space and original data space<jupyter_code>import umap

umap_embedding_2_4kb = umap.UMAP(
    n_neighbors=10,
    min_dist=0.01,
    metric='l2',
).fit_transform(data_2_4kb)<jupyter_output><empty_output><jupyter_text>2.4kb Search<jupyter_code># Some genomic windows with interesting peaks
targets = [20196, 26570, 26576, 26790, 26791, 31670, 25513]
k = 5

knn_ae = []
knn_eq = []
knn_dtw = []
knn_umap = []
top_xcross = []

for target in targets:
    print('Search for window #{}'.format(target), end='', flush=True)
    knn_ae.append(knn(latent, target, k))
    print('.', end='', flush=True)
    knn_eq.append(knn(data_test, target, k))
    print('.', end='', flush=True)
    knn_dtw.append(pooled_dtw(data_test, target, k))
    print('.', end='', flush=True)
    knn_umap.append(knn(umap_embedding, target, k))
    print('.', end='', flush=True)
    top_xcross.append(xcorrelation(data_test, target, k, normalize=True, zero_normalize=True))
    print('. done!')<jupyter_output>Search for window #20196..Computing the DTW distance for 42753 windows took 287.4487359523773
... done!
Search for window #26570..Computing the DTW distance for 42753 windows took 279.80624294281006
... done!
Search for window #26576..Computing the DTW distance for 42753 windows took 282.0815329551697
... done!
Search for window #26790..Computing the DTW distance for 42753 windows took 279.5405559539795
... done!
Search for window #26791..Computing the DTW distance for 42753 windows took 276.1450231075287
... done!
Search for window #31670..Computing the DTW distance for 42753 windows took 276.2667832374573
... done!
Search for window #25513..Computing the DTW distance for 42753 windows took 271.9153709411621
... done!<jupyter_text>2.4kb visual comparison of search strategiesCompare encoding-based knn similarity search strategies against normalize cross correlation.| Rows  | Description                                      ||-------|--------------------------------------------------|| 1     | Search target. Manually picked.                  || 2-6   | k-NN in the CAE-based latent space.              || 7-11  | k-NN in the Euclide space.                       || 11-16 | k-NN on the dynamic time warping space.          || 17-21 | k-NN on the UMAP embedded space.                 || 21-26 | Zero-normalized cross correlation.               |<jupyter_code>import matplotlib.pyplot as plt 
import numpy as np

N = (k + 1) * 6

T = len(targets)
sz = data_test[0].size

plt.figure(figsize=(5 * T, N))

ymax = 1.0

show_predictions = False

for i, target in enumerate(targets):
    ax = plt.subplot(N, T, (i + 1))
        
    ax.set_facecolor("#D9FAFF")
    
    if show_predictions:
        plt.bar(np.arange(sz), data_test[target], color='black', alpha=0.33)
        plt.bar(np.arange(sz), data_test[target], color='black', alpha=0.33)
    else:
        plt.bar(np.arange(sz), data_test[target], color='#004754')

    plt.ylim(0, ymax)
    plt.xticks([], [])
    plt.yticks([], [])

    for j, hit in enumerate(knn_ae[i]):
        plt.subplot(N, T, ((j + 1) * T) + (i + 1))
        if show_predictions:
            plt.bar(np.arange(sz), predicted[hit], color='green', alpha=0.33)
            plt.bar(np.arange(sz), data_test[hit], color='black', alpha=0.33)
        else:
            plt.bar(np.arange(sz), data_test[hit], color='#d24f00')
        plt.ylim(0, ymax)
        plt.xticks([], [])
        plt.yticks([], [])
        plt.subplots_adjust(top=0.9)

    for j, hit in enumerate(knn_eq[i]):
        plt.subplot(N, T, ((j + 6) * T) + (i + 1))
        if show_predictions:
            plt.bar(np.arange(sz), predicted[hit], color='green', alpha=0.33)
            plt.bar(np.arange(sz), data_test[hit], color='black', alpha=0.33)
        else:
            plt.bar(np.arange(sz), data_test[hit], color='#008ca8')
        plt.ylim(0, ymax)
        plt.xticks([], [])
        plt.yticks([], [])

    for j, hit in enumerate(knn_dtw[i]):
        plt.subplot(N, T, ((j + 11) * T) + (i + 1))
        if show_predictions:
            plt.bar(np.arange(sz), predicted[hit], color='red', alpha=0.33)
            plt.bar(np.arange(sz), data_test[hit], color='black', alpha=0.33)
        else:
            plt.bar(np.arange(sz), data_test[hit], color='#93209e')
        plt.ylim(0, ymax)
        plt.xticks([], [])
        plt.yticks([], [])

    for j, hit in enumerate(knn_umap[i]):
        plt.subplot(N, T, ((j + 16) * T) + (i + 1))
        if show_predictions:
            plt.bar(np.arange(sz), predicted[hit], color='yellow', alpha=0.33)
            plt.bar(np.arange(sz), data_test[hit], color='black')
        else:
            plt.bar(np.arange(sz), data_test[hit], color='#209e4e')
        plt.ylim(0, ymax)
        plt.xticks([], [])
        plt.yticks([], [])

    for j, hit in enumerate(top_xcross[i]):
        plt.subplot(N, T, ((j + 21) * T) + (i + 1))
        plt.bar(np.arange(sz), data_test[hit], color='#cc9900')
        plt.ylim(0, ymax)
        plt.xticks([], [])
        plt.yticks([], [])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13966633239128454, 'fraction_numerical': 0.04658583280153159, 'mean_word_length': 3.125611132004513, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '56076', 'license': 'Apache-2.0', 'n_tokens_mistral': 4236, 'n_tokens_neox': 3849, 'n_words': 961, 'path': 'experiments/notebooks/Compare similarity search strategies (old).ipynb', 'repo_name': 'flekschas/peax', 'size': '10969'}