8207753731

Committed 08 Mar 2024 06:56PM UTC coverage: 93.468% (+0.006%) from 93.462%

Build # 8207753731

Build Type

Pull #1124

github

Committed by

web-flow

Commit Message

Merge d3cd7a3a9 into 887721348

Pull Request Pull Request #1124: Don't require user to have all pixel SOM clusters represented during assignment

Run Details

2122 of 2381 branches covered (89.12%)

Branch coverage included in aggregate %.

2 of 3 new or added lines in 1 file covered. (66.67%)

18 existing lines in 2 files now uncovered.

4088 of 4263 relevant lines covered (95.89%)

5.72 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.61

/src/ark/phenotyping/pixel_cluster_utils.py

import os
import random
import warnings
from typing import List

import feather
import numpy as np
import pandas as pd
import scipy.ndimage as ndimage
from alpineer import image_utils, io_utils, load_utils, misc_utils
from pyarrow.lib import ArrowInvalid
from skimage.io import imread
import natsort as ns


def calculate_channel_percentiles(tiff_dir, fovs, channels, img_sub_folder,
                                  percentile):
    """Calculates average percentile for each channel in the dataset

    Args:
        tiff_dir (str):
            Name of the directory containing the tiff files
        fovs (list):
            List of fovs to include
        channels (list):
            List of channels to include
        img_sub_folder (str):
            Sub folder within each FOV containing image data
        percentile (float):
            The specific percentile to compute

    Returns:
        pd.DataFrame:
            The mapping between each channel and its normalization value
    """

    # create list to hold percentiles
    percentile_means = []

    # loop over channels and FOVs
    for channel in channels:
        percentile_list = []
        for fov in fovs:
            # load image data and remove 0 valued pixels
            img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,
                                                 channels=[channel], fovs=[fov]).values[0, :, :, 0]
            img = img[img > 0]

            # record and store percentile, skip if no non-zero pixels
            if len(img) > 0:
                img_percentile = np.quantile(img, percentile)
                percentile_list.append(img_percentile)

        # save channel-wide average
        percentile_means.append(np.mean(percentile_list))

    percentile_df = pd.DataFrame(np.expand_dims(percentile_means, axis=0), columns=channels)
    percentile_df.sort_index(axis="columns", key=ns.natsort_key, inplace=True)

    return percentile_df


def calculate_pixel_intensity_percentile(tiff_dir, fovs, channels, img_sub_folder,
                                         channel_percentiles, percentile=0.05):
    """Calculates average percentile per FOV for total signal in each pixel

    Args:
        tiff_dir (str):
            Name of the directory containing the tiff files
        fovs (list):
            List of fovs to include
        channels (list):
            List of channels to include
        img_sub_folder (str):
            Sub folder within each FOV containing image data
        channel_percentiles (pd.DataFrame):
            The mapping between each channel and its normalization value
            Computed by `calculate_channel_percentiles`
        percentile (float):
            The pixel intensity percentile per FOV to average over


    Returns:
        float:
            The average percentile per FOV for total signal in each pixel
    """

    # create vector of channel percentiles to enable broadcasting
    norm_vect = channel_percentiles.iloc[0].values
    norm_vect = norm_vect.reshape([1, 1, len(norm_vect)])

    intensity_percentile_list = []

    for fov in fovs:
        # load image data
        img_data = load_utils.load_imgs_from_tree(data_dir=tiff_dir, fovs=[fov],
                                                  channels=channels, img_sub_folder=img_sub_folder)

        # normalize each channel by its percentile value
        norm_data = img_data[0].values / norm_vect

        # sum channels together to determine total intensity
        summed_data = np.sum(norm_data, axis=-1)
        intensity_percentile_list.append(np.quantile(summed_data, percentile))

    return np.mean(intensity_percentile_list)


def normalize_rows(pixel_data, channels, include_seg_label=True):
    """Normalizes the rows of a pixel matrix by their sum

    Args:
        pixel_data (pandas.DataFrame):
            The dataframe containing the pixel data for a given fov
            Includes channel and meta (`fov`, `label`, etc.) columns
        channels (list):
            List of channels to subset over
        include_seg_label (bool):
            Whether to include `'label'` as a metadata column

    Returns:
        pandas.DataFrame:
            The pixel data with rows normalized and 0-sum rows removed
    """

    # subset the fov data by the channels the user trained the pixel SOM on
    pixel_data_sub = pixel_data[channels]

    # divide each row by their sum
    pixel_data_sub = pixel_data_sub.div(pixel_data_sub.sum(axis=1), axis=0)

    # define the meta columns to add back
    meta_cols = ['fov', 'row_index', 'column_index']

    # add the label column if it should be kept
    if include_seg_label:
        meta_cols.append('label')

    # add back meta columns, making sure to remove 0-row indices
    pixel_data_sub[meta_cols] = pixel_data.loc[pixel_data_sub.index.values, meta_cols]

    return pixel_data_sub


