2014-09-23 1 views
1

배경 : MEME file format에서MEME 파일 형식 시퀀스 PWM의 엔트로피?

내가 가지고있는 DNA의 주제는 position-weight-matrices (PWMs)로 표현. 다음은 Here에서 가져온 예제입니다 (자세한 내용은 Here 참조).

MOTIF c585_n005_AGCWTAATC_d_0.034 

letter-probability matrix: alength= 4 w= 10 nsites= 1000 E= 0.05 
0.4000 0.2000 0.2000 0.2000 
0.0000 0.0000 1.0000 0.0000 
1.0000 0.0000 0.0000 0.0000 
0.0000 0.0000 0.0000 1.0000 
0.0000 0.0000 0.1130 0.8870 
1.0000 0.0000 0.0000 0.0000 
0.9217 0.0000 0.0000 0.0783 
0.0000 0.0000 1.0000 0.0000 
0.0000 1.0000 0.0000 0.0000 
0.0500 0.0500 0.0500 0.8500 

MOTIF c101_n004_ATCARTACA_d_0.042 

letter-probability matrix: alength= 4 w= 9 nsites= 1000 E= 0.05 

1.0000 0.0000 0.0000 0.0000 
0.0000 0.0000 0.0402 0.9598 
0.0000 1.0000 0.0000 0.0000 
1.0000 0.0000 0.0000 0.0000 
0.2696 0.0000 0.7304 0.0000 
0.0000 0.0000 0.0000 1.0000 
1.0000 0.0000 0.0000 0.0000 
0.0000 1.0000 0.0000 0.0000 
1.0000 0.0000 0.0000 0.0000 

질문 : 각 주제에 대한 Shannon entropy을 계산할 수있는 방법을 간단한 스크립트로

?

저는 (Python, Perl, MATLAB)에 익숙하지만 이해할 수없는 언어로 스크립트를 검색했습니다. 프로그래밍 방식으로 간단한 방법으로 또는 Biopython/Bioperl/MATLAB-toolbox에서 사용할 수있는 도구를 사용하여이를 수행하는 방법에 대한 조언을 보내 주시면 감사하겠습니다.

답변

2

이 시도 :

def entropy(X, Y): 
    import numpy as np 
    probs = [] 
    for c1 in set(X): 
     for c2 in set(Y): 
      probs.append(np.mean(np.logical_and(X == c1, Y == c2))) 
    return np.sum(-p * np.log2(p) for p in probs) 

def calcShannonEnt(dataSet): 
    import math 
    numEntries = len(dataSet) 
    labelCounts = {} 
    for featVec in dataSet: #the the number of unique elements and their occurance 
     currentLabel = featVec[-1] 
     if currentLabel not in labelCounts.keys(): labelCounts[currentLabel] = 0 
     labelCounts[currentLabel] += 1 
    shannonEnt = 0.0 
    for key in labelCounts: 
     prob = float(labelCounts[key])/numEntries 
     shannonEnt -= prob * math.log(prob,3) 
    return shannonEnt 

dataSet = [[0, 1, 1, 'a'], 
    [0, 1, 1, 'b'], 
    [0, 1, 1, 'c'], 
    [0, 1, 1, 'c'], 
    [1, 0, 1, 'c']] 

a = calcShannonEnt(dataSet) 

from cogent import LoadSeqs 
from cogent.parse.meme import MemeParser 

print a