angelolab / toffy / 405

pending completion

Build # 405

Build Type

Pull #50

travis-ci-com

Committed by

web-flow

Commit Message

<a href="https://github.com/angelolab/toffy/commit/<a class=hub.com/angelolab/toffy/commit/7b5626f3e8a9ce5b59c2891b8b7a8e612cac320d">7b5626f3e<a href="https://github.com/angelolab/toffy/commit/7b5626f3e8a9ce5b59c2891b8b7a8e612cac320d">">Merge </a><a class="double-link" href="https://github.com/angelolab/toffy/commit/<a class="double-link" href="https://github.com/angelolab/toffy/commit/c580e70600663cbb239a44fdcf19e85895444ff0">c580e7060</a>">c580e7060</a><a href="https://github.com/angelolab/toffy/commit/7b5626f3e8a9ce5b59c2891b8b7a8e612cac320d"> into 0cebac9e6">0cebac9e6</a>

Pull Request Pull Request #50: Add watcher compatible QC notebook

Run Details

47 of 47 new or added lines in 1 file covered. (100.0%)

832 of 950 relevant lines covered (87.58%)

0.88 hits per line

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

95.59

/toffy/qc_comp.py

import copy
import matplotlib.pyplot as plt
import os
import numpy as np
import pandas as pd
import re
from requests.exceptions import HTTPError
from scipy.ndimage import gaussian_filter
import seaborn as sns
from shutil import rmtree
from skimage.io import imsave

from toffy.mibitracker_utils import MibiTrackerError
from toffy.mibitracker_utils import MibiRequests

import ark.settings as settings
import ark.utils.io_utils as io_utils
import ark.utils.load_utils as load_utils
import ark.utils.misc_utils as misc_utils
import mibi_bin_tools.bin_files as bin_files

# needed to prevent UserWarning: low contrast image barf when saving images
import warnings
warnings.filterwarnings('ignore')


def create_mibitracker_request_helper(email, password):
    """Create a mibitracker request helper to access a user's MIBItracker info on Ionpath

    Args:
        email (str):
            The user's MIBItracker email address
        password (str):
            The user's MIBItracker password

    Returns:
        ark.mibi.mibitracker_utils.MibiRequests:
            A request helper module instance to access a user's MIBItracker info
    """

    try:
        return MibiRequests(settings.MIBITRACKER_BACKEND, email, password)
    except HTTPError:
        print("Invalid MIBItracker email or password provided")