def check_for_modified_channels(tiff_dir, test_fov, img_sub_folder, channels):
    """Checks to make sure the user selected newly modified channels

    Args:
        tiff_dir (str):
            Name of the directory containing the tiff files
        test_fov (str):
            example fov used to check channel names
        img_sub_folder (str):
            sub-folder within each FOV containing image data
        channels (list):
            list of channels to use for analysis
    """

    # convert to path-compatible format
    if img_sub_folder is None:
        img_sub_folder = ''

    # get all channels within example FOV
    all_channels = io_utils.list_files(os.path.join(tiff_dir, test_fov, img_sub_folder))
    all_channels = io_utils.remove_file_extensions(all_channels)
    # define potential modifications to channel names
    mods = ['_smoothed', '_nuc_include', '_nuc_exclude']

    # loop over each user-provided channel
    for channel in channels:
        for mod in mods:
            # check for substring matching
            chan_mod = channel + mod
            if chan_mod in all_channels:
                warnings.warn('You selected {} as the channel to analyze, but there were potential'
                              ' modified channels found: {}. Make sure you selected the correct '
                              'version of the channel for inclusion in '
                              'clustering'.format(channel, chan_mod))
            else:
                pass


def smooth_channels(fovs, tiff_dir, img_sub_folder, channels, smooth_vals):
    """Adds additional smoothing for selected channels as a preprocessing step

    Args:
        fovs (list):
            List of fovs to process
        tiff_dir (str):
            Name of the directory containing the tiff files
        img_sub_folder (str):
            sub-folder within each FOV containing image data
        channels (list):
            list of channels to apply smoothing to
        smooth_vals (list or int):
            amount to smooth channels. If a single int, applies
            to all channels. Otherwise, a custom value per channel can be supplied

    """

    # no output if no channels specified
    if channels is None or len(channels) == 0:
        return

    # convert to path-compatible format
    if img_sub_folder is None:
        img_sub_folder = ''

    # convert int to list of same length
    if type(smooth_vals) is int:
        smooth_vals = [smooth_vals for _ in range(len(channels))]
    elif type(smooth_vals) is list:
        if len(smooth_vals) != len(channels):
            raise ValueError("A list was provided for variable smooth_vals, but it does not "
                             "have the same length as the list of channels provided")
    else:
        raise ValueError("Variable smooth_vals must be either a single integer or a list")

    for fov in fovs:
        for idx, chan in enumerate(channels):
            img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,
                                                 fovs=[fov], channels=[chan]).values[0, :, :, 0]
            chan_out = ndimage.gaussian_filter(img, sigma=smooth_vals[idx])
            image_utils.save_image(
                os.path.join(tiff_dir, fov, img_sub_folder, chan + '_smoothed.tiff'),
                chan_out
            )


def filter_with_nuclear_mask(fovs: List, tiff_dir: str, seg_dir: str, channel: str,
                             nuc_seg_suffix: str = "_nuclear.tiff", img_sub_folder: str = None,
                             exclude: bool = True):
    """Filters out background staining using subcellular marker localization.

    Non-nuclear signal is removed from nuclear markers and vice-versa for membrane markers.

    Args:
        fovs (list):
            The list of fovs to filter
        tiff_dir (str):
            Name of the directory containing the tiff files
        seg_dir (str):
            Name of the directory containing the segmented files
        channel (str):
            Channel to apply filtering to
        nuc_seg_suffix (str):
            The suffix for the nuclear channel.
            (i.e. for "fov1", a suffix of "_nuclear.tiff" would make a file named
            "fov1_nuclear.tiff")
        img_sub_folder (str):
            Name of the subdirectory inside `tiff_dir` containing the tiff files.
            Set to `None` if there isn't any.
        exclude (bool):
            Whether to filter out nuclear or membrane signal
    """

    # if seg_dir is None, the user cannot run filtering
    if seg_dir is None:
        print('No seg_dir provided, you must provide one to run nuclear filtering')
        return

    # raise an error if the provided seg_dir does not exist
    io_utils.validate_paths(seg_dir)

    # convert to path-compatible format
    if img_sub_folder is None:
        img_sub_folder = ''

    for fov in fovs:
        # load the channel image in
        img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,
                                             fovs=[fov], channels=[channel]).values[0, :, :, 0]

        # load the segmented image in
        seg_img_name: str = f"{fov}{nuc_seg_suffix}"
        seg_img = imread(os.path.join(seg_dir, seg_img_name))[0, ...]

        # mask out the nucleus
        if exclude:
            suffix = "_nuc_exclude.tiff"
            seg_mask = seg_img > 0
        # mask out the membrane
        else:
            suffix = "_nuc_include.tiff"
            seg_mask = seg_img == 0

        # filter out the nucleus or membrane depending on exclude parameter
        img[seg_mask] = 0

        # save filtered image
        image_utils.save_image(os.path.join(tiff_dir, fov, img_sub_folder, channel + suffix), img)


