-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathgene_pseudotime.m
69 lines (61 loc) · 2.19 KB
/
gene_pseudotime.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
63
64
65
66
67
68
function [pseudotime1] = gene_pseudotime(cell_location,ordered_cell1,cluster_order1,gene_index,gene_name,data,C,mm,folder1,folder)
zzz = get(gca,'colororder');
mycolor = zeros(11,3);
mycolor(1:7,:) = zzz;
mycolor(8,:) = [0 0 0];
mycolor(9,:) = [0 0 0.803922];
mycolor(10,:) = [1 0 1];
mycolor(11,:) = [0.5 1 0];
%based on lineage and cell location to get pseudotime
%computer the Euclidean dist of neighboring cells
dist2 = zeros(1,length(ordered_cell1)-1);
for i = 1:length(dist2)
dist2(i) = norm(cell_location(ordered_cell1(i+1),:)-cell_location(ordered_cell1(i),:));
end
pseudotime1 = [0,cumsum(dist2)];
data_used = data(ordered_cell1,gene_index);
x = 0:1/(length(ordered_cell1)-1):1; % evenly distributed
%x = pseudotime1;
y = data_used';
z = C(ordered_cell1);
idx_b1 = [];
for i = 1:length(cluster_order1)
idx_b1 = [idx_b1;find(z==cluster_order1(i))];
end
y_b1 = y(idx_b1);
p_b1 = polyfit(x(sort(idx_b1)),y_b1,mm);
f_b1 = polyval(p_b1,x(sort(idx_b1)));
% plot gene expression along ptime (label)
for ik = 1:length(cluster_order1)
if length(cluster_order1)~=3 || ik~=3
scatter(x(find(z==cluster_order1(ik))),y(find(z==cluster_order1(ik))),20,mycolor(ik,:),'filled','MarkerEdgeAlpha',0.6,'MarkerFaceAlpha',0.6);
hold on;
else
scatter(x(find(z==cluster_order1(ik))),y(find(z==cluster_order1(ik))),20,mycolor(4,:),'filled','MarkerEdgeAlpha',0.6,'MarkerFaceAlpha',0.6);
hold on;
end
end
hold on;
%plot(x(sort(idx_b1)),f_b1,'-','LineWidth',3,'Color','k','MarkerSize',10);
%plot(x(sort(idx_b1)),f_b1,':','LineWidth',3,'Color','k','MarkerSize',10);
plot(x(sort(idx_b1)),f_b1,'LineWidth',3,'Color','k','MarkerSize',10);
%xticks([])
L = get(gca,'YLim');
set(gca,'YTick',round(linspace(L(1),L(2),2),1))
%ylim([floor(min(y)),ceil(max(y))])
yticks([ceil(min(y)),floor(max(y))])
%yticks([round(min(y),1),round(max(y),1)])
xticks([]);
hold off
fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
% fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 0 4 3];
%
fig.Units = 'Inches';
fig.Position = [0 0 3 2.3];
fig.OuterPosition=fig.InnerPosition;
print([folder1 '\' folder '\' gene_name],'-dpdf','-r600','-fillpage');
end