<jupyter_start><jupyter_text>First we need to create the mapping between read ids and Fast5 files so that we can pair full length reads from the bam files with the signal that they originated from. This information is present in the sequencing summaries file:<jupyter_code>read_id_filemap = get_fast5_read_id_mapping(
    '/cluster/ggs_lab/mtparker/analysis_notebooks/fiveprime_adapter_identification/data/201804_col0_5adapterLIG/5prime_adapter_rebasecalled/',
    '/cluster/ggs_lab/mtparker/analysis_notebooks/fiveprime_adapter_identification/data/201804_col0_5adapterLIG/5prime_adapter_fast5s/'
)
positive_bam_fn = '/cluster/gjb_lab/nschurch/git/NS_Nanopore_paper/manuscript/supplementary/fulllength_data/201902_col0_2916_5adapter_exp2_fulllength.bam'
negative_bam_fn = '/cluster/gjb_lab/nschurch/git/NS_Nanopore_paper/manuscript/supplementary/fulllength_data/201902_col0_2916_5adapter_exp2_short_relaxed.bam'

with pysam.AlignmentFile(positive_bam_fn) as p:
    print(p.mapped)
with pysam.AlignmentFile(negative_bam_fn) as n:
    print(n.mapped)
keras<jupyter_output><empty_output><jupyter_text>Now we can extract the 5' end signals (last 3000 signal measurements) from each Fast5 file. We also extract a bunch of internal signals to use as negative examples as well<jupyter_code>def get_fiveprime_for_reads_in_bam(bam_fn, read_id_filemap, size=3000, internal=False):
    with pysam.AlignmentFile(bam_fn) as bam:
        signals = Parallel(n_jobs=16)(
            delayed(get_fast5_fiveprime)(
                r.query_name, read_id_filemap[r.query_name], size, internal)
            for r in bam.fetch()
        )
    if internal:
        read_ids, sig_lens, fiveprime_signals, internal_signals = zip(*signals)
        internal_signals = [x for x in internal_signals if x is not None and len(x) == size]
        return np.asarray(fiveprime_signals), np.asarray(internal_signals)
    else:
        read_ids, sig_lens, fiveprime_signals = zip(*signals)
        return np.asarray(fiveprime_signals)

pos_signals = get_fiveprime_for_reads_in_bam(
    positive_bam_fn, read_id_filemap, size=3000, internal=False)
neg_signals, neg_signals_internal = get_fiveprime_for_reads_in_bam(
    negative_bam_fn, read_id_filemap, size=3000, internal=True)
def write_to_hdf5(h5_fn, pos_signals, neg_signals, neg_signals_internal):
    with h5py.File(h5_fn, 'w') as f:
        f.create_dataset('pos_signals', data=pos_signals)
        f.create_dataset('neg_signals', data=neg_signals)
        f.create_dataset('neg_internal_signals', data=neg_signals_internal)
write_to_hdf5(
    'data/training_data_round1.h5',
    pos_signals, neg_signals, neg_signals_internal
)
EPOCHS = 100
BATCH_SIZE = 128
STEPS_PER_EPOCH = 100
VAL_STEPS_PER_EPOCH = 20

train_gen, val_gen, test_data, _ = generate_training_data(
    'data/training_data_round1.h5',
    batch_size=BATCH_SIZE,
    steps_per_epoch=STEPS_PER_EPOCH,
    val_steps_per_epoch=VAL_STEPS_PER_EPOCH
)
model, history = train_model(
    model, train_gen, val_gen,
    epochs=EPOCHS,
    steps_per_epoch=STEPS_PER_EPOCH,
    val_steps_per_epoch=VAL_STEPS_PER_EPOCH
)
X_test, y_test = test_data
y_pred = model.predict(X_test)
fpr, tpr, thresh = roc_curve(y_test, y_pred)
i = len(thresh) - np.searchsorted(thresh[::-1], 0.9)
print(thresh[i], fpr[i], tpr[i])
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(fpr, tpr, color=pal[0], lw=2, label='AUC={:.3f}'.format(roc_auc_score(y_test, y_pred)))
ax.plot([0, 1], [0, 1], color='#555555', ls='--')
plt.legend()
plt.show()
fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(y_pred[y_test == 1], bins=np.linspace(0, 1, 101), color=pal[0], alpha=0.5)
ax.hist(y_pred[y_test == 0], bins=np.linspace(0, 1, 101), color=pal[1], alpha=0.5)
plt.show()
# use the existing model to remove false negatives from BLAST

with h5py.File('data/training_data_round1.h5') as f:
    pos_signals = f['pos_signals'][:]
    neg_signals = f['neg_signals'][:]
    neg_signals_internal = f['neg_internal_signals'][:]
    
preds = model.predict(neg_signals.reshape(-1, 3000, 1)).squeeze()
neg_signals_filt = neg_signals[preds < 0.5]
write_to_hdf5(
    'data/training_data_round2.h5',
    pos_signals, neg_signals_filt, neg_signals_internal
)
with h5py.File('data/training_data_round2.h5') as f:
    print(f['pos_signals'].shape)
    print(f['neg_signals'].shape)
