拟合分布,拟合优度,p值.用Scipy(Python)可以做到这一点吗？

Question

简介:我是生物信息学家.在我对所有人类基因(约20 000)进行的分析中,我搜索特定的短序列基序,以检查每个基因中出现这个基序的次数.

基因以四个字母(A,T,G,C)的线性序列"书写".例如:CGTAGGGGGTTTAC ......这是遗传密码的四个字母的字母表,就像每个细胞的秘密语言一样,它就是DNA实际存储信息的方式.

我怀疑在一些基因中频繁重复特定的短基序列(AGTGGAC)在细胞的特定生化过程中是至关重要的.由于基序本身非常短,因此用计算工具很难区分基因中的真实功能性实例和偶然看起来相似的实例.为了避免这个问题,我得到了所有基因的序列并连接成一个字符串并进行了改组.存储每个原始基因的长度.然后,对于每个原始序列长度,通过从连接序列中随机重复地挑选A或T或G或C并将其转移到随机序列来构建随机序列.以这种方式,得到的随机序列组具有相同的长度分布,以及总体A,T,G,C组成.然后我在这些随机序列中搜索主题.我将此程序置于1000次并对结果取平均值.

15000个不含给定基序的基因5000个基因含有1个基序3000个基因,含有2个基序1000个含有3个基序的基因... 1个含有6个基序的基因

因此,即使经过1000次真正遗传密码的随机化,也没有任何基因具有超过6个基序.但是在真正的遗传密码中,有一些基因含有超过20个基序的出现,这表明这些重复可能是有效的,并且它不可能通过纯粹的机会找到它们如此丰富.

问题:我想知道找到一个基因的可能性,假设我的分布中出现了20个基序.所以我想知道偶然发现这样一个基因的可能性.我想在Python中实现它,但我不知道如何.

我可以在Python中进行这样的分析吗？

任何帮助,将不胜感激.

Answer 1

在SciPy文档中,您将找到所有已实现的连续分发函数的列表.每个都有一个fit()方法,返回相应的形状参数.

即使您不知道使用哪种发行版,您也可以同时尝试许多发行版,并选择更适合您数据的发行版,如下面的代码所示.请注意,如果您不了解分布,则可能难以拟合样本.

import matplotlib.pyplot as plt
import scipy
import scipy.stats
size = 20000
x = scipy.arange(size)
# creating the dummy sample (using beta distribution)
y = scipy.int_(scipy.round_(scipy.stats.beta.rvs(6,2,size=size)*47))
# creating the histogram
h = plt.hist(y, bins=range(48))

dist_names = ['alpha', 'beta', 'arcsine',
              'weibull_min', 'weibull_max', 'rayleigh']

for dist_name in dist_names:
    dist = getattr(scipy.stats, dist_name)
    param = dist.fit(y)
    pdf_fitted = dist.pdf(x, *param[:-2], loc=param[-2], scale=param[-1]) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc='upper left')
plt.show()

参考文献:

- 与Scipy配电

- 使用Scipy(Python)将经验分布拟合到理论分布？