def download_mibitracker_data(email, password, run_name, run_label, base_dir, tiff_dir,
                              overwrite_tiff_dir=False, img_sub_folder=None,
                              fovs=None, channels=None):
    """Download a specific run's image data off of MIBITracker
    in an `ark` compatible directory structure

    Args:
        email (str):
            The user's MIBItracker email address
        password (str):
            The user's MIBItracker password
        run_name (str):
            The name of the run (specified on the user's MIBItracker run page)
        run_label (str):
            The label of the run (specified on the user's MIBItracker run page)
        base_dir (str):
            Where to place the created `tiff_dir`
        overwrite_tiff_dir (bool):
            Whether to overwrite the data already in `tiff_dir`
        tiff_dir (str):
            The name of the data directory in `base_dir` to write the run's image data to
        img_sub_folder (str):
            If specified, the subdirectory inside each FOV folder in `data_dir` to place
            the image data into
        fovs (list):
            A list of FOVs to subset over. If `None`, uses all FOVs.
        channels (lsit):
            A list of channels to subset over. If `None`, uses all channels.

    Returns:
        list:
            A list of tuples containing (point name, point id), sorted by point id.
            This defines the run acquisition order needed for plotting the QC graphs.
    """

    # verify that base_dir provided exists
    if not os.path.exists(base_dir):
        raise FileNotFoundError("base_dir %s does not exist" % base_dir)

    # create the MIBItracker request helper
    mr = create_mibitracker_request_helper(email, password)

    # get the run info using the run_name and the run_label
    # NOTE: there will only be one entry in the 'results' key with run_name and run_label specified
    run_info = mr.search_runs(run_name, run_label)

    # if no results are returned, invalid run_name and/or run_label provided
    if len(run_info['results']) == 0:
        raise ValueError('No data found for run_name %s and run_label %s' % (run_name, run_label))

    # extract the name of the FOVs and their associated internal IDs
    run_fov_names = [img['number'] for img in run_info['results'][0]['imageset']['images']]
    run_fov_ids = [img['id'] for img in run_info['results'][0]['imageset']['images']]

    # if fovs is None, ensure all of the fovs in run_fov_names_ids are chosen
    if fovs is None:
        fovs = run_fov_names

    # ensure all of the fovs are valid (important if the user explicitly specifies fovs)
    misc_utils.verify_in_list(
        provided_fovs=fovs,
        mibitracker_run_fovs=run_fov_names
    )

    # iterate through each fov and get the longest set of channels
    # this is to prevent getting no channels if the first FOV is a Moly point or is incomplete
    image_data = run_info['results'][0]['imageset']['images']
    _, run_channels = max(
        {i: image_data[i]['pngs'] for i in np.arange(len(image_data))}.items(),
        key=lambda x: len(set(x[1]))
    )

    # if channels is None, ensure all of the channels in run_channels are chosen
    if channels is None:
        channels = run_channels

    # ensure all of the channels are valid (important if the user explicitly specifies channels)
    misc_utils.verify_in_list(
        provided_chans=channels,
        mibitracker_run_chans=run_channels
    )

    # if the desired tiff_dir exists, remove it if overwrite_tiff_dir is True
    # otherwise, throw an error
    if os.path.exists(os.path.join(base_dir, tiff_dir)):
        if overwrite_tiff_dir:
            print("Overwriting existing data in tiff_dir %s" % tiff_dir)
            rmtree(os.path.join(base_dir, tiff_dir))
        else:
            raise ValueError("tiff_dir %s already exists in %s" % (tiff_dir, base_dir))

    # make the image directory
    os.mkdir(os.path.join(base_dir, tiff_dir))

    # ensure sub_folder gets set to "" if img_sub_folder is None (for os.path.join convenience)
    if not img_sub_folder:
        img_sub_folder = ""

    # define the run order list to return
    run_order = []

    # iterate over each FOV of the run
    for img in run_info['results'][0]['imageset']['images']:
        # if the image fov name is not specified, move on
        if img['number'] not in fovs:
            continue

        print("Creating data for fov %s" % img['number'])

        # make the fov directory
        os.mkdir(os.path.join(base_dir, tiff_dir, img['number']))

        # make the img_sub_folder inside the fov directory if specified
        if len(img_sub_folder) > 0:
            os.mkdir(os.path.join(base_dir, tiff_dir, img['number'], img_sub_folder))

        # iterate over each provided channel
        for chan in channels:
            # extract the channel data from MIBItracker as a numpy array
            # fail on the whole FOV if the channel is not found (most likely a Moly point)
            try:
                chan_data = mr.get_channel_data(img['id'], chan)
            except MibiTrackerError as mte:
                print("On FOV %s, failed to download channel %s, moving on to the next FOV. "
                      "If FOV %s is a Moly point, ignore. "
                      "Otherwise, please ensure that channel %s exists on MibiTracker for FOV %s"
                      % (img['number'], chan, img['number'], chan, img['number']))

                # clean the FOV: we will not have a folder for it (in case of Moly point)
                rmtree(os.path.join(base_dir, tiff_dir, img['number']))

                # do not attempt to download any more channels
                break

            # define the name of the channel file
            chan_file = '%s.tiff' % chan

            # write the data to a .tiff file in the FOV directory structure
            imsave(
                os.path.join(base_dir, tiff_dir, img['number'], img_sub_folder, chan_file),
                chan_data
            )

        # append the run name and run id to the list
        run_order.append((img['number'], img['id']))

    return run_order


def compute_nonzero_mean_intensity(image_data):
    """Compute the nonzero mean of a specific fov/chan pair

    Args:
        image_data (numpy.ndarray):
            the image data for a specific fov/chan pair

    Returns:
        float:
            The nonzero mean intensity of the fov/chan pair (`np.nan` if channel contains all 0s)
    """

    # take just the non-zero pixels
    image_data_nonzero = image_data[image_data != 0]

    # take the mean of the non-zero pixels and assign to (fov, channel) in array
    # unless there are no non-zero pixels, in which case default to 0
    if len(image_data_nonzero) > 0:
        nonzero_mean_intensity = image_data_nonzero.mean()
    else:
        nonzero_mean_intensity = 0

    return nonzero_mean_intensity


def compute_total_intensity(image_data):
    """Compute the sum of all pixels of a specific fov/chan pair

    Args:
        image_data (numpy.ndarray):
            the image data for a specific fov/chan pair

    Returns:
        float:
            The total intensity of the fov/chan pair (`np.nan` if channel contains all 0s)
    """

    return np.sum(image_data)