train_gen, val_gen, test_data, _ = generate_training_data(
    'data/training_data_round2.h5',
    batch_size=BATCH_SIZE,
    steps_per_epoch=STEPS_PER_EPOCH,
    val_steps_per_epoch=VAL_STEPS_PER_EPOCH
)
model = build_model(3000)
model, history = train_model(
    model, train_gen, val_gen,
    epochs=EPOCHS,
    steps_per_epoch=STEPS_PER_EPOCH,
    val_steps_per_epoch=VAL_STEPS_PER_EPOCH
)
X_test, y_test = test_data
y_pred = model.predict(X_test)
fpr, tpr, thresh = roc_curve(y_test, y_pred)
i = len(thresh) - np.searchsorted(thresh[::-1], 0.9)
print(thresh[i], fpr[i], tpr[i])
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(fpr, tpr, color=pal[0], lw=2, label='AUC={:.3f}'.format(roc_auc_score(y_test, y_pred)))
ax.plot([0, 1], [0, 1], color='#555555', ls='--')
plt.legend()
plt.show()
fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(y_pred[y_test == 1], bins=np.linspace(0, 1, 101), color=pal[0], alpha=0.5)
ax.hist(y_pred[y_test == 0], bins=np.linspace(0, 1, 101), color=pal[1], alpha=0.5)
plt.show()
model.save('data/model_weights_070319.h5')<jupyter_output><empty_output><jupyter_text>Use K Fold cross validation to do proper model evaluation:<jupyter_code>def get_data(h5_fn):
    with h5py.File(h5_fn) as h5_file:
        pos = h5_file['pos_signals'][:]
        neg = np.concatenate([
            h5_file['neg_signals'][:],
            h5_file['neg_internal_signals'][:]
        ])
    return pos, neg

pos_data, neg_data = get_data('data/training_data_round2.h5')
pos_kf = KFold(n_splits=5, shuffle=False).split(pos_data)
neg_kf = KFold(n_splits=5, shuffle=False).split(neg_data)
pos_pred = []
neg_pred = []

for (pos_train_idx, pos_test_idx), (neg_train_idx, neg_test_idx) in zip(pos_kf, neg_kf):
    pos_train = pos_data[pos_train_idx]
    pos_test = pos_data[pos_test_idx]
    neg_train = neg_data[neg_train_idx]
    neg_test = neg_data[neg_test_idx]
    pos_train, pos_val = train_test_split(pos_train, test_size=0.1)
    neg_train, neg_val = train_test_split(neg_train, test_size=0.1)
    train_gen = SquiggleSequence(
        pos_train, neg_train,
        batch_size=128, steps_per_epoch=100,
    )
    val_gen = SquiggleSequence(
        pos_val, neg_val,
        batch_size=128, steps_per_epoch=20,
    )
    model = build_model(3000)
    model, history = train_model(
        model, train_gen, val_gen,
        epochs=100,
        steps_per_epoch=100,
        val_steps_per_epoch=20
    )
    pos_pred.append(model.predict(pos_test.reshape(-1, 3000, 1)).squeeze())
    neg_pred.append(model.predict(neg_test.reshape(-1, 3000, 1)).squeeze())

pos_pred = np.concatenate(pos_pred)
neg_pred = np.concatenate(neg_pred)
y_pred = np.concatenate([pos_pred, neg_pred])
y_true = np.concatenate([[1] * len(pos_pred), [0] * len(neg_pred)])

fpr, tpr, thresh = roc_curve(y_true, y_pred)
prec, rec, pr_thresh = precision_recall_curve(y_true, y_pred)
i = len(thresh) - np.searchsorted(thresh[::-1], 0.9)
print(thresh[i], fpr[i], tpr[i])
fig, axes = plt.subplots(figsize=(10, 5), ncols=2)
axes[0].plot(fpr, tpr, color=pal[0], lw=2)
axes[0].plot([0, 1], [0, 1], color='#555555', ls='--')
axes[0].set_xlabel('False Positive Rate')
axes[0].set_ylabel('True Positive Rate')
axes[1].plot(rec, prec)
axes[1].plot([0, 1], [prec[1], prec[1]], color='#555555', ls='--')
axes[1].set_xlabel('Recall')
axes[1].set_ylabel('Precision')
axes[1].set_ylim(0, 1.05)
plt.tight_layout()
plt.savefig('figures/roc_pr_curve.svg')
plt.show()<jupyter_output>0.8999793 0.0055848859025262884 0.9163880526615152
{'fraction_non_alphanumeric': 0.13374406366320113, 'fraction_numerical': 0.04274162495186754, 'mean_word_length': 5.329813160032494, 'pattern_counts': {'":': 0, '<': 13, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '25835', 'license': 'MIT', 'n_tokens_mistral': 3388, 'n_tokens_neox': 3161, 'n_words': 518, 'path': 'pipeline/adapter_detection_pipeline/scripts/adapter_detector/notebooks/train_model.ipynb', 'repo_name': 'elifesciences-publications/Simpson_Barton_Nanopore_1', 'size': '7791'}