拟合分布、拟合优度、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中进行这种分析吗？

任何帮助都将不胜感激。

- s_sherly

---

请注意，您已经大幅修改了您的问题。是否可能将此问题恢复为原始问题，并添加一个清晰的“更新”部分以涵盖所有新细节？或者只是提出一个新问题？谢谢。 - eat

你可以考虑在BioStar上提问。 - ars

我提出一个新问题：https://dev59.com/qGw15IYBdhLWcg3wbLLU - s_sherly

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)
    params = dist.fit(y)
    arg = params[:-2]
    loc = params[-2]
    scale = params[-1]
    if arg:
        pdf_fitted = dist.pdf(x, *arg, loc=loc, scale=scale) * size
    else:
        pdf_fitted = dist.pdf(x, loc=loc, scale=loc) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc='upper left')
plt.show()

参考资料：

- 使用Scipy进行分布拟合

- 如何使用Scipy（Python）将经验分布拟合到理论分布？

- Saullo G. P. Castro

---

上面的代码在“dist = getattr(scipy.stats, dist_name)”语句处抛出以下错误消息：“AttributeError: 'module' object has no attribute 'arcsineweibull_min'”。我的版本是：scipy为0.13.3，numpy为1.8.0，matplotlib为1.3.1。 - srodriguex

@srodriguex 谢谢您！有一个小错误，我刚刚修正了。 - Saullo G. P. Castro

@SaulloCastro，我如何将fit()函数应用于三维曲面拟合，我刚刚使用scipy.linalg.lstsq实现了这一点？我该如何确认我对数据进行了良好的拟合。谢谢。 - diffracteD

这个例子现在不能正常工作了，你能更新一下吗？ - Khris

@SaulloG.P.Castro 好的，我一直在尝试这个。代码仍然可以工作，但是 arcsine 拟合出现了问题，可能是函数调用发生了变化。此外，其他函数看起来与您的图像有所不同，但它们的拟合效果现在要好得多。 - Khris

Answer 2

import numpy as np
import pandas as pd
import scipy.stats as st
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
import random

mpl.style.use("ggplot")

def danoes_formula(data):
    """
    DANOE'S FORMULA
    https://en.wikipedia.org/wiki/Histogram#Doane's_formula
    """
    N = len(data)
    skewness = st.skew(data)
    sigma_g1 = math.sqrt((6*(N-2))/((N+1)*(N+3)))
    num_bins = 1 + math.log(N,2) + math.log(1+abs(skewness)/sigma_g1,2)
    num_bins = round(num_bins)
    return num_bins

def plot_histogram(data, results, n):
    ## n first distribution of the ranking
    N_DISTRIBUTIONS = {k: results[k] for k in list(results)[:n]}

    ## Histogram of data
    plt.figure(figsize=(10, 5))
    plt.hist(data, density=True, ec='white', color=(63/235, 149/235, 170/235))
    plt.title('HISTOGRAM')
    plt.xlabel('Values')
    plt.ylabel('Frequencies')

    ## Plot n distributions
    for distribution, result in N_DISTRIBUTIONS.items():
        # print(i, distribution)
        sse = result[0]
        arg = result[1]
        loc = result[2]
        scale = result[3]
        x_plot = np.linspace(min(data), max(data), 1000)
        y_plot = distribution.pdf(x_plot, loc=loc, scale=scale, *arg)
        plt.plot(x_plot, y_plot, label=str(distribution)[32:-34] + ": " + str(sse)[0:6], color=(random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)))
    
    plt.legend(title='DISTRIBUTIONS', bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.show()

def fit_data(data):
    ## st.frechet_r,st.frechet_l: are disbled in current SciPy version
    ## st.levy_stable: a lot of time of estimation parameters
    ALL_DISTRIBUTIONS = [        
        st.alpha,st.anglit,st.arcsine,st.beta,st.betaprime,st.bradford,st.burr,st.cauchy,st.chi,st.chi2,st.cosine,
        st.dgamma,st.dweibull,st.erlang,st.expon,st.exponnorm,st.exponweib,st.exponpow,st.f,st.fatiguelife,st.fisk,
        st.foldcauchy,st.foldnorm, st.genlogistic,st.genpareto,st.gennorm,st.genexpon,
        st.genextreme,st.gausshyper,st.gamma,st.gengamma,st.genhalflogistic,st.gilbrat,st.gompertz,st.gumbel_r,
        st.gumbel_l,st.halfcauchy,st.halflogistic,st.halfnorm,st.halfgennorm,st.hypsecant,st.invgamma,st.invgauss,
        st.invweibull,st.johnsonsb,st.johnsonsu,st.ksone,st.kstwobign,st.laplace,st.levy,st.levy_l,
        st.logistic,st.loggamma,st.loglaplace,st.lognorm,st.lomax,st.maxwell,st.mielke,st.nakagami,st.ncx2,st.ncf,
        st.nct,st.norm,st.pareto,st.pearson3,st.powerlaw,st.powerlognorm,st.powernorm,st.rdist,st.reciprocal,
        st.rayleigh,st.rice,st.recipinvgauss,st.semicircular,st.t,st.triang,st.truncexpon,st.truncnorm,st.tukeylambda,
        st.uniform,st.vonmises,st.vonmises_line,st.wald,st.weibull_min,st.weibull_max,st.wrapcauchy
    ]
    
    MY_DISTRIBUTIONS = [st.beta, st.expon, st.norm, st.uniform, st.johnsonsb, st.gennorm, st.gausshyper]

    ## Calculae Histogram
    num_bins = danoes_formula(data)
    frequencies, bin_edges = np.histogram(data, num_bins, density=True)
    central_values = [(bin_edges[i] + bin_edges[i+1])/2 for i in range(len(bin_edges)-1)]

    results = {}
    for distribution in MY_DISTRIBUTIONS:
        ## Get parameters of distribution
        params = distribution.fit(data)
        
        ## Separate parts of parameters
        arg = params[:-2]
        loc = params[-2]
        scale = params[-1]
    
        ## Calculate fitted PDF and error with fit in distribution
        pdf_values = [distribution.pdf(c, loc=loc, scale=scale, *arg) for c in central_values]
        
        ## Calculate SSE (sum of squared estimate of errors)
        sse = np.sum(np.power(frequencies - pdf_values, 2.0))
        
        ## Build results and sort by sse
        results[distribution] = [sse, arg, loc, scale]
        
    results = {k: results[k] for k in sorted(results, key=results.get)}
    return results
        
def main():
    ## Import data
    data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())
    results = fit_data(data)
    plot_histogram(data, results, 5)

if __name__ == "__main__":
    main()
    
    

- Sebastian Jose