Merge pull request #8 from wanglab-broad/dev

27. fix environment.yaml

Author: jiahaoh

code-base/src/FilterReads.m

@@ -23,7 +23,7 @@
         numel(barcodes));
     fprintf(s);
 
-    obj.signal.scores = [score_1]; 
+    obj.signal.scores = [obj.signal.scores score_1]; 
 
     for i=1:numel(end_base)
 

code-base/src/FilterReadsMultiSegment.m

@@ -37,6 +37,19 @@
 
     fprintf('Filtration Statistics:\n');
 
+    % filter reads based on codebook
+    codebook_barcodes = obj.codebook.seqToGene.keys;
+    barcodes_in_codebook = contains(color_seq, codebook_barcodes);
+    
+    score_1 = sum(barcodes_in_codebook)/numel(color_seq);
+    s = sprintf('%f [%d / %d] percent of good reads are in codebook\n',...
+        sum(barcodes_in_codebook)/numel(color_seq),...
+        sum(barcodes_in_codebook),...
+        numel(color_seq));
+    fprintf(s);
+
+    obj.signal.scores = [obj.signal.scores score_1]; 
+
     for i=1:numel(end_base)
 
         current_end_base_char = end_base_char(i, :);
@@ -50,20 +63,9 @@
             current_end_base_char(1),...
             current_end_base_char(2));
         fprintf(s);
-    end
 
-    % filter reads based on codebook
-    codebook_barcodes = obj.codebook.seqToGene.keys;
-    barcodes_in_codebook = contains(color_seq, codebook_barcodes);
-    
-    score_1 = sum(barcodes_in_codebook)/numel(color_seq);
-    s = sprintf('%f [%d / %d] percent of good reads are in codebook\n',...
-        sum(barcodes_in_codebook)/numel(color_seq),...
-        sum(barcodes_in_codebook),...
-        numel(color_seq));
-    fprintf(s);
-
-    obj.signal.scores = [score_1]; 
+        obj.signal.scores = [obj.signal.scores current_score]; 
+    end
 
     obj.signal.goodSpots = obj.signal.allSpots(barcodes_in_codebook, :);
     % obj.signal.goodSpots{:, "barcode"} = barcodes(barcodes_in_codebook);

code-base/src/MakeSubtileTable.m

@@ -7,6 +7,7 @@
 
     output_table = table([],[],[],[],[],[],[],[],[],[],[],...
                         'VariableNames', column_headers);
+
     sub_order = [];
     for i = 0:(sqrt_pieces-1)
         for j = 0:(sqrt_pieces-1)
@@ -25,15 +26,15 @@
         start_coords_y = tile_idx(2) * tile_size - overlap_half + 1;
         end_coords_y = (tile_idx(2) + 1) * tile_size + overlap_half;
 
-        %% compensate in edge
+        % compensate in edge
         if tile_idx(1) == 0
             start_coords_x = start_coords_x + overlap_half;
         end
         if tile_idx(2) == 0
             start_coords_y = start_coords_y + overlap_half;
 
         end
-        %% compensate in edge
+        % compensate in edge
         if tile_idx(1) == sqrt_pieces - 1
             end_coords_x = dims(1);
         end
@@ -46,7 +47,7 @@
 
         dims_t = dims;
         dims_t(1:2) = [end_coords_x - start_coords_x + 1,end_coords_y - start_coords_y + 1];
