nanocompore_peaks.py

import pandas as pd
import numpy as np
from scipy.signal import find_peaks
from collections import OrderedDict

test = "GMM_logit_pvalue"
df = pd.read_csv("out_nanocompore_results.tsv", sep="\t")
df["Peak"] = 0 
df = df[["pos", "chr", "genomicPos", "ref_id", "strand", "ref_kmer", "Peak", test]]
transcripts = set(df["ref_id"])
p_val_lim = 0.01
sig_lim = -np.log10(p_val_lim)

i=1
for tx in transcripts:
    if(not i%50): print(i)
    i+=1
    res = df[df.ref_id==tx]
    x = -np.log10(res[test])
    x = x.fillna(0)
    threshold = sig_lim
    peaks, extra = find_peaks(x, height=threshold, distance=9)
    peaks_indexes = res.iloc[peaks].index
    df.loc[peaks_indexes, "Peak"] = extra["peak_heights"]
df.to_csv("out_nanocompore_results_peaks.txt", index=False, sep="\t")