def compute_99_9_intensity(image_data):
    """Compute the 99.9% pixel intensity value of a specific fov/chan pair

    Args:
        image_data (numpy.ndarray):
            the image data for a specific fov/chan pair

    Returns:
        float:
            The 99.9% pixel intensity value of a specific fov/chan pair
    """

    return np.percentile(image_data, q=99.9)


def compute_qc_metrics(bin_file_path, fov_name, panel_path,
                       gaussian_blur=False, blur_factor=1):
    """Compute the QC metric matrices for the image data provided

    Args:
        bin_file_path (str):
            the directory to the MIBI bin files for extraction,
            also where the fov-level QC metric files will be written
        fov_name (str):
            the name of the FOV to extract from `bin_file_path`, needs to correspond with JSON name
        panel_path (str):
            the path to the file defining the panel info for bin file extraction
        gaussian_blur (bool):
            whether or not to add Gaussian blurring
        blur_factor (int):
            the sigma (standard deviation) to use for Gaussian blurring
            set to 0 to use raw inputs without Gaussian blurring
            ignored if `gaussian_blur` set to `False`
    """

    # path validation checks
    if not os.path.exists(bin_file_path):
        raise FileNotFoundError("bin_file_path %s does not exist" % bin_file_path)

    if not os.path.exists(panel_path):
        raise FileNotFoundError("panel_path %s does not exist" % panel_path)

    if not os.path.exists(os.path.join(bin_file_path, fov_name + '.json')):
        raise FileNotFoundError("fov file %s.json not found in bin_file_path" % fov_name)

    # run the bin file extraction, we'll extract all FOVs
    # the image coords should be: ['fov', 'type', 'x', 'y', 'channel']
    image_data = bin_files.extract_bin_files(
        data_dir=bin_file_path,
        out_dir=None,
        include_fovs=[fov_name],
        panel=pd.read_csv(panel_path)
    )

    # there's only 1 FOV and 1 type ('pulse'), so subset on that
    image_data = image_data.loc[fov_name, 'pulse', :, :, :]

    # define the list of channels to use
    chans = image_data.channel.values

    # define numpy arrays for all the metrics to extract, more efficient indexing than pandas
    blank_arr = np.zeros(image_data.shape[2], dtype='float32')
    nonzero_mean_intensity = copy.deepcopy(blank_arr)
    total_intensity = copy.deepcopy(blank_arr)
    intensity_99_9 = copy.deepcopy(blank_arr)

    # it's faster to loop through the individual channels rather than broadcasting
    for i, chan in enumerate(chans):
        # subset on the channel, cast to float32 to prevent truncation
        image_data_np = image_data.loc[:, :, chan].values.astype(np.float32)

        # STEP 1: gaussian blur (if specified)
        if gaussian_blur:
            image_data_np = gaussian_filter(
                image_data_np, sigma=blur_factor, mode='nearest', truncate=2.0
            )

        # STEP 2: extract non-zero mean intensity
        nonzero_mean_intensity[i] = compute_nonzero_mean_intensity(image_data_np)

        # STEP 3: extract total intensity
        total_intensity[i] = compute_total_intensity(image_data_np)

        # STEP 4: take 99.9% value of the data and assign
        intensity_99_9[i] = compute_99_9_intensity(image_data_np)

    # define a pandas DataFrame for each metric
    df_nonzero_mean = pd.DataFrame(
        columns=['fov', 'channel', 'Non-zero mean intensity'], dtype=object
    )
    df_total_intensity = pd.DataFrame(
        columns=['fov', 'channel', 'Total intensity'], dtype=object
    )
    df_99_9_intensity = pd.DataFrame(
        columns=['fov', 'channel', '99.9 intensity value'], dtype=object
    )

    # append the stats
    df_nonzero_mean['Non-zero mean intensity'] = nonzero_mean_intensity
    df_total_intensity['Total intensity'] = total_intensity
    df_99_9_intensity['99.9 intensity value'] = intensity_99_9

    # append the FOV name and channel name
    df_nonzero_mean['fov'] = fov_name
    df_total_intensity['fov'] = fov_name
    df_99_9_intensity['fov'] = fov_name

    df_nonzero_mean['channel'] = chans
    df_total_intensity['channel'] = chans
    df_99_9_intensity['channel'] = chans

    # save the data to the qc_dir
    df_nonzero_mean.to_csv(
        os.path.join(bin_file_path, '%s_nonzero_mean_stats.csv') % fov_name, index=False
    )
    df_total_intensity.to_csv(
        os.path.join(bin_file_path, '%s_total_intensity_stats.csv') % fov_name, index=False
    )
    df_99_9_intensity.to_csv(
        os.path.join(bin_file_path, '%s_percentile_99_9_stats.csv') % fov_name, index=False
    )


