smFSmodels/kinetics.py at master · daviddesancho/smFSmodels · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
import numpy as np

def calc_life(trajs, ub=5, lb=-5):
    """
    Identifies transition paths and returns lifetimes of states.

    Parameters
    ----------
    trajs : list of lists
        Set of trajectories.

    ub, lb : float
        Cutoff value for upper and lower states.
    """
    try:
        assert ub > lb
    except AssertionError:
        print (" Upper bound is lower than lower bound")
        return

    lifeA = []
    lifeB = []
    time = 0
    for tr in trajs:
        state = None
        ntp = 0
        time_prev = 0
        for t,q in enumerate(tr):
            # assign state when beyond boundaries
            if q > ub: # state "B"
                if state == 'A':
                    ntp +=1
                    lifeA.append(time - time_prev)
                    time_prev = time
                state = 'B'
            elif q < lb: # state "A"
                if state == 'B':
                    ntp +=1
                    lifeB.append(time - time_prev)
                    time_prev = time
                state = 'A'
            else:
                if state == 'A' and q < ub:
                    time = t
                elif state == 'B' and q > lb:
                    time = t
    return lifeA, lifeB


def calc_life_multi(trajs, bounds=[[-3,-1], [1,3], [6,8]]):
    """
    Identifies transition paths and returns lifetimes of states.

    Parameters
    ----------
    trajs : list of lists
        Set of trajectories.

    bounds : list
        Limits for states.
    """

    life = [[],[],[]]
    tau = {}
    time = 0
    for tr in trajs:
        state = None
        ntp = 0
        for t,q in enumerate(tr):
            # assign state when beyond boundaries
            for i,b in enumerate(bounds):
                if b[0] < q < b[1]:
                    if state != i:
                        try:
                            life[state].append(time - time_prev)
                            tau[i, state].append(time - time_prev)
                        except TypeError:
                            pass
                        except KeyError:
                            tau[(i, state)] = [time - time_prev]
                        state = i
                        time_prev = t
                        break
                    state = i
                time = t

    return life, tau

def lifetimes(data, f_bound=-5, u_bound=5):
    """
    Estimates lifetimes using a transition path analysis. Transitions are
    only assigned from one state to the other when the core of the other
    state is reached.

    Parameters
    ----------
    data : np.array
        Time series data containing times and corrected extensions.

    Returns
    -------
    tau_f : list
        Waiting times in the unfolded state.

    tau_u : list
        Waiting times in the folded state.

    data_f : list
        Stretches of data corresponding to the folded segments.

    data_u : list
        Stretches of data corresponding to the unfolded segments.

    tp_f : list
        Transition paths for folding.

    tp_u : list
        Transition paths for unfolding.
    """
    folded = False
    unfolded = False
    maybetp = []
    recrossings = []
    t = data[:,0]
    dist = data[:,1]
    data_f = []
    data_u = []
    tau_u = []
    tau_f = []
    tp_f = []
    tp_u = []
    time = 0
    data_cum = []
    for t, d in zip(t,dist):
        if d <= f_bound:
            folded = True
            if unfolded:
                #print ' Refolding event: %g'%t,
                tp_u.append(np.array(maybetp))
                tau_f.append(t-time)
                data_u.append(np.array(data_cum))
                unfolded = False
                time = t
                data_cum = []
            else:
                recrossings.append(np.array(maybetp))
                for tt,dd in maybetp:
                    data_cum.append([tt,dd])
            data_cum.append([t,d])

            maybetp = []

        elif u_bound <= d :
            unfolded = True
            if folded:
                tp_f.append(np.array(maybetp))
                #print ' Unfolding event: %g'%t,
                tau_u.append(t-time)
                data_f.append(np.array(data_cum))
                folded = False
                time = t
                data_cum = []
            else:
                recrossings.append(np.array(maybetp))
                for tt,dd in maybetp:
                    data_cum.append([tt,dd])
            data_cum.append([t,d])
            maybetp = []
        else:
            maybetp.append([t,d])

    if unfolded:
        #print ' Refolding event: %g'%t,
        tau_f.append(t-time)
        data_u.append(np.array(data_cum))

    if folded:
#        tp_f.append(np.array(maybetp))
        #print ' Unfolding event: %g'%t,
        tau_u.append(t-time)
        data_f.append(np.array(data_cum))
    return tau_f, tau_u, data_f, data_u, tp_f, tp_u, recrossings