2011-07-07 44 views
17

簡介:我是一名生物信息學家。在我對所有人類基因(約20 000個)進行的分析中,我搜索了一個特定的短序列模體,以檢查這個模體在每個基因中發生了多少次。擬合分佈,擬合度,p值。 Scipy(Python)可以做到這一點嗎?

基因以四個字母(A,T,G,C)以線性順序「書寫」。例如:CGTAGGGGGTTTAC ...這是遺傳密碼的四字母字母表,就像每個細胞的祕密語言,它是DNA如何存儲信息。

我懷疑某些基因中特定短基序(AGTGGAC)的頻繁重複在細胞的特定生化過程中至關重要。由於圖案本身非常短,因此用計算工具很難區分基因中真正的功能性例子和偶然看起來相似的例子。爲了避免這個問題,我得到了所有基因的序列並連接成一個單一的字符串,然後洗牌。存儲每個原始基因的長度。然後,對於每個原始序列長度,通過從級聯序列中隨機地重複選擇A或T或G或C並將其轉移到隨機序列來構建隨機序列。通過這種方式,得到的一組隨機序列具有相同的長度分佈,以及總的A,T,G,C組成。然後我在這些隨機序列中搜索基序。我對這個程序進行了1000次修改並取得了平均結果。

不包含給定的主題是包含含有2個圖案 1000包含3個圖案 基因1個基 3000基因 5000個基因...包含6個圖案

所以 1基因

15000個基因即使經過1000次真遺傳密碼隨機化,也沒有任何基因具有6個以上的基序。但是在真正的遺傳密碼中,有幾個基因含有20個以上的主題,這表明這些重複可能是功能性的,並且不可能通過純粹的機會在如此豐富的情況下找到它們。

問題: 我想知道找到一個基因的概率,讓我們說我的分佈中的20個主題的發生。所以我想知道偶然發現這種基因的可能性。我想用Python實現它,但我不知道如何。

我可以在Python中做這樣的分析嗎?

任何幫助,將不勝感激。

+2

請注意,您已大幅修改了你的問題。是否有可能將這個問題還原爲您的原始問題以及所有新細節的明確「更新」部分?或者也許只是一個新問題?謝謝 – eat

+0

你可能會考慮詢問[BioStar](http://biostar.stackexchange.com/questions) – ars

+1

我問一個新的問題:http://stackoverflow.com/questions/6620471/fitting-empirical-distribution-to -theoretical-ones-with-scipy-python –

回答

28

In SciPy documentation您可以找到所有實施的連續分配功能的列表。每一個都有a fit() method,它返回相應的形狀參數。

即使您不知道要使用哪種分配,也可以同時嘗試許多分配,並選擇更適合您數據的分配,如下面的代碼所示。請注意,如果您對分配不瞭解,則可能難以擬合樣本。

enter image description here

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() 

參考文獻:

- Distribution fitting with Scipy

- Fitting empirical distribution to theoretical ones with Scipy (Python)?

+0

上面的代碼會引發以下消息錯誤:「Dist = getattr(scipy.stats,dist_name)」語句中的「AttributeError:'模塊'對象沒有屬性'arcsineweibull_min'」 。我的版本是:scipy是0.13.3,numpy是1.8.0,matplotlib是1.3.1。 – srodriguex

+1

@srodriguex謝謝!有一個小錯字,我只是修復它 –

+0

@SaulloCastro我怎樣才能將這個'fit()'函數用於3D表面擬合,我剛剛使用'scipy.linalg.lstsq'完成了?我如何確認我對數據非常合適。謝謝。 – diffracteD