diff --git a/README.md b/README.md index 87ee2ed..70b4da6 100644 --- a/README.md +++ b/README.md @@ -28,13 +28,27 @@ Module for single-cell data extraction given a segmentation mask and multi-chann See https://en.wikipedia.org/wiki/Gini_coefficient for more information. * `--intensity_props intensity_median` : Will calculate the median of intensity values per labeled object in the mask. * `--intensity_props intensity_sum` : Will calculate the sum of intensity values per labelled object in the mask. This can be useful if you want to count RNA molecules from FISH based images for example. + * `--intensity_props intensity_std` : Will calculate the standard deviation of intensity values per labeled object in the mask. + + Further available are GLCM-derived graycoprops (see https://scikit-image.org/docs/stable/api/skimage.feature.html). Currently these metrices can only be calculated for one angle and distance at a time. + * `--intensity_props contrast` : Will calculate the glcm derived symmetric and normalized contrast of intensity values per labeled object in the mask. + * `--intensity_props dissimilarity` : Will calculate the glcm derived symmetric and normalized dissimilarity of intensity values per labeled object in the mask. + * `--intensity_props homogeneity` : Will calculate the glcm derived symmetric and normalized homogeneity of intensity values per labeled object in the mask. + * `--intensity_props energy` : Will calculate the glcm derived symmetric and normalized energy of intensity values per labeled object in the mask. + * `--intensity_props correlation` : Will calculate the glcm derived symmetric and normalized correlation of intensity values per labeled object in the mask. + * `--intensity_props ASM` : Will calculate the glcm derived symmetric and normalized ASM of intensity values per labeled object in the mask. + + Parameters for GLCM calculation + * `--glcm_angle` Angle in radians used for calculating the GLCM per label. Default is 0 radians. + * `--glcm_distance` Distance in pixels used for calculating the GLCM per label. Default is 1 pixel. + # Run script -`python CommandSingleCellExtraction.py --masks ./segmentation/cellMask.tif ./segmentation/membraneMask.tif --image ./registration/Exemplar_001.h5 --output ./feature_extraction --channel_names ./my_channels.csv` +`python CLI.py --masks ./segmentation/cellMask.tif ./segmentation/membraneMask.tif --image ./registration/Exemplar_001.h5 --output ./feature_extraction --channel_names ./my_channels.csv` # Main developer Denis Schapiro (https://github.com/DenisSch) Joshua Hess (https://github.com/JoshuaHess12) -Jeremy Muhlich (https://github.com/jmuhlich) \ No newline at end of file +Jeremy Muhlich (https://github.com/jmuhlich) diff --git a/mcquant/ParseInput.py b/mcquant/ParseInput.py index 1b98f5c..9318bc5 100644 --- a/mcquant/ParseInput.py +++ b/mcquant/ParseInput.py @@ -1,6 +1,6 @@ #Functions for parsing command line arguments for ome ilastik prep import argparse -from . import __version__ +from . import __version__ # This still has to be adjusted with __init__.py function def ParseInputDataExtract(): @@ -29,11 +29,29 @@ def ParseInputDataExtract(): By default only mean intensity is calculated. If the metric doesn't depend on signal intensity, use --mask-props instead. See list at https://scikit-image.org/docs/dev/api/skimage.measure.html#regionprops - Additionally available is gini_index, which calculates a single number - between 0 and 1, representing how unequal the signal is distributed in each region. - See https://en.wikipedia.org/wiki/Gini_coefficient + + Additionally available is: gini_index, intensity_median, intensity_sum, intensity_std, + contrast, dissimilarity, homogeneity, energy, correlation, ASM. + Further information on the parameters: + gini_index: + which calculates a single number between 0 and 1, representing how unequal the signal is distributed in each region. + Will be calculated for every marker + See https://en.wikipedia.org/wiki/Gini_coefficient + contrast, dissimilarity, homogeneity, energy, correlation, ASM: + glcm/graycoprops features from scikit-image features. The default distance is set to 1 pixel and the default angle is set to 0 rad. + Both parameters can be controlled via CLI inputs. However, both parameters are limited to 1 input each. + Will be calculated for every marker + See https://scikit-image.org/docs/stable/api/skimage.feature.html """ ) + parser.add_argument( + '--glcm_angle', type=float, default=0, + help="Angle in radians for GLCM calculation. Default is 0 radians. Currently limited to 1 angle" + ) + parser.add_argument( + '--glcm_distance', type=int, default=1, + help="Distance in pixels for GLCM calculation. Default is 1 pixel." + ) #parser.add_argument('--suffix') parser.add_argument('--version', action='version', version=f'mcquant {__version__}') args = parser.parse_args() @@ -42,6 +60,9 @@ def ParseInputDataExtract(): 'channel_names': args.channel_names,'output':args.output, 'intensity_props': set(args.intensity_props if args.intensity_props is not None else []).union(["intensity_mean"]), 'mask_props': args.mask_props, + 'glcm_angle': args.glcm_angle, + 'glcm_distance': args.glcm_distance, + } #Print the dictionary object print(dict) diff --git a/mcquant/SingleCellDataExtraction.py b/mcquant/SingleCellDataExtraction.py index 2f97351..1541831 100644 --- a/mcquant/SingleCellDataExtraction.py +++ b/mcquant/SingleCellDataExtraction.py @@ -9,8 +9,8 @@ import os import skimage.measure import skimage.measure._regionprops +from skimage.feature import graycomatrix, graycoprops import tifffile - from pathlib import Path #### Additional functions that can be specified by the user via intensity_props @@ -23,6 +23,10 @@ def intensity_median(mask, intensity): def intensity_sum(mask, intensity): return np.sum(intensity[mask]) +## Function to calculate standard deviation of intensity values +def intensity_std(mask, intensity): + return np.std(intensity[mask]) + ## Function to calculate the gini index: https://en.wikipedia.org/wiki/Gini_coefficient def gini_index(mask, intensity): x = intensity[mask] @@ -31,7 +35,8 @@ def gini_index(mask, intensity): cumx = np.cumsum(sorted_x, dtype=float) return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n -def MaskChannel(mask_loaded, image_loaded_z, intensity_props=["intensity_mean"]): + +def MaskChannel(mask_loaded, image_loaded_z, intensity_props=["intensity_mean"], glcm_angle=0, glcm_distance=1): """Function for quantifying a single channel image Returns a table with CellID according to the mask and the mean pixel intensity @@ -41,11 +46,37 @@ def MaskChannel(mask_loaded, image_loaded_z, intensity_props=["intensity_mean"]) builtin_props = set(intensity_props).intersection(standard_props) # Otherwise look for them in this module extra_props = set(intensity_props).difference(standard_props) + + # All possible scikit-image gracoprops features from glcm + glcm_features_set = {"contrast", "dissimilarity", "homogeneity", "energy", "correlation", "ASM"} + glcm_features = {} + + # Function to calculate the glcm and associated graycoprops per label only if user inputs a glcm metric into --intensity_props + if glcm_features_set.intersection(intensity_props): + for region in skimage.measure.regionprops(mask_loaded, image_loaded_z): + # Get the label/cell + label = region.label + # Rescale the image to uint8, which is needed for glcm calculation if level attribute is not set. + image_uint8 = ((region.intensity_image - np.min(region.intensity_image)) * (255 / (np.max(region.intensity_image) - np.min(region.intensity_image)))).astype(np.uint8) + # Calculate the glcm once per label/cell + glcm = graycomatrix(image_uint8, distances = [glcm_distance], angles = [glcm_angle], symmetric = True, normed = True) + glcm_props = {} + # Calculate the user-specified feature(s) per label/cell + for prop in glcm_features_set.intersection(intensity_props): + glcm_props[prop] = graycoprops(glcm, prop)[0, 0] + glcm_features[label] = glcm_props + dat = skimage.measure.regionprops_table( mask_loaded, image_loaded_z, properties = tuple(builtin_props), - extra_properties = [globals()[n] for n in extra_props] + extra_properties = [globals()[n] for n in extra_props if n not in glcm_features_set] ) + + # Add the glcm dict to the regionproperties + if glcm_features: + for label, features in glcm_features.items(): + for prop, value in features.items(): + dat[prop] = dat.get(prop, []) + [value] return dat @@ -88,7 +119,7 @@ def n_channels(image): image_path = Path(image) - if image_path.suffix in ['.tiff', '.tif', '.btf']: + if image_path.suffix in ['.tiff', '.tif', '.btf', 'qptiff']: s = tifffile.TiffFile(image).series[0] ndim = len(s.shape) if ndim == 2: return 1 @@ -105,12 +136,12 @@ def n_channels(image): def PrepareData(image,z): """Function for preparing input for maskzstack function. Connecting function - to use with mc micro ilastik pipeline""" + to use with mcmicro ilastik pipeline""" image_path = Path(image) #Check to see if image tif(f) - if image_path.suffix in ['.tiff', '.tif', '.btf']: + if image_path.suffix in ['.tiff', '.tif', '.btf', 'qptiff']: image_loaded_z = tifffile.imread(image, key=z) #Check to see if image is hdf5 @@ -129,7 +160,7 @@ def PrepareData(image,z): return image_loaded_z -def MaskZstack(masks_loaded,image,channel_names_loaded, mask_props=None, intensity_props=["intensity_mean"]): +def MaskZstack(masks_loaded,image,channel_names_loaded, mask_props=None, intensity_props=["intensity_mean"], glcm_angle = 0, glcm_distance = 1): """This function will extract the stats for each cell mask through each channel in the input image @@ -151,7 +182,7 @@ def MaskZstack(masks_loaded,image,channel_names_loaded, mask_props=None, intensi for nm in range(len(mask_names)): #Use the above information to mask z stack dict_of_chan[mask_names[nm]].append( - MaskChannel(masks_loaded[mask_names[nm]],image_loaded_z, intensity_props=intensity_props) + MaskChannel(masks_loaded[mask_names[nm]],image_loaded_z, intensity_props=intensity_props, glcm_angle=glcm_angle, glcm_distance=glcm_distance) ) #Print progress print("Finished "+str(z)) @@ -198,7 +229,7 @@ def col_sort(x): # Return the dict of dataframes for each mask return dict_of_chan -def ExtractSingleCells(masks,image,channel_names,output, mask_props=None, intensity_props=["intensity_mean"]): +def ExtractSingleCells(masks,image,channel_names,output, mask_props=None, intensity_props=["intensity_mean"], glcm_angle = 0, glcm_distance = 1): """Function for extracting single cell information from input path containing single-cell masks, z_stack path, and channel_names path.""" @@ -242,7 +273,7 @@ def ExtractSingleCells(masks,image,channel_names,output, mask_props=None, intens m_name = m_full_name.split('.')[0] masks_loaded.update({str(m_name):skimage.io.imread(m,plugin='tifffile')}) - scdata_z = MaskZstack(masks_loaded,image,channel_names_loaded_checked, mask_props=mask_props, intensity_props=intensity_props) + scdata_z = MaskZstack(masks_loaded,image,channel_names_loaded_checked, mask_props=mask_props, intensity_props=intensity_props, glcm_angle=glcm_angle, glcm_distance=glcm_distance) #Write the singe cell data to a csv file using the image name # Determine the image name by cutting off its extension @@ -262,14 +293,14 @@ def ExtractSingleCells(masks,image,channel_names,output, mask_props=None, intens ) -def MultiExtractSingleCells(masks,image,channel_names,output, mask_props=None, intensity_props=["intensity_mean"]): +def MultiExtractSingleCells(masks,image,channel_names,output, mask_props=None, intensity_props=["intensity_mean"], glcm_angle = 0, glcm_distance = 1): """Function for iterating over a list of z_stacks and output locations to export single-cell data from image masks""" print("Extracting single-cell data for "+str(image)+'...') #Run the ExtractSingleCells function for this image - ExtractSingleCells(masks,image,channel_names,output, mask_props=mask_props, intensity_props=intensity_props) + ExtractSingleCells(masks,image,channel_names,output, mask_props=mask_props, intensity_props=intensity_props, glcm_angle=glcm_angle, glcm_distance=glcm_distance) #Print update im_full_name = os.path.basename(image) diff --git a/mcquant/__init__.py b/mcquant/__init__.py index 7dd1069..23bdad8 100644 --- a/mcquant/__init__.py +++ b/mcquant/__init__.py @@ -1,4 +1,4 @@ try: - from ._version import version as __version__ + from ._version import version as __version__ # Version file? except ImportError: __version__ = "unknown"