A kernel density estimate (KDE) models a probability density by placing one localized kernel at every observed data point and averaging their contributions. With a Gaussian kernel, the estimate averages Gaussians centered on the observations; a general kernel κ(x) provides the broader formulation. section 9.2.1
Use KDE when choosing a parametric density model or mixture count is undesirable:
- Avoid model fitting apart from fine-tuning the bandwidth parameter h. section 9.2.1
- Avoid selecting the number of mixture components K, because every observation supplies a cluster center. section 9.2.1 section 9.8
- Trade memory and evaluation speed for accuracy, making KDE most suitable for relatively small datasets. section 9.2.1
- Adjust smoothing through the bandwidth; the resulting density can be compared with a histogram, whose smoothing depends on bin size. section 9.2.1 section 9.6
RNA-seq promoter-flux estimation
Estimate the flux of a T7 promoter from RNA-seq data, then compare a histogram-based density estimate with a Gaussian KDE. The Gaussian estimate appears smoother, while its smoothing can be further tuned through the bandwidth. section 9.2.1
Listing 9.2 Kernel density estimate
import numpy as np
import matplotlib.pyplot as plt
class KDE():
def __init__(self):
#Histogram and Gaussian Kernel Estimator used to
#analyze RNA-seq data for flux estimation of a T7 promoter
self.G = 1e9 #1
self.C = 1e3 #2
self.L = 100 #3
self.N = 1e6 #4
self.M = 1e4 #5
self.LN = 1000 #6
self.FDR = 0.05 #7
#uniform sampling (poisson model)
self.lmbda = (self.N * self.L) / self.G #8
self.C_est = self.M/(1-np.exp(-self.lmbda)) #9
self.C_cvrg = self.G - self.G *
np.exp(-self.lmbda) #10
self.N_gaps = self.N * np.exp(-self.lmbda) #11
#gamma prior sampling (negative binomial model)
#X = "number of failures before rth success"
self.k = 0.5 #12
self.p = self.lmbda/(self.lmbda + 1/self.k) #13
self.r = 1/self.k #14
#RNAP binding data (RNA-seq)
self.data = np.random.negative_binomial(self.r, self.p, size=self.LN)
def histogram(self):
self.bin_delta = 1 #15
self.bin_range = np.arange(1, np.max(self.data), self.bin_delta)
self.bin_counts, _ = np.histogram(self.data, bins=self.bin_range)
#histogram density estimation
#P = integral_R p(x) dx, where X is in R^3
#p(x) = K/(NxV), where K=number of points in region R
#N=total number of points, V=volume of region R
rnap_density_est = self.bin_counts/(sum(self.bin_counts) *
self.bin_delta)
return rnap_density_est
def kernel(self):
#Gaussian kernel density estimator with smoothing parameter h
#sum N Guassians centered at each data point, parameterized by
common std dev h
x_dim = 1 #16
h = 10 #17
rnap_density_support = np.arange(np.max(self.data))
rnap_density_est = 0
for i in range(np.sum(self.bin_counts)):
rnap_density_est += (1/(2*np.pi*h**2)**(x_dim/2.0))*np.exp(-
(rnap_density_support - self.data[i])**2 / (2.0*h**2))
#end for
rnap_density_est = rnap_density_est / np.sum(rnap_density_est)
return rnap_density_est
if __name__ == "__main__":
kde = KDE()
est1 = kde.histogram()
est2 = kde.kernel()
plt.figure()
plt.plot(est1, c='b', label='histogram')
plt.plot(est2, c='r', label='gaussian kernel')
plt.title("RNA-seq density estimate based on negative binomial sampling
model")
plt.xlabel("read length, [base pairs]"); plt.ylabel("density");
plt.legend()
plt.show()
A Python implementation defines a KDE class and sets computational and probabilistic parameters for the RNA-seq analysis, including a Poisson sampling model and a negative binomial model. section 9.2.1