def compute_pixel_cluster_channel_avg(fovs, channels, base_dir, pixel_cluster_col,
                                      num_pixel_clusters=None,
                                      pixel_data_dir='pixel_mat_data',
                                      num_fovs_subset=100, seed=42, keep_count=False):
    """Compute the average channel values across each pixel SOM cluster.

    To improve performance, number of FOVs is downsampled by `fov_subset_proportion`

    Args:
        fovs (list):
            The list of fovs to subset on
        channels (list):
            The list of channels to subset on
        base_dir (str):
            The path to the data directories
        pixel_cluster_col (str):
            Name of the column to group by
        num_pixel_clusters (int):
            The number of pixel clusters that are desired, if None then no fixed amount required
        pixel_data_dir (str):
            Name of the directory containing the pixel data with cluster labels
        num_fovs_subset (float):
            The number of FOVs to subset on. Note that if `len(fovs) < num_fovs_subset`, all of
            the FOVs will still be selected
        seed (int):
            The random seed to use for subsetting FOVs
        keep_count (bool):
            Whether to keep the count column when aggregating or not
            This should only be set to `True` for visualization purposes

    Returns:
        pandas.DataFrame:
            Contains the average channel values for each pixel SOM/meta cluster
    """

    # verify the pixel cluster col specified is valid
    misc_utils.verify_in_list(
        provided_cluster_col=[pixel_cluster_col],
        valid_cluster_cols=['pixel_som_cluster', 'pixel_meta_cluster']
    )

    # verify num_pixel_clusters is valid if set
    if num_pixel_clusters is not None and num_pixel_clusters <= 0:
        raise ValueError("If set, number of pixel clusters desired must be a positive integer")

    # verify fovs subset value is valid
    if num_fovs_subset <= 0:
        raise ValueError("Number of fovs to subset must be a positive integer")

    # define a list to hold the cluster averages for each FOV
    fov_cluster_avgs = []

    # warn the user that we can only select so many FOVs if len(fovs) < num_fovs_subset
    if len(fovs) < num_fovs_subset:
        warnings.warn(
            'Provided num_fovs_subset=%d but only %d FOVs in dataset, '
            'subsetting just the %d FOVs' %
            (num_fovs_subset, len(fovs), len(fovs))
        )

    # subset number of FOVs based on num_fovs_subset if less than total number of fovs
    random.seed(seed)
    fovs_sub = random.sample(fovs, num_fovs_subset) if num_fovs_subset < len(fovs) else fovs

    for fov in fovs_sub:
        # read in the fovs data
        try:
            fov_pixel_data = feather.read_dataframe(
                os.path.join(base_dir, pixel_data_dir, fov + '.feather')
            )
        except (ArrowInvalid, OSError, IOError):
            print("The data for FOV %s has been corrupted, skipping" % fov)
            continue

        # aggregate the sums and counts
        sum_by_cluster = fov_pixel_data.groupby(
            pixel_cluster_col
        )[channels].sum()
        count_by_cluster = fov_pixel_data.groupby(
            pixel_cluster_col
        )[channels].size().to_frame('count')

        # merge the results by column
        agg_results = pd.merge(
            sum_by_cluster, count_by_cluster, left_index=True, right_index=True
        ).reset_index()

        # append the result to cluster_avgs
        fov_cluster_avgs.append(agg_results)

    cluster_avgs = pd.concat(fov_cluster_avgs)

    # reset the index of cluster_avgs for consistency
    cluster_avgs = cluster_avgs.reset_index(drop=True)

    # sum the counts and the channel sums
    sum_count_totals = cluster_avgs.groupby(
        pixel_cluster_col
    )[channels + ['count']].sum().reset_index()

    # error out if any clusters were lost during the averaging process
    if num_pixel_clusters is not None and sum_count_totals.shape[0] < num_pixel_clusters:
        raise ValueError(
            'Averaged data contains just %d clusters out of %d. '
            'Average expression file not written. '
            'Consider increasing your num_fovs_subset value.' %
            (sum_count_totals.shape[0], num_pixel_clusters)
        )

    # now compute the means using the count column
    sum_count_totals[channels] = sum_count_totals[channels].div(sum_count_totals['count'], axis=0)

    # convert cluster column to integer type
    sum_count_totals[pixel_cluster_col] = sum_count_totals[pixel_cluster_col].astype(int)

    # sort cluster col in ascending order
    sum_count_totals = sum_count_totals.sort_values(by=pixel_cluster_col)

    # drop the count column if specified
    if not keep_count:
        sum_count_totals = sum_count_totals.drop('count', axis=1)

    return sum_count_totals


