-
Notifications
You must be signed in to change notification settings - Fork 9
/
plotGeneProfiles.m
62 lines (41 loc) · 1.59 KB
/
plotGeneProfiles.m
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
function plotGeneProfiles(dataset, selectedGenes)
% function to plot profiles of selected genes ordered by pseudotime
[~, processDataMat, ~, ~, ~, ~, ~, ~] = initialization(dataset);
load(processDataMat)
[~, locs] = ismember(selectedGenes, pro.gname);
pro.pseudotime = (pro.pseudotime- min(pro.pseudotime))/(max(pro.pseudotime)-min(pro.pseudotime));
[pseudotimeSorted, pseudotimeSortIndices] = sort(pro.pseudotime);
expr = pro.expr(pseudotimeSortIndices,locs);
nWindows = 100;
t0 = pseudotimeSorted(1);
t1 = pseudotimeSorted(end);
stepSize = (t1-t0)/2/nWindows;
boundariesAndCenters = t0:stepSize:t1;
windowCenters = boundariesAndCenters(2:2:end);
windowRadius = 0.08*(t1-t0);
smoothExpr = calculateWindowMedians(windowCenters, windowRadius, pseudotimeSorted, expr);
%%
figure
colorList = {'b', 'r', 'g', 'k', 'm'};
for index = 1:size(smoothExpr,2)
plot(windowCenters, smoothExpr(:,index),colorList{index},'linewidth',2)
hold on
end
set(gca,'fontsize',18)
legend(selectedGenes, 'location','SouthEast', 'fontsize',12)
xlabel('Pseudotime')
ylabel('Normalized Intensity')
ylim([0 1.05])
xlim([0 max(pro.pseudotime)])
title(dataset)
%%
end
function smoothExpr = calculateWindowMedians(windowCenters, windowWidths, pseudotime, expr)
smoothExpr = zeros(length(windowCenters),size(expr,2));
for index = 1:length(windowCenters)
windowLocations = pseudotime > windowCenters(index) - windowWidths & ...
pseudotime < windowCenters(index) + windowWidths;
smoothExpr(index,:) = median(expr(windowLocations,:));
end
smoothExpr = smoothExpr./repmat(max(smoothExpr),length(windowCenters),1);
end