forked from neuhoffmj/SingleMoleculeAnalysisScripts
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathPlotTraceMarginal.m
More file actions
107 lines (94 loc) · 3.62 KB
/
PlotTraceMarginal.m
File metadata and controls
107 lines (94 loc) · 3.62 KB
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
%% Import Data from ebFRET GUI export as time series (.mat)
clear;
close all;
clc;
load('23_2_9_6bp_eb_export.mat')
bp = 6;
frameRate = 10;
exposureTime = 1/frameRate;
for i = 1:traces(end,1)
donor{i} = traces(traces(:,1)==i,2);
acceptor{i} = traces(traces(:,1)==i,3);
FRET{i} = traces(traces(:,1)==i, 4);
viterbiState{i} = traces(traces(:,1)==i, 5);
viterbiMean{i} = traces(traces(:,1)==i, 6);
end
[thigh, tlow,thigh_Ind, tlow_Ind, histTlow, histThigh, cumtlow, cumthigh] = eb_dwelltimes_traceInd(viterbiState', exposureTime);
%% Compile FRET States
highStates = [];
lowStates = [];
highRawStates = [];
lowRawStates = [];
for i = 1:size(viterbiMean, 2)
highStates = cat(1, highStates, viterbiMean{i}(viterbiMean{i}==max(viterbiMean{i})));
lowStates = cat(1, lowStates, viterbiMean{i}(viterbiMean{i}==min(viterbiMean{i})));
highRawStates = cat(1, highRawStates, FRET{i}(viterbiMean{i}==max(viterbiMean{i})));
lowRawStates = cat(1, lowRawStates, FRET{i}(viterbiMean{i}==min(viterbiMean{i})));
end
%% Make FRET Histograms
%
% histos = figure();
% histogram(highStates, 20, "Normalization","probability", "FaceColor", 'red', 'DisplayName', 'High state')
% hold on
% % figure;
% histogram(lowStates, 20, "Normalization","probability", "FaceColor", 'blue', 'DisplayName', 'Low State')
% title(sprintf('%d bp HMM FRET Value Histograms', bp))
% ylabel('Frequency')
% xlabel('FRET (arbitrary units)')
% legend('Location', 'best')
% saveas(histos, sprintf('%d bp HMM FRET Histograms', bp), 'png')
%
%% Make FRET Histograms Normalization fix
%% Old Version
% [countHigh, xHigh] = hist(highRawStates, 20, "Normalization","count", "FaceColor", 'red', 'DisplayName', 'High state');
% [countLow, xLow] = hist(lowRawStates, 20, "Normalization","probability", "FaceColor", 'blue', 'DisplayName', 'Low State');
% histos = figure();
% totalCounts = sum(countHigh)+sum(countLow);
% bar(xHigh, countHigh/totalCounts, 'b', 'DisplayName', 'High State');
% hold on;
% bar(xLow, countLow/totalCounts, 'r', 'DisplayName', 'Low State');
% title(sprintf('%d bp HMM FRET Value Histograms', bp))
% ylabel('Frequency')
% xlabel('FRET (arbitrary units)')
% legend('Location', 'best')
% saveas(histos, sprintf('%d bp HMM FRET Histograms', bp), 'png')
%% New Version
histos = figure();
histogram(highRawStates, 'BinWidth', .01, 'DisplayName', 'HMM High FRET States')
hold on
histogram(lowRawStates, 'BinWidth', .01, 'DisplayName', 'HMM Low FRET States')
framesHighFRET = size(highRawStates,1);
framesLowFRET = size(lowRawStates,1);
timeAvgFRET = framesHighFRET/(framesHighFRET + framesLowFRET);
title(sprintf('%d bp HMM FRET Value Histograms', bp))
ylabel('Frames Spent at FRET Efficiency')
xlabel('FRET (arbitrary units)')
lgd = legend('Location', 'best');
title(lgd, sprintf('Normalized Average FRET: %.2f', timeAvgFRET))
saveas(histos, sprintf('%d bp Raw FRET Histograms', bp), 'png')
%% Plot Traces
tracesfig = figure('Position', [100 100 1200 200]);
plotNum = 1;
for i = 24:24
% figure('Position', [500 500 1200 400])
ax1 = subplot(1,8,[1:7]);
ax1.PositionConstraint = 'outerposition';
length = size(viterbiMean{i},1);
t = 0:exposureTime:(exposureTime*length-exposureTime);
plot(t, viterbiMean{i}, "LineWidth", 1.5, 'Color', 'r')
hold on;
plot(t, FRET{i}, "LineWidth", 0.5, 'Color', 'k')
plotNum = plotNum + 1;
xlim([0 max(t)])
ylim([0,1])
hold on;
end
ylabel('FRET (arbitrary units)')
xlabel('Time (s)')
title(sprintf('%d bp Sample Traces', bp))
ax2 = subplot(1,8,[8:8]);
histogram(FRET{i}, 20);
xlim([0,1])
set(ax2,'xticklabel',[], 'yticklabel',[])
view([90 -90]);
saveas(tracesfig, sprintf('%d bp Sample Traces Figure', bp), 'png')