Changeset 66 for liacs/dbdm/dbdm_4


Ignore:
Timestamp:
Jan 3, 2010, 8:20:21 PM (15 years ago)
Author:
Rick van der Zwet
Message:

Annotations and readablity fixes

Location:
liacs/dbdm/dbdm_4
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • liacs/dbdm/dbdm_4/ecoli_hmm.py

    r64 r66  
    2929    }
    3030
     31# Transform the state to something human readable
    3132dna_ascii_translate = {
    32     '0' : '*',
     33    '0' : '.',
    3334    '1' : '<',
    3435    '2' : '<',
    3536    '3' : '<',
    36     '4' : '-',
     37    '4' : '|',
    3738    '5' : '>',
    3839    '6' : '>',
    … …  
    4243dna_ascii = MultiReplace(dna_ascii_translate)
    4344
    44 def pretty_print(test_seq, ans_seq, v, length=70, parts=10, seperator=''):
     45def pretty_print(test_seq, ans_seq, v, length=100, parts=10, seperator=''):
    4546    """ Pretty printing of output for verification purposes """
    4647    for i in range(0,len(v[0]),length):
    … …  
    5455            result.append(dna_ascii.replace(''.join(map(str,v[0][t:t+parts]))))
    5556
    56         print seperator.join(seq)
    57         print seperator.join(ans)
    58         print seperator.join(result)
     57        print "ORG: " + seperator.join(seq)
     58        print "CSV: " + seperator.join(ans)
     59        print "HMM: " + seperator.join(result)
    5960        print ''
    6061    print "fairness of test_seq: ", v[1]
    … …  
    7980    # 7) Stop-codon  : end of gene - part 3
    8081    A = [
    81             [0.8, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    82             [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    83             [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
    84             [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
    85             [0.0, 0.0, 0.0, 0.0, 0.7, 0.3, 0.0, 0.0],
    86             [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
    87             [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
    88             [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
     82            [0.8, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],   # 0
     83            [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],   # 1
     84            [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],   # 2
     85            [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],   # 3
     86            [0.0, 0.0, 0.0, 0.0, 0.9, 0.1, 0.0, 0.0],   # 4
     87            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],   # 5
     88            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],   # 6
     89            [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],   # 7
    8990        ]
    9091
    … …  
    9495    B = [
    9596            # e.g. state 0 -> emission probability
    96             [0.2, 0.2, 0.2, 0.2, 0.2] ,
    97             [0.9, 0.0, 0.1, 0.0, 0.0] ,
    98             [0.0, 0.0, 0.0, 1.0, 0.0] ,
    99             [0.0, 0.0, 1.0, 0.0, 0.0] ,
    100             [0.2, 0.2, 0.2, 0.2, 0.2] ,
    101             [0.0, 0.0, 0.0, 1.0, 0.0] ,
    102             [0.7, 0.0, 0.3, 0.0, 0.0] ,
    103             [0.7, 0.0, 0.3, 0.0, 0.0] ,
     97            [0.2, 0.2, 0.2, 0.2, 0.2] ,     # 0
     98            [0.9, 0.0, 0.1, 0.0, 0.0] ,     # 1
     99            [0.0, 0.0, 0.0, 1.0, 0.0] ,     # 2
     100            [0.0, 0.0, 1.0, 0.0, 0.0] ,     # 3
     101            [0.2, 0.2, 0.2, 0.2, 0.2] ,     # 4
     102            [0.0, 0.0, 0.0, 1.0, 0.0] ,     # 5
     103            [0.7, 0.0, 0.3, 0.0, 0.0] ,     # 6
     104            [0.7, 0.0, 0.3, 0.0, 0.0] ,     # 7
    104105        ]
    105106   
    … …  
    109110    m = ghmm.HMMFromMatrices(sigma,ghmm.DiscreteDistribution(sigma),A ,B, pi)
    110111    print "Initial HMM"
    111     print m
     112    print m.verboseStr()
    112113
    113114    obs_seq = m.sampleSingle(20)
    … …  
    130131    handle.close()
    131132
    132     test_seq = contig_seq['AE005174v2-2'][0:490]
    133     ans_seq = answer['AE005174v2-2'][0:490]
     133    test_size = 1000;
     134    test_seq = contig_seq['AE005174v2-2'][0:test_size]
     135    ans_seq = answer['AE005174v2-2'][0:test_size]
    134136    test_eseq=ghmm.EmissionSequence(sigma,list(test_seq))
    135137    v = m.viterbi(test_eseq)
    … …  
    139141    # Train sequence
    140142    print "Training baumWelch"
    141     train_seq = ghmm.EmissionSequence(sigma,list(contig_seq['AE005174v2-1']))
     143    train_seq = ghmm.EmissionSequence(sigma,list(contig_seq['AE005174v2-1'][0:10000]))
    142144    v = m.baumWelch(train_seq)
    143     print m
     145    print m.verboseStr()
    144146   
    145147

  • liacs/dbdm/dbdm_4/gene_detect.py

    r64 r66  
    5858                        handle = open(contig + ".raw","rU")
    5959                        contig_seq[contig] = handle.read()
    60                         abstract_contig[contig] = list('*' * len(contig_seq[contig]))
     60                        abstract_contig[contig] = list('.' * len(contig_seq[contig]))
    6161                        handle.close()
    6262                    except IOError:
    … …  
    9494                    if start_pos < stop_pos:
    9595                        abstract_contig[contig][start_pos:start_pos+3] = ['<'] * 3
    96                         abstract_contig[contig][start_pos+3:stop_pos] = ['-'] * (stop_pos - (start_pos + 3))
     96                        abstract_contig[contig][start_pos+3:stop_pos] = ['|'] * (stop_pos - (start_pos + 3))
    9797                        abstract_contig[contig][stop_pos:stop_pos+3] = ['>'] * 3
    9898                    break
Note: See TracChangeset for help on using the changeset viewer.