def combine_qc_metrics(bin_file_path):
    """Aggregates the QC results of each FOV into one `.csv`

    Args:
        bin_file_path (str):
            the name of the folder containing the QC metric files
    """

    # path validation check
    if not os.path.exists(bin_file_path):
        raise FileNotFoundError('bin_file_path %s does not exist' % bin_file_path)

    # iterate over each metric substr
    metric_substrs = [
        'nonzero_mean_stats.csv', 'total_intensity_stats.csv', 'percentile_99_9_stats.csv'
    ]
    for ms in metric_substrs:
        # define an aggregated metric DataFrame
        metric_df = pd.DataFrame()

        # list all the files corresponding to this metric
        metric_files = io_utils.list_files(bin_file_path, substrs=ms)

        # this prevents an already-existing combined .csv file from sneaking in
        metric_files = [
            mf_match[0] + '_%s' % ms for mf in metric_files if
            len(mf_match := re.findall(r'fov-\d+-scan-\d+', mf)) == 1
        ]  # go walrus operators

        # sort the files to ensure consistency
        # NOTE: all FOVs are named fov-m-scan-n, 'm' determines run acquisition order
        # TODO: need an additional sort by scan?
        metric_files = sorted(metric_files, key=lambda x: re.findall(r'\d+', x)[0])

        # iterate over each metric file and append the data to metric_df
        for mf in metric_files:
            metric_df = pd.concat([metric_df, pd.read_csv(os.path.join(bin_file_path, mf))])

        # write the aggregated metric data
        # NOTE: if this combined metric file already exists, it will be overwritten
        metric_df.to_csv(os.path.join(bin_file_path, 'combined_%s' % ms), index=False)


def visualize_qc_metrics(qc_metric_df, metric_name, axes_size=16, wrap=6, dpi=None, save_dir=None):
    """Visualize a barplot of a specific QC metric

    Args:
        qc_metric_df (pandas.DataFrame):
            A QC metric matrix as returned by `compute_qc_metrics`, melted
        metric_name (str):
            The name of the QC metric, used as the y-axis label
        axes_size (int):
            The font size of the axes labels
        wrap (int):
            How many plots to display per row
        dpi (int):
            If saving, the resolution of the image to use
            Ignored if save_dir is None
        save_dir (str):
            If saving, the name of the directory to save visualization to
    """

    # catplot allows for easy facets on a barplot
    g = sns.catplot(
        x='fov',
        y=metric_name,
        col='channel',
        col_wrap=wrap,
        data=qc_metric_df,
        kind='bar',
        color='black',
        sharex=True,
        sharey=False
    )

    # per Erin's visualization, don't show the hundreds of fov labels on the x-axis
    _ = plt.xticks([])

    # remove the 'channel =' in each subplot title
    _ = g.set_titles(template='{col_name}')

    # per Erin's visualization remove the default axis title on the y-axis
    # and instead show 'fov' along x-axis and the metric name along the y-axis (overarching)
    _ = g.set_axis_labels('', '')
    _ = g.fig.text(
        x=0.5,
        y=0,
        horizontalalignment='center',
        s='fov',
        size=axes_size
    )
    _ = g.fig.text(
        x=0,
        y=0.5,
        verticalalignment='center',
        s=metric_name,
        size=axes_size,
        rotation=90
    )

    # save the figure if specified
    if save_dir is not None:
        misc_utils.save_figure(save_dir, '%s_barplot_stats.png' % metric_name, dpi=dpi)