-        disp([tile_idx,start_coords_x,end_coords_x,start_coords_y,end_coords_y,upper_left(1:2),dims_t(1:2)]);
+        % disp([tile_idx,start_coords_x,end_coords_x,start_coords_y,end_coords_y,upper_left(1:2),dims_t(1:2)]);
         output_table_t = table(t, tile_idx(1), tile_idx(2),...
                                 start_coords_x, start_coords_y, end_coords_x, end_coords_y,...
                                 upper_left(1), upper_left(2), dims_t(1), dims_t(2),...

code-base/src/STARMapDataset.m

@@ -73,11 +73,13 @@
             obj.images = dictionary();
             obj.projections = dictionary();
             obj.metadata = dictionary();
+            obj.subtile = struct();
             obj.layers = struct();
             obj.registration = dictionary();
             obj.signal = struct();
             obj.codebook = struct();
 
+            obj.subtile.index = 0;
             obj.layers.seq = [];
             obj.layers.other = [];
             obj.layers.ref = [];
@@ -720,12 +722,16 @@
             defaultSplitIndex = [5];
             addParameter(p, 'split_index', defaultSplitIndex);
 
+            defaultSaveScores = false;
+            addParameter(p, 'save_scores', defaultSaveScores);
+
             parse(p, varargin{:});
 
             fprintf('====Reads Filtration====\n');
             end_base = string(p.Results.end_base);
 
             % remove reads with incorrect colors
+            n_spots = size(obj.signal.allSpots, 1);
             no_color_spots = contains(obj.signal.allSpots.color_seq, "N");
             multi_color_spots = contains(obj.signal.allSpots.color_seq, "M");
             to_keep = ~or(no_color_spots, multi_color_spots);
@@ -734,6 +740,7 @@
             fprintf(sprintf('Number of spots were dropped because of multi-max color: %d\n', sum(multi_color_spots)));
             fprintf('Comparing with codebook...\n');
             fprintf(sprintf('Number of barcode segments: %d\n', p.Results.n_barcode_segments));
+            obj.signal.scores = [n_spots sum(no_color_spots) sum(multi_color_spots)]; 
 
             if numel(end_base) > 1
                 end_base_msg = strjoin(end_base, " or ");
@@ -748,6 +755,23 @@
                 obj = FilterReadsMultiSegment(obj, end_base, p.Results.split_index);
             end
 
+            if p.Results.save_scores
+                current_fname = fullfile(obj.outputPath, "log", "rsf_scores.csv");
+                current_output_folder_msg = strrep(current_fname, '\', '\\');
+                fprintf(sprintf('Saving scores to %s\n', current_output_folder_msg));
+                if obj.subtile.index > 0
+                    scores_to_save = [string(obj.fovID) string(obj.subtile.index) string(obj.signal.scores)];
+                else
+                    scores_to_save = [string(obj.fovID) string(obj.signal.scores)];
+                end
+
+                if ~exist(current_fname, 'file')
+                    writematrix(scores_to_save, current_fname, 'Delimiter', ',');
+                else
+                    writematrix(scores_to_save, current_fname, 'Delimiter', ',', 'WriteMode', 'append');
+                end
+            end
+
             % change metadata
             obj.jobFinished.ReadsFiltration = true;
 
@@ -847,13 +871,21 @@
 
             switch p.Results.signal_slot
                 case "goodSpots"
-                    current_fname = fullfile(p.Results.output_path, sprintf("%s_goodSpots.csv", obj.fovID));
+                    if obj.subtile.index > 0
+                        current_fname = fullfile(obj.outputPath, "output", "subtile", obj.fovID, sprintf("subtile_goodSpots_%d.csv", obj.subtile.index));
+                    else
+                        current_fname = fullfile(p.Results.output_path, sprintf("%s_goodSpots.csv", obj.fovID));
+                    end
                     current_output_folder_msg = strrep(current_fname, '\', '\\');
                     fprintf(sprintf('Saving goodSpots to %s\n', current_output_folder_msg));
                     output_table = obj.signal.goodSpots(:, p.Results.field_to_keep);
                     writetable(output_table, current_fname);
                 case "allSpots"
-                    current_fname = fullfile(p.Results.output_path, sprintf("%s_allSpots.csv", obj.fovID));
+                    if obj.subtile.index > 0
+                        current_fname = fullfile(obj.outputPath, "output", "subtile", obj.fovID, sprintf("subtile_allSpots_%d.csv", obj.subtile.index));
+                    else
+                        current_fname = fullfile(p.Results.output_path, sprintf("%s_allSpots.csv", obj.fovID));
+                    end
                     current_output_folder_msg = strrep(current_fname, '\', '\\');
                     fprintf(sprintf('Saving allSpots to %s\n', current_output_folder_msg));
                     output_table = obj.signal.allSpots(:, p.Results.field_to_keep);
@@ -879,7 +911,7 @@
             defaultSave = false;
             addOptional(p, 'save', defaultSave);
 
-            defaultOutputFolder = fullfile(obj.outputPath, "temp");
+            defaultOutputFolder = fullfile(obj.outputPath, "output");
             addOptional(p, 'output_path', defaultOutputFolder);
 
             defaultRefLayer = obj.layers.ref;
@@ -891,12 +923,41 @@
                 mkdir(p.Results.output_path)
             end
 
+            current_subtile_folder = fullfile(p.Results.output_path, "subtile");
+            if ~exist(current_subtile_folder, 'dir')
+                mkdir(current_subtile_folder)
+            end
+
             current_metadata = obj.metadata{p.Results.ref_layer};
-            current_metadata = current_metadata{1};
-            obj.subtile = MakeSubtileTable(current_metadata.dims, p.Results.sqrt_pieces);     
+            obj.subtile.coords = MakeSubtileTable(current_metadata.dims, p.Results.sqrt_pieces);     
 
             if p.Results.save
-                writetable(obj.subtile, fullfile(p.Results.output_path, sprintf("%s_subtile.csv", obj.fovID)), 'Delimiter',',', 'QuoteStrings',false);
+                current_output_folder = fullfile(current_subtile_folder, obj.fovID);
+                if ~exist(current_output_folder, 'dir')
+                    mkdir(current_output_folder)
+                end
+                % save the coords of subtiles
+                writetable(obj.subtile.coords, fullfile(current_output_folder, "subtile_coords.csv"), 'Delimiter',',', 'QuoteStrings',false);
+
+                % save the object of subtiles
+                for i=1:size(obj.subtile.coords)
+                    tile_idx = obj.subtile.coords.t(i);
+                    start_coords_x = obj.subtile.coords.scoords_x(i);
+                    start_coords_y = obj.subtile.coords.scoords_y(i);
+
+                    end_coords_x = obj.subtile.coords.ecoords_x(i);
+                    end_coords_y = obj.subtile.coords.ecoords_y(i);
+
+                    subtile_data = obj;
+                    for l=1:length(obj.layers.seq)
+                        current_layer = obj.layers.seq{l};
+                        subtile_data.images{current_layer} = subtile_data.images{current_layer}(start_coords_y:end_coords_y, start_coords_x:end_coords_x, :, :);
+                    end
+                    subtile_data.subtile.index = tile_idx;
+
+                    save(fullfile(current_output_folder, sprintf('subtile_data_%d.mat', tile_idx)), "subtile_data");
+                end
+
             end
 
             % change metadata

config/environment.yaml

@@ -375,7 +375,6 @@ dependencies:
     - tifffile==2023.4.12
     - tqdm==4.65.0
     - typeguard==4.1.5
-    - typing-extensions==4.5.0
     - umap-learn==0.5.3
     - visions==0.7.5
     - werkzeug==2.2.3