def find_fovs_missing_col(base_dir, data_dir, missing_col):
    """Identify FOV names in `data_dir` without `missing_col`

    Args:
        base_dir (str):
            The path to the data directories
        data_dir (str):
            Name of the directory which contains the full preprocessed pixel data
        missing_col (str):
            Name of the column to identify

    Returns:
        list:
            List of FOVs without `missing_col`
    """

    # define the main data path
    data_path = os.path.join(base_dir, data_dir)

    # define the temp data path
    temp_path = os.path.join(base_dir, data_dir + '_temp')

    # verify the data path exists
    io_utils.validate_paths(data_path)

    # if the temp path does not exist, either all the FOVs need to be run or none of them do
    if not os.path.exists(temp_path):
        # read in one of the FOV files in data_path
        fov_files = io_utils.list_files(data_path, substrs='.feather')

        # read in a sample FOV, do it this way to avoid potentially corrupted files
        # NOTE: this assumes only a few files get corrupted each run per MIBIAN runs
        # NOTE: handling of corrupted files gets propagated to respective R script
        i = 0
        while i < len(fov_files):
            try:
                fov_data = feather.read_dataframe(os.path.join(data_path, fov_files[i]))
            except (ArrowInvalid, OSError, IOError):
                i += 1
                continue
            break

        # if the missing_col is not found in fov_data, we need to run all the FOVs
        if missing_col not in fov_data.columns.values:
            # we will also make the temp directory to store the new files
            os.mkdir(temp_path)
            return io_utils.remove_file_extensions(fov_files)
        # otherwise, the column has already been assigned for this cohort, no need to run anything
        else:
            return []
    # if the temp path does exist, we have FOVs that need further processing
    else:
        # retrieve the FOV file names from both data_path and temp_path
        data_files = set(io_utils.list_files(data_path, substrs='.feather'))
        temp_files = set(io_utils.list_files(temp_path, substrs='.feather'))

        # get the difference between the two of these to determine which ones are valid
        leftover_files = list(data_files.difference(temp_files))

        return io_utils.remove_file_extensions(leftover_files)