1	import copy	1✔
2	import matplotlib.pyplot as plt	1✔
3	import os	1✔
4	import numpy as np	1✔
5	import pandas as pd	1✔
6	import re	1✔
7	from requests.exceptions import HTTPError	1✔
8	from scipy.ndimage import gaussian_filter	1✔
9	import seaborn as sns	1✔
10	from shutil import rmtree	1✔
11	from skimage.io import imsave	1✔
12
13	from toffy.mibitracker_utils import MibiTrackerError	1✔
14	from toffy.mibitracker_utils import MibiRequests	1✔
15
16	import ark.settings as settings	1✔
17	import ark.utils.io_utils as io_utils	1✔
18	import ark.utils.load_utils as load_utils	1✔
19	import ark.utils.misc_utils as misc_utils	1✔
20	import mibi_bin_tools.bin_files as bin_files	1✔
21
22	# needed to prevent UserWarning: low contrast image barf when saving images
23	import warnings	1✔
24	warnings.filterwarnings('ignore')	1✔
25
26
27	def create_mibitracker_request_helper(email, password):	1✔
28	"""Create a mibitracker request helper to access a user's MIBItracker info on Ionpath
29
30	Args:
31	email (str):
32	The user's MIBItracker email address
33	password (str):
34	The user's MIBItracker password
35
36	Returns:
37	ark.mibi.mibitracker_utils.MibiRequests:
38	A request helper module instance to access a user's MIBItracker info
39	"""
40
41	try:	1✔
42	return MibiRequests(settings.MIBITRACKER_BACKEND, email, password)	1✔
43	except HTTPError:	1✔
44	print("Invalid MIBItracker email or password provided")	1✔
45
46
47	def download_mibitracker_data(email, password, run_name, run_label, base_dir, tiff_dir,	1✔
48	overwrite_tiff_dir=False, img_sub_folder=None,
49	fovs=None, channels=None):
50	"""Download a specific run's image data off of MIBITracker
51	in an `ark` compatible directory structure
52
53	Args:
54	email (str):
55	The user's MIBItracker email address
56	password (str):
57	The user's MIBItracker password
58	run_name (str):
59	The name of the run (specified on the user's MIBItracker run page)
60	run_label (str):
61	The label of the run (specified on the user's MIBItracker run page)
62	base_dir (str):
63	Where to place the created `tiff_dir`
64	overwrite_tiff_dir (bool):
65	Whether to overwrite the data already in `tiff_dir`
66	tiff_dir (str):
67	The name of the data directory in `base_dir` to write the run's image data to
68	img_sub_folder (str):
69	If specified, the subdirectory inside each FOV folder in `data_dir` to place
70	the image data into
71	fovs (list):
72	A list of FOVs to subset over. If `None`, uses all FOVs.
73	channels (lsit):
74	A list of channels to subset over. If `None`, uses all channels.
75
76	Returns:
77	list:
78	A list of tuples containing (point name, point id), sorted by point id.
79	This defines the run acquisition order needed for plotting the QC graphs.
80	"""
81
82	# verify that base_dir provided exists
83	if not os.path.exists(base_dir):	1✔
84	raise FileNotFoundError("base_dir %s does not exist" % base_dir)	1✔
85
86	# create the MIBItracker request helper
87	mr = create_mibitracker_request_helper(email, password)	1✔
88
89	# get the run info using the run_name and the run_label
90	# NOTE: there will only be one entry in the 'results' key with run_name and run_label specified
91	run_info = mr.search_runs(run_name, run_label)	1✔
92
93	# if no results are returned, invalid run_name and/or run_label provided
94	if len(run_info['results']) == 0:	1✔
95	raise ValueError('No data found for run_name %s and run_label %s' % (run_name, run_label))	1✔
96
97	# extract the name of the FOVs and their associated internal IDs
98	run_fov_names = [img['number'] for img in run_info['results'][0]['imageset']['images']]	1✔
99	run_fov_ids = [img['id'] for img in run_info['results'][0]['imageset']['images']]	1✔
100
101	# if fovs is None, ensure all of the fovs in run_fov_names_ids are chosen
102	if fovs is None:	1✔
103	fovs = run_fov_names	1✔
104
105	# ensure all of the fovs are valid (important if the user explicitly specifies fovs)
106	misc_utils.verify_in_list(	1✔
107	provided_fovs=fovs,
108	mibitracker_run_fovs=run_fov_names
109	)
110
111	# iterate through each fov and get the longest set of channels
112	# this is to prevent getting no channels if the first FOV is a Moly point or is incomplete
113	image_data = run_info['results'][0]['imageset']['images']	1✔
114	_, run_channels = max(	1✔
115	{i: image_data[i]['pngs'] for i in np.arange(len(image_data))}.items(),
116	key=lambda x: len(set(x[1]))
117	)
118
119	# if channels is None, ensure all of the channels in run_channels are chosen
120	if channels is None:	1✔
121	channels = run_channels	1✔
122
123	# ensure all of the channels are valid (important if the user explicitly specifies channels)
124	misc_utils.verify_in_list(	1✔
125	provided_chans=channels,
126	mibitracker_run_chans=run_channels
127	)
128
129	# if the desired tiff_dir exists, remove it if overwrite_tiff_dir is True
130	# otherwise, throw an error
131	if os.path.exists(os.path.join(base_dir, tiff_dir)):	1✔
132	if overwrite_tiff_dir:	1✔
133	print("Overwriting existing data in tiff_dir %s" % tiff_dir)	1✔
134	rmtree(os.path.join(base_dir, tiff_dir))	1✔
135	else:
136	raise ValueError("tiff_dir %s already exists in %s" % (tiff_dir, base_dir))	1✔
137
138	# make the image directory
139	os.mkdir(os.path.join(base_dir, tiff_dir))	1✔
140
141	# ensure sub_folder gets set to "" if img_sub_folder is None (for os.path.join convenience)
142	if not img_sub_folder:	1✔
143	img_sub_folder = ""	1✔
144
145	# define the run order list to return
146	run_order = []	1✔
147
148	# iterate over each FOV of the run
149	for img in run_info['results'][0]['imageset']['images']:	1✔
150	# if the image fov name is not specified, move on
151	if img['number'] not in fovs:	1✔
152	continue	1✔
153
154	print("Creating data for fov %s" % img['number'])	1✔
155
156	# make the fov directory
157	os.mkdir(os.path.join(base_dir, tiff_dir, img['number']))	1✔
158
159	# make the img_sub_folder inside the fov directory if specified
160	if len(img_sub_folder) > 0:	1✔
161	os.mkdir(os.path.join(base_dir, tiff_dir, img['number'], img_sub_folder))	1✔
162
163	# iterate over each provided channel
164	for chan in channels:	1✔
165	# extract the channel data from MIBItracker as a numpy array
166	# fail on the whole FOV if the channel is not found (most likely a Moly point)
167	try:	1✔
168	chan_data = mr.get_channel_data(img['id'], chan)	1✔
169	except MibiTrackerError as mte:	×
170	print("On FOV %s, failed to download channel %s, moving on to the next FOV. "	×
171	"If FOV %s is a Moly point, ignore. "
172	"Otherwise, please ensure that channel %s exists on MibiTracker for FOV %s"
173	% (img['number'], chan, img['number'], chan, img['number']))
174
175	# clean the FOV: we will not have a folder for it (in case of Moly point)
176	rmtree(os.path.join(base_dir, tiff_dir, img['number']))	×
177
178	# do not attempt to download any more channels
179	break	×
180
181	# define the name of the channel file
182	chan_file = '%s.tiff' % chan	1✔
183
184	# write the data to a .tiff file in the FOV directory structure
185	imsave(	1✔
186	os.path.join(base_dir, tiff_dir, img['number'], img_sub_folder, chan_file),
187	chan_data
188	)
189
190	# append the run name and run id to the list
191	run_order.append((img['number'], img['id']))	1✔
192
193	return run_order	1✔
194
195
196	def compute_nonzero_mean_intensity(image_data):	1✔
197	"""Compute the nonzero mean of a specific fov/chan pair
198
199	Args:
200	image_data (numpy.ndarray):
201	the image data for a specific fov/chan pair
202
203	Returns:
204	float:
205	The nonzero mean intensity of the fov/chan pair (`np.nan` if channel contains all 0s)
206	"""
207
208	# take just the non-zero pixels
209	image_data_nonzero = image_data[image_data != 0]	1✔
210
211	# take the mean of the non-zero pixels and assign to (fov, channel) in array
212	# unless there are no non-zero pixels, in which case default to 0
213	if len(image_data_nonzero) > 0:	1✔
214	nonzero_mean_intensity = image_data_nonzero.mean()	1✔
215	else:
216	nonzero_mean_intensity = 0	×
217
218	return nonzero_mean_intensity	1✔
219
220
221	def compute_total_intensity(image_data):	1✔
222	"""Compute the sum of all pixels of a specific fov/chan pair
223
224	Args:
225	image_data (numpy.ndarray):
226	the image data for a specific fov/chan pair
227
228	Returns:
229	float:
230	The total intensity of the fov/chan pair (`np.nan` if channel contains all 0s)
231	"""
232
233	return np.sum(image_data)	1✔
234
235
236	def compute_99_9_intensity(image_data):	1✔
237	"""Compute the 99.9% pixel intensity value of a specific fov/chan pair
238
239	Args:
240	image_data (numpy.ndarray):
241	the image data for a specific fov/chan pair
242
243	Returns:
244	float:
245	The 99.9% pixel intensity value of a specific fov/chan pair
246	"""
247
248	return np.percentile(image_data, q=99.9)	1✔
249
250
251	def compute_qc_metrics(bin_file_path, fov_name, panel_path,	1✔
252	gaussian_blur=False, blur_factor=1):
253	"""Compute the QC metric matrices for the image data provided
254
255	Args:
256	bin_file_path (str):
257	the directory to the MIBI bin files for extraction,
258	also where the fov-level QC metric files will be written
259	fov_name (str):
260	the name of the FOV to extract from `bin_file_path`, needs to correspond with JSON name
261	panel_path (str):
262	the path to the file defining the panel info for bin file extraction
263	gaussian_blur (bool):
264	whether or not to add Gaussian blurring
265	blur_factor (int):
266	the sigma (standard deviation) to use for Gaussian blurring
267	set to 0 to use raw inputs without Gaussian blurring
268	ignored if `gaussian_blur` set to `False`
269	"""
270
271	# path validation checks
272	if not os.path.exists(bin_file_path):	1✔
273	raise FileNotFoundError("bin_file_path %s does not exist" % bin_file_path)	1✔
274
275	if not os.path.exists(panel_path):	1✔
276	raise FileNotFoundError("panel_path %s does not exist" % panel_path)	1✔
277
278	if not os.path.exists(os.path.join(bin_file_path, fov_name + '.json')):	1✔
279	raise FileNotFoundError("fov file %s.json not found in bin_file_path" % fov_name)	1✔
280
281	# run the bin file extraction, we'll extract all FOVs
282	# the image coords should be: ['fov', 'type', 'x', 'y', 'channel']
283	image_data = bin_files.extract_bin_files(	1✔
284	data_dir=bin_file_path,
285	out_dir=None,
286	include_fovs=[fov_name],
287	panel=pd.read_csv(panel_path)
288	)
289
290	# there's only 1 FOV and 1 type ('pulse'), so subset on that
291	image_data = image_data.loc[fov_name, 'pulse', :, :, :]	1✔
292
293	# define the list of channels to use
294	chans = image_data.channel.values	1✔
295
296	# define numpy arrays for all the metrics to extract, more efficient indexing than pandas
297	blank_arr = np.zeros(image_data.shape[2], dtype='float32')	1✔
298	nonzero_mean_intensity = copy.deepcopy(blank_arr)	1✔
299	total_intensity = copy.deepcopy(blank_arr)	1✔
300	intensity_99_9 = copy.deepcopy(blank_arr)	1✔
301
302	# it's faster to loop through the individual channels rather than broadcasting
303	for i, chan in enumerate(chans):	1✔
304	# subset on the channel, cast to float32 to prevent truncation
305	image_data_np = image_data.loc[:, :, chan].values.astype(np.float32)	1✔
306
307	# STEP 1: gaussian blur (if specified)
308	if gaussian_blur:	1✔
309	image_data_np = gaussian_filter(	1✔
310	image_data_np, sigma=blur_factor, mode='nearest', truncate=2.0
311	)
312
313	# STEP 2: extract non-zero mean intensity
314	nonzero_mean_intensity[i] = compute_nonzero_mean_intensity(image_data_np)	1✔
315
316	# STEP 3: extract total intensity
317	total_intensity[i] = compute_total_intensity(image_data_np)	1✔
318
319	# STEP 4: take 99.9% value of the data and assign
320	intensity_99_9[i] = compute_99_9_intensity(image_data_np)	1✔
321
322	# define a pandas DataFrame for each metric
323	df_nonzero_mean = pd.DataFrame(	1✔
324	columns=['fov', 'channel', 'Non-zero mean intensity'], dtype=object
325	)
326	df_total_intensity = pd.DataFrame(	1✔
327	columns=['fov', 'channel', 'Total intensity'], dtype=object
328	)
329	df_99_9_intensity = pd.DataFrame(	1✔
330	columns=['fov', 'channel', '99.9 intensity value'], dtype=object
331	)
332
333	# append the stats
334	df_nonzero_mean['Non-zero mean intensity'] = nonzero_mean_intensity	1✔
335	df_total_intensity['Total intensity'] = total_intensity	1✔
336	df_99_9_intensity['99.9 intensity value'] = intensity_99_9	1✔
337
338	# append the FOV name and channel name
339	df_nonzero_mean['fov'] = fov_name	1✔
340	df_total_intensity['fov'] = fov_name	1✔
341	df_99_9_intensity['fov'] = fov_name	1✔
342
343	df_nonzero_mean['channel'] = chans	1✔
344	df_total_intensity['channel'] = chans	1✔
345	df_99_9_intensity['channel'] = chans	1✔
346
347	# save the data to the qc_dir
348	df_nonzero_mean.to_csv(	1✔
349	os.path.join(bin_file_path, '%s_nonzero_mean_stats.csv') % fov_name, index=False
350	)
351	df_total_intensity.to_csv(	1✔
352	os.path.join(bin_file_path, '%s_total_intensity_stats.csv') % fov_name, index=False
353	)
354	df_99_9_intensity.to_csv(	1✔
355	os.path.join(bin_file_path, '%s_percentile_99_9_stats.csv') % fov_name, index=False
356	)
357
358
359	def combine_qc_metrics(bin_file_path):	1✔
360	"""Aggregates the QC results of each FOV into one `.csv`
361
362	Args:
363	bin_file_path (str):
364	the name of the folder containing the QC metric files
365	"""
366
367	# path validation check
368	if not os.path.exists(bin_file_path):	1✔
369	raise FileNotFoundError('bin_file_path %s does not exist' % bin_file_path)	×
370
371	# iterate over each metric substr
372	metric_substrs = [	1✔
373	'nonzero_mean_stats.csv', 'total_intensity_stats.csv', 'percentile_99_9_stats.csv'
374	]
375	for ms in metric_substrs:	1✔
376	# define an aggregated metric DataFrame
377	metric_df = pd.DataFrame()	1✔
378
379	# list all the files corresponding to this metric
380	metric_files = io_utils.list_files(bin_file_path, substrs=ms)	1✔
381
382	# this prevents an already-existing combined .csv file from sneaking in
383	metric_files = [	1✔
384	mf_match[0] + '_%s' % ms for mf in metric_files if
385	len(mf_match := re.findall(r'fov-\d+-scan-\d+', mf)) == 1
386	] # go walrus operators
387
388	# sort the files to ensure consistency
389	# NOTE: all FOVs are named fov-m-scan-n, 'm' determines run acquisition order
390	# TODO: need an additional sort by scan?
391	metric_files = sorted(metric_files, key=lambda x: re.findall(r'\d+', x)[0])	1✔
392
393	# iterate over each metric file and append the data to metric_df
394	for mf in metric_files:	1✔
395	metric_df = pd.concat([metric_df, pd.read_csv(os.path.join(bin_file_path, mf))])	1✔
396
397	# write the aggregated metric data
398	# NOTE: if this combined metric file already exists, it will be overwritten
399	metric_df.to_csv(os.path.join(bin_file_path, 'combined_%s' % ms), index=False)	1✔
400
401
402	def visualize_qc_metrics(qc_metric_df, metric_name, axes_size=16, wrap=6, dpi=None, save_dir=None):	1✔
403	"""Visualize a barplot of a specific QC metric
404
405	Args:
406	qc_metric_df (pandas.DataFrame):
407	A QC metric matrix as returned by `compute_qc_metrics`, melted
408	metric_name (str):
409	The name of the QC metric, used as the y-axis label
410	axes_size (int):
411	The font size of the axes labels
412	wrap (int):
413	How many plots to display per row
414	dpi (int):
415	If saving, the resolution of the image to use
416	Ignored if save_dir is None
417	save_dir (str):
418	If saving, the name of the directory to save visualization to
419	"""
420
421	# catplot allows for easy facets on a barplot
422	g = sns.catplot(	1✔
423	x='fov',
424	y=metric_name,
425	col='channel',
426	col_wrap=wrap,
427	data=qc_metric_df,
428	kind='bar',
429	color='black',
430	sharex=True,
431	sharey=False
432	)
433
434	# per Erin's visualization, don't show the hundreds of fov labels on the x-axis
435	_ = plt.xticks([])	1✔
436
437	# remove the 'channel =' in each subplot title
438	_ = g.set_titles(template='{col_name}')	1✔
439
440	# per Erin's visualization remove the default axis title on the y-axis
441	# and instead show 'fov' along x-axis and the metric name along the y-axis (overarching)
442	_ = g.set_axis_labels('', '')	1✔
443	_ = g.fig.text(	1✔
444	x=0.5,
445	y=0,
446	horizontalalignment='center',
447	s='fov',
448	size=axes_size
449	)
450	_ = g.fig.text(	1✔
451	x=0,
452	y=0.5,
453	verticalalignment='center',
454	s=metric_name,
455	size=axes_size,
456	rotation=90
457	)
458
459	# save the figure if specified
460	if save_dir is not None:	1✔
461	misc_utils.save_figure(save_dir, '%s_barplot_stats.png' % metric_name, dpi=dpi)	1✔

angelolab / toffy / 405

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