使用Scipy（Python）拟合分布、优度检验和p值是否可行？
介绍：我是一名生物信息学家。在我的分析中，我研究所有人类基因（大约有20000个），寻找一个特定的短序列模式，以检查这个模式在每个基因中出现的次数。
基因是由四个字母（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中进行这样的分析吗？
任何帮助都将不胜感激。
2 个回答
3
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()
如果你想进行更详细的分析，我推荐使用 Phitter库
import phitter
data: list[int | float] = [...]
phitter_cont = phitter.PHITTER(
data=data,
fit_type="continuous",
num_bins=15,
confidence_level=0.95,
minimum_sse=1e-2,
distributions_to_fit=["beta", "normal", "fatigue_life", "triangular"],
)
phitter_cont.fit(n_workers=6)
37
在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()
参考资料：