1	import os	6✔
2	import random	6✔
3	import warnings	6✔
4	from typing import List	6✔
5
6	import feather	6✔
7	import numpy as np	6✔
8	import pandas as pd	6✔
9	import scipy.ndimage as ndimage	6✔
10	from alpineer import image_utils, io_utils, load_utils, misc_utils	6✔
11	from pyarrow.lib import ArrowInvalid	6✔
12	from skimage.io import imread	6✔
13	import natsort as ns	6✔
14
15
16	def calculate_channel_percentiles(tiff_dir, fovs, channels, img_sub_folder,	6✔
17	percentile):
18	"""Calculates average percentile for each channel in the dataset
19
20	Args:
21	tiff_dir (str):
22	Name of the directory containing the tiff files
23	fovs (list):
24	List of fovs to include
25	channels (list):
26	List of channels to include
27	img_sub_folder (str):
28	Sub folder within each FOV containing image data
29	percentile (float):
30	The specific percentile to compute
31
32	Returns:
33	pd.DataFrame:
34	The mapping between each channel and its normalization value
35	"""
36
37	# create list to hold percentiles
38	percentile_means = []	6✔
39
40	# loop over channels and FOVs
41	for channel in channels:	6✔
42	percentile_list = []	6✔
43	for fov in fovs:	6!
44	# load image data and remove 0 valued pixels
45	img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,	6✔
46	channels=[channel], fovs=[fov]).values[0, :, :, 0]
47	img = img[img > 0]	6✔
48
49	# record and store percentile, skip if no non-zero pixels
50	if len(img) > 0:	6✔
51	img_percentile = np.quantile(img, percentile)	6✔
52	percentile_list.append(img_percentile)	6✔
53
54	# save channel-wide average
55	percentile_means.append(np.mean(percentile_list))	6✔
56
57	percentile_df = pd.DataFrame(np.expand_dims(percentile_means, axis=0), columns=channels)	6✔
58	percentile_df.sort_index(axis="columns", key=ns.natsort_key, inplace=True)	6✔
59
60	return percentile_df	6✔
61
62
63	def calculate_pixel_intensity_percentile(tiff_dir, fovs, channels, img_sub_folder,	6✔
64	channel_percentiles, percentile=0.05):
65	"""Calculates average percentile per FOV for total signal in each pixel
66
67	Args:
68	tiff_dir (str):
69	Name of the directory containing the tiff files
70	fovs (list):
71	List of fovs to include
72	channels (list):
73	List of channels to include
74	img_sub_folder (str):
75	Sub folder within each FOV containing image data
76	channel_percentiles (pd.DataFrame):
77	The mapping between each channel and its normalization value
78	Computed by `calculate_channel_percentiles`
79	percentile (float):
80	The pixel intensity percentile per FOV to average over
81
82
83	Returns:
84	float:
85	The average percentile per FOV for total signal in each pixel
86	"""
87
88	# create vector of channel percentiles to enable broadcasting
89	norm_vect = channel_percentiles.iloc[0].values	6✔
90	norm_vect = norm_vect.reshape([1, 1, len(norm_vect)])	6✔
91
92	intensity_percentile_list = []	6✔
93
94	for fov in fovs:	6✔
95	# load image data
96	img_data = load_utils.load_imgs_from_tree(data_dir=tiff_dir, fovs=[fov],	6✔
97	channels=channels, img_sub_folder=img_sub_folder)
98
99	# normalize each channel by its percentile value
100	norm_data = img_data[0].values / norm_vect	6✔
101
102	# sum channels together to determine total intensity
103	summed_data = np.sum(norm_data, axis=-1)	6✔
104	intensity_percentile_list.append(np.quantile(summed_data, percentile))	6✔
105
106	return np.mean(intensity_percentile_list)	6✔
107
108
109	def normalize_rows(pixel_data, channels, include_seg_label=True):	6✔
110	"""Normalizes the rows of a pixel matrix by their sum
111
112	Args:
113	pixel_data (pandas.DataFrame):
114	The dataframe containing the pixel data for a given fov
115	Includes channel and meta (`fov`, `label`, etc.) columns
116	channels (list):
117	List of channels to subset over
118	include_seg_label (bool):
119	Whether to include `'label'` as a metadata column
120
121	Returns:
122	pandas.DataFrame:
123	The pixel data with rows normalized and 0-sum rows removed
124	"""
125
126	# subset the fov data by the channels the user trained the pixel SOM on
127	pixel_data_sub = pixel_data[channels]	6✔
128
129	# divide each row by their sum
130	pixel_data_sub = pixel_data_sub.div(pixel_data_sub.sum(axis=1), axis=0)	6✔
131
132	# define the meta columns to add back
133	meta_cols = ['fov', 'row_index', 'column_index']	6✔
134
135	# add the label column if it should be kept
136	if include_seg_label:	6✔
137	meta_cols.append('label')	6✔
138
139	# add back meta columns, making sure to remove 0-row indices
140	pixel_data_sub[meta_cols] = pixel_data.loc[pixel_data_sub.index.values, meta_cols]	6✔
141
142	return pixel_data_sub	6✔
143
144
145	def check_for_modified_channels(tiff_dir, test_fov, img_sub_folder, channels):	6✔
146	"""Checks to make sure the user selected newly modified channels
147
148	Args:
149	tiff_dir (str):
150	Name of the directory containing the tiff files
151	test_fov (str):
152	example fov used to check channel names
153	img_sub_folder (str):
154	sub-folder within each FOV containing image data
155	channels (list):
156	list of channels to use for analysis
157	"""
158
159	# convert to path-compatible format
160	if img_sub_folder is None:	6✔
161	img_sub_folder = ''	6✔
162
163	# get all channels within example FOV
164	all_channels = io_utils.list_files(os.path.join(tiff_dir, test_fov, img_sub_folder))	6✔
165	all_channels = io_utils.remove_file_extensions(all_channels)	6✔
166	# define potential modifications to channel names
167	mods = ['_smoothed', '_nuc_include', '_nuc_exclude']	6✔
168
169	# loop over each user-provided channel
170	for channel in channels:	6✔
171	for mod in mods:	6✔
172	# check for substring matching
173	chan_mod = channel + mod	6✔
174	if chan_mod in all_channels:	6✔
175	warnings.warn('You selected {} as the channel to analyze, but there were potential'	6✔
176	' modified channels found: {}. Make sure you selected the correct '
177	'version of the channel for inclusion in '
178	'clustering'.format(channel, chan_mod))
179	else:
180	pass	4!
181
182
183	def smooth_channels(fovs, tiff_dir, img_sub_folder, channels, smooth_vals):	6✔
184	"""Adds additional smoothing for selected channels as a preprocessing step
185
186	Args:
187	fovs (list):
188	List of fovs to process
189	tiff_dir (str):
190	Name of the directory containing the tiff files
191	img_sub_folder (str):
192	sub-folder within each FOV containing image data
193	channels (list):
194	list of channels to apply smoothing to
195	smooth_vals (list or int):
196	amount to smooth channels. If a single int, applies
197	to all channels. Otherwise, a custom value per channel can be supplied
198
199	"""
200
201	# no output if no channels specified
202	if channels is None or len(channels) == 0:	6✔
203	return	6✔
204
205	# convert to path-compatible format
206	if img_sub_folder is None:	6✔
207	img_sub_folder = ''	×
208
209	# convert int to list of same length
210	if type(smooth_vals) is int:	6✔
211	smooth_vals = [smooth_vals for _ in range(len(channels))]	6✔
212	elif type(smooth_vals) is list:	6✔
213	if len(smooth_vals) != len(channels):	6✔
214	raise ValueError("A list was provided for variable smooth_vals, but it does not "	6✔
215	"have the same length as the list of channels provided")
216	else:
217	raise ValueError("Variable smooth_vals must be either a single integer or a list")	6✔
218
219	for fov in fovs:	6✔
220	for idx, chan in enumerate(channels):	6✔
221	img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,	6✔
222	fovs=[fov], channels=[chan]).values[0, :, :, 0]
223	chan_out = ndimage.gaussian_filter(img, sigma=smooth_vals[idx])	6✔
224	image_utils.save_image(	6✔
225	os.path.join(tiff_dir, fov, img_sub_folder, chan + '_smoothed.tiff'),
226	chan_out
227	)
228
229
230	def filter_with_nuclear_mask(fovs: List, tiff_dir: str, seg_dir: str, channel: str,	6✔
231	nuc_seg_suffix: str = "_nuclear.tiff", img_sub_folder: str = None,
232	exclude: bool = True):
233	"""Filters out background staining using subcellular marker localization.
234
235	Non-nuclear signal is removed from nuclear markers and vice-versa for membrane markers.
236
237	Args:
238	fovs (list):
239	The list of fovs to filter
240	tiff_dir (str):
241	Name of the directory containing the tiff files
242	seg_dir (str):
243	Name of the directory containing the segmented files
244	channel (str):
245	Channel to apply filtering to
246	nuc_seg_suffix (str):
247	The suffix for the nuclear channel.
248	(i.e. for "fov1", a suffix of "_nuclear.tiff" would make a file named
249	"fov1_nuclear.tiff")
250	img_sub_folder (str):
251	Name of the subdirectory inside `tiff_dir` containing the tiff files.
252	Set to `None` if there isn't any.
253	exclude (bool):
254	Whether to filter out nuclear or membrane signal
255	"""
256
257	# if seg_dir is None, the user cannot run filtering
258	if seg_dir is None:	6✔
259	print('No seg_dir provided, you must provide one to run nuclear filtering')	6✔
260	return	6✔
261
262	# raise an error if the provided seg_dir does not exist
263	io_utils.validate_paths(seg_dir)	6✔
264
265	# convert to path-compatible format
266	if img_sub_folder is None:	6✔
267	img_sub_folder = ''	6✔
268
269	for fov in fovs:	6✔
270	# load the channel image in
271	img = load_utils.load_imgs_from_tree(data_dir=tiff_dir, img_sub_folder=img_sub_folder,	6✔
272	fovs=[fov], channels=[channel]).values[0, :, :, 0]
273
274	# load the segmented image in
275	seg_img_name: str = f"{fov}{nuc_seg_suffix}"	6✔
276	seg_img = imread(os.path.join(seg_dir, seg_img_name))[0, ...]	6✔
277
278	# mask out the nucleus
279	if exclude:	6✔
280	suffix = "_nuc_exclude.tiff"	6✔
281	seg_mask = seg_img > 0	6✔
282	# mask out the membrane
283	else:
284	suffix = "_nuc_include.tiff"	6✔
285	seg_mask = seg_img == 0	6✔
286
287	# filter out the nucleus or membrane depending on exclude parameter
288	img[seg_mask] = 0	6✔
289
290	# save filtered image
291	image_utils.save_image(os.path.join(tiff_dir, fov, img_sub_folder, channel + suffix), img)	6✔
292
293
294	def compute_pixel_cluster_channel_avg(fovs, channels, base_dir, pixel_cluster_col,	6✔
295	num_pixel_clusters=None,
296	pixel_data_dir='pixel_mat_data',
297	num_fovs_subset=100, seed=42, keep_count=False):
298	"""Compute the average channel values across each pixel SOM cluster.
299
300	To improve performance, number of FOVs is downsampled by `fov_subset_proportion`
301
302	Args:
303	fovs (list):
304	The list of fovs to subset on
305	channels (list):
306	The list of channels to subset on
307	base_dir (str):
308	The path to the data directories
309	pixel_cluster_col (str):
310	Name of the column to group by
311	num_pixel_clusters (int):
312	The number of pixel clusters that are desired, if None then no fixed amount required
313	pixel_data_dir (str):
314	Name of the directory containing the pixel data with cluster labels
315	num_fovs_subset (float):
316	The number of FOVs to subset on. Note that if `len(fovs) < num_fovs_subset`, all of
317	the FOVs will still be selected
318	seed (int):
319	The random seed to use for subsetting FOVs
320	keep_count (bool):
321	Whether to keep the count column when aggregating or not
322	This should only be set to `True` for visualization purposes
323
324	Returns:
325	pandas.DataFrame:
326	Contains the average channel values for each pixel SOM/meta cluster
327	"""
328
329	# verify the pixel cluster col specified is valid
330	misc_utils.verify_in_list(	6✔
331	provided_cluster_col=[pixel_cluster_col],
332	valid_cluster_cols=['pixel_som_cluster', 'pixel_meta_cluster']
333	)
334
335	# verify num_pixel_clusters is valid if set
336	if num_pixel_clusters is not None and num_pixel_clusters <= 0:	6✔
NEW 337	raise ValueError("If set, number of pixel clusters desired must be a positive integer")	×
338
339	# verify fovs subset value is valid
340	if num_fovs_subset <= 0:	6✔
341	raise ValueError("Number of fovs to subset must be a positive integer")	×
342
343	# define a list to hold the cluster averages for each FOV
344	fov_cluster_avgs = []	6✔
345
346	# warn the user that we can only select so many FOVs if len(fovs) < num_fovs_subset
347	if len(fovs) < num_fovs_subset:	6✔
348	warnings.warn(	6✔
349	'Provided num_fovs_subset=%d but only %d FOVs in dataset, '
350	'subsetting just the %d FOVs' %
351	(num_fovs_subset, len(fovs), len(fovs))
352	)
353
354	# subset number of FOVs based on num_fovs_subset if less than total number of fovs
355	random.seed(seed)	6✔
356	fovs_sub = random.sample(fovs, num_fovs_subset) if num_fovs_subset < len(fovs) else fovs	6✔
357
358	for fov in fovs_sub:	6✔
359	# read in the fovs data
360	try:	6✔
361	fov_pixel_data = feather.read_dataframe(	6✔
362	os.path.join(base_dir, pixel_data_dir, fov + '.feather')
363	)
364	except (ArrowInvalid, OSError, IOError):
365	print("The data for FOV %s has been corrupted, skipping" % fov)
366	continue
367
368	# aggregate the sums and counts
369	sum_by_cluster = fov_pixel_data.groupby(	6✔
370	pixel_cluster_col
371	)[channels].sum()
372	count_by_cluster = fov_pixel_data.groupby(	6✔
373	pixel_cluster_col
374	)[channels].size().to_frame('count')
375
376	# merge the results by column
377	agg_results = pd.merge(	6✔
378	sum_by_cluster, count_by_cluster, left_index=True, right_index=True
379	).reset_index()
380
381	# append the result to cluster_avgs
382	fov_cluster_avgs.append(agg_results)	6✔
383
384	cluster_avgs = pd.concat(fov_cluster_avgs)	6✔
385
386	# reset the index of cluster_avgs for consistency
387	cluster_avgs = cluster_avgs.reset_index(drop=True)	6✔
388
389	# sum the counts and the channel sums
390	sum_count_totals = cluster_avgs.groupby(	6✔
391	pixel_cluster_col
392	)[channels + ['count']].sum().reset_index()
393
394	# error out if any clusters were lost during the averaging process
395	if num_pixel_clusters is not None and sum_count_totals.shape[0] < num_pixel_clusters:	6✔
396	raise ValueError(	6✔
397	'Averaged data contains just %d clusters out of %d. '
398	'Average expression file not written. '
399	'Consider increasing your num_fovs_subset value.' %
400	(sum_count_totals.shape[0], num_pixel_clusters)
401	)
402
403	# now compute the means using the count column
404	sum_count_totals[channels] = sum_count_totals[channels].div(sum_count_totals['count'], axis=0)	6✔
405
406	# convert cluster column to integer type
407	sum_count_totals[pixel_cluster_col] = sum_count_totals[pixel_cluster_col].astype(int)	6✔
408
409	# sort cluster col in ascending order
410	sum_count_totals = sum_count_totals.sort_values(by=pixel_cluster_col)	6✔
411
412	# drop the count column if specified
413	if not keep_count:	6✔
414	sum_count_totals = sum_count_totals.drop('count', axis=1)	6✔
415
416	return sum_count_totals	6✔
417
418
419	def find_fovs_missing_col(base_dir, data_dir, missing_col):	6✔
420	"""Identify FOV names in `data_dir` without `missing_col`
421
422	Args:
423	base_dir (str):
424	The path to the data directories
425	data_dir (str):
426	Name of the directory which contains the full preprocessed pixel data
427	missing_col (str):
428	Name of the column to identify
429
430	Returns:
431	list:
432	List of FOVs without `missing_col`
433	"""
434
435	# define the main data path
436	data_path = os.path.join(base_dir, data_dir)	6✔
437
438	# define the temp data path
439	temp_path = os.path.join(base_dir, data_dir + '_temp')	6✔
440
441	# verify the data path exists
442	io_utils.validate_paths(data_path)	6✔
443
444	# if the temp path does not exist, either all the FOVs need to be run or none of them do
445	if not os.path.exists(temp_path):	6✔
446	# read in one of the FOV files in data_path
447	fov_files = io_utils.list_files(data_path, substrs='.feather')	6✔
448
449	# read in a sample FOV, do it this way to avoid potentially corrupted files
450	# NOTE: this assumes only a few files get corrupted each run per MIBIAN runs
451	# NOTE: handling of corrupted files gets propagated to respective R script
452	i = 0	6✔
453	while i < len(fov_files):	6✔
454	try:	6✔
455	fov_data = feather.read_dataframe(os.path.join(data_path, fov_files[i]))	6✔
456	except (ArrowInvalid, OSError, IOError):	3✔
457	i += 1	3✔
458	continue	3✔
459	break	4✔
460
461	# if the missing_col is not found in fov_data, we need to run all the FOVs
462	if missing_col not in fov_data.columns.values:	6!
463	# we will also make the temp directory to store the new files
464	os.mkdir(temp_path)	6✔
465	return io_utils.remove_file_extensions(fov_files)	6✔
466	# otherwise, the column has already been assigned for this cohort, no need to run anything
467	else:
468	return []	6✔
469	# if the temp path does exist, we have FOVs that need further processing
470	else:
471	# retrieve the FOV file names from both data_path and temp_path
472	data_files = set(io_utils.list_files(data_path, substrs='.feather'))	6✔
473	temp_files = set(io_utils.list_files(temp_path, substrs='.feather'))	6✔
474
475	# get the difference between the two of these to determine which ones are valid
476	leftover_files = list(data_files.difference(temp_files))	6✔
477
478	return io_utils.remove_file_extensions(leftover_files)	6✔

angelolab / ark-analysis / 8207753731

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous