How to add Gaussian fitting to a genome coverage plot?
0
0
Entering edit mode
16 months ago

Hi,

I am making a genome coverage plot using python and trying to add a gaussian fitting. What I want is like this: expected gaussian fitting But when I use my code to plot, it gives me something like this: actual result

Please only use the red and black curve and ignore the x,y label and the green&blue curves in the first image. Here's my code:

cov_roi = cov[4569000:4666000] #select the genome region of interests (roi) using pandas dataframe

x = cov_roi.loci
y = cov_roi.coverage

n = len(x)
mean = sum(y)/n
sigma = (sum((y-mean)**2)/n) **(0.5)
mu = np.median(x)

def gaussian(x, mu, sig):
return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

plt.figure(figsize=(6,4))

plt.plot(x,y, "r." ,label = 'MAPLE', ms=0.3,color = 'r')

plt.plot(x, gaussian(x, mu, sigma), 'k-', alpha=0.5,lw=2, label='Gaussian fitting')
plt.ylim((-500, 10000))

plt.xlabel('Base Position (Kbp)')
plt.ylabel('Coverage')

plt.title('Coverage Map')
plt.show()


Could anybody help me?

Thank you so much!

genome python gaussian fitting • 210 views