-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathst_read10xdir.m
113 lines (94 loc) · 3.4 KB
/
st_read10xdir.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
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
function [ste] = st_read10xdir(selpath)
%Read 10x folder
%Read files from a 10x Visium cellranger output folder
%readVisium
if nargin < 1, selpath = uigetdir; end
fprintf('Processing %s...\n', selpath);
[~, aff] = i_guessmtxfile(selpath);
if ~isempty(aff)
h5fname = fullfile(selpath, sprintf('%sfiltered_feature_bc_matrix.h5', aff));
zh5fname = fullfile(selpath, sprintf('%sfiltered_feature_bc_matrix.h5.gz', aff));
else
h5fname = fullfile(selpath, 'filtered_feature_bc_matrix.h5');
zh5fname = fullfile(selpath, 'filtered_feature_bc_matrix.h5.gz');
end
% h5fname=fullfile(selpath,'filtered_feature_bc_matrix.h5');
% h5fname=fullfile(selpath,'raw_feature_bc_matrix.h5');
if ~exist(h5fname, 'file')
if ~exist(zh5fname, 'file')
error('[st_read10xdir] No matrix.h5 file.');
else
[~, nametxt] = fileparts(zh5fname);
fprintf('Unzipping %s.gz...\n', nametxt);
gunzip(zh5fname);
end
end
imgfolder = fullfile(selpath, 'spatial');
if ~isempty(aff)
image_file = fullfile(imgfolder, sprintf('%stissue_hires_image.png', aff));
else
image_file = fullfile(imgfolder, 'tissue_hires_image.png');
end
img = imread(image_file);
json_paths = fullfile(imgfolder, sprintf('%sscalefactors_json.json', aff));
txt = fileread(json_paths);
scalef = jsondecode(txt);
% scalef.tissue_hires_scalef=value.tissue_hires_scalef;
% scalef.spot_diameter_fullres=value.spot_diameter_fullres;
% scalef.fiducial_diameter_fullres=value.fiducial_diameter_fullres;
% scalef.tissue_lowres_scalef=value.tissue_lowres_scalef;
position_file = fullfile(imgfolder, sprintf('%stissue_positions_list.csv', aff));
if ~exist(position_file, 'file')
position_file = fullfile(imgfolder, sprintf('%stissue_positions.csv', aff));
T = readtable(position_file, 'ReadVariableNames', true);
else
T = readtable(position_file, 'ReadVariableNames', false);
end
[X, genelist, celllist] = sc_readhdf5file(h5fname);
sce = SingleCellExperiment(X, genelist);
sce.c_cell_id = celllist;
metainfo = sprintf("Source: %s", selpath);
sce = sce.appendmetainfo(metainfo);
if ismember('barcode', T.Properties.VariableNames)
assert(all(ismember(sce.c_cell_id, string(T.barcode))))
[~, idx] = ismember(sce.c_cell_id, string(T.barcode));
else
assert(all(ismember(sce.c_cell_id, string(T.Var1))))
[~, idx] = ismember(sce.c_cell_id, string(T.Var1));
end
t = T(idx, :);
if ismember('pxl_row_in_fullres', T.Properties.VariableNames) && ...
ismember('pxl_col_in_fullres', T.Properties.VariableNames)
x_pixel = t.pxl_row_in_fullres;
y_pixel = t.pxl_col_in_fullres;
else
%s=[t.Var3,t.Var4];
%ste.s=s;
x_pixel = t.Var5;
y_pixel = t.Var6;
end
xy = [x_pixel, y_pixel];
% sce.s=[x_pixel y_pixel];
xy = xy - mean(xy);
r = [range(xy(:, 1)), range(xy(:, 2))];
a = size(img, 1:2);
xy = xy .* (0.65 * a ./ r);
xy = xy + a / 2;
ste = SpatialTranscriptomicsExperiment(sce, xy, img);
ste.tissue_positions_list = T;
ste.scalefactors_json = scalef;
metainfo = sprintf("Source: %s", selpath);
ste = ste.appendmetainfo(metainfo);
function [out, aff] = i_guessmtxfile(selpath)
out = [];
aff = [];
a = dir(selpath);
for k = 1:length(a)
if contains(a(k).name, 'filtered_feature_bc_matrix.h5')
out = a(k).name;
aff = extractBefore(out, 'filtered_feature_bc_matrix.h5');
continue;
end
end
end % guess_mtx
end