SingleMoleculeAnalysisScripts/Neuhoff_Generate_save_data.m at main · NilsSomm/SingleMoleculeAnalysisScripts · 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
 clear all
import ebfret.analysis.*

%%change directory to where your trace files are at
% cd = 'Z:\Darcy\DNA Origami\2019 smHinge\dsDNA 3.3.3 9bp\AllConcForAnalysis'
cd 'C:\Users\Michael\Documents\Grad School\Research\Data\Single Molecule Data\ebFRET Workspaces\Type1 Oscillator\7bp'

%% change this to match your acquisition rate
secperFrame = 1;

%%get all filenames that you want to analyze
files = dir('aggregatecatcut.mat')
filecounter = 1;

%%panel is a package that allows dynamic resizing of multiple plots side by
%%side
pan = panel();

%%these arrays are for storing the parameter outputs of the fitting
%%procedure: the amplitude of the decay A and the corresponding rates k1
%%and k2 for both transitions from high states and transitions from low
%%states

Ahigh = [];
Alow = [];
k1high = [];
k1low = [];
k2high = [];
k2low = [];

%%I thought it looked best with figures packed into a 3x5 grid
pan.pack(3, 5);


%%now we want to analyze each file separately
file = files';

%file.namedat = strcat(char(file.name),'.dat')
%formatEB(char(file.name) ,file.namedat);%% ebFRET part


%%these are for saving the data (aka fit traces) and plots (aka parameters and dwell
%%times)
save_name = strcat(file.name, '_save.dat')
mat_save_name = strcat('plotsave', file.name)

%%These parameters are for ebFRET. K_values is how many states to use
%%for each run and can be a list [2 3 4] if desired
K_values = [2];
%%restarts is how many times should we attempt the entire fitting and
%%then take the best attempt
restarts = 3;
%%this is how many cores we have in the parallel computing pool
num_cpu = 6;


% options for loading and stripping traces, see load_fret_defaults and
% eb_defaults for more info


%% write more comments here
opts.load_fret = load_fret_defaults();

opts.load_fret.min_length = 50;

opts.load_fret.max_outliers = 1;

opts.load_fret.strip_first = true;

opts.load_fret.remove_bleaching = true;


% options for vbem algorithm

opts.eb = eb_defaults();
[data, raw] = load_fret(file.name, 'variable', 'data');


%% here we actually do the eb_fret fitting.
runs = eb_fret_waitbar(data.fret, K_values, restarts...
    , 'eb', opts.eb, 'num_cpu', num_cpu);


% w = runs.vb.w;
% u = runs.u;
% M_values =2 ;
% test = eb_fret_pmm_dir(data.fret, w , u, 2);

%%This line ends the parallel pool. If you stop the program while it is
%%fitting traces, you will need to close the parallel pool with this
%%command in console or the script wil
% l crash when it tries to open a
%%new pool on top of the old one.
delete(gcp('nocreate'))


%%runs.vit contains the Viterbi idealized traces.
dummy = struct2table(runs.vit);

%%now we need to put it into the right structure to extract dwell times
for i = 1:length(dummy.z)
    cellpaths{i,1} = i*ones(length(cell2mat(dummy.x(i))),1);
    cellpaths{i,2} = cell2mat(dummy.x(i));
end

paths = cell2mat(cellpaths);

%% here we save the fret traces and idealized fits to the filename we made before
save('cellpaths')
save('cellpaths.dat', 'cellpaths', '-ascii')