jimg_int.intensity
1import json 2import os 3import random 4import re 5from itertools import combinations, permutations 6 7import matplotlib.pyplot as plt 8import numpy as np 9import pandas as pd 10import pingouin as pg 11from scipy import stats 12from scipy.stats import wasserstein_distance 13from tqdm import tqdm 14 15import jimg_int.config as cfg 16 17from .utils import * 18 19random.seed(42) 20 21 22class FeatureIntensity(ImageTools): 23 r""" 24 Class for quantitative analysis of pixel intensity and size measurements 25 in 2D/3D biological images. Supports projection of 3D stacks, mask-based 26 intensity normalization, region size estimation and metadata extraction. 27 28 Parameters 29 ---------- 30 input_image : ndarray, optional 31 Input image or 3D stack for analysis. If 3D, projection must be applied. 32 33 image : ndarray, optional 34 2D projected image (internal use). 35 36 normalized_image_values : dict, optional 37 Dictionary storing normalized intensity statistics. 38 39 mask : ndarray, optional 40 Binary mask of region of interest (ROI). 41 42 background_mask : ndarray, optional 43 Binary mask used for background estimation. If not provided, `mask` is used. 44 45 typ : {"avg", "median", "std", "var", "max", "min"}, optional 46 Projection type for 3D stacks. Default is `"avg"`. 47 48 size_info : dict, optional 49 Dictionary storing ROI size measurements. 50 51 correction_factor : float, optional 52 Normalization correction factor applied to background intensity. 53 Must satisfy 0 < factor < 1. Default is 0.1. 54 55 img_type : str, optional 56 Image type metadata. 57 58 scale : float, optional 59 Pixel resolution in physical units (e.g. µm/px). Used in size calculations. 60 61 stack_selection : list of int, optional 62 List of Z-indices to remove when projecting a 3D image. 63 64 Attributes 65 ---------- 66 input_image : ndarray or None 67 Loaded input image. 68 69 image : ndarray or None 70 Projected 2D image. 71 72 mask : ndarray or None 73 Region of interest mask. 74 75 background_mask : ndarray or None 76 Background normalization mask. 77 78 scale : float or None 79 Scale value for size estimation. 80 81 normalized_image_values : dict or None 82 Dictionary containing intensity metrics. 83 84 size_info : dict or None 85 Dictionary with ROI size information. 86 87 typ : str 88 Selected projection type for 3D images. 89 90 stack_selection : list of int 91 Z-levels excluded from projection. 92 93 Notes 94 ----- 95 The intensity normalization formula applied per pixel is: 96 97 .. math:: 98 99 R_{i,j} = T_{i,j} - \\left( \\mu_B (1 + c) \\right) 100 101 where 102 * ``T_{i,j}`` – pixel intensity in ROI 103 * ``μ_B`` – mean intensity in background region 104 * ``c`` – correction factor 105 * ``R_{i,j}`` – normalized pixel intensity 106 107 Examples 108 -------- 109 Load a 3D image, mask and compute statistics: 110 111 >>> fi = FeatureIntensity() 112 >>> fi.load_image_3D("stack.tiff") 113 >>> fi.load_mask_("mask.png") 114 >>> fi.set_projection("median") 115 >>> fi.run_calculations() 116 >>> results = fi.get_results() 117 >>> results["intensity"]["norm_mean"] 118 """ 119 120 def __init__( 121 self, 122 input_image=None, 123 image=None, 124 normalized_image_values=None, 125 mask=None, 126 background_mask=None, 127 typ=None, 128 size_info=None, 129 correction_factor=None, 130 img_type=None, 131 scale=None, 132 stack_selection=None, 133 ): 134 """ 135 Initialize a FeatureIntensity analysis instance. 136 137 Parameters 138 ---------- 139 input_image : ndarray, optional 140 Input image or 3D stack used for analysis. If the image is 3D, a 141 projection will be computed depending on the `typ` parameter. 142 143 image : ndarray, optional 144 2D image buffer used internally after projection of the input image. 145 Should not be set manually. 146 147 normalized_image_values : dict, optional 148 Dictionary containing normalized intensity statistics. Usually filled 149 automatically after running `run_calculations()`. 150 151 mask : ndarray, optional 152 Binary mask of the target region of interest (ROI). Required for 153 intensity and size calculations. 154 155 background_mask : ndarray, optional 156 Binary mask specifying the background region used to compute the 157 normalization threshold. If not provided, the ROI mask is also used 158 as the background reference. 159 160 typ : {"avg", "median", "std", "var", "max", "min"}, optional 161 Projection method for 3D images. Determines how the z-stack is 162 collapsed into a 2D image. Default is `"avg"`. 163 164 size_info : dict, optional 165 Dictionary storing computed size metrics of the ROI. Populated after 166 invoking `size_calculations()`. 167 168 correction_factor : float, optional 169 Correction term used during intensity normalization. Must satisfy 170 0 < correction_factor < 1. Default is 0.1. 171 172 img_type : str, optional 173 Optional metadata about the image type (e.g., "tiff", "png"). 174 175 scale : float, optional 176 Pixel resolution in physical units (e.g., µm/px). Required for 177 real-size estimation in `size_calculations()`. 178 179 stack_selection : list of int, optional 180 Indices of z-planes to exclude during projection of a 3D stack. 181 182 Notes 183 ----- 184 Values not provided are initialized to `None`, except for `typ`, which 185 defaults to `"avg"`, and `correction_factor`, which defaults to 0.1. 186 187 The class is designed to be populated by loading functions: 188 `load_image_()`, `load_image_3D()`, `load_mask_()`, 189 and optionally `load_normalization_mask_()` and `load_JIMG_project_()`. 190 """ 191 192 self.input_image = input_image or None 193 """ Input image or 3D stack used for analysis. If the image is 3D, a 194 projection will be computed depending on the `typ` parameter.""" 195 196 self.image = image or None 197 """ 2D image buffer used internally after projection of the input image. 198 Should not be set manually.""" 199 200 self.normalized_image_values = normalized_image_values or None 201 """Dictionary containing normalized intensity statistics. Usually filled 202 automatically after running `run_calculations()`.""" 203 204 self.mask = mask or None 205 """Binary mask of the target region of interest (ROI). Required for 206 intensity and size calculations.""" 207 208 self.background_mask = background_mask or None 209 """ Binary mask specifying the background region used to compute the 210 normalization threshold. If not provided, the ROI mask is also used 211 as the background reference.""" 212 213 self.typ = typ or "avg" 214 """Projection method for 3D images. Determines how the z-stack is 215 collapsed into a 2D image. Default is `"avg"`.""" 216 217 self.size_info = size_info or None 218 """Dictionary storing computed size metrics of the ROI. Populated after 219 invoking `size_calculations()`.""" 220 221 self.correction_factor = correction_factor or 0.1 222 """ Correction term used during intensity normalization. Must satisfy 223 0 < correction_factor < 1. Default is 0.1.""" 224 225 self.scale = scale or None 226 """ Pixel resolution in physical units (e.g., µm/px). Required for 227 real-size estimation in `size_calculations()`.""" 228 229 self.stack_selection = stack_selection or [] 230 """Indices of z-planes to exclude during projection of a 3D stack.""" 231 232 @property 233 def current_metadata(self): 234 r""" 235 Return current metadata parameters used in image processing and normalization. 236 237 Returns 238 ------- 239 tuple 240 A tuple containing: 241 242 projection_type : str 243 Projection method used for 3D image reduction (e.g., "avg", "median"). 244 245 correction_factor : float 246 Correction factor used for background subtraction during intensity 247 normalization. The applied formula is: 248 249 .. math:: 250 251 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 252 253 where 254 * ``R_{i,j}`` — normalized pixel intensity 255 * ``T_{i,j}`` — original pixel intensity 256 * ``μ_B`` — mean background intensity 257 * ``c`` — correction factor 258 scale : float or None 259 Pixel resolution (unit/px), loaded via `load_JIMG_project_()` or set manually 260 using `set_scale()`. 261 262 stack_selection : list of int 263 Indices of z-slices excluded from projection of a 3D image. 264 265 Notes 266 ----- 267 This property also prints the metadata values to the console for quick inspection. 268 """ 269 270 print(f"Projection type: {self.typ}") 271 print(f"Correction factor: {self.correction_factor}") 272 print(f"Scale (unit/px): {self.scale}") 273 print(f"Selected stac to remove: {self.stack_selection}") 274 275 return self.typ, self.correction_factor, self.scale, self.stack_selection 276 277 def set_projection(self, projection: str): 278 """ 279 Set the projection method for 3D image stack reduction. 280 281 Parameters 282 ---------- 283 projection : {"avg", "median", "std", "var", "max", "min"} 284 Projection method to reduce a 3D image stack to a 2D image. Default is `"avg"`. 285 286 Notes 287 ----- 288 This method updates the `typ` attribute of the class. The selected projection 289 determines how the z-stack is collapsed: 290 - `"avg"` : average intensity across slices 291 - `"median"` : median intensity across slices 292 - `"std"` : standard deviation across slices 293 - `"var"` : variance across slices 294 - `"max"` : maximum intensity across slices 295 - `"min"` : minimum intensity across slices 296 """ 297 298 t = ["avg", "median", "std", "var", "max", "min"] 299 if projection in t: 300 self.typ = projection 301 else: 302 print(f"\nProvided parameter is incorrect. Avaiable projection types: {t}") 303 304 def set_correction_factorn(self, factor: float): 305 r""" 306 Set the correction factor for background subtraction during intensity normalization. 307 308 Parameters 309 ---------- 310 factor : float 311 Correction factor to adjust background subtraction. Must satisfy 0 < factor < 1. 312 Default is 0.1. 313 314 Notes 315 ----- 316 The correction is applied per pixel in the target mask using the formula: 317 318 .. math:: 319 320 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 321 322 where 323 * ``R_{i,j}`` — normalized pixel intensity 324 * ``T_{i,j}`` — original pixel intensity 325 * ``μ_B`` — mean intensity in the background mask 326 * ``c`` — correction factor 327 """ 328 329 if factor < 1 and factor > 0: 330 self.correction_factor = factor 331 else: 332 print( 333 "\nProvided parameter is incorrect. The factor should be a floating-point value within the range of 0 to 1." 334 ) 335 336 def set_scale(self, scale): 337 """ 338 Set the scale for converting pixel measurements to physical units. 339 340 Parameters 341 ---------- 342 scale : float 343 Pixel resolution in physical units (e.g., µm/px). Used to calculate the 344 actual size of the tissue or organ. 345 346 Notes 347 ----- 348 The scale can also be automatically loaded from a JIMG project using 349 `load_JIMG_project_()`. This value is required for size calculations in 350 `size_calculations()`. 351 """ 352 353 self.scale = scale 354 355 def set_selection_list(self, rm_list: list): 356 """ 357 Set the list of z-slices to exclude when projecting a 3D image stack. 358 359 Parameters 360 ---------- 361 rm_list : list of int 362 List of indices corresponding to z-slices that should be removed from 363 the full 3D image stack before projection. 364 365 Notes 366 ----- 367 This updates the `stack_selection` attribute, which is used by the 368 `stack_selection_()` method during projection. 369 """ 370 371 self.stack_selection = rm_list 372 373 def load_JIMG_project_(self, path): 374 """ 375 Load a JIMG project from a `.pjm` file. 376 377 Parameters 378 ---------- 379 file_path : str 380 Path to the JIMG project file. The file must have a `.pjm` extension. 381 382 Returns 383 ------- 384 project : object 385 Loaded project object containing images and metadata. 386 387 Raises 388 ------ 389 ValueError 390 If the provided file path does not point to a `.pjm` file. 391 392 Notes 393 ----- 394 The method attempts to automatically set the `scale` and `stack_selection` 395 attributes from the project metadata if available. 396 """ 397 398 path = os.path.abspath(path) 399 400 if ".pjm" in path: 401 metadata = self.load_JIMG_project(path) 402 403 try: 404 self.scale = metadata.metadata["X_resolution[um/px]"] 405 except: 406 407 try: 408 self.scale = metadata.images_dict["metadata"][0][ 409 "X_resolution[um/px]" 410 ] 411 412 except: 413 print( 414 '\nUnable to set scale on this project! Set scale using "set_scale()"' 415 ) 416 417 self.stack_selection = metadata.removal_list 418 419 else: 420 print( 421 "\nWrong path. The provided path does not point to a JIMG project (*.pjm)." 422 ) 423 424 def stack_selection_(self): 425 """ 426 Remove selected z-slices from a 3D image stack based on `stack_selection`. 427 428 Notes 429 ----- 430 Only works if `input_image` is a 3D ndarray. The slices with indices listed 431 in `stack_selection` are excluded from the stack. Updates `input_image` 432 in-place. 433 434 Prints a warning if `stack_selection` is empty. 435 """ 436 437 if len(self.input_image.shape) == 3: 438 if len(self.stack_selection) > 0: 439 self.input_image = self.input_image[ 440 [ 441 x 442 for x in range(self.input_image.shape[0]) 443 if x not in self.stack_selection 444 ] 445 ] 446 else: 447 print("\nImages to remove from the stack were not selected!") 448 449 def projection(self): 450 """ 451 Project a 3D image stack into a 2D image using the method defined by `typ`. 452 453 Notes 454 ----- 455 Updates the `image` attribute with the projected 2D result. 456 457 Supported projection types (`typ`): 458 - "avg" : mean intensity across slices 459 - "median" : median intensity across slices 460 - "std" : standard deviation across slices 461 - "var" : variance across slices 462 - "max" : maximum intensity across slices 463 - "min" : minimum intensity across slices 464 465 Raises 466 ------ 467 AttributeError 468 If `input_image` is not defined. 469 """ 470 471 if self.typ == "avg": 472 img = np.mean(self.input_image, axis=0) 473 474 elif self.typ == "std": 475 img = np.std(self.input_image, axis=0) 476 477 elif self.typ == "median": 478 img = np.median(self.input_image, axis=0) 479 480 elif self.typ == "var": 481 img = np.var(self.input_image, axis=0) 482 483 elif self.typ == "max": 484 img = np.max(self.input_image, axis=0) 485 486 elif self.typ == "min": 487 img = np.min(self.input_image, axis=0) 488 489 self.image = img 490 491 def detect_img(self): 492 """ 493 Detect whether the input image is 2D or 3D and perform appropriate preprocessing. 494 495 Notes 496 ----- 497 - For 3D images, applies `stack_selection_()` and then `projection()`. 498 - For 2D images, no projection is applied. 499 - Prints status messages indicating the type of image and applied operations. 500 501 Raises 502 ------ 503 AttributeError 504 If `input_image` is not defined. 505 """ 506 check = len(self.input_image.shape) 507 508 if check == 3: 509 print("\n3D image detected! Starting processing for 3D image...") 510 print(f"Projection - {self.typ}...") 511 512 self.stack_selection_() 513 self.projection() 514 515 elif check == 2: 516 print("\n2D image detected! Starting processing for 2D image...") 517 518 else: 519 print("\nData does not match any image type!") 520 521 def load_image_3D(self, path): 522 """ 523 Load a 3D image stack from a TIFF file. 524 525 Parameters 526 ---------- 527 path : str 528 Path to the 3D image file (*.tiff) to be loaded. 529 530 Notes 531 ----- 532 The loaded image is stored in the `input_image` attribute as a 3D ndarray. 533 """ 534 path = os.path.abspath(path) 535 536 self.input_image = self.load_3D_tiff(path) 537 538 def load_image_(self, path): 539 """ 540 Load a 2D image into the class. 541 542 Parameters 543 ---------- 544 path : str 545 Path to the image file to be loaded. 546 547 Notes 548 ----- 549 The loaded image is stored in the `input_image` attribute as a 2D ndarray. 550 """ 551 path = os.path.abspath(path) 552 553 self.input_image = self.load_image(path) 554 555 def load_mask_(self, path): 556 r""" 557 Load a binary mask into the class and optionally set it as the normalization mask. 558 559 Parameters 560 ---------- 561 path : str 562 Path to the mask image file. Supported formats include 8-bit or 16-bit images 563 with extensions such as `.png` or `.jpeg`. The mask must be binary 564 (e.g., 0/255, 0/2**16-1, 0/1). 565 566 Notes 567 ----- 568 - If `load_normalization_mask_()` is not called, this mask is also used as the 569 background mask for intensity normalization. 570 - Normalization is applied per pixel using the formula: 571 572 .. math:: 573 574 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 575 576 where 577 * ``R_{i,j}`` — normalized pixel intensity 578 * ``T_{i,j}`` — pixel intensity in the target mask 579 * ``μ_B`` — mean intensity of the background (reversed mask) 580 * ``c`` — correction factor 581 """ 582 583 path = os.path.abspath(path) 584 585 self.mask = self.load_mask(path) 586 587 print( 588 "\nThis mask was also set as the reverse background mask. If you want a different background mask for normalization, use 'load_normalization_mask()'." 589 ) 590 self.background_mask = self.load_mask(path) 591 592 def load_normalization_mask_(self, path): 593 r""" 594 Load a binary mask for normalization into the class. 595 596 Parameters 597 ---------- 598 path : str 599 Path to the mask image file. Supported formats include 8-bit or 16-bit 600 images (e.g., `.png`, `.jpeg`). The mask must be binary (0/255, 0/2**16-1, 0/1). 601 602 Notes 603 ----- 604 - The mask defines the area of interest. Normalization is applied to the inverse 605 of this area (reversed mask). 606 - Normalization formula applied per pixel: 607 608 .. math:: 609 610 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 611 612 where 613 * ``R_{i,j}`` — normalized pixel intensity 614 * ``T_{i,j}`` — pixel intensity in the target mask 615 * ``μ_B`` — mean intensity of the background (reversed mask) 616 * ``c`` — correction factor 617 """ 618 619 path = os.path.abspath(path) 620 621 self.background_mask = self.load_mask(path) 622 623 def intensity_calculations(self): 624 """ 625 Calculate normalized and raw intensity statistics from the image based on masks. 626 627 This method performs intensity calculations using the main mask (`self.mask`) 628 and the background mask (`self.background_mask`). The pixel intensities within 629 the mask of interest are normalized by subtracting a threshold derived from the 630 background region and applying a correction factor (`self.correction_factor`). 631 Negative values after normalization are clipped to zero. 632 633 The following statistics are computed for both normalized and raw values: 634 - Minimum 635 - Maximum 636 - Mean 637 - Median 638 - Standard deviation 639 - Variance 640 - List of all normalized values (only for normalized data) 641 642 Notes 643 ----- 644 - The method updates the instance attribute `self.normalized_image_values` 645 with a dictionary containing both normalized and raw statistics. 646 - Normalization formula applied for each pixel in the selected mask: 647 final_val = selected_value - (threshold + threshold * correction_factor) 648 where threshold is the mean intensity in the background mask. 649 - Negative values after normalization are set to zero. 650 """ 651 652 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 653 tmp_bmask = self.ajd_mask_size(image=self.image, mask=self.background_mask) 654 655 selected_values = self.image[tmp_mask == np.max(tmp_mask)] 656 657 threshold = np.mean(self.image[tmp_bmask == np.min(tmp_bmask)]) 658 659 # normalization 660 final_val = selected_values - (threshold + (threshold * self.correction_factor)) 661 662 final_val[final_val < 0] = 0 663 664 tmp_dict = { 665 "norm_min": np.min(final_val), 666 "norm_max": np.max(final_val), 667 "norm_mean": np.mean(final_val), 668 "norm_median": np.median(final_val), 669 "norm_std": np.std(final_val), 670 "norm_var": np.var(final_val), 671 "norm_values": final_val.tolist(), 672 "min": np.min(selected_values), 673 "max": np.max(selected_values), 674 "mean": np.mean(selected_values), 675 "median": np.median(selected_values), 676 "std": np.std(selected_values), 677 "var": np.var(selected_values), 678 } 679 680 self.normalized_image_values = tmp_dict 681 682 def size_calculations(self): 683 """ 684 Calculates the size and bounding dimensions of the masked region in the image. 685 686 This method computes the following metrics based on the current mask: 687 - Total number of pixels in the mask (`px_size`) 688 - Real-world size if a scale is provided (`size`) 689 - Maximum lengths along x and y axes (`max_length_x_axis`, `max_length_y_axis`) 690 691 If `self.scale` is defined (unit per pixel), the real-world size is calculated. 692 If not, `size` will be `None` and a warning message is printed. 693 694 Returns: 695 Updates the following attributes in the class: 696 - self.size_info (dict) containing: 697 - 'size' (float or None): real-world size of the mask 698 - 'px_size' (int): number of pixels in the masked region 699 - 'max_length_x_axis' (int): length of the bounding box along the x-axis 700 - 'max_length_y_axis' (int): length of the bounding box along the y-axis 701 702 Example: 703 analysis.size_calculations() 704 print(analysis.size_info) 705 """ 706 707 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 708 709 size_px = int(len(tmp_mask[tmp_mask > np.min(tmp_mask)])) 710 711 if self.scale is not None: 712 size = float(size_px * self.scale) 713 else: 714 size = None 715 print( 716 '\nUnable to calculate real size, scale (unit/px) not provided, use "set_scale()" or load JIMG project .pjm metadata "load_pjm()" to set scale for calculations!' 717 ) 718 719 non_zero_indices = np.where(tmp_mask == np.max(tmp_mask)) 720 721 min_y, max_y = np.min(non_zero_indices[0]), np.max(non_zero_indices[0]) 722 min_x, max_x = np.min(non_zero_indices[1]), np.max(non_zero_indices[1]) 723 724 max_length_x_axis = int(max_x - min_x + 1) 725 max_length_y_axis = int(max_y - min_y + 1) 726 727 tmp_val = { 728 "size": size, 729 "px_size": size_px, 730 "max_length_x_axis": max_length_x_axis, 731 "max_length_y_axis": max_length_y_axis, 732 } 733 734 self.size_info = tmp_val 735 736 def run_calculations(self): 737 """ 738 Run the full analysis pipeline on the loaded image using the provided masks. 739 740 Notes 741 ----- 742 - The input image must be loaded via `load_image_()` or `load_image_3D()`. 743 - The ROI mask must be loaded via `load_mask_()`. Optionally, a normalization 744 mask can be loaded via `load_normalization_mask_()`. 745 - Parameters such as projection type and correction factor can be set with 746 `set_projection()` and `set_correction_factor()`. 747 - Scale and stack selection can also influence calculations if defined. 748 - To view current parameters, use the `current_metadata` property. 749 750 Returns 751 ------- 752 None 753 The results are stored internally and can be retrieved using 754 `get_results()`. 755 """ 756 757 if self.input_image is not None: 758 759 if self.mask is not None: 760 761 print("\nStart...") 762 self.detect_img() 763 self.intensity_calculations() 764 self.size_calculations() 765 print("\nCompleted!") 766 767 def get_results(self): 768 """ 769 Return the results from the analysis performed by `run_calculations()`. 770 771 Returns 772 ------- 773 results_dict : dict or None 774 Dictionary containing intensity and size results. Structure: 775 - 'intensity' : dict with normalized and raw intensity statistics 776 - 'size' : dict with ROI size metrics 777 778 Notes 779 ----- 780 If analysis has not been run yet, prints a message and returns None. 781 """ 782 783 if self.normalized_image_values is not None and self.size_info is not None: 784 785 results = { 786 "intensity": self.normalized_image_values, 787 "size": self.size_info, 788 } 789 790 return results 791 792 else: 793 print('\nAnalysis were not conducted. Run analysis "run_calculations()"') 794 795 def save_results( 796 self, 797 path="", 798 mask_region: str = "", 799 feature_name: str = "", 800 individual_number: int = 0, 801 individual_name: str = "", 802 ): 803 """ 804 Save the analysis results to a `.int` (JSON) file. 805 806 Parameters 807 ---------- 808 path : str, optional 809 Directory path where the file will be saved. Defaults to the current working directory. 810 811 mask_region : str 812 Name or identifier of the mask region (e.g., tissue, part of tissue). 813 814 feature_name : str 815 Name of the feature being analyzed. Underscores or spaces are replaced with periods. 816 817 individual_number : int 818 Unique identifier for the individual in the analysis (e.g., 1, 2, 3). 819 820 individual_name : str 821 Name of the individual (e.g., species name, tissue, organoid). 822 823 Notes 824 ----- 825 - The method validates that all required parameters are provided and that 826 analysis results exist (`normalized_image_values` and `size_info`). 827 - Creates the directory if it does not exist. 828 - File name format: 829 '<individual_name>_<individual_number>_<mask_region>_<feature_name>.int' 830 831 Raises 832 ------ 833 FileNotFoundError 834 If the specified path cannot be created or accessed. 835 836 ValueError 837 If any of `mask_region`, `feature_name`, `individual_number`, or 838 `individual_name` are missing or invalid. 839 """ 840 841 path = os.path.abspath(path) 842 843 if ( 844 len(mask_region) > 1 845 and len(feature_name) > 1 846 and individual_number != 0 847 and len(individual_name) > 1 848 ): 849 850 if self.normalized_image_values is not None and self.size_info is not None: 851 852 results = { 853 "intensity": self.normalized_image_values, 854 "size": self.size_info, 855 } 856 857 mask_region = re.sub(r"[_\s]+", ".", mask_region) 858 feature_name = re.sub(r"[_\s]+", ".", feature_name) 859 individual_number = re.sub(r"[_\s]+", ".", str(individual_number)) 860 individual_name = re.sub(r"[_\s]+", ".", individual_name) 861 862 full_name = f"{individual_name}_{individual_number}_{mask_region}_{feature_name}" 863 864 isExist = os.path.exists(path) 865 if not isExist: 866 os.makedirs(path, exist_ok=True) 867 868 full_path = os.path.join( 869 path, re.sub("\\.json", "", full_name) + ".int" 870 ) 871 872 with open(full_path, "w") as file: 873 json.dump(results, file, indent=4) 874 875 else: 876 print( 877 '\nAnalysis were not conducted. Run analysis "run_calculations()"' 878 ) 879 880 else: 881 print( 882 "\nAny of 'mask_region', 'feature_name', 'individual_number', 'individual_name' parameters were provided wrong!" 883 ) 884 885 def concatenate_intensity_data(self, directory: str = "", name: str = ""): 886 """ 887 Concatenate intensity data from multiple `.int` files and save as CSV. 888 889 Parameters 890 ---------- 891 directory : str, optional 892 Path to the directory containing `.int` files. Defaults to the current working directory. 893 894 name : str 895 Prefix for the output CSV file names. CSV files are saved in the format 896 '<name>_<gene>_<region>.csv'. 897 898 Raises 899 ------ 900 FileNotFoundError 901 If the directory cannot be accessed or no `.int` files are found. 902 903 ValueError 904 If an `.int` file is missing expected data or has an incorrect format. 905 906 Notes 907 ----- 908 - The method groups intensity data by gene (feature) and mask region. 909 - Outputs one CSV file per unique gene-region combination, saved in the specified directory. 910 """ 911 912 directory = os.path.abspath(directory) 913 914 files_list = [f for f in os.listdir(directory) if f.endswith(".int")] 915 916 genes_set = set([re.sub("\\.int", "", x.split("_")[3]) for x in files_list]) 917 regions_set = set([re.sub("\\.int", "", x.split("_")[2]) for x in files_list]) 918 919 for g in genes_set: 920 for r in regions_set: 921 json_to_save = { 922 "individual_name": [], 923 "individual_number": [], 924 "norm_intensity": [], 925 "size": [], 926 } 927 928 for f in tqdm(files_list): 929 if g in f and r in f: 930 with open(os.path.join(directory, f), "r") as file: 931 data = json.load(file) 932 933 json_to_save["norm_intensity"] = ( 934 json_to_save["norm_intensity"] 935 + data["intensity"]["norm_values"] 936 ) 937 json_to_save["individual_name"] = json_to_save[ 938 "individual_name" 939 ] + [f.split("_")[0]] * len( 940 data["intensity"]["norm_values"] 941 ) 942 json_to_save["individual_number"] = json_to_save[ 943 "individual_number" 944 ] + [f.split("_")[1]] * len( 945 data["intensity"]["norm_values"] 946 ) 947 json_to_save["size"] = json_to_save["size"] + [ 948 data["size"]["px_size"] 949 ] * len(data["intensity"]["norm_values"]) 950 951 pd.DataFrame(json_to_save).to_csv(f"{name}_{g}_{r}.csv", index=False) 952 953 954class IntensityAnalysis: 955 """ 956 Class for performing percentile-based statistical analysis on grouped data. 957 958 This class provides methods to calculate percentiles, remove outliers, aggregate 959 data into percentile bins, perform Welch's ANOVA, Kolmogorov-Smirnov and Fisher's exact tests, 960 evaluate Wasserstein distances, calculate Fold Change, and visualize results via comparative 961 histograms. It is designed to handle both single-column and multi-column combinations 962 of values for group-based analysis. 963 964 Methods 965 ------- 966 drop_up_df(data, group_col, values_col) 967 Removes upper outliers from a DataFrame based on a grouping column. 968 969 percentiles_calculation(values, sep_perc=1) 970 Calculates percentiles and creates loopable percentile ranges. 971 972 to_percentil(values, percentiles, percentiles_loop, values_col, replication_col) 973 Aggregates statistics based on percentile ranges. 974 975 df_to_percentiles(data, group_col, values_col, replication_col, sep_perc=1, drop_outlires=True) 976 Computes percentile statistics for grouped DataFrame data. 977 978 round_to_scientific_notation(num) 979 Formats a number in scientific notation or standard format. 980 981 aov(data, testes_col, comb="*") 982 Performs Welch's ANOVA on percentile-based group data. 983 984 post_aov(data, testes_col, comb="*") 985 Performs Welch's ANOVA with pairwise t-tests. 986 987 ks_percentiles(input_hist) 988 Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups. 989 990 fisher_percentiles(input_hist) 991 Perform pairwise Fisher's exact tests on percentile data across all groups. 992 993 to_wasserstein_distance(data) 994 Calculates scaled pairwise Wasserstein distances for grouped distributions. 995 996 to_fold_change(data, tested_value) 997 Calculates the Fold Change (FC) between all directed permutations of groups. 998 999 get_stats(data, tested_value) 1000 Calculates and aggregates overall statistical metrics (ANOVA, Fisher's exact, Kolmogorov-Smirnov, FC, Wasserstein distance). 1001 1002 hist_compare_plot(data, queue=None, p_adj=True, txt_size=20) 1003 Generates comparative histograms with statistical test results and metrics. 1004 """ 1005 1006 def drop_up_df(self, data: pd.DataFrame, group_col: str, values_col: str): 1007 """ 1008 Remove upper outliers from a DataFrame based on a specified value column, grouped by a grouping column. 1009 1010 Outliers are calculated and removed separately for each group defined by `group_col`. 1011 The upper outliers are defined using the interquartile range (IQR) method: 1012 values greater than Q3 + 1.5 * IQR are considered outliers. 1013 1014 Parameters 1015 ---------- 1016 data : pd.DataFrame 1017 The input DataFrame containing the data. 1018 1019 group_col : str 1020 The name of the column used for grouping the data. 1021 1022 values_col : str 1023 The column containing the values from which upper outliers will be removed. 1024 1025 Returns 1026 ------- 1027 filtered_data : pd.DataFrame 1028 A filtered DataFrame with the upper outliers removed for each group. 1029 1030 Notes 1031 ----- 1032 - Outliers are removed separately within each group. 1033 - The original DataFrame is not modified; a new filtered DataFrame is returned. 1034 """ 1035 1036 def iqr_filter(group): 1037 q75 = np.quantile(group[values_col], 0.75) 1038 q25 = np.quantile(group[values_col], 0.25) 1039 itq = q75 - q25 1040 return group[group[values_col] <= (q75 + 1.5 * itq)] 1041 1042 filtered_data = data.groupby(group_col).apply(iqr_filter).reset_index(drop=True) 1043 1044 return filtered_data 1045 1046 def percentiles_calculation(self, values, sep_perc: int = 1): 1047 """ 1048 Calculate percentiles for a set of values and generate consecutive percentile ranges. 1049 1050 This function computes percentiles from 0 to 100 at intervals defined by `sep_perc`. 1051 It also generates a list of consecutive percentile ranges that can be used for further analysis or binning. 1052 1053 Parameters 1054 ---------- 1055 values : array-like 1056 The input data values for which the percentiles are calculated. 1057 1058 sep_perc : int, optional 1059 Separation interval between percentiles (default is 1, meaning percentiles are calculated every 1%). 1060 1061 Returns 1062 ------- 1063 percentiles : np.ndarray 1064 Array of calculated percentile values. 1065 1066 percentiles_loop : list of tuple 1067 List of consecutive percentile ranges as tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1068 1069 Notes 1070 ----- 1071 - The first percentile is set to 0 to avoid issues with zero values. 1072 - `percentiles_loop` is useful for iterating through percentile ranges when aggregating statistics. 1073 """ 1074 1075 per_vector = values.copy() 1076 1077 percentiles = np.percentile(per_vector, np.arange(0, 101, sep_perc)) 1078 percentiles[0] = 0 1079 1080 percentiles_loop = [(i, i + 1) for i in range(int(100 / sep_perc))] 1081 1082 return percentiles, percentiles_loop 1083 1084 def to_percentil( 1085 self, values, percentiles, percentiles_loop, values_col, replication_col 1086 ): 1087 """ 1088 Aggregate statistics for a set of values based on percentile ranges, including replications. 1089 1090 This function calculates summary statistics (count, proportion, mean, median, 1091 standard deviation, variance) for each percentile range defined in `percentiles_loop`. 1092 It computes these statistics both for the combined data ('mutual') and separately 1093 for each individual replication. It also calculates overall metrics per replication. 1094 1095 Parameters 1096 ---------- 1097 values : pd.DataFrame 1098 Input DataFrame containing the data to be analyzed. 1099 percentiles : np.ndarray 1100 Array of percentile values used to define the boundaries of each range. 1101 percentiles_loop : list of tuple 1102 List of consecutive percentile ranges as index tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1103 values_col : str 1104 The column name in `values` containing the numeric data to aggregate. 1105 replication_col : str 1106 The column name in `values` used to identify distinct replications or samples. 1107 1108 Returns 1109 ------- 1110 full_data : dict 1111 A nested dictionary containing the calculated statistics with the following structure: 1112 - 'percentiles' : dict 1113 - 'mutual' : dict 1114 Lists of statistics ('n', 'n_standarized', 'avg', 'median', 'std', 'var') 1115 aggregated across all replications for each percentile bin. 1116 - 'replications' : dict 1117 Keys are replication names. Values are dictionaries of statistics (same as above) 1118 calculated specifically for that replication within each bin. 1119 - 'values' : dict 1120 Lists of overall statistics ('avg', 'median', 'std', 'var', 'replication') 1121 calculated for each replication as a whole (ignoring bins). 1122 1123 Notes 1124 ----- 1125 - If a percentile range contains no elements, statistics are set to 0 and count is set to 1 to avoid empty lists. 1126 """ 1127 1128 full_data = {} 1129 per_vector = values[values_col] 1130 amount = len(per_vector) 1131 1132 data_mutual = { 1133 "n": [], 1134 "n_standarized": [], 1135 "avg": [], 1136 "median": [], 1137 "std": [], 1138 "var": [], 1139 } 1140 1141 for x in percentiles_loop: 1142 subset = per_vector[ 1143 (per_vector > percentiles[x[0]]) & (per_vector <= percentiles[x[1]]) 1144 ] 1145 n_subset = len(subset) 1146 1147 if n_subset > 0: 1148 data_mutual["n"].append(n_subset) 1149 data_mutual["n_standarized"].append(n_subset / amount) 1150 data_mutual["avg"].append(np.mean(subset)) 1151 data_mutual["median"].append(np.median(subset)) 1152 data_mutual["std"].append(np.std(subset)) 1153 data_mutual["var"].append(np.var(subset)) 1154 else: 1155 data_mutual["n"].append(0) 1156 data_mutual["n_standarized"].append(0) 1157 data_mutual["avg"].append(0) 1158 data_mutual["median"].append(0) 1159 data_mutual["std"].append(0) 1160 data_mutual["var"].append(0) 1161 1162 full_data["percentiles"] = {"mutual": data_mutual, "replications": {}} 1163 1164 unique_names = set(values[replication_col]) 1165 1166 for nam in unique_names: 1167 per_vector_rep = values[values_col][values[replication_col] == nam] 1168 1169 data_rep = { 1170 "n": [], 1171 "n_standarized": [], 1172 "avg": [], 1173 "median": [], 1174 "std": [], 1175 "var": [], 1176 } 1177 1178 for x in percentiles_loop: 1179 subset = per_vector_rep[ 1180 (per_vector_rep > percentiles[x[0]]) 1181 & (per_vector_rep <= percentiles[x[1]]) 1182 ] 1183 n_subset = len(subset) 1184 1185 if n_subset > 0: 1186 data_rep["n"].append(n_subset) 1187 data_rep["n_standarized"].append(n_subset / amount) 1188 data_rep["avg"].append(np.mean(subset)) 1189 data_rep["median"].append(np.median(subset)) 1190 data_rep["std"].append(np.std(subset)) 1191 data_rep["var"].append(np.var(subset)) 1192 else: 1193 data_rep["n"].append(0) 1194 data_rep["n_standarized"].append(0) 1195 data_rep["avg"].append(0) 1196 data_rep["median"].append(0) 1197 data_rep["std"].append(0) 1198 data_rep["var"].append(0) 1199 1200 full_data["percentiles"]["replications"][nam] = data_rep 1201 1202 unique_names = set(values[replication_col]) 1203 1204 data_rep = {"avg": [], "median": [], "std": [], "var": [], "replication": []} 1205 1206 for nam in unique_names: 1207 per_vector_rep = values[values_col][values[replication_col] == nam] 1208 1209 data_rep["avg"].append(np.mean(per_vector_rep)) 1210 data_rep["median"].append(np.median(per_vector_rep)) 1211 data_rep["std"].append(np.std(per_vector_rep)) 1212 data_rep["var"].append(np.var(per_vector_rep)) 1213 data_rep["replication"].append(nam) 1214 1215 full_data["values"] = data_rep 1216 1217 return full_data 1218 1219 def df_to_percentiles( 1220 self, 1221 data: pd.DataFrame, 1222 group_col: str = "individual_name", 1223 values_col: str = "norm_intensity", 1224 replication_col: str = "individual_number", 1225 sep_perc: int = 1, 1226 drop_outlires: bool = True, 1227 ): 1228 """ 1229 Calculate summary statistics based on percentile ranges for each group in a DataFrame. 1230 1231 This method groups the input DataFrame by `group_col`, computes global percentile ranges 1232 based on `values_col`, and then aggregates statistics for each percentile bin. The aggregation 1233 is performed both mutually for the group and individually per replication. Optionally, 1234 upper outliers can be removed before the calculations. 1235 1236 Parameters 1237 ---------- 1238 data : pd.DataFrame 1239 Input DataFrame containing the grouped data. 1240 group_col : str, optional 1241 Column name used to define groups (default is 'individual_name'). 1242 values_col : str, optional 1243 Column name containing the numeric values for percentile calculations 1244 (default is 'norm_intensity'). 1245 replication_col : str, optional 1246 Column name used to identify separate replications within the groups 1247 (default is 'individual_number'). 1248 sep_perc : int, optional 1249 Separation interval for percentiles (default is 1, meaning 1% steps). 1250 drop_outlires : bool, optional 1251 If True, removes upper outliers from the data using the IQR method before 1252 performing calculations (default is True). 1253 1254 Returns 1255 ------- 1256 full_data : dict 1257 A dictionary where each key is a unique group name (from `group_col`). 1258 The corresponding value is the nested dictionary returned by `to_percentil()`, 1259 which includes bin-wise statistics ('mutual' and 'replications') and overall 1260 metrics ('values'). 1261 1262 Notes 1263 ----- 1264 - Outlier removal uses the IQR method within each group if `drop_outlires` is True. 1265 """ 1266 1267 full_data = {} 1268 1269 if drop_outlires == True: 1270 data = self.drop_up_df( 1271 data=data, group_col=group_col, values_col=values_col 1272 ) 1273 1274 groups = set(data[group_col]) 1275 val_dat = [x for x in data[values_col] if x > 0] 1276 1277 percentiles, percentiles_loop = self.percentiles_calculation( 1278 val_dat, sep_perc=sep_perc 1279 ) 1280 1281 for g in groups: 1282 1283 print(f"Group: {g} ...") 1284 1285 tmp_values = data[data[group_col] == g] 1286 1287 per_dat = self.to_percentil( 1288 tmp_values, percentiles, percentiles_loop, values_col, replication_col 1289 ) 1290 1291 full_data[g] = per_dat 1292 1293 return full_data 1294 1295 def round_to_scientific_notation(self, num): 1296 """ 1297 Round a number to scientific notation if very small, otherwise to one decimal place. 1298 1299 Parameters 1300 ---------- 1301 num : float 1302 The number to round. 1303 1304 Returns 1305 ------- 1306 str 1307 The rounded number as a string. 1308 - If `num` is 0, returns "0.0". 1309 - If `abs(num) < 1e-4`, returns scientific notation with 1 decimal and 1-digit exponent. 1310 - Otherwise, returns the number rounded to one decimal place. 1311 """ 1312 1313 if num == 0: 1314 return "0.0" 1315 1316 if abs(num) < 0.0001: 1317 rounded_num = np.format_float_scientific(num, precision=1, exp_digits=1) 1318 return rounded_num 1319 else: 1320 return f"{num:.1f}" 1321 1322 def aov(self, data, testes_col, comb: str = "*"): 1323 """ 1324 Perform a Welch's ANOVA analysis. 1325 1326 This function calculates group values by aggregating specified columns (testes_col) 1327 via the comb method and then conducts a Welch's ANOVA. This approach is ideal for 1328 comparing group means when data exhibits unequal variances across groups. 1329 1330 Parameters 1331 ---------- 1332 data : dict of pd.DataFrame 1333 Dictionary where keys are group names and values are DataFrames containing the data. 1334 1335 testes_col : str or list of str 1336 Column name(s) from which the group values are derived. If a list is provided, columns 1337 will be combined based on the `comb` operation. 1338 1339 comb : str, optional 1340 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1341 '*' : multiplication 1342 '+' : addition 1343 '**': exponentiation 1344 '-' : subtraction 1345 '/' : division 1346 Default is '*'. 1347 1348 Returns 1349 ------- 1350 F : float 1351 F-statistic from Welch's ANOVA. 1352 1353 p_val : float 1354 Uncorrected p-value from Welch's ANOVA, testing for significant differences between groups. 1355 1356 Notes 1357 ----- 1358 - If `testes_col` is a single string, no combination is performed, and the group values 1359 are taken directly from that column. 1360 - Welch's ANOVA is used as it accounts for unequal variances between groups. 1361 - The `df.melt()` method is used to reshape the data, allowing the ANOVA to be applied to all groups. 1362 1363 Examples 1364 -------- 1365 >>> welch_F, welch_p = self.aov(data, testes_col=['col1', 'col2'], comb='+') 1366 >>> print(f"Welch's ANOVA F-statistic: {welch_F}, p-value: {welch_p}") 1367 """ 1368 1369 groups = [] 1370 1371 for d in data.keys(): 1372 1373 if isinstance(testes_col, str): 1374 g = data[d]["values"][testes_col] 1375 elif isinstance(testes_col, list): 1376 g = [1] * len(data[d]["values"][testes_col[0]]) 1377 for t in testes_col: 1378 if comb == "*": 1379 g = [a * b for a, b in zip(g, data[d]["values"][t])] 1380 elif comb == "+": 1381 g = [a + b for a, b in zip(g, data[d]["values"][t])] 1382 elif comb == "**": 1383 g = [a**b for a, b in zip(g, data[d]["values"][t])] 1384 elif comb == "-": 1385 g = [a - b for a, b in zip(g, data[d]["values"][t])] 1386 elif comb == "/": 1387 g = [a / b for a, b in zip(g, data[d]["values"][t])] 1388 1389 groups.append(g) 1390 1391 df = pd.DataFrame({f"group_{i}": group for i, group in enumerate(groups)}) 1392 1393 df_melted = df.melt(var_name="group", value_name="value") 1394 1395 welch_results = pg.welch_anova(data=df_melted, dv="value", between="group") 1396 1397 return welch_results["F"].values[0], welch_results["p-unc"].values[0] 1398 1399 def post_aov(self, data, testes_col, comb: str = "*"): 1400 """ 1401 Perform Welch's ANOVA and pairwise Welch's t-tests on grouped data. 1402 1403 This method first conducts a Welch's ANOVA to detect significant differences 1404 in group means. It then performs pairwise Welch's t-tests across all group 1405 combinations to identify specific differences. All p-values are adjusted using 1406 the Bonferroni correction to account for multiple comparisons. 1407 1408 Parameters 1409 ---------- 1410 data : dict of pd.DataFrame 1411 Dictionary where keys are group names and values are DataFrames containing the data. 1412 1413 testes_col : str or list of str 1414 Column name(s) from which the group values are derived. If a list is provided, 1415 columns will be combined according to the `comb` operation. 1416 1417 comb : str, optional 1418 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1419 '*' : multiplication 1420 '+' : addition 1421 '**': exponentiation 1422 '-' : subtraction 1423 '/' : division 1424 Default is '*'. 1425 1426 Returns 1427 ------- 1428 p_val : float 1429 Uncorrected p-value from the Welch's ANOVA. 1430 1431 final_results : dict 1432 Dictionary containing results of pairwise Welch's t-tests with keys: 1433 'group1' : list of first group names in each comparison 1434 'group2' : list of second group names in each comparison 1435 'stat' : list of t-statistics for each comparison 1436 'p_val' : list of uncorrected p-values for each comparison 1437 'adj_p_val' : list of Bonferroni-adjusted p-values for multiple comparisons 1438 """ 1439 1440 p_val = self.aov(data=data, testes_col=testes_col, comb=comb)[1] 1441 1442 pairs = list(combinations(data, 2)) 1443 final_results = { 1444 "group1": [], 1445 "group2": [], 1446 "stat": [], 1447 "p_val": [], 1448 "adj_p_val": [], 1449 } 1450 1451 for group1, group2 in pairs: 1452 if isinstance(testes_col, str): 1453 g1 = data[group1]["values"][testes_col] 1454 elif isinstance(testes_col, list): 1455 g1 = [1] * len(data[group1]["values"][testes_col[0]]) 1456 for t in testes_col: 1457 if comb == "*": 1458 g1 = [a * b for a, b in zip(g1, data[group1]["values"][t])] 1459 elif comb == "+": 1460 g1 = [a + b for a, b in zip(g1, data[group1]["values"][t])] 1461 elif comb == "**": 1462 g1 = [a**b for a, b in zip(g1, data[group1]["values"][t])] 1463 elif comb == "-": 1464 g1 = [a - b for a, b in zip(g1, data[group1]["values"][t])] 1465 elif comb == "/": 1466 g1 = [a / b for a, b in zip(g1, data[group1]["values"][t])] 1467 1468 if isinstance(testes_col, str): 1469 g2 = data[group2]["values"][testes_col] 1470 elif isinstance(testes_col, list): 1471 g2 = [1] * len(data[group2]["values"][testes_col[0]]) 1472 for t in testes_col: 1473 if comb == "*": 1474 g2 = [a * b for a, b in zip(g2, data[group2]["values"][t])] 1475 elif comb == "+": 1476 g2 = [a + b for a, b in zip(g2, data[group2]["values"][t])] 1477 elif comb == "**": 1478 g2 = [a**b for a, b in zip(g2, data[group2]["values"][t])] 1479 elif comb == "-": 1480 g2 = [a - b for a, b in zip(g2, data[group2]["values"][t])] 1481 elif comb == "/": 1482 g2 = [a / b for a, b in zip(g2, data[group2]["values"][t])] 1483 1484 stat, p_val = stats.ttest_ind( 1485 g1, g2, alternative="two-sided", equal_var=False 1486 ) 1487 g = sorted([group1, group2]) 1488 final_results["group1"].append(g[0]) 1489 final_results["group2"].append(g[1]) 1490 final_results["stat"].append(stat) 1491 final_results["p_val"].append(p_val) 1492 adj = p_val * len(pairs) 1493 if adj > 1: 1494 final_results["adj_p_val"].append(1) 1495 else: 1496 final_results["adj_p_val"].append(adj) 1497 1498 return p_val, final_results 1499 1500 def ks_percentiles(self, input_hist): 1501 """ 1502 Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups. 1503 1504 This method extracts the percentile levels and computes the average value for 1505 each percentile to obtain a lower-dimensional representation of the data, thereby 1506 reducing the Big Data scale problem for each group. Using these metrics, it reconstructs 1507 the underlying empirical distributions to evaluate both structural proportions and scale. 1508 1509 To further mitigate the large sample size problem ("curse of Big Data") where inflating 1510 pixel counts yields artificially significant results, a controlled downsampling (resampling) 1511 is applied to standardize the sample sizes across groups. 1512 1513 A two-sample Kolmogorov-Smirnov (KS) test is then performed for every possible pair 1514 of groups to detect differences in distribution shapes. Finally, p-values are adjusted 1515 using the Bonferroni correction method to account for multiple comparisons and control 1516 the family-wise error rate. 1517 1518 Parameters 1519 ---------- 1520 input_hist : dict 1521 A nested dictionary where keys are group names. Each group must contain 1522 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1523 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1524 1525 Returns 1526 ------- 1527 final_results : dict 1528 A dictionary containing the results of the pairwise comparisons with keys: 1529 - 'group1': list of the first group names in the pairs. 1530 - 'group2': list of the second group names in the pairs. 1531 - 'K-S': list of Kolmogorov-Smirnov test statistics. 1532 - 'p_val': list of unadjusted p-values. 1533 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1534 1535 Example 1536 ------- 1537 >>> results = self.ks_percentiles(input_hist) 1538 """ 1539 1540 ks_data = {} 1541 1542 for d in input_hist.keys(): 1543 tmp_dic = {} 1544 1545 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["avg"]): 1546 tmp_dic[f"p{n+1}"] = c 1547 1548 ks_data[d] = tmp_dic 1549 1550 df_cleaned = pd.DataFrame(ks_data).T 1551 1552 pairs = list(combinations(df_cleaned.index, 2)) 1553 1554 final_results = { 1555 "group1": [], 1556 "group2": [], 1557 "K-S": [], 1558 "p_val": [], 1559 "adj_p_val": [], 1560 } 1561 1562 for group1, group2 in pairs: 1563 1564 g = sorted([group1, group2]) 1565 1566 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1567 1568 res = stats.ks_2samp(table_pair.iloc[:, 0], table_pair.iloc[:, 1]) 1569 1570 final_results["group1"].append(g[0]) 1571 final_results["group2"].append(g[1]) 1572 final_results["K-S"].append(res.statistic) 1573 final_results["p_val"].append(res.pvalue) 1574 adj = res.pvalue * len(pairs) 1575 if adj > 1: 1576 final_results["adj_p_val"].append(1) 1577 else: 1578 final_results["adj_p_val"].append(adj) 1579 1580 return final_results 1581 1582 def fisher_percentiles(self, input_hist): 1583 """ 1584 Perform pairwise Fisher's exact tests on percentile data across all groups. 1585 1586 This method extracts the raw counts (N) for each percentile bin across all 1587 groups to construct a contingency table representation of the data. By utilizing 1588 the discrete frequency counts per bin rather than continuous average values, it 1589 evaluates both structural distribution proportions and sample size scaling 1590 differences simultaneously. 1591 1592 An exact testing approach is applied to every unique pair of groups by extracting 1593 their corresponding sub-tables. For each pair, a Fisher's exact test (or its 1594 extension for larger contingency tables) is performed to detect statistically 1595 significant deviations in distribution profiles. 1596 1597 Finally, p-values are manually adjusted using the Bonferroni correction method 1598 by multiplying the raw p-values by the total number of comparisons to control 1599 the family-wise error rate across multiple pair-wise tests. 1600 the family-wise error rate. 1601 1602 Parameters 1603 ---------- 1604 input_hist : dict 1605 A nested dictionary where keys are group names. Each group must contain 1606 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1607 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1608 1609 Returns 1610 ------- 1611 final_results : dict 1612 A dictionary containing the results of the pairwise comparisons with keys: 1613 - 'group1': list of the first group names in the pairs. 1614 - 'group2': list of the second group names in the pairs. 1615 - 'fish': list of Fisher's exact test statistics. 1616 - 'p_val': list of unadjusted p-values. 1617 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1618 1619 Example 1620 ------- 1621 >>> results = self.fisher_percentiles(input_hist) 1622 """ 1623 1624 fish_data = {} 1625 1626 for d in input_hist.keys(): 1627 tmp_dic = {} 1628 1629 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["n"]): 1630 tmp_dic[f"p{n+1}"] = c 1631 1632 fish_data[d] = tmp_dic 1633 1634 df_cleaned = pd.DataFrame(fish_data).T 1635 1636 pairs = list(combinations(df_cleaned.index, 2)) 1637 1638 final_results = { 1639 "group1": [], 1640 "group2": [], 1641 "fish": [], 1642 "p_val": [], 1643 "adj_p_val": [], 1644 } 1645 1646 for group1, group2 in pairs: 1647 1648 g = sorted([group1, group2]) 1649 1650 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1651 1652 res = stats.fisher_exact(table_pair) 1653 1654 final_results["group1"].append(g[0]) 1655 final_results["group2"].append(g[1]) 1656 final_results["fish"].append(res.statistic) 1657 final_results["p_val"].append(res.pvalue) 1658 adj = res.pvalue * len(pairs) 1659 if adj > 1: 1660 final_results["adj_p_val"].append(1) 1661 else: 1662 final_results["adj_p_val"].append(adj) 1663 1664 return final_results 1665 1666 def to_wasserstein_distance(self, data): 1667 """ 1668 Calculate scaled pairwise Wasserstein distances for grouped distributions. 1669 1670 This method computes the 1D Wasserstein distance (Earth Mover's Distance) 1671 between all possible combinations of groups in the provided dataset. 1672 Before calculating the distance, the standardized frequencies are scaled 1673 by a factor representing the average total count (sample size) of the 1674 two compared groups. 1675 1676 Parameters 1677 ---------- 1678 data : dict 1679 A nested dictionary where keys are group names. For each group, the 1680 method expects the following internal data structure: 1681 - `data[group_name]['percentiles']['mutual']['n']` : list-like 1682 Absolute counts or sample sizes for the distribution. 1683 - `data[group_name]['percentiles']['mutual']['n_standarized']` : list-like 1684 Standardized frequencies or probabilities to be compared. 1685 1686 Returns 1687 ------- 1688 final_results : dict 1689 A dictionary containing the results of the pairwise distance calculations: 1690 - 'group1' : list of str 1691 The name of the first group in the comparison. 1692 - 'group2' : list of str 1693 The name of the second group in the comparison. 1694 - 'wasserstein_distance' : list of float 1695 The computed scaled Wasserstein distance for each pair. 1696 """ 1697 1698 pairs = list(combinations(data.keys(), 2)) 1699 1700 final_results = {"group1": [], "group2": [], "wasserstein_distance": []} 1701 1702 for group1, group2 in pairs: 1703 1704 factor = ( 1705 sum(data[group1]["percentiles"]["mutual"]["n"]) 1706 + sum(data[group2]["percentiles"]["mutual"]["n"]) 1707 ) / 2 1708 1709 dist = wasserstein_distance( 1710 [ 1711 x * factor 1712 for x in data[group1]["percentiles"]["mutual"]["n_standarized"] 1713 ], 1714 [ 1715 x * factor 1716 for x in data[group2]["percentiles"]["mutual"]["n_standarized"] 1717 ], 1718 ) 1719 1720 g = sorted([group1, group2]) 1721 final_results["group1"].append(g[0]) 1722 final_results["group2"].append(g[1]) 1723 final_results["wasserstein_distance"].append(dist) 1724 1725 return final_results 1726 1727 def to_fold_change(self, data, tested_value): 1728 """ 1729 Calculate the Fold Change (FC) between all permutations of groups. 1730 1731 This method computes the ratio of the mean values of a specified feature 1732 (`tested_value`) for every directed pair of groups. Because permutations 1733 are used, the calculation is directional (i.e., both Group A / Group B 1734 and Group B / Group A are computed). 1735 1736 Parameters 1737 ---------- 1738 data : dict 1739 A nested dictionary where keys are group names. For each group, the 1740 method expects the following internal structure: 1741 - `data[group_name]['values'][tested_value]` : array-like 1742 Numeric values used to compute the mean for the group. 1743 1744 tested_value : str 1745 The specific key or column name within the 'values' dictionary 1746 indicating which feature's fold change should be calculated. 1747 1748 Returns 1749 ------- 1750 final_results : dict 1751 A dictionary containing the results of the pairwise fold change calculations: 1752 - 'group1' : list of str 1753 The name of the numerator group in the comparison. 1754 - 'group2' : list of str 1755 The name of the denominator group in the comparison. 1756 - 'FC' : list of float 1757 The calculated fold change (mean of group1 / mean of group2). 1758 """ 1759 1760 pairs = list(permutations(data.keys(), 2)) 1761 1762 final_results = {"group1": [], "group2": [], "FC": []} 1763 1764 values = [] 1765 for group1, group2 in pairs: 1766 1767 values = values + data[group1]["values"][tested_value] 1768 values = values + data[group2]["values"][tested_value] 1769 1770 values_min = min([x for x in values if x > 0]) 1771 values_min = values_min / 2 1772 1773 for group1, group2 in pairs: 1774 1775 g1 = np.mean(data[group1]["values"][tested_value]) 1776 g2 = np.mean(data[group2]["values"][tested_value]) 1777 1778 if g1 == 0: 1779 g1 = g1 + values_min 1780 1781 if g2 == 0: 1782 g2 = g2 + values_min 1783 1784 fc = g1 / g2 1785 1786 final_results["group1"].append(group1) 1787 final_results["group2"].append(group2) 1788 final_results["FC"].append(fc) 1789 1790 return final_results 1791 1792 def get_stats(self, data, tested_value): 1793 """ 1794 Calculate and aggregate statistical metrics (ANOVA, Fisher's Exact, 1795 Kolmogorov-Smirnov, Fold Change, Wasserstein distance). 1796 1797 This method computes overall statistics and pairwise comparisons for grouped data. 1798 To properly capture both structural proportions and total count variations across 1799 percentiles while avoiding the curse of Big Data, it runs two distinct tests: 1800 1. Fisher's exact test on discrete percentile counts to evaluate absolute scale 1801 and profile differences. 1802 2. Two-sample Kolmogorov-Smirnov (KS) test on reconstructed empirical 1803 distributions to evaluate discrepancies in distribution shapes. 1804 1805 Additionally, it calculates the Fold Change (FC) and evaluates Wasserstein 1806 distances. If the average number of replicates per group is at least 3, 1807 it conducts Welch's ANOVA. The input dictionary is modified in-place to 1808 include a new 'statistics' key containing all results. 1809 1810 Parameters 1811 ---------- 1812 data : dict 1813 A nested dictionary where keys are group names. Each group's dictionary 1814 must contain the structure `['values']['replication']` to verify sample sizes, 1815 along with the necessary data structures required by downstream statistical methods. 1816 1817 tested_value : str 1818 The key or column name representing the specific variable to evaluate 1819 (e.g., used for ANOVA and Fold Change calculations). 1820 1821 Returns 1822 ------- 1823 data : dict 1824 The original input dictionary, extended with a new `data['statistics']` key 1825 that houses the computed statistical results, including `percintiles_fish` 1826 and `percintiles_ks`. 1827 1828 Example 1829 ------- 1830 stats = self.get_stats( 1831 data, 1832 tested_value='n', 1833 ) 1834 """ 1835 1836 # parametric selected value 1837 sum_k = 0 1838 n = 0 1839 for k in data.keys(): 1840 if k != "statistics": 1841 n += 1 1842 sum_k += len(data[k]["values"]["replication"]) 1843 1844 sum_k = sum_k / n 1845 1846 if sum_k >= 3: 1847 pk, dfk = self.post_aov(data, testes_col=tested_value) 1848 1849 # fish 1850 fish = self.fisher_percentiles(data) 1851 1852 # K_S 1853 ks = self.ks_percentiles(data) 1854 1855 dw = self.to_wasserstein_distance(data) 1856 1857 fc = self.to_fold_change(data, tested_value) 1858 1859 data["statistics"] = {} 1860 1861 data["statistics"]["percintiles_fish"] = fish 1862 1863 data["statistics"]["percintiles_ks"] = ks 1864 1865 if sum_k >= 3: 1866 data["statistics"]["ANOVA"] = {} 1867 1868 data["statistics"]["ANOVA"]["p_value"] = pk 1869 data["statistics"]["ANOVA"]["pair-comparison"] = dfk 1870 else: 1871 import warnings 1872 1873 warnings.warn( 1874 f"Insufficient replicates for statistical analysis. " 1875 f"At least 3 replicates per group (3 vs 3) are required. " 1876 f"The average number of samples per probe in this dataset was {n}.", 1877 RuntimeWarning, 1878 ) 1879 1880 data["statistics"]["FC"] = fc 1881 1882 data["statistics"]["wasserstein_distance"] = dw 1883 1884 data["statistics"]["tested_value"] = tested_value 1885 1886 return data 1887 1888 def hist_compare_plot( 1889 self, data, queue=None, p_adj: bool = True, txt_size: int = 20 1890 ): 1891 """ 1892 Generate comparative histograms and display results of statistical tests 1893 (ANOVA / T-test [if min 3 vs. 3 comparison available], Kolmogorov-Smirnov, Fisher percentiles) 1894 and statistics (FC, Wasserstein distance). 1895 1896 1897 Parameters 1898 ---------- 1899 data : dict 1900 Dictionary where keys are group names and values are containing histogram data. 1901 Each DataFrame should include the column specified by `tested_value`. 1902 1903 queue : list of str or None 1904 Defines the order of groups to be plotted. 1905 1906 p_adj : bool, optional 1907 If True, applies Bonferroni correction for multiple comparisons (default is True). 1908 1909 txt_size : int, optional 1910 Font size for text annotations in the plot (default is 20). 1911 1912 Returns 1913 ------- 1914 fig : matplotlib.figure.Figure 1915 Matplotlib figure object containing the generated histograms and statistical test results. 1916 1917 Example 1918 ------- 1919 fig = self.hist_compare_plot( 1920 data, 1921 queue=['group1', 'group2', 'group3'], 1922 p_adj=True, 1923 txt_size=18 1924 ) 1925 plt.show() 1926 """ 1927 1928 if queue is None: 1929 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1930 1931 if sorted(queue) != [x for x in sorted(data.keys()) if x != "statistics"]: 1932 print( 1933 "\n Wrong queue provided! The queue will be sorted with default settings!" 1934 ) 1935 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1936 1937 # parametric selected value 1938 tested_value = data["statistics"]["tested_value"] 1939 1940 ############################################################################## 1941 1942 standarized_max, standarized_min, value_max, value_min = [], [], [], [] 1943 for d in queue: 1944 standarized_max.append( 1945 max(data[d]["percentiles"]["mutual"]["n_standarized"]) 1946 ) 1947 standarized_min.append( 1948 min(data[d]["percentiles"]["mutual"]["n_standarized"]) 1949 ) 1950 value_max.append(max(data[d]["percentiles"]["mutual"][tested_value])) 1951 value_min.append(min(data[d]["percentiles"]["mutual"][tested_value])) 1952 1953 num_columns = len(queue) + 1 1954 1955 fig, axs = plt.subplots( 1956 3, 1957 num_columns, 1958 figsize=(8 * num_columns, 10), 1959 gridspec_kw={"width_ratios": [1] * len(queue) + [0.5], "wspace": 0.05}, 1960 ) 1961 1962 for i, d in enumerate(queue): 1963 tmp_data = data[d]["percentiles"]["mutual"] 1964 1965 axs[0, i].bar( 1966 [str(n) for n in range(len(tmp_data["n_standarized"]))], 1967 tmp_data["n_standarized"], 1968 width=0.95, 1969 color="gold", 1970 ) 1971 1972 # line 1973 n_groups = len(data[d]["percentiles"]["replications"].keys()) 1974 colors = plt.cm.OrRd(np.linspace(0.3, 0.9, n_groups)) 1975 1976 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 1977 1978 color = colors[ix] 1979 1980 y = data[d]["percentiles"]["replications"][dn]["n_standarized"] 1981 x = np.arange(len(y)) 1982 1983 axs[0, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 1984 1985 axs[0, i].plot( 1986 x, 1987 y, 1988 color=color, 1989 linewidth=1, 1990 marker="o", 1991 ) 1992 1993 axs[0, i].set_ylim( 1994 min(standarized_min) * 0.9995, max(standarized_max) * 1.0005 1995 ) 1996 1997 if i == 0: 1998 axs[0, i].set_ylabel("Standarized\nfrequency", fontsize=txt_size) 1999 else: 2000 axs[0, i].set_yticks([]) 2001 2002 axs[0, i].set_xticks([]) 2003 axs[0, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2004 2005 axs[1, i].bar( 2006 [str(n) for n in range(len(tmp_data[tested_value]))], 2007 tmp_data[tested_value], 2008 width=0.95, 2009 color="orange", 2010 ) 2011 2012 # line 2013 2014 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2015 2016 color = colors[ix] 2017 2018 y = data[d]["percentiles"]["replications"][dn][tested_value] 2019 x = np.arange(len(y)) 2020 2021 axs[1, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2022 2023 axs[1, i].plot( 2024 x, 2025 y, 2026 color=color, 2027 linewidth=1, 2028 marker="o", 2029 ) 2030 2031 mean_value = np.mean(data[d]["values"][tested_value]) 2032 axs[1, i].axhline(y=mean_value, color="red", linestyle="--") 2033 2034 axs[1, i].set_ylim(min(value_min) * 0.9995, max(value_max) * 1.0005) 2035 2036 if i == 0: 2037 axs[1, i].set_ylabel(f"Normalized\n{tested_value}", fontsize=txt_size) 2038 else: 2039 axs[1, i].set_yticks([]) 2040 2041 axs[1, i].set_xticks([]) 2042 axs[1, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2043 2044 axs[2, i].bar( 2045 [str(n) for n in range(len(tmp_data["n_standarized"]))], 2046 [ 2047 a * b 2048 for a, b in zip(tmp_data[tested_value], tmp_data["n_standarized"]) 2049 ], 2050 width=0.95, 2051 color="goldenrod", 2052 ) 2053 2054 # line 2055 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2056 2057 color = colors[ix] 2058 2059 y = [ 2060 a * b 2061 for a, b in zip( 2062 data[d]["percentiles"]["replications"][dn][tested_value], 2063 data[d]["percentiles"]["replications"][dn]["n_standarized"], 2064 ) 2065 ] 2066 x = np.arange(len(y)) 2067 2068 axs[2, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2069 2070 axs[2, i].plot( 2071 x, 2072 y, 2073 color=color, 2074 linewidth=1, 2075 marker="o", 2076 ) 2077 2078 mean_value = np.mean( 2079 data[d]["values"][data["statistics"]["tested_value"]] 2080 ) * np.mean(tmp_data["n_standarized"]) 2081 2082 axs[2, i].axhline(y=mean_value, color="red", linestyle="--") 2083 2084 axs[2, i].set_ylim( 2085 (min(standarized_min) * min(value_min)) * 0.9995, 2086 (max(standarized_max) * max(value_max) * 1.0005), 2087 ) 2088 axs[2, i].set_xlabel(d, fontsize=txt_size) 2089 2090 if i == 0: 2091 axs[2, i].set_ylabel( 2092 f"Standarized\nnorm_{tested_value}", fontsize=txt_size 2093 ) 2094 else: 2095 axs[2, i].set_yticks([]) 2096 2097 axs[2, i].set_xticks([]) 2098 axs[2, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2099 2100 # statistics 2101 2102 # ANOVA / t-test 2103 2104 if "ANOVA" in data["statistics"].keys(): 2105 pk = data["statistics"]["ANOVA"]["p_value"] 2106 dfk = data["statistics"]["ANOVA"]["pair-comparison"] 2107 dfk = pd.DataFrame(dfk) 2108 2109 dfk = dfk.sort_values( 2110 by=["group1", "group2"], 2111 key=lambda col: [ 2112 queue.index(val) if val in queue else -1 for val in col 2113 ], 2114 ).reset_index(drop=True) 2115 2116 sign = "ns" 2117 if float(self.round_to_scientific_notation(pk)) < 0.001: 2118 sign = "***" 2119 elif float(self.round_to_scientific_notation(pk)) < 0.01: 2120 sign = "**" 2121 elif float(self.round_to_scientific_notation(pk)) < 0.05: 2122 sign = "*" 2123 2124 text = f"Test Welch's ANOVA\non '{tested_value}' values\np-value: {self.round_to_scientific_notation(pk)} - {sign}\n" 2125 2126 if p_adj == True: 2127 for i in range(len(dfk["group1"])): 2128 sign = "ns" 2129 if dfk["adj_p_val"][i] < 0.001: 2130 sign = "***" 2131 elif dfk["adj_p_val"][i] < 0.01: 2132 sign = "**" 2133 elif dfk["adj_p_val"][i] < 0.05: 2134 sign = "*" 2135 2136 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['adj_p_val'][i])} - {sign}\n" 2137 else: 2138 for i in range(len(dfk["group1"])): 2139 sign = "ns" 2140 if dfk["p_val"][i] < 0.001: 2141 sign = "***" 2142 elif dfk["p_val"][i] < 0.01: 2143 sign = "**" 2144 elif dfk["p_val"][i] < 0.05: 2145 sign = "*" 2146 2147 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['p_val'][i])} - {sign}\n" 2148 2149 axs[2, -1].text( 2150 0.5, 2151 0.5, 2152 text, 2153 ha="center", 2154 va="center", 2155 fontsize=txt_size * 0.7, 2156 wrap=True, 2157 ) 2158 axs[2, -1].set_axis_off() 2159 else: 2160 axs[2, -1].set_axis_off() 2161 2162 # FC / Distance 2163 2164 ranking_FC = pd.DataFrame(data["statistics"]["FC"]) 2165 2166 ranking_dw = pd.DataFrame(data["statistics"]["wasserstein_distance"]) 2167 2168 ranking_combined = pd.merge( 2169 ranking_FC, ranking_dw, on=["group1", "group2"], how="right" 2170 ) 2171 2172 ranking_combined = ranking_combined.sort_values( 2173 by=["group1", "group2"], 2174 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2175 ).reset_index(drop=True) 2176 2177 text = "FC / Wasserstein distance\n" 2178 for i in range(len(ranking_combined)): 2179 group1 = ranking_combined["group1"][i] 2180 group2 = ranking_combined["group2"][i] 2181 fc_val = ranking_combined["FC"][i] 2182 wasserstein_val = ranking_combined["wasserstein_distance"][i] 2183 2184 text += f"{group1} vs. {group2}\n {fc_val:.2f} | {wasserstein_val:.2f}\n" 2185 2186 axs[1, -1].text( 2187 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2188 ) 2189 axs[1, -1].set_axis_off() 2190 2191 # fish 2192 2193 fish = pd.DataFrame(data["statistics"]["percintiles_fish"]) 2194 2195 # K-S 2196 2197 ks = pd.DataFrame(data["statistics"]["percintiles_ks"]) 2198 2199 fish = fish.sort_values( 2200 by=["group1", "group2"], 2201 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2202 ).reset_index(drop=True) 2203 2204 ks = ks.sort_values( 2205 by=["group1", "group2"], 2206 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2207 ).reset_index(drop=True) 2208 2209 text = f"Tests\nFisher's exact / Kolmogorov-Smirnov\n" 2210 2211 if p_adj == True: 2212 for i in range(len(fish["group1"])): 2213 sign1 = "ns" 2214 if fish["adj_p_val"][i] < 0.001: 2215 sign1 = "***" 2216 elif fish["adj_p_val"][i] < 0.01: 2217 sign1 = "**" 2218 elif fish["adj_p_val"][i] < 0.05: 2219 sign1 = "*" 2220 2221 sign2 = "ns" 2222 if ks["adj_p_val"][i] < 0.001: 2223 sign2 = "***" 2224 elif ks["adj_p_val"][i] < 0.01: 2225 sign2 = "**" 2226 elif ks["adj_p_val"][i] < 0.05: 2227 sign2 = "*" 2228 2229 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2230 2231 else: 2232 for i in range(len(fish["group1"])): 2233 sign1 = "ns" 2234 if fish["p_val"][i] < 0.001: 2235 sign1 = "***" 2236 elif fish["p_val"][i] < 0.01: 2237 sign1 = "**" 2238 elif fish["p_val"][i] < 0.05: 2239 sign1 = "*" 2240 2241 sign2 = "ns" 2242 if ks["p_val"][i] < 0.001: 2243 sign2 = "***" 2244 elif ks["p_val"][i] < 0.01: 2245 sign2 = "**" 2246 elif ks["p_val"][i] < 0.05: 2247 sign2 = "*" 2248 2249 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2250 2251 axs[0, -1].text( 2252 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2253 ) 2254 axs[0, -1].set_axis_off() 2255 2256 plt.tight_layout() 2257 2258 if cfg._DISPLAY_MODE: 2259 plt.show() 2260 2261 return fig
23class FeatureIntensity(ImageTools): 24 r""" 25 Class for quantitative analysis of pixel intensity and size measurements 26 in 2D/3D biological images. Supports projection of 3D stacks, mask-based 27 intensity normalization, region size estimation and metadata extraction. 28 29 Parameters 30 ---------- 31 input_image : ndarray, optional 32 Input image or 3D stack for analysis. If 3D, projection must be applied. 33 34 image : ndarray, optional 35 2D projected image (internal use). 36 37 normalized_image_values : dict, optional 38 Dictionary storing normalized intensity statistics. 39 40 mask : ndarray, optional 41 Binary mask of region of interest (ROI). 42 43 background_mask : ndarray, optional 44 Binary mask used for background estimation. If not provided, `mask` is used. 45 46 typ : {"avg", "median", "std", "var", "max", "min"}, optional 47 Projection type for 3D stacks. Default is `"avg"`. 48 49 size_info : dict, optional 50 Dictionary storing ROI size measurements. 51 52 correction_factor : float, optional 53 Normalization correction factor applied to background intensity. 54 Must satisfy 0 < factor < 1. Default is 0.1. 55 56 img_type : str, optional 57 Image type metadata. 58 59 scale : float, optional 60 Pixel resolution in physical units (e.g. µm/px). Used in size calculations. 61 62 stack_selection : list of int, optional 63 List of Z-indices to remove when projecting a 3D image. 64 65 Attributes 66 ---------- 67 input_image : ndarray or None 68 Loaded input image. 69 70 image : ndarray or None 71 Projected 2D image. 72 73 mask : ndarray or None 74 Region of interest mask. 75 76 background_mask : ndarray or None 77 Background normalization mask. 78 79 scale : float or None 80 Scale value for size estimation. 81 82 normalized_image_values : dict or None 83 Dictionary containing intensity metrics. 84 85 size_info : dict or None 86 Dictionary with ROI size information. 87 88 typ : str 89 Selected projection type for 3D images. 90 91 stack_selection : list of int 92 Z-levels excluded from projection. 93 94 Notes 95 ----- 96 The intensity normalization formula applied per pixel is: 97 98 .. math:: 99 100 R_{i,j} = T_{i,j} - \\left( \\mu_B (1 + c) \\right) 101 102 where 103 * ``T_{i,j}`` – pixel intensity in ROI 104 * ``μ_B`` – mean intensity in background region 105 * ``c`` – correction factor 106 * ``R_{i,j}`` – normalized pixel intensity 107 108 Examples 109 -------- 110 Load a 3D image, mask and compute statistics: 111 112 >>> fi = FeatureIntensity() 113 >>> fi.load_image_3D("stack.tiff") 114 >>> fi.load_mask_("mask.png") 115 >>> fi.set_projection("median") 116 >>> fi.run_calculations() 117 >>> results = fi.get_results() 118 >>> results["intensity"]["norm_mean"] 119 """ 120 121 def __init__( 122 self, 123 input_image=None, 124 image=None, 125 normalized_image_values=None, 126 mask=None, 127 background_mask=None, 128 typ=None, 129 size_info=None, 130 correction_factor=None, 131 img_type=None, 132 scale=None, 133 stack_selection=None, 134 ): 135 """ 136 Initialize a FeatureIntensity analysis instance. 137 138 Parameters 139 ---------- 140 input_image : ndarray, optional 141 Input image or 3D stack used for analysis. If the image is 3D, a 142 projection will be computed depending on the `typ` parameter. 143 144 image : ndarray, optional 145 2D image buffer used internally after projection of the input image. 146 Should not be set manually. 147 148 normalized_image_values : dict, optional 149 Dictionary containing normalized intensity statistics. Usually filled 150 automatically after running `run_calculations()`. 151 152 mask : ndarray, optional 153 Binary mask of the target region of interest (ROI). Required for 154 intensity and size calculations. 155 156 background_mask : ndarray, optional 157 Binary mask specifying the background region used to compute the 158 normalization threshold. If not provided, the ROI mask is also used 159 as the background reference. 160 161 typ : {"avg", "median", "std", "var", "max", "min"}, optional 162 Projection method for 3D images. Determines how the z-stack is 163 collapsed into a 2D image. Default is `"avg"`. 164 165 size_info : dict, optional 166 Dictionary storing computed size metrics of the ROI. Populated after 167 invoking `size_calculations()`. 168 169 correction_factor : float, optional 170 Correction term used during intensity normalization. Must satisfy 171 0 < correction_factor < 1. Default is 0.1. 172 173 img_type : str, optional 174 Optional metadata about the image type (e.g., "tiff", "png"). 175 176 scale : float, optional 177 Pixel resolution in physical units (e.g., µm/px). Required for 178 real-size estimation in `size_calculations()`. 179 180 stack_selection : list of int, optional 181 Indices of z-planes to exclude during projection of a 3D stack. 182 183 Notes 184 ----- 185 Values not provided are initialized to `None`, except for `typ`, which 186 defaults to `"avg"`, and `correction_factor`, which defaults to 0.1. 187 188 The class is designed to be populated by loading functions: 189 `load_image_()`, `load_image_3D()`, `load_mask_()`, 190 and optionally `load_normalization_mask_()` and `load_JIMG_project_()`. 191 """ 192 193 self.input_image = input_image or None 194 """ Input image or 3D stack used for analysis. If the image is 3D, a 195 projection will be computed depending on the `typ` parameter.""" 196 197 self.image = image or None 198 """ 2D image buffer used internally after projection of the input image. 199 Should not be set manually.""" 200 201 self.normalized_image_values = normalized_image_values or None 202 """Dictionary containing normalized intensity statistics. Usually filled 203 automatically after running `run_calculations()`.""" 204 205 self.mask = mask or None 206 """Binary mask of the target region of interest (ROI). Required for 207 intensity and size calculations.""" 208 209 self.background_mask = background_mask or None 210 """ Binary mask specifying the background region used to compute the 211 normalization threshold. If not provided, the ROI mask is also used 212 as the background reference.""" 213 214 self.typ = typ or "avg" 215 """Projection method for 3D images. Determines how the z-stack is 216 collapsed into a 2D image. Default is `"avg"`.""" 217 218 self.size_info = size_info or None 219 """Dictionary storing computed size metrics of the ROI. Populated after 220 invoking `size_calculations()`.""" 221 222 self.correction_factor = correction_factor or 0.1 223 """ Correction term used during intensity normalization. Must satisfy 224 0 < correction_factor < 1. Default is 0.1.""" 225 226 self.scale = scale or None 227 """ Pixel resolution in physical units (e.g., µm/px). Required for 228 real-size estimation in `size_calculations()`.""" 229 230 self.stack_selection = stack_selection or [] 231 """Indices of z-planes to exclude during projection of a 3D stack.""" 232 233 @property 234 def current_metadata(self): 235 r""" 236 Return current metadata parameters used in image processing and normalization. 237 238 Returns 239 ------- 240 tuple 241 A tuple containing: 242 243 projection_type : str 244 Projection method used for 3D image reduction (e.g., "avg", "median"). 245 246 correction_factor : float 247 Correction factor used for background subtraction during intensity 248 normalization. The applied formula is: 249 250 .. math:: 251 252 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 253 254 where 255 * ``R_{i,j}`` — normalized pixel intensity 256 * ``T_{i,j}`` — original pixel intensity 257 * ``μ_B`` — mean background intensity 258 * ``c`` — correction factor 259 scale : float or None 260 Pixel resolution (unit/px), loaded via `load_JIMG_project_()` or set manually 261 using `set_scale()`. 262 263 stack_selection : list of int 264 Indices of z-slices excluded from projection of a 3D image. 265 266 Notes 267 ----- 268 This property also prints the metadata values to the console for quick inspection. 269 """ 270 271 print(f"Projection type: {self.typ}") 272 print(f"Correction factor: {self.correction_factor}") 273 print(f"Scale (unit/px): {self.scale}") 274 print(f"Selected stac to remove: {self.stack_selection}") 275 276 return self.typ, self.correction_factor, self.scale, self.stack_selection 277 278 def set_projection(self, projection: str): 279 """ 280 Set the projection method for 3D image stack reduction. 281 282 Parameters 283 ---------- 284 projection : {"avg", "median", "std", "var", "max", "min"} 285 Projection method to reduce a 3D image stack to a 2D image. Default is `"avg"`. 286 287 Notes 288 ----- 289 This method updates the `typ` attribute of the class. The selected projection 290 determines how the z-stack is collapsed: 291 - `"avg"` : average intensity across slices 292 - `"median"` : median intensity across slices 293 - `"std"` : standard deviation across slices 294 - `"var"` : variance across slices 295 - `"max"` : maximum intensity across slices 296 - `"min"` : minimum intensity across slices 297 """ 298 299 t = ["avg", "median", "std", "var", "max", "min"] 300 if projection in t: 301 self.typ = projection 302 else: 303 print(f"\nProvided parameter is incorrect. Avaiable projection types: {t}") 304 305 def set_correction_factorn(self, factor: float): 306 r""" 307 Set the correction factor for background subtraction during intensity normalization. 308 309 Parameters 310 ---------- 311 factor : float 312 Correction factor to adjust background subtraction. Must satisfy 0 < factor < 1. 313 Default is 0.1. 314 315 Notes 316 ----- 317 The correction is applied per pixel in the target mask using the formula: 318 319 .. math:: 320 321 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 322 323 where 324 * ``R_{i,j}`` — normalized pixel intensity 325 * ``T_{i,j}`` — original pixel intensity 326 * ``μ_B`` — mean intensity in the background mask 327 * ``c`` — correction factor 328 """ 329 330 if factor < 1 and factor > 0: 331 self.correction_factor = factor 332 else: 333 print( 334 "\nProvided parameter is incorrect. The factor should be a floating-point value within the range of 0 to 1." 335 ) 336 337 def set_scale(self, scale): 338 """ 339 Set the scale for converting pixel measurements to physical units. 340 341 Parameters 342 ---------- 343 scale : float 344 Pixel resolution in physical units (e.g., µm/px). Used to calculate the 345 actual size of the tissue or organ. 346 347 Notes 348 ----- 349 The scale can also be automatically loaded from a JIMG project using 350 `load_JIMG_project_()`. This value is required for size calculations in 351 `size_calculations()`. 352 """ 353 354 self.scale = scale 355 356 def set_selection_list(self, rm_list: list): 357 """ 358 Set the list of z-slices to exclude when projecting a 3D image stack. 359 360 Parameters 361 ---------- 362 rm_list : list of int 363 List of indices corresponding to z-slices that should be removed from 364 the full 3D image stack before projection. 365 366 Notes 367 ----- 368 This updates the `stack_selection` attribute, which is used by the 369 `stack_selection_()` method during projection. 370 """ 371 372 self.stack_selection = rm_list 373 374 def load_JIMG_project_(self, path): 375 """ 376 Load a JIMG project from a `.pjm` file. 377 378 Parameters 379 ---------- 380 file_path : str 381 Path to the JIMG project file. The file must have a `.pjm` extension. 382 383 Returns 384 ------- 385 project : object 386 Loaded project object containing images and metadata. 387 388 Raises 389 ------ 390 ValueError 391 If the provided file path does not point to a `.pjm` file. 392 393 Notes 394 ----- 395 The method attempts to automatically set the `scale` and `stack_selection` 396 attributes from the project metadata if available. 397 """ 398 399 path = os.path.abspath(path) 400 401 if ".pjm" in path: 402 metadata = self.load_JIMG_project(path) 403 404 try: 405 self.scale = metadata.metadata["X_resolution[um/px]"] 406 except: 407 408 try: 409 self.scale = metadata.images_dict["metadata"][0][ 410 "X_resolution[um/px]" 411 ] 412 413 except: 414 print( 415 '\nUnable to set scale on this project! Set scale using "set_scale()"' 416 ) 417 418 self.stack_selection = metadata.removal_list 419 420 else: 421 print( 422 "\nWrong path. The provided path does not point to a JIMG project (*.pjm)." 423 ) 424 425 def stack_selection_(self): 426 """ 427 Remove selected z-slices from a 3D image stack based on `stack_selection`. 428 429 Notes 430 ----- 431 Only works if `input_image` is a 3D ndarray. The slices with indices listed 432 in `stack_selection` are excluded from the stack. Updates `input_image` 433 in-place. 434 435 Prints a warning if `stack_selection` is empty. 436 """ 437 438 if len(self.input_image.shape) == 3: 439 if len(self.stack_selection) > 0: 440 self.input_image = self.input_image[ 441 [ 442 x 443 for x in range(self.input_image.shape[0]) 444 if x not in self.stack_selection 445 ] 446 ] 447 else: 448 print("\nImages to remove from the stack were not selected!") 449 450 def projection(self): 451 """ 452 Project a 3D image stack into a 2D image using the method defined by `typ`. 453 454 Notes 455 ----- 456 Updates the `image` attribute with the projected 2D result. 457 458 Supported projection types (`typ`): 459 - "avg" : mean intensity across slices 460 - "median" : median intensity across slices 461 - "std" : standard deviation across slices 462 - "var" : variance across slices 463 - "max" : maximum intensity across slices 464 - "min" : minimum intensity across slices 465 466 Raises 467 ------ 468 AttributeError 469 If `input_image` is not defined. 470 """ 471 472 if self.typ == "avg": 473 img = np.mean(self.input_image, axis=0) 474 475 elif self.typ == "std": 476 img = np.std(self.input_image, axis=0) 477 478 elif self.typ == "median": 479 img = np.median(self.input_image, axis=0) 480 481 elif self.typ == "var": 482 img = np.var(self.input_image, axis=0) 483 484 elif self.typ == "max": 485 img = np.max(self.input_image, axis=0) 486 487 elif self.typ == "min": 488 img = np.min(self.input_image, axis=0) 489 490 self.image = img 491 492 def detect_img(self): 493 """ 494 Detect whether the input image is 2D or 3D and perform appropriate preprocessing. 495 496 Notes 497 ----- 498 - For 3D images, applies `stack_selection_()` and then `projection()`. 499 - For 2D images, no projection is applied. 500 - Prints status messages indicating the type of image and applied operations. 501 502 Raises 503 ------ 504 AttributeError 505 If `input_image` is not defined. 506 """ 507 check = len(self.input_image.shape) 508 509 if check == 3: 510 print("\n3D image detected! Starting processing for 3D image...") 511 print(f"Projection - {self.typ}...") 512 513 self.stack_selection_() 514 self.projection() 515 516 elif check == 2: 517 print("\n2D image detected! Starting processing for 2D image...") 518 519 else: 520 print("\nData does not match any image type!") 521 522 def load_image_3D(self, path): 523 """ 524 Load a 3D image stack from a TIFF file. 525 526 Parameters 527 ---------- 528 path : str 529 Path to the 3D image file (*.tiff) to be loaded. 530 531 Notes 532 ----- 533 The loaded image is stored in the `input_image` attribute as a 3D ndarray. 534 """ 535 path = os.path.abspath(path) 536 537 self.input_image = self.load_3D_tiff(path) 538 539 def load_image_(self, path): 540 """ 541 Load a 2D image into the class. 542 543 Parameters 544 ---------- 545 path : str 546 Path to the image file to be loaded. 547 548 Notes 549 ----- 550 The loaded image is stored in the `input_image` attribute as a 2D ndarray. 551 """ 552 path = os.path.abspath(path) 553 554 self.input_image = self.load_image(path) 555 556 def load_mask_(self, path): 557 r""" 558 Load a binary mask into the class and optionally set it as the normalization mask. 559 560 Parameters 561 ---------- 562 path : str 563 Path to the mask image file. Supported formats include 8-bit or 16-bit images 564 with extensions such as `.png` or `.jpeg`. The mask must be binary 565 (e.g., 0/255, 0/2**16-1, 0/1). 566 567 Notes 568 ----- 569 - If `load_normalization_mask_()` is not called, this mask is also used as the 570 background mask for intensity normalization. 571 - Normalization is applied per pixel using the formula: 572 573 .. math:: 574 575 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 576 577 where 578 * ``R_{i,j}`` — normalized pixel intensity 579 * ``T_{i,j}`` — pixel intensity in the target mask 580 * ``μ_B`` — mean intensity of the background (reversed mask) 581 * ``c`` — correction factor 582 """ 583 584 path = os.path.abspath(path) 585 586 self.mask = self.load_mask(path) 587 588 print( 589 "\nThis mask was also set as the reverse background mask. If you want a different background mask for normalization, use 'load_normalization_mask()'." 590 ) 591 self.background_mask = self.load_mask(path) 592 593 def load_normalization_mask_(self, path): 594 r""" 595 Load a binary mask for normalization into the class. 596 597 Parameters 598 ---------- 599 path : str 600 Path to the mask image file. Supported formats include 8-bit or 16-bit 601 images (e.g., `.png`, `.jpeg`). The mask must be binary (0/255, 0/2**16-1, 0/1). 602 603 Notes 604 ----- 605 - The mask defines the area of interest. Normalization is applied to the inverse 606 of this area (reversed mask). 607 - Normalization formula applied per pixel: 608 609 .. math:: 610 611 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 612 613 where 614 * ``R_{i,j}`` — normalized pixel intensity 615 * ``T_{i,j}`` — pixel intensity in the target mask 616 * ``μ_B`` — mean intensity of the background (reversed mask) 617 * ``c`` — correction factor 618 """ 619 620 path = os.path.abspath(path) 621 622 self.background_mask = self.load_mask(path) 623 624 def intensity_calculations(self): 625 """ 626 Calculate normalized and raw intensity statistics from the image based on masks. 627 628 This method performs intensity calculations using the main mask (`self.mask`) 629 and the background mask (`self.background_mask`). The pixel intensities within 630 the mask of interest are normalized by subtracting a threshold derived from the 631 background region and applying a correction factor (`self.correction_factor`). 632 Negative values after normalization are clipped to zero. 633 634 The following statistics are computed for both normalized and raw values: 635 - Minimum 636 - Maximum 637 - Mean 638 - Median 639 - Standard deviation 640 - Variance 641 - List of all normalized values (only for normalized data) 642 643 Notes 644 ----- 645 - The method updates the instance attribute `self.normalized_image_values` 646 with a dictionary containing both normalized and raw statistics. 647 - Normalization formula applied for each pixel in the selected mask: 648 final_val = selected_value - (threshold + threshold * correction_factor) 649 where threshold is the mean intensity in the background mask. 650 - Negative values after normalization are set to zero. 651 """ 652 653 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 654 tmp_bmask = self.ajd_mask_size(image=self.image, mask=self.background_mask) 655 656 selected_values = self.image[tmp_mask == np.max(tmp_mask)] 657 658 threshold = np.mean(self.image[tmp_bmask == np.min(tmp_bmask)]) 659 660 # normalization 661 final_val = selected_values - (threshold + (threshold * self.correction_factor)) 662 663 final_val[final_val < 0] = 0 664 665 tmp_dict = { 666 "norm_min": np.min(final_val), 667 "norm_max": np.max(final_val), 668 "norm_mean": np.mean(final_val), 669 "norm_median": np.median(final_val), 670 "norm_std": np.std(final_val), 671 "norm_var": np.var(final_val), 672 "norm_values": final_val.tolist(), 673 "min": np.min(selected_values), 674 "max": np.max(selected_values), 675 "mean": np.mean(selected_values), 676 "median": np.median(selected_values), 677 "std": np.std(selected_values), 678 "var": np.var(selected_values), 679 } 680 681 self.normalized_image_values = tmp_dict 682 683 def size_calculations(self): 684 """ 685 Calculates the size and bounding dimensions of the masked region in the image. 686 687 This method computes the following metrics based on the current mask: 688 - Total number of pixels in the mask (`px_size`) 689 - Real-world size if a scale is provided (`size`) 690 - Maximum lengths along x and y axes (`max_length_x_axis`, `max_length_y_axis`) 691 692 If `self.scale` is defined (unit per pixel), the real-world size is calculated. 693 If not, `size` will be `None` and a warning message is printed. 694 695 Returns: 696 Updates the following attributes in the class: 697 - self.size_info (dict) containing: 698 - 'size' (float or None): real-world size of the mask 699 - 'px_size' (int): number of pixels in the masked region 700 - 'max_length_x_axis' (int): length of the bounding box along the x-axis 701 - 'max_length_y_axis' (int): length of the bounding box along the y-axis 702 703 Example: 704 analysis.size_calculations() 705 print(analysis.size_info) 706 """ 707 708 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 709 710 size_px = int(len(tmp_mask[tmp_mask > np.min(tmp_mask)])) 711 712 if self.scale is not None: 713 size = float(size_px * self.scale) 714 else: 715 size = None 716 print( 717 '\nUnable to calculate real size, scale (unit/px) not provided, use "set_scale()" or load JIMG project .pjm metadata "load_pjm()" to set scale for calculations!' 718 ) 719 720 non_zero_indices = np.where(tmp_mask == np.max(tmp_mask)) 721 722 min_y, max_y = np.min(non_zero_indices[0]), np.max(non_zero_indices[0]) 723 min_x, max_x = np.min(non_zero_indices[1]), np.max(non_zero_indices[1]) 724 725 max_length_x_axis = int(max_x - min_x + 1) 726 max_length_y_axis = int(max_y - min_y + 1) 727 728 tmp_val = { 729 "size": size, 730 "px_size": size_px, 731 "max_length_x_axis": max_length_x_axis, 732 "max_length_y_axis": max_length_y_axis, 733 } 734 735 self.size_info = tmp_val 736 737 def run_calculations(self): 738 """ 739 Run the full analysis pipeline on the loaded image using the provided masks. 740 741 Notes 742 ----- 743 - The input image must be loaded via `load_image_()` or `load_image_3D()`. 744 - The ROI mask must be loaded via `load_mask_()`. Optionally, a normalization 745 mask can be loaded via `load_normalization_mask_()`. 746 - Parameters such as projection type and correction factor can be set with 747 `set_projection()` and `set_correction_factor()`. 748 - Scale and stack selection can also influence calculations if defined. 749 - To view current parameters, use the `current_metadata` property. 750 751 Returns 752 ------- 753 None 754 The results are stored internally and can be retrieved using 755 `get_results()`. 756 """ 757 758 if self.input_image is not None: 759 760 if self.mask is not None: 761 762 print("\nStart...") 763 self.detect_img() 764 self.intensity_calculations() 765 self.size_calculations() 766 print("\nCompleted!") 767 768 def get_results(self): 769 """ 770 Return the results from the analysis performed by `run_calculations()`. 771 772 Returns 773 ------- 774 results_dict : dict or None 775 Dictionary containing intensity and size results. Structure: 776 - 'intensity' : dict with normalized and raw intensity statistics 777 - 'size' : dict with ROI size metrics 778 779 Notes 780 ----- 781 If analysis has not been run yet, prints a message and returns None. 782 """ 783 784 if self.normalized_image_values is not None and self.size_info is not None: 785 786 results = { 787 "intensity": self.normalized_image_values, 788 "size": self.size_info, 789 } 790 791 return results 792 793 else: 794 print('\nAnalysis were not conducted. Run analysis "run_calculations()"') 795 796 def save_results( 797 self, 798 path="", 799 mask_region: str = "", 800 feature_name: str = "", 801 individual_number: int = 0, 802 individual_name: str = "", 803 ): 804 """ 805 Save the analysis results to a `.int` (JSON) file. 806 807 Parameters 808 ---------- 809 path : str, optional 810 Directory path where the file will be saved. Defaults to the current working directory. 811 812 mask_region : str 813 Name or identifier of the mask region (e.g., tissue, part of tissue). 814 815 feature_name : str 816 Name of the feature being analyzed. Underscores or spaces are replaced with periods. 817 818 individual_number : int 819 Unique identifier for the individual in the analysis (e.g., 1, 2, 3). 820 821 individual_name : str 822 Name of the individual (e.g., species name, tissue, organoid). 823 824 Notes 825 ----- 826 - The method validates that all required parameters are provided and that 827 analysis results exist (`normalized_image_values` and `size_info`). 828 - Creates the directory if it does not exist. 829 - File name format: 830 '<individual_name>_<individual_number>_<mask_region>_<feature_name>.int' 831 832 Raises 833 ------ 834 FileNotFoundError 835 If the specified path cannot be created or accessed. 836 837 ValueError 838 If any of `mask_region`, `feature_name`, `individual_number`, or 839 `individual_name` are missing or invalid. 840 """ 841 842 path = os.path.abspath(path) 843 844 if ( 845 len(mask_region) > 1 846 and len(feature_name) > 1 847 and individual_number != 0 848 and len(individual_name) > 1 849 ): 850 851 if self.normalized_image_values is not None and self.size_info is not None: 852 853 results = { 854 "intensity": self.normalized_image_values, 855 "size": self.size_info, 856 } 857 858 mask_region = re.sub(r"[_\s]+", ".", mask_region) 859 feature_name = re.sub(r"[_\s]+", ".", feature_name) 860 individual_number = re.sub(r"[_\s]+", ".", str(individual_number)) 861 individual_name = re.sub(r"[_\s]+", ".", individual_name) 862 863 full_name = f"{individual_name}_{individual_number}_{mask_region}_{feature_name}" 864 865 isExist = os.path.exists(path) 866 if not isExist: 867 os.makedirs(path, exist_ok=True) 868 869 full_path = os.path.join( 870 path, re.sub("\\.json", "", full_name) + ".int" 871 ) 872 873 with open(full_path, "w") as file: 874 json.dump(results, file, indent=4) 875 876 else: 877 print( 878 '\nAnalysis were not conducted. Run analysis "run_calculations()"' 879 ) 880 881 else: 882 print( 883 "\nAny of 'mask_region', 'feature_name', 'individual_number', 'individual_name' parameters were provided wrong!" 884 ) 885 886 def concatenate_intensity_data(self, directory: str = "", name: str = ""): 887 """ 888 Concatenate intensity data from multiple `.int` files and save as CSV. 889 890 Parameters 891 ---------- 892 directory : str, optional 893 Path to the directory containing `.int` files. Defaults to the current working directory. 894 895 name : str 896 Prefix for the output CSV file names. CSV files are saved in the format 897 '<name>_<gene>_<region>.csv'. 898 899 Raises 900 ------ 901 FileNotFoundError 902 If the directory cannot be accessed or no `.int` files are found. 903 904 ValueError 905 If an `.int` file is missing expected data or has an incorrect format. 906 907 Notes 908 ----- 909 - The method groups intensity data by gene (feature) and mask region. 910 - Outputs one CSV file per unique gene-region combination, saved in the specified directory. 911 """ 912 913 directory = os.path.abspath(directory) 914 915 files_list = [f for f in os.listdir(directory) if f.endswith(".int")] 916 917 genes_set = set([re.sub("\\.int", "", x.split("_")[3]) for x in files_list]) 918 regions_set = set([re.sub("\\.int", "", x.split("_")[2]) for x in files_list]) 919 920 for g in genes_set: 921 for r in regions_set: 922 json_to_save = { 923 "individual_name": [], 924 "individual_number": [], 925 "norm_intensity": [], 926 "size": [], 927 } 928 929 for f in tqdm(files_list): 930 if g in f and r in f: 931 with open(os.path.join(directory, f), "r") as file: 932 data = json.load(file) 933 934 json_to_save["norm_intensity"] = ( 935 json_to_save["norm_intensity"] 936 + data["intensity"]["norm_values"] 937 ) 938 json_to_save["individual_name"] = json_to_save[ 939 "individual_name" 940 ] + [f.split("_")[0]] * len( 941 data["intensity"]["norm_values"] 942 ) 943 json_to_save["individual_number"] = json_to_save[ 944 "individual_number" 945 ] + [f.split("_")[1]] * len( 946 data["intensity"]["norm_values"] 947 ) 948 json_to_save["size"] = json_to_save["size"] + [ 949 data["size"]["px_size"] 950 ] * len(data["intensity"]["norm_values"]) 951 952 pd.DataFrame(json_to_save).to_csv(f"{name}_{g}_{r}.csv", index=False)
Class for quantitative analysis of pixel intensity and size measurements in 2D/3D biological images. Supports projection of 3D stacks, mask-based intensity normalization, region size estimation and metadata extraction.
Parameters
input_image : ndarray, optional Input image or 3D stack for analysis. If 3D, projection must be applied.
image : ndarray, optional 2D projected image (internal use).
normalized_image_values : dict, optional Dictionary storing normalized intensity statistics.
mask : ndarray, optional Binary mask of region of interest (ROI).
background_mask : ndarray, optional
Binary mask used for background estimation. If not provided, mask is used.
typ : {"avg", "median", "std", "var", "max", "min"}, optional
Projection type for 3D stacks. Default is "avg".
size_info : dict, optional Dictionary storing ROI size measurements.
correction_factor : float, optional Normalization correction factor applied to background intensity. Must satisfy 0 < factor < 1. Default is 0.1.
img_type : str, optional Image type metadata.
scale : float, optional Pixel resolution in physical units (e.g. µm/px). Used in size calculations.
stack_selection : list of int, optional List of Z-indices to remove when projecting a 3D image.
Attributes
input_image : ndarray or None Loaded input image.
image : ndarray or None Projected 2D image.
mask : ndarray or None Region of interest mask.
background_mask : ndarray or None Background normalization mask.
scale : float or None Scale value for size estimation.
normalized_image_values : dict or None Dictionary containing intensity metrics.
size_info : dict or None Dictionary with ROI size information.
typ : str Selected projection type for 3D images.
stack_selection : list of int Z-levels excluded from projection.
Notes
The intensity normalization formula applied per pixel is:
$$R_{i,j} = T_{i,j} - \left( \mu_B (1 + c) \right)$$
where
T_{i,j}– pixel intensity in ROIμ_B– mean intensity in background regionc– correction factorR_{i,j}– normalized pixel intensity
Examples
Load a 3D image, mask and compute statistics:
>>> fi = FeatureIntensity()
>>> fi.load_image_3D("stack.tiff")
>>> fi.load_mask_("mask.png")
>>> fi.set_projection("median")
>>> fi.run_calculations()
>>> results = fi.get_results()
>>> results["intensity"]["norm_mean"]
FeatureIntensity( input_image=None, image=None, normalized_image_values=None, mask=None, background_mask=None, typ=None, size_info=None, correction_factor=None, img_type=None, scale=None, stack_selection=None)
121 def __init__( 122 self, 123 input_image=None, 124 image=None, 125 normalized_image_values=None, 126 mask=None, 127 background_mask=None, 128 typ=None, 129 size_info=None, 130 correction_factor=None, 131 img_type=None, 132 scale=None, 133 stack_selection=None, 134 ): 135 """ 136 Initialize a FeatureIntensity analysis instance. 137 138 Parameters 139 ---------- 140 input_image : ndarray, optional 141 Input image or 3D stack used for analysis. If the image is 3D, a 142 projection will be computed depending on the `typ` parameter. 143 144 image : ndarray, optional 145 2D image buffer used internally after projection of the input image. 146 Should not be set manually. 147 148 normalized_image_values : dict, optional 149 Dictionary containing normalized intensity statistics. Usually filled 150 automatically after running `run_calculations()`. 151 152 mask : ndarray, optional 153 Binary mask of the target region of interest (ROI). Required for 154 intensity and size calculations. 155 156 background_mask : ndarray, optional 157 Binary mask specifying the background region used to compute the 158 normalization threshold. If not provided, the ROI mask is also used 159 as the background reference. 160 161 typ : {"avg", "median", "std", "var", "max", "min"}, optional 162 Projection method for 3D images. Determines how the z-stack is 163 collapsed into a 2D image. Default is `"avg"`. 164 165 size_info : dict, optional 166 Dictionary storing computed size metrics of the ROI. Populated after 167 invoking `size_calculations()`. 168 169 correction_factor : float, optional 170 Correction term used during intensity normalization. Must satisfy 171 0 < correction_factor < 1. Default is 0.1. 172 173 img_type : str, optional 174 Optional metadata about the image type (e.g., "tiff", "png"). 175 176 scale : float, optional 177 Pixel resolution in physical units (e.g., µm/px). Required for 178 real-size estimation in `size_calculations()`. 179 180 stack_selection : list of int, optional 181 Indices of z-planes to exclude during projection of a 3D stack. 182 183 Notes 184 ----- 185 Values not provided are initialized to `None`, except for `typ`, which 186 defaults to `"avg"`, and `correction_factor`, which defaults to 0.1. 187 188 The class is designed to be populated by loading functions: 189 `load_image_()`, `load_image_3D()`, `load_mask_()`, 190 and optionally `load_normalization_mask_()` and `load_JIMG_project_()`. 191 """ 192 193 self.input_image = input_image or None 194 """ Input image or 3D stack used for analysis. If the image is 3D, a 195 projection will be computed depending on the `typ` parameter.""" 196 197 self.image = image or None 198 """ 2D image buffer used internally after projection of the input image. 199 Should not be set manually.""" 200 201 self.normalized_image_values = normalized_image_values or None 202 """Dictionary containing normalized intensity statistics. Usually filled 203 automatically after running `run_calculations()`.""" 204 205 self.mask = mask or None 206 """Binary mask of the target region of interest (ROI). Required for 207 intensity and size calculations.""" 208 209 self.background_mask = background_mask or None 210 """ Binary mask specifying the background region used to compute the 211 normalization threshold. If not provided, the ROI mask is also used 212 as the background reference.""" 213 214 self.typ = typ or "avg" 215 """Projection method for 3D images. Determines how the z-stack is 216 collapsed into a 2D image. Default is `"avg"`.""" 217 218 self.size_info = size_info or None 219 """Dictionary storing computed size metrics of the ROI. Populated after 220 invoking `size_calculations()`.""" 221 222 self.correction_factor = correction_factor or 0.1 223 """ Correction term used during intensity normalization. Must satisfy 224 0 < correction_factor < 1. Default is 0.1.""" 225 226 self.scale = scale or None 227 """ Pixel resolution in physical units (e.g., µm/px). Required for 228 real-size estimation in `size_calculations()`.""" 229 230 self.stack_selection = stack_selection or [] 231 """Indices of z-planes to exclude during projection of a 3D stack."""
Initialize a FeatureIntensity analysis instance.
Parameters
input_image : ndarray, optional
Input image or 3D stack used for analysis. If the image is 3D, a
projection will be computed depending on the typ parameter.
image : ndarray, optional 2D image buffer used internally after projection of the input image. Should not be set manually.
normalized_image_values : dict, optional
Dictionary containing normalized intensity statistics. Usually filled
automatically after running run_calculations().
mask : ndarray, optional Binary mask of the target region of interest (ROI). Required for intensity and size calculations.
background_mask : ndarray, optional Binary mask specifying the background region used to compute the normalization threshold. If not provided, the ROI mask is also used as the background reference.
typ : {"avg", "median", "std", "var", "max", "min"}, optional
Projection method for 3D images. Determines how the z-stack is
collapsed into a 2D image. Default is "avg".
size_info : dict, optional
Dictionary storing computed size metrics of the ROI. Populated after
invoking size_calculations().
correction_factor : float, optional Correction term used during intensity normalization. Must satisfy 0 < correction_factor < 1. Default is 0.1.
img_type : str, optional Optional metadata about the image type (e.g., "tiff", "png").
scale : float, optional
Pixel resolution in physical units (e.g., µm/px). Required for
real-size estimation in size_calculations().
stack_selection : list of int, optional Indices of z-planes to exclude during projection of a 3D stack.
Notes
Values not provided are initialized to None, except for typ, which
defaults to "avg", and correction_factor, which defaults to 0.1.
The class is designed to be populated by loading functions:
load_image_(), load_image_3D(), load_mask_(),
and optionally load_normalization_mask_() and load_JIMG_project_().
input_image
Input image or 3D stack used for analysis. If the image is 3D, a
projection will be computed depending on the typ parameter.
image
2D image buffer used internally after projection of the input image. Should not be set manually.
normalized_image_values
Dictionary containing normalized intensity statistics. Usually filled
automatically after running run_calculations().
mask
Binary mask of the target region of interest (ROI). Required for intensity and size calculations.
background_mask
Binary mask specifying the background region used to compute the normalization threshold. If not provided, the ROI mask is also used as the background reference.
typ
Projection method for 3D images. Determines how the z-stack is
collapsed into a 2D image. Default is "avg".
size_info
Dictionary storing computed size metrics of the ROI. Populated after
invoking size_calculations().
correction_factor
Correction term used during intensity normalization. Must satisfy 0 < correction_factor < 1. Default is 0.1.
scale
Pixel resolution in physical units (e.g., µm/px). Required for
real-size estimation in size_calculations().
current_metadata
233 @property 234 def current_metadata(self): 235 r""" 236 Return current metadata parameters used in image processing and normalization. 237 238 Returns 239 ------- 240 tuple 241 A tuple containing: 242 243 projection_type : str 244 Projection method used for 3D image reduction (e.g., "avg", "median"). 245 246 correction_factor : float 247 Correction factor used for background subtraction during intensity 248 normalization. The applied formula is: 249 250 .. math:: 251 252 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 253 254 where 255 * ``R_{i,j}`` — normalized pixel intensity 256 * ``T_{i,j}`` — original pixel intensity 257 * ``μ_B`` — mean background intensity 258 * ``c`` — correction factor 259 scale : float or None 260 Pixel resolution (unit/px), loaded via `load_JIMG_project_()` or set manually 261 using `set_scale()`. 262 263 stack_selection : list of int 264 Indices of z-slices excluded from projection of a 3D image. 265 266 Notes 267 ----- 268 This property also prints the metadata values to the console for quick inspection. 269 """ 270 271 print(f"Projection type: {self.typ}") 272 print(f"Correction factor: {self.correction_factor}") 273 print(f"Scale (unit/px): {self.scale}") 274 print(f"Selected stac to remove: {self.stack_selection}") 275 276 return self.typ, self.correction_factor, self.scale, self.stack_selection
Return current metadata parameters used in image processing and normalization.
Returns
tuple A tuple containing:
projection_type : str
Projection method used for 3D image reduction (e.g., "avg", "median").
correction_factor : float
Correction factor used for background subtraction during intensity
normalization. The applied formula is:
$$R_{i,j} = T_{i,j} - ( \mu_B (1 + c) )$$
where
* ``R_{i,j}`` — normalized pixel intensity
* ``T_{i,j}`` — original pixel intensity
* ``μ_B`` — mean background intensity
* ``c`` — correction factor
scale : float or None
Pixel resolution (unit/px), loaded via `load_JIMG_project_()` or set manually
using `set_scale()`.
stack_selection : list of int
Indices of z-slices excluded from projection of a 3D image.
Notes
This property also prints the metadata values to the console for quick inspection.
def
set_projection(self, projection: str):
278 def set_projection(self, projection: str): 279 """ 280 Set the projection method for 3D image stack reduction. 281 282 Parameters 283 ---------- 284 projection : {"avg", "median", "std", "var", "max", "min"} 285 Projection method to reduce a 3D image stack to a 2D image. Default is `"avg"`. 286 287 Notes 288 ----- 289 This method updates the `typ` attribute of the class. The selected projection 290 determines how the z-stack is collapsed: 291 - `"avg"` : average intensity across slices 292 - `"median"` : median intensity across slices 293 - `"std"` : standard deviation across slices 294 - `"var"` : variance across slices 295 - `"max"` : maximum intensity across slices 296 - `"min"` : minimum intensity across slices 297 """ 298 299 t = ["avg", "median", "std", "var", "max", "min"] 300 if projection in t: 301 self.typ = projection 302 else: 303 print(f"\nProvided parameter is incorrect. Avaiable projection types: {t}")
Set the projection method for 3D image stack reduction.
Parameters
projection : {"avg", "median", "std", "var", "max", "min"}
Projection method to reduce a 3D image stack to a 2D image. Default is "avg".
Notes
This method updates the typ attribute of the class. The selected projection
determines how the z-stack is collapsed:
"avg": average intensity across slices"median": median intensity across slices"std": standard deviation across slices"var": variance across slices"max": maximum intensity across slices"min": minimum intensity across slices
def
set_correction_factorn(self, factor: float):
305 def set_correction_factorn(self, factor: float): 306 r""" 307 Set the correction factor for background subtraction during intensity normalization. 308 309 Parameters 310 ---------- 311 factor : float 312 Correction factor to adjust background subtraction. Must satisfy 0 < factor < 1. 313 Default is 0.1. 314 315 Notes 316 ----- 317 The correction is applied per pixel in the target mask using the formula: 318 319 .. math:: 320 321 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 322 323 where 324 * ``R_{i,j}`` — normalized pixel intensity 325 * ``T_{i,j}`` — original pixel intensity 326 * ``μ_B`` — mean intensity in the background mask 327 * ``c`` — correction factor 328 """ 329 330 if factor < 1 and factor > 0: 331 self.correction_factor = factor 332 else: 333 print( 334 "\nProvided parameter is incorrect. The factor should be a floating-point value within the range of 0 to 1." 335 )
Set the correction factor for background subtraction during intensity normalization.
Parameters
factor : float Correction factor to adjust background subtraction. Must satisfy 0 < factor < 1. Default is 0.1.
Notes
The correction is applied per pixel in the target mask using the formula:
$$R_{i,j} = T_{i,j} - ( \mu_B (1 + c) )$$
where
R_{i,j}— normalized pixel intensityT_{i,j}— original pixel intensityμ_B— mean intensity in the background maskc— correction factor
def
set_scale(self, scale):
337 def set_scale(self, scale): 338 """ 339 Set the scale for converting pixel measurements to physical units. 340 341 Parameters 342 ---------- 343 scale : float 344 Pixel resolution in physical units (e.g., µm/px). Used to calculate the 345 actual size of the tissue or organ. 346 347 Notes 348 ----- 349 The scale can also be automatically loaded from a JIMG project using 350 `load_JIMG_project_()`. This value is required for size calculations in 351 `size_calculations()`. 352 """ 353 354 self.scale = scale
Set the scale for converting pixel measurements to physical units.
Parameters
scale : float Pixel resolution in physical units (e.g., µm/px). Used to calculate the actual size of the tissue or organ.
Notes
The scale can also be automatically loaded from a JIMG project using
load_JIMG_project_(). This value is required for size calculations in
size_calculations().
def
set_selection_list(self, rm_list: list):
356 def set_selection_list(self, rm_list: list): 357 """ 358 Set the list of z-slices to exclude when projecting a 3D image stack. 359 360 Parameters 361 ---------- 362 rm_list : list of int 363 List of indices corresponding to z-slices that should be removed from 364 the full 3D image stack before projection. 365 366 Notes 367 ----- 368 This updates the `stack_selection` attribute, which is used by the 369 `stack_selection_()` method during projection. 370 """ 371 372 self.stack_selection = rm_list
Set the list of z-slices to exclude when projecting a 3D image stack.
Parameters
rm_list : list of int List of indices corresponding to z-slices that should be removed from the full 3D image stack before projection.
Notes
This updates the stack_selection attribute, which is used by the
stack_selection_() method during projection.
def
load_JIMG_project_(self, path):
374 def load_JIMG_project_(self, path): 375 """ 376 Load a JIMG project from a `.pjm` file. 377 378 Parameters 379 ---------- 380 file_path : str 381 Path to the JIMG project file. The file must have a `.pjm` extension. 382 383 Returns 384 ------- 385 project : object 386 Loaded project object containing images and metadata. 387 388 Raises 389 ------ 390 ValueError 391 If the provided file path does not point to a `.pjm` file. 392 393 Notes 394 ----- 395 The method attempts to automatically set the `scale` and `stack_selection` 396 attributes from the project metadata if available. 397 """ 398 399 path = os.path.abspath(path) 400 401 if ".pjm" in path: 402 metadata = self.load_JIMG_project(path) 403 404 try: 405 self.scale = metadata.metadata["X_resolution[um/px]"] 406 except: 407 408 try: 409 self.scale = metadata.images_dict["metadata"][0][ 410 "X_resolution[um/px]" 411 ] 412 413 except: 414 print( 415 '\nUnable to set scale on this project! Set scale using "set_scale()"' 416 ) 417 418 self.stack_selection = metadata.removal_list 419 420 else: 421 print( 422 "\nWrong path. The provided path does not point to a JIMG project (*.pjm)." 423 )
Load a JIMG project from a .pjm file.
Parameters
file_path : str
Path to the JIMG project file. The file must have a .pjm extension.
Returns
project : object Loaded project object containing images and metadata.
Raises
ValueError
If the provided file path does not point to a .pjm file.
Notes
The method attempts to automatically set the scale and stack_selection
attributes from the project metadata if available.
def
stack_selection_(self):
425 def stack_selection_(self): 426 """ 427 Remove selected z-slices from a 3D image stack based on `stack_selection`. 428 429 Notes 430 ----- 431 Only works if `input_image` is a 3D ndarray. The slices with indices listed 432 in `stack_selection` are excluded from the stack. Updates `input_image` 433 in-place. 434 435 Prints a warning if `stack_selection` is empty. 436 """ 437 438 if len(self.input_image.shape) == 3: 439 if len(self.stack_selection) > 0: 440 self.input_image = self.input_image[ 441 [ 442 x 443 for x in range(self.input_image.shape[0]) 444 if x not in self.stack_selection 445 ] 446 ] 447 else: 448 print("\nImages to remove from the stack were not selected!")
Remove selected z-slices from a 3D image stack based on stack_selection.
Notes
Only works if input_image is a 3D ndarray. The slices with indices listed
in stack_selection are excluded from the stack. Updates input_image
in-place.
Prints a warning if stack_selection is empty.
def
projection(self):
450 def projection(self): 451 """ 452 Project a 3D image stack into a 2D image using the method defined by `typ`. 453 454 Notes 455 ----- 456 Updates the `image` attribute with the projected 2D result. 457 458 Supported projection types (`typ`): 459 - "avg" : mean intensity across slices 460 - "median" : median intensity across slices 461 - "std" : standard deviation across slices 462 - "var" : variance across slices 463 - "max" : maximum intensity across slices 464 - "min" : minimum intensity across slices 465 466 Raises 467 ------ 468 AttributeError 469 If `input_image` is not defined. 470 """ 471 472 if self.typ == "avg": 473 img = np.mean(self.input_image, axis=0) 474 475 elif self.typ == "std": 476 img = np.std(self.input_image, axis=0) 477 478 elif self.typ == "median": 479 img = np.median(self.input_image, axis=0) 480 481 elif self.typ == "var": 482 img = np.var(self.input_image, axis=0) 483 484 elif self.typ == "max": 485 img = np.max(self.input_image, axis=0) 486 487 elif self.typ == "min": 488 img = np.min(self.input_image, axis=0) 489 490 self.image = img
Project a 3D image stack into a 2D image using the method defined by typ.
Notes
Updates the image attribute with the projected 2D result.
Supported projection types (typ):
- "avg" : mean intensity across slices
- "median" : median intensity across slices
- "std" : standard deviation across slices
- "var" : variance across slices
- "max" : maximum intensity across slices
- "min" : minimum intensity across slices
Raises
AttributeError
If input_image is not defined.
def
detect_img(self):
492 def detect_img(self): 493 """ 494 Detect whether the input image is 2D or 3D and perform appropriate preprocessing. 495 496 Notes 497 ----- 498 - For 3D images, applies `stack_selection_()` and then `projection()`. 499 - For 2D images, no projection is applied. 500 - Prints status messages indicating the type of image and applied operations. 501 502 Raises 503 ------ 504 AttributeError 505 If `input_image` is not defined. 506 """ 507 check = len(self.input_image.shape) 508 509 if check == 3: 510 print("\n3D image detected! Starting processing for 3D image...") 511 print(f"Projection - {self.typ}...") 512 513 self.stack_selection_() 514 self.projection() 515 516 elif check == 2: 517 print("\n2D image detected! Starting processing for 2D image...") 518 519 else: 520 print("\nData does not match any image type!")
Detect whether the input image is 2D or 3D and perform appropriate preprocessing.
Notes
- For 3D images, applies
stack_selection_()and thenprojection(). - For 2D images, no projection is applied.
- Prints status messages indicating the type of image and applied operations.
Raises
AttributeError
If input_image is not defined.
def
load_image_3D(self, path):
522 def load_image_3D(self, path): 523 """ 524 Load a 3D image stack from a TIFF file. 525 526 Parameters 527 ---------- 528 path : str 529 Path to the 3D image file (*.tiff) to be loaded. 530 531 Notes 532 ----- 533 The loaded image is stored in the `input_image` attribute as a 3D ndarray. 534 """ 535 path = os.path.abspath(path) 536 537 self.input_image = self.load_3D_tiff(path)
Load a 3D image stack from a TIFF file.
Parameters
path : str Path to the 3D image file (*.tiff) to be loaded.
Notes
The loaded image is stored in the input_image attribute as a 3D ndarray.
def
load_image_(self, path):
539 def load_image_(self, path): 540 """ 541 Load a 2D image into the class. 542 543 Parameters 544 ---------- 545 path : str 546 Path to the image file to be loaded. 547 548 Notes 549 ----- 550 The loaded image is stored in the `input_image` attribute as a 2D ndarray. 551 """ 552 path = os.path.abspath(path) 553 554 self.input_image = self.load_image(path)
Load a 2D image into the class.
Parameters
path : str Path to the image file to be loaded.
Notes
The loaded image is stored in the input_image attribute as a 2D ndarray.
def
load_mask_(self, path):
556 def load_mask_(self, path): 557 r""" 558 Load a binary mask into the class and optionally set it as the normalization mask. 559 560 Parameters 561 ---------- 562 path : str 563 Path to the mask image file. Supported formats include 8-bit or 16-bit images 564 with extensions such as `.png` or `.jpeg`. The mask must be binary 565 (e.g., 0/255, 0/2**16-1, 0/1). 566 567 Notes 568 ----- 569 - If `load_normalization_mask_()` is not called, this mask is also used as the 570 background mask for intensity normalization. 571 - Normalization is applied per pixel using the formula: 572 573 .. math:: 574 575 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 576 577 where 578 * ``R_{i,j}`` — normalized pixel intensity 579 * ``T_{i,j}`` — pixel intensity in the target mask 580 * ``μ_B`` — mean intensity of the background (reversed mask) 581 * ``c`` — correction factor 582 """ 583 584 path = os.path.abspath(path) 585 586 self.mask = self.load_mask(path) 587 588 print( 589 "\nThis mask was also set as the reverse background mask. If you want a different background mask for normalization, use 'load_normalization_mask()'." 590 ) 591 self.background_mask = self.load_mask(path)
Load a binary mask into the class and optionally set it as the normalization mask.
Parameters
path : str
Path to the mask image file. Supported formats include 8-bit or 16-bit images
with extensions such as .png or .jpeg. The mask must be binary
(e.g., 0/255, 0/2**16-1, 0/1).
Notes
- If
load_normalization_mask_()is not called, this mask is also used as the background mask for intensity normalization. Normalization is applied per pixel using the formula:
$$R_{i,j} = T_{i,j} - ( \mu_B (1 + c) )$$
where
R_{i,j}— normalized pixel intensityT_{i,j}— pixel intensity in the target maskμ_B— mean intensity of the background (reversed mask)c— correction factor
def
load_normalization_mask_(self, path):
593 def load_normalization_mask_(self, path): 594 r""" 595 Load a binary mask for normalization into the class. 596 597 Parameters 598 ---------- 599 path : str 600 Path to the mask image file. Supported formats include 8-bit or 16-bit 601 images (e.g., `.png`, `.jpeg`). The mask must be binary (0/255, 0/2**16-1, 0/1). 602 603 Notes 604 ----- 605 - The mask defines the area of interest. Normalization is applied to the inverse 606 of this area (reversed mask). 607 - Normalization formula applied per pixel: 608 609 .. math:: 610 611 R_{i,j} = T_{i,j} - ( \mu_B (1 + c) ) 612 613 where 614 * ``R_{i,j}`` — normalized pixel intensity 615 * ``T_{i,j}`` — pixel intensity in the target mask 616 * ``μ_B`` — mean intensity of the background (reversed mask) 617 * ``c`` — correction factor 618 """ 619 620 path = os.path.abspath(path) 621 622 self.background_mask = self.load_mask(path)
Load a binary mask for normalization into the class.
Parameters
path : str
Path to the mask image file. Supported formats include 8-bit or 16-bit
images (e.g., .png, .jpeg). The mask must be binary (0/255, 0/2**16-1, 0/1).
Notes
- The mask defines the area of interest. Normalization is applied to the inverse of this area (reversed mask).
Normalization formula applied per pixel:
$$R_{i,j} = T_{i,j} - ( \mu_B (1 + c) )$$
where
R_{i,j}— normalized pixel intensityT_{i,j}— pixel intensity in the target maskμ_B— mean intensity of the background (reversed mask)c— correction factor
def
intensity_calculations(self):
624 def intensity_calculations(self): 625 """ 626 Calculate normalized and raw intensity statistics from the image based on masks. 627 628 This method performs intensity calculations using the main mask (`self.mask`) 629 and the background mask (`self.background_mask`). The pixel intensities within 630 the mask of interest are normalized by subtracting a threshold derived from the 631 background region and applying a correction factor (`self.correction_factor`). 632 Negative values after normalization are clipped to zero. 633 634 The following statistics are computed for both normalized and raw values: 635 - Minimum 636 - Maximum 637 - Mean 638 - Median 639 - Standard deviation 640 - Variance 641 - List of all normalized values (only for normalized data) 642 643 Notes 644 ----- 645 - The method updates the instance attribute `self.normalized_image_values` 646 with a dictionary containing both normalized and raw statistics. 647 - Normalization formula applied for each pixel in the selected mask: 648 final_val = selected_value - (threshold + threshold * correction_factor) 649 where threshold is the mean intensity in the background mask. 650 - Negative values after normalization are set to zero. 651 """ 652 653 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 654 tmp_bmask = self.ajd_mask_size(image=self.image, mask=self.background_mask) 655 656 selected_values = self.image[tmp_mask == np.max(tmp_mask)] 657 658 threshold = np.mean(self.image[tmp_bmask == np.min(tmp_bmask)]) 659 660 # normalization 661 final_val = selected_values - (threshold + (threshold * self.correction_factor)) 662 663 final_val[final_val < 0] = 0 664 665 tmp_dict = { 666 "norm_min": np.min(final_val), 667 "norm_max": np.max(final_val), 668 "norm_mean": np.mean(final_val), 669 "norm_median": np.median(final_val), 670 "norm_std": np.std(final_val), 671 "norm_var": np.var(final_val), 672 "norm_values": final_val.tolist(), 673 "min": np.min(selected_values), 674 "max": np.max(selected_values), 675 "mean": np.mean(selected_values), 676 "median": np.median(selected_values), 677 "std": np.std(selected_values), 678 "var": np.var(selected_values), 679 } 680 681 self.normalized_image_values = tmp_dict
Calculate normalized and raw intensity statistics from the image based on masks.
This method performs intensity calculations using the main mask (self.mask)
and the background mask (self.background_mask). The pixel intensities within
the mask of interest are normalized by subtracting a threshold derived from the
background region and applying a correction factor (self.correction_factor).
Negative values after normalization are clipped to zero.
The following statistics are computed for both normalized and raw values:
- Minimum
- Maximum
- Mean
- Median
- Standard deviation
- Variance
- List of all normalized values (only for normalized data)
Notes
- The method updates the instance attribute
self.normalized_image_valueswith a dictionary containing both normalized and raw statistics. - Normalization formula applied for each pixel in the selected mask: final_val = selected_value - (threshold + threshold * correction_factor) where threshold is the mean intensity in the background mask.
- Negative values after normalization are set to zero.
def
size_calculations(self):
683 def size_calculations(self): 684 """ 685 Calculates the size and bounding dimensions of the masked region in the image. 686 687 This method computes the following metrics based on the current mask: 688 - Total number of pixels in the mask (`px_size`) 689 - Real-world size if a scale is provided (`size`) 690 - Maximum lengths along x and y axes (`max_length_x_axis`, `max_length_y_axis`) 691 692 If `self.scale` is defined (unit per pixel), the real-world size is calculated. 693 If not, `size` will be `None` and a warning message is printed. 694 695 Returns: 696 Updates the following attributes in the class: 697 - self.size_info (dict) containing: 698 - 'size' (float or None): real-world size of the mask 699 - 'px_size' (int): number of pixels in the masked region 700 - 'max_length_x_axis' (int): length of the bounding box along the x-axis 701 - 'max_length_y_axis' (int): length of the bounding box along the y-axis 702 703 Example: 704 analysis.size_calculations() 705 print(analysis.size_info) 706 """ 707 708 tmp_mask = self.ajd_mask_size(image=self.image, mask=self.mask) 709 710 size_px = int(len(tmp_mask[tmp_mask > np.min(tmp_mask)])) 711 712 if self.scale is not None: 713 size = float(size_px * self.scale) 714 else: 715 size = None 716 print( 717 '\nUnable to calculate real size, scale (unit/px) not provided, use "set_scale()" or load JIMG project .pjm metadata "load_pjm()" to set scale for calculations!' 718 ) 719 720 non_zero_indices = np.where(tmp_mask == np.max(tmp_mask)) 721 722 min_y, max_y = np.min(non_zero_indices[0]), np.max(non_zero_indices[0]) 723 min_x, max_x = np.min(non_zero_indices[1]), np.max(non_zero_indices[1]) 724 725 max_length_x_axis = int(max_x - min_x + 1) 726 max_length_y_axis = int(max_y - min_y + 1) 727 728 tmp_val = { 729 "size": size, 730 "px_size": size_px, 731 "max_length_x_axis": max_length_x_axis, 732 "max_length_y_axis": max_length_y_axis, 733 } 734 735 self.size_info = tmp_val
Calculates the size and bounding dimensions of the masked region in the image.
This method computes the following metrics based on the current mask:
- Total number of pixels in the mask (px_size)
- Real-world size if a scale is provided (size)
- Maximum lengths along x and y axes (max_length_x_axis, max_length_y_axis)
If self.scale is defined (unit per pixel), the real-world size is calculated.
If not, size will be None and a warning message is printed.
Returns: Updates the following attributes in the class: - self.size_info (dict) containing: - 'size' (float or None): real-world size of the mask - 'px_size' (int): number of pixels in the masked region - 'max_length_x_axis' (int): length of the bounding box along the x-axis - 'max_length_y_axis' (int): length of the bounding box along the y-axis
Example: analysis.size_calculations() print(analysis.size_info)
def
run_calculations(self):
737 def run_calculations(self): 738 """ 739 Run the full analysis pipeline on the loaded image using the provided masks. 740 741 Notes 742 ----- 743 - The input image must be loaded via `load_image_()` or `load_image_3D()`. 744 - The ROI mask must be loaded via `load_mask_()`. Optionally, a normalization 745 mask can be loaded via `load_normalization_mask_()`. 746 - Parameters such as projection type and correction factor can be set with 747 `set_projection()` and `set_correction_factor()`. 748 - Scale and stack selection can also influence calculations if defined. 749 - To view current parameters, use the `current_metadata` property. 750 751 Returns 752 ------- 753 None 754 The results are stored internally and can be retrieved using 755 `get_results()`. 756 """ 757 758 if self.input_image is not None: 759 760 if self.mask is not None: 761 762 print("\nStart...") 763 self.detect_img() 764 self.intensity_calculations() 765 self.size_calculations() 766 print("\nCompleted!")
Run the full analysis pipeline on the loaded image using the provided masks.
Notes
- The input image must be loaded via
load_image_()orload_image_3D(). - The ROI mask must be loaded via
load_mask_(). Optionally, a normalization mask can be loaded viaload_normalization_mask_(). - Parameters such as projection type and correction factor can be set with
set_projection()andset_correction_factor(). - Scale and stack selection can also influence calculations if defined.
- To view current parameters, use the
current_metadataproperty.
Returns
None
The results are stored internally and can be retrieved using
get_results().
def
get_results(self):
768 def get_results(self): 769 """ 770 Return the results from the analysis performed by `run_calculations()`. 771 772 Returns 773 ------- 774 results_dict : dict or None 775 Dictionary containing intensity and size results. Structure: 776 - 'intensity' : dict with normalized and raw intensity statistics 777 - 'size' : dict with ROI size metrics 778 779 Notes 780 ----- 781 If analysis has not been run yet, prints a message and returns None. 782 """ 783 784 if self.normalized_image_values is not None and self.size_info is not None: 785 786 results = { 787 "intensity": self.normalized_image_values, 788 "size": self.size_info, 789 } 790 791 return results 792 793 else: 794 print('\nAnalysis were not conducted. Run analysis "run_calculations()"')
Return the results from the analysis performed by run_calculations().
Returns
results_dict : dict or None Dictionary containing intensity and size results. Structure: - 'intensity' : dict with normalized and raw intensity statistics - 'size' : dict with ROI size metrics
Notes
If analysis has not been run yet, prints a message and returns None.
def
save_results( self, path='', mask_region: str = '', feature_name: str = '', individual_number: int = 0, individual_name: str = ''):
796 def save_results( 797 self, 798 path="", 799 mask_region: str = "", 800 feature_name: str = "", 801 individual_number: int = 0, 802 individual_name: str = "", 803 ): 804 """ 805 Save the analysis results to a `.int` (JSON) file. 806 807 Parameters 808 ---------- 809 path : str, optional 810 Directory path where the file will be saved. Defaults to the current working directory. 811 812 mask_region : str 813 Name or identifier of the mask region (e.g., tissue, part of tissue). 814 815 feature_name : str 816 Name of the feature being analyzed. Underscores or spaces are replaced with periods. 817 818 individual_number : int 819 Unique identifier for the individual in the analysis (e.g., 1, 2, 3). 820 821 individual_name : str 822 Name of the individual (e.g., species name, tissue, organoid). 823 824 Notes 825 ----- 826 - The method validates that all required parameters are provided and that 827 analysis results exist (`normalized_image_values` and `size_info`). 828 - Creates the directory if it does not exist. 829 - File name format: 830 '<individual_name>_<individual_number>_<mask_region>_<feature_name>.int' 831 832 Raises 833 ------ 834 FileNotFoundError 835 If the specified path cannot be created or accessed. 836 837 ValueError 838 If any of `mask_region`, `feature_name`, `individual_number`, or 839 `individual_name` are missing or invalid. 840 """ 841 842 path = os.path.abspath(path) 843 844 if ( 845 len(mask_region) > 1 846 and len(feature_name) > 1 847 and individual_number != 0 848 and len(individual_name) > 1 849 ): 850 851 if self.normalized_image_values is not None and self.size_info is not None: 852 853 results = { 854 "intensity": self.normalized_image_values, 855 "size": self.size_info, 856 } 857 858 mask_region = re.sub(r"[_\s]+", ".", mask_region) 859 feature_name = re.sub(r"[_\s]+", ".", feature_name) 860 individual_number = re.sub(r"[_\s]+", ".", str(individual_number)) 861 individual_name = re.sub(r"[_\s]+", ".", individual_name) 862 863 full_name = f"{individual_name}_{individual_number}_{mask_region}_{feature_name}" 864 865 isExist = os.path.exists(path) 866 if not isExist: 867 os.makedirs(path, exist_ok=True) 868 869 full_path = os.path.join( 870 path, re.sub("\\.json", "", full_name) + ".int" 871 ) 872 873 with open(full_path, "w") as file: 874 json.dump(results, file, indent=4) 875 876 else: 877 print( 878 '\nAnalysis were not conducted. Run analysis "run_calculations()"' 879 ) 880 881 else: 882 print( 883 "\nAny of 'mask_region', 'feature_name', 'individual_number', 'individual_name' parameters were provided wrong!" 884 )
Save the analysis results to a .int (JSON) file.
Parameters
path : str, optional Directory path where the file will be saved. Defaults to the current working directory.
mask_region : str Name or identifier of the mask region (e.g., tissue, part of tissue).
feature_name : str Name of the feature being analyzed. Underscores or spaces are replaced with periods.
individual_number : int Unique identifier for the individual in the analysis (e.g., 1, 2, 3).
individual_name : str Name of the individual (e.g., species name, tissue, organoid).
Notes
- The method validates that all required parameters are provided and that
analysis results exist (
normalized_image_valuesandsize_info). - Creates the directory if it does not exist.
- File name format:
'
_ _ _ .int'
Raises
FileNotFoundError If the specified path cannot be created or accessed.
ValueError
If any of mask_region, feature_name, individual_number, or
individual_name are missing or invalid.
def
concatenate_intensity_data(self, directory: str = '', name: str = ''):
886 def concatenate_intensity_data(self, directory: str = "", name: str = ""): 887 """ 888 Concatenate intensity data from multiple `.int` files and save as CSV. 889 890 Parameters 891 ---------- 892 directory : str, optional 893 Path to the directory containing `.int` files. Defaults to the current working directory. 894 895 name : str 896 Prefix for the output CSV file names. CSV files are saved in the format 897 '<name>_<gene>_<region>.csv'. 898 899 Raises 900 ------ 901 FileNotFoundError 902 If the directory cannot be accessed or no `.int` files are found. 903 904 ValueError 905 If an `.int` file is missing expected data or has an incorrect format. 906 907 Notes 908 ----- 909 - The method groups intensity data by gene (feature) and mask region. 910 - Outputs one CSV file per unique gene-region combination, saved in the specified directory. 911 """ 912 913 directory = os.path.abspath(directory) 914 915 files_list = [f for f in os.listdir(directory) if f.endswith(".int")] 916 917 genes_set = set([re.sub("\\.int", "", x.split("_")[3]) for x in files_list]) 918 regions_set = set([re.sub("\\.int", "", x.split("_")[2]) for x in files_list]) 919 920 for g in genes_set: 921 for r in regions_set: 922 json_to_save = { 923 "individual_name": [], 924 "individual_number": [], 925 "norm_intensity": [], 926 "size": [], 927 } 928 929 for f in tqdm(files_list): 930 if g in f and r in f: 931 with open(os.path.join(directory, f), "r") as file: 932 data = json.load(file) 933 934 json_to_save["norm_intensity"] = ( 935 json_to_save["norm_intensity"] 936 + data["intensity"]["norm_values"] 937 ) 938 json_to_save["individual_name"] = json_to_save[ 939 "individual_name" 940 ] + [f.split("_")[0]] * len( 941 data["intensity"]["norm_values"] 942 ) 943 json_to_save["individual_number"] = json_to_save[ 944 "individual_number" 945 ] + [f.split("_")[1]] * len( 946 data["intensity"]["norm_values"] 947 ) 948 json_to_save["size"] = json_to_save["size"] + [ 949 data["size"]["px_size"] 950 ] * len(data["intensity"]["norm_values"]) 951 952 pd.DataFrame(json_to_save).to_csv(f"{name}_{g}_{r}.csv", index=False)
Concatenate intensity data from multiple .int files and save as CSV.
Parameters
directory : str, optional
Path to the directory containing .int files. Defaults to the current working directory.
name : str
Prefix for the output CSV file names. CSV files are saved in the format
'
Raises
FileNotFoundError
If the directory cannot be accessed or no .int files are found.
ValueError
If an .int file is missing expected data or has an incorrect format.
Notes
- The method groups intensity data by gene (feature) and mask region.
- Outputs one CSV file per unique gene-region combination, saved in the specified directory.
class
IntensityAnalysis:
955class IntensityAnalysis: 956 """ 957 Class for performing percentile-based statistical analysis on grouped data. 958 959 This class provides methods to calculate percentiles, remove outliers, aggregate 960 data into percentile bins, perform Welch's ANOVA, Kolmogorov-Smirnov and Fisher's exact tests, 961 evaluate Wasserstein distances, calculate Fold Change, and visualize results via comparative 962 histograms. It is designed to handle both single-column and multi-column combinations 963 of values for group-based analysis. 964 965 Methods 966 ------- 967 drop_up_df(data, group_col, values_col) 968 Removes upper outliers from a DataFrame based on a grouping column. 969 970 percentiles_calculation(values, sep_perc=1) 971 Calculates percentiles and creates loopable percentile ranges. 972 973 to_percentil(values, percentiles, percentiles_loop, values_col, replication_col) 974 Aggregates statistics based on percentile ranges. 975 976 df_to_percentiles(data, group_col, values_col, replication_col, sep_perc=1, drop_outlires=True) 977 Computes percentile statistics for grouped DataFrame data. 978 979 round_to_scientific_notation(num) 980 Formats a number in scientific notation or standard format. 981 982 aov(data, testes_col, comb="*") 983 Performs Welch's ANOVA on percentile-based group data. 984 985 post_aov(data, testes_col, comb="*") 986 Performs Welch's ANOVA with pairwise t-tests. 987 988 ks_percentiles(input_hist) 989 Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups. 990 991 fisher_percentiles(input_hist) 992 Perform pairwise Fisher's exact tests on percentile data across all groups. 993 994 to_wasserstein_distance(data) 995 Calculates scaled pairwise Wasserstein distances for grouped distributions. 996 997 to_fold_change(data, tested_value) 998 Calculates the Fold Change (FC) between all directed permutations of groups. 999 1000 get_stats(data, tested_value) 1001 Calculates and aggregates overall statistical metrics (ANOVA, Fisher's exact, Kolmogorov-Smirnov, FC, Wasserstein distance). 1002 1003 hist_compare_plot(data, queue=None, p_adj=True, txt_size=20) 1004 Generates comparative histograms with statistical test results and metrics. 1005 """ 1006 1007 def drop_up_df(self, data: pd.DataFrame, group_col: str, values_col: str): 1008 """ 1009 Remove upper outliers from a DataFrame based on a specified value column, grouped by a grouping column. 1010 1011 Outliers are calculated and removed separately for each group defined by `group_col`. 1012 The upper outliers are defined using the interquartile range (IQR) method: 1013 values greater than Q3 + 1.5 * IQR are considered outliers. 1014 1015 Parameters 1016 ---------- 1017 data : pd.DataFrame 1018 The input DataFrame containing the data. 1019 1020 group_col : str 1021 The name of the column used for grouping the data. 1022 1023 values_col : str 1024 The column containing the values from which upper outliers will be removed. 1025 1026 Returns 1027 ------- 1028 filtered_data : pd.DataFrame 1029 A filtered DataFrame with the upper outliers removed for each group. 1030 1031 Notes 1032 ----- 1033 - Outliers are removed separately within each group. 1034 - The original DataFrame is not modified; a new filtered DataFrame is returned. 1035 """ 1036 1037 def iqr_filter(group): 1038 q75 = np.quantile(group[values_col], 0.75) 1039 q25 = np.quantile(group[values_col], 0.25) 1040 itq = q75 - q25 1041 return group[group[values_col] <= (q75 + 1.5 * itq)] 1042 1043 filtered_data = data.groupby(group_col).apply(iqr_filter).reset_index(drop=True) 1044 1045 return filtered_data 1046 1047 def percentiles_calculation(self, values, sep_perc: int = 1): 1048 """ 1049 Calculate percentiles for a set of values and generate consecutive percentile ranges. 1050 1051 This function computes percentiles from 0 to 100 at intervals defined by `sep_perc`. 1052 It also generates a list of consecutive percentile ranges that can be used for further analysis or binning. 1053 1054 Parameters 1055 ---------- 1056 values : array-like 1057 The input data values for which the percentiles are calculated. 1058 1059 sep_perc : int, optional 1060 Separation interval between percentiles (default is 1, meaning percentiles are calculated every 1%). 1061 1062 Returns 1063 ------- 1064 percentiles : np.ndarray 1065 Array of calculated percentile values. 1066 1067 percentiles_loop : list of tuple 1068 List of consecutive percentile ranges as tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1069 1070 Notes 1071 ----- 1072 - The first percentile is set to 0 to avoid issues with zero values. 1073 - `percentiles_loop` is useful for iterating through percentile ranges when aggregating statistics. 1074 """ 1075 1076 per_vector = values.copy() 1077 1078 percentiles = np.percentile(per_vector, np.arange(0, 101, sep_perc)) 1079 percentiles[0] = 0 1080 1081 percentiles_loop = [(i, i + 1) for i in range(int(100 / sep_perc))] 1082 1083 return percentiles, percentiles_loop 1084 1085 def to_percentil( 1086 self, values, percentiles, percentiles_loop, values_col, replication_col 1087 ): 1088 """ 1089 Aggregate statistics for a set of values based on percentile ranges, including replications. 1090 1091 This function calculates summary statistics (count, proportion, mean, median, 1092 standard deviation, variance) for each percentile range defined in `percentiles_loop`. 1093 It computes these statistics both for the combined data ('mutual') and separately 1094 for each individual replication. It also calculates overall metrics per replication. 1095 1096 Parameters 1097 ---------- 1098 values : pd.DataFrame 1099 Input DataFrame containing the data to be analyzed. 1100 percentiles : np.ndarray 1101 Array of percentile values used to define the boundaries of each range. 1102 percentiles_loop : list of tuple 1103 List of consecutive percentile ranges as index tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1104 values_col : str 1105 The column name in `values` containing the numeric data to aggregate. 1106 replication_col : str 1107 The column name in `values` used to identify distinct replications or samples. 1108 1109 Returns 1110 ------- 1111 full_data : dict 1112 A nested dictionary containing the calculated statistics with the following structure: 1113 - 'percentiles' : dict 1114 - 'mutual' : dict 1115 Lists of statistics ('n', 'n_standarized', 'avg', 'median', 'std', 'var') 1116 aggregated across all replications for each percentile bin. 1117 - 'replications' : dict 1118 Keys are replication names. Values are dictionaries of statistics (same as above) 1119 calculated specifically for that replication within each bin. 1120 - 'values' : dict 1121 Lists of overall statistics ('avg', 'median', 'std', 'var', 'replication') 1122 calculated for each replication as a whole (ignoring bins). 1123 1124 Notes 1125 ----- 1126 - If a percentile range contains no elements, statistics are set to 0 and count is set to 1 to avoid empty lists. 1127 """ 1128 1129 full_data = {} 1130 per_vector = values[values_col] 1131 amount = len(per_vector) 1132 1133 data_mutual = { 1134 "n": [], 1135 "n_standarized": [], 1136 "avg": [], 1137 "median": [], 1138 "std": [], 1139 "var": [], 1140 } 1141 1142 for x in percentiles_loop: 1143 subset = per_vector[ 1144 (per_vector > percentiles[x[0]]) & (per_vector <= percentiles[x[1]]) 1145 ] 1146 n_subset = len(subset) 1147 1148 if n_subset > 0: 1149 data_mutual["n"].append(n_subset) 1150 data_mutual["n_standarized"].append(n_subset / amount) 1151 data_mutual["avg"].append(np.mean(subset)) 1152 data_mutual["median"].append(np.median(subset)) 1153 data_mutual["std"].append(np.std(subset)) 1154 data_mutual["var"].append(np.var(subset)) 1155 else: 1156 data_mutual["n"].append(0) 1157 data_mutual["n_standarized"].append(0) 1158 data_mutual["avg"].append(0) 1159 data_mutual["median"].append(0) 1160 data_mutual["std"].append(0) 1161 data_mutual["var"].append(0) 1162 1163 full_data["percentiles"] = {"mutual": data_mutual, "replications": {}} 1164 1165 unique_names = set(values[replication_col]) 1166 1167 for nam in unique_names: 1168 per_vector_rep = values[values_col][values[replication_col] == nam] 1169 1170 data_rep = { 1171 "n": [], 1172 "n_standarized": [], 1173 "avg": [], 1174 "median": [], 1175 "std": [], 1176 "var": [], 1177 } 1178 1179 for x in percentiles_loop: 1180 subset = per_vector_rep[ 1181 (per_vector_rep > percentiles[x[0]]) 1182 & (per_vector_rep <= percentiles[x[1]]) 1183 ] 1184 n_subset = len(subset) 1185 1186 if n_subset > 0: 1187 data_rep["n"].append(n_subset) 1188 data_rep["n_standarized"].append(n_subset / amount) 1189 data_rep["avg"].append(np.mean(subset)) 1190 data_rep["median"].append(np.median(subset)) 1191 data_rep["std"].append(np.std(subset)) 1192 data_rep["var"].append(np.var(subset)) 1193 else: 1194 data_rep["n"].append(0) 1195 data_rep["n_standarized"].append(0) 1196 data_rep["avg"].append(0) 1197 data_rep["median"].append(0) 1198 data_rep["std"].append(0) 1199 data_rep["var"].append(0) 1200 1201 full_data["percentiles"]["replications"][nam] = data_rep 1202 1203 unique_names = set(values[replication_col]) 1204 1205 data_rep = {"avg": [], "median": [], "std": [], "var": [], "replication": []} 1206 1207 for nam in unique_names: 1208 per_vector_rep = values[values_col][values[replication_col] == nam] 1209 1210 data_rep["avg"].append(np.mean(per_vector_rep)) 1211 data_rep["median"].append(np.median(per_vector_rep)) 1212 data_rep["std"].append(np.std(per_vector_rep)) 1213 data_rep["var"].append(np.var(per_vector_rep)) 1214 data_rep["replication"].append(nam) 1215 1216 full_data["values"] = data_rep 1217 1218 return full_data 1219 1220 def df_to_percentiles( 1221 self, 1222 data: pd.DataFrame, 1223 group_col: str = "individual_name", 1224 values_col: str = "norm_intensity", 1225 replication_col: str = "individual_number", 1226 sep_perc: int = 1, 1227 drop_outlires: bool = True, 1228 ): 1229 """ 1230 Calculate summary statistics based on percentile ranges for each group in a DataFrame. 1231 1232 This method groups the input DataFrame by `group_col`, computes global percentile ranges 1233 based on `values_col`, and then aggregates statistics for each percentile bin. The aggregation 1234 is performed both mutually for the group and individually per replication. Optionally, 1235 upper outliers can be removed before the calculations. 1236 1237 Parameters 1238 ---------- 1239 data : pd.DataFrame 1240 Input DataFrame containing the grouped data. 1241 group_col : str, optional 1242 Column name used to define groups (default is 'individual_name'). 1243 values_col : str, optional 1244 Column name containing the numeric values for percentile calculations 1245 (default is 'norm_intensity'). 1246 replication_col : str, optional 1247 Column name used to identify separate replications within the groups 1248 (default is 'individual_number'). 1249 sep_perc : int, optional 1250 Separation interval for percentiles (default is 1, meaning 1% steps). 1251 drop_outlires : bool, optional 1252 If True, removes upper outliers from the data using the IQR method before 1253 performing calculations (default is True). 1254 1255 Returns 1256 ------- 1257 full_data : dict 1258 A dictionary where each key is a unique group name (from `group_col`). 1259 The corresponding value is the nested dictionary returned by `to_percentil()`, 1260 which includes bin-wise statistics ('mutual' and 'replications') and overall 1261 metrics ('values'). 1262 1263 Notes 1264 ----- 1265 - Outlier removal uses the IQR method within each group if `drop_outlires` is True. 1266 """ 1267 1268 full_data = {} 1269 1270 if drop_outlires == True: 1271 data = self.drop_up_df( 1272 data=data, group_col=group_col, values_col=values_col 1273 ) 1274 1275 groups = set(data[group_col]) 1276 val_dat = [x for x in data[values_col] if x > 0] 1277 1278 percentiles, percentiles_loop = self.percentiles_calculation( 1279 val_dat, sep_perc=sep_perc 1280 ) 1281 1282 for g in groups: 1283 1284 print(f"Group: {g} ...") 1285 1286 tmp_values = data[data[group_col] == g] 1287 1288 per_dat = self.to_percentil( 1289 tmp_values, percentiles, percentiles_loop, values_col, replication_col 1290 ) 1291 1292 full_data[g] = per_dat 1293 1294 return full_data 1295 1296 def round_to_scientific_notation(self, num): 1297 """ 1298 Round a number to scientific notation if very small, otherwise to one decimal place. 1299 1300 Parameters 1301 ---------- 1302 num : float 1303 The number to round. 1304 1305 Returns 1306 ------- 1307 str 1308 The rounded number as a string. 1309 - If `num` is 0, returns "0.0". 1310 - If `abs(num) < 1e-4`, returns scientific notation with 1 decimal and 1-digit exponent. 1311 - Otherwise, returns the number rounded to one decimal place. 1312 """ 1313 1314 if num == 0: 1315 return "0.0" 1316 1317 if abs(num) < 0.0001: 1318 rounded_num = np.format_float_scientific(num, precision=1, exp_digits=1) 1319 return rounded_num 1320 else: 1321 return f"{num:.1f}" 1322 1323 def aov(self, data, testes_col, comb: str = "*"): 1324 """ 1325 Perform a Welch's ANOVA analysis. 1326 1327 This function calculates group values by aggregating specified columns (testes_col) 1328 via the comb method and then conducts a Welch's ANOVA. This approach is ideal for 1329 comparing group means when data exhibits unequal variances across groups. 1330 1331 Parameters 1332 ---------- 1333 data : dict of pd.DataFrame 1334 Dictionary where keys are group names and values are DataFrames containing the data. 1335 1336 testes_col : str or list of str 1337 Column name(s) from which the group values are derived. If a list is provided, columns 1338 will be combined based on the `comb` operation. 1339 1340 comb : str, optional 1341 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1342 '*' : multiplication 1343 '+' : addition 1344 '**': exponentiation 1345 '-' : subtraction 1346 '/' : division 1347 Default is '*'. 1348 1349 Returns 1350 ------- 1351 F : float 1352 F-statistic from Welch's ANOVA. 1353 1354 p_val : float 1355 Uncorrected p-value from Welch's ANOVA, testing for significant differences between groups. 1356 1357 Notes 1358 ----- 1359 - If `testes_col` is a single string, no combination is performed, and the group values 1360 are taken directly from that column. 1361 - Welch's ANOVA is used as it accounts for unequal variances between groups. 1362 - The `df.melt()` method is used to reshape the data, allowing the ANOVA to be applied to all groups. 1363 1364 Examples 1365 -------- 1366 >>> welch_F, welch_p = self.aov(data, testes_col=['col1', 'col2'], comb='+') 1367 >>> print(f"Welch's ANOVA F-statistic: {welch_F}, p-value: {welch_p}") 1368 """ 1369 1370 groups = [] 1371 1372 for d in data.keys(): 1373 1374 if isinstance(testes_col, str): 1375 g = data[d]["values"][testes_col] 1376 elif isinstance(testes_col, list): 1377 g = [1] * len(data[d]["values"][testes_col[0]]) 1378 for t in testes_col: 1379 if comb == "*": 1380 g = [a * b for a, b in zip(g, data[d]["values"][t])] 1381 elif comb == "+": 1382 g = [a + b for a, b in zip(g, data[d]["values"][t])] 1383 elif comb == "**": 1384 g = [a**b for a, b in zip(g, data[d]["values"][t])] 1385 elif comb == "-": 1386 g = [a - b for a, b in zip(g, data[d]["values"][t])] 1387 elif comb == "/": 1388 g = [a / b for a, b in zip(g, data[d]["values"][t])] 1389 1390 groups.append(g) 1391 1392 df = pd.DataFrame({f"group_{i}": group for i, group in enumerate(groups)}) 1393 1394 df_melted = df.melt(var_name="group", value_name="value") 1395 1396 welch_results = pg.welch_anova(data=df_melted, dv="value", between="group") 1397 1398 return welch_results["F"].values[0], welch_results["p-unc"].values[0] 1399 1400 def post_aov(self, data, testes_col, comb: str = "*"): 1401 """ 1402 Perform Welch's ANOVA and pairwise Welch's t-tests on grouped data. 1403 1404 This method first conducts a Welch's ANOVA to detect significant differences 1405 in group means. It then performs pairwise Welch's t-tests across all group 1406 combinations to identify specific differences. All p-values are adjusted using 1407 the Bonferroni correction to account for multiple comparisons. 1408 1409 Parameters 1410 ---------- 1411 data : dict of pd.DataFrame 1412 Dictionary where keys are group names and values are DataFrames containing the data. 1413 1414 testes_col : str or list of str 1415 Column name(s) from which the group values are derived. If a list is provided, 1416 columns will be combined according to the `comb` operation. 1417 1418 comb : str, optional 1419 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1420 '*' : multiplication 1421 '+' : addition 1422 '**': exponentiation 1423 '-' : subtraction 1424 '/' : division 1425 Default is '*'. 1426 1427 Returns 1428 ------- 1429 p_val : float 1430 Uncorrected p-value from the Welch's ANOVA. 1431 1432 final_results : dict 1433 Dictionary containing results of pairwise Welch's t-tests with keys: 1434 'group1' : list of first group names in each comparison 1435 'group2' : list of second group names in each comparison 1436 'stat' : list of t-statistics for each comparison 1437 'p_val' : list of uncorrected p-values for each comparison 1438 'adj_p_val' : list of Bonferroni-adjusted p-values for multiple comparisons 1439 """ 1440 1441 p_val = self.aov(data=data, testes_col=testes_col, comb=comb)[1] 1442 1443 pairs = list(combinations(data, 2)) 1444 final_results = { 1445 "group1": [], 1446 "group2": [], 1447 "stat": [], 1448 "p_val": [], 1449 "adj_p_val": [], 1450 } 1451 1452 for group1, group2 in pairs: 1453 if isinstance(testes_col, str): 1454 g1 = data[group1]["values"][testes_col] 1455 elif isinstance(testes_col, list): 1456 g1 = [1] * len(data[group1]["values"][testes_col[0]]) 1457 for t in testes_col: 1458 if comb == "*": 1459 g1 = [a * b for a, b in zip(g1, data[group1]["values"][t])] 1460 elif comb == "+": 1461 g1 = [a + b for a, b in zip(g1, data[group1]["values"][t])] 1462 elif comb == "**": 1463 g1 = [a**b for a, b in zip(g1, data[group1]["values"][t])] 1464 elif comb == "-": 1465 g1 = [a - b for a, b in zip(g1, data[group1]["values"][t])] 1466 elif comb == "/": 1467 g1 = [a / b for a, b in zip(g1, data[group1]["values"][t])] 1468 1469 if isinstance(testes_col, str): 1470 g2 = data[group2]["values"][testes_col] 1471 elif isinstance(testes_col, list): 1472 g2 = [1] * len(data[group2]["values"][testes_col[0]]) 1473 for t in testes_col: 1474 if comb == "*": 1475 g2 = [a * b for a, b in zip(g2, data[group2]["values"][t])] 1476 elif comb == "+": 1477 g2 = [a + b for a, b in zip(g2, data[group2]["values"][t])] 1478 elif comb == "**": 1479 g2 = [a**b for a, b in zip(g2, data[group2]["values"][t])] 1480 elif comb == "-": 1481 g2 = [a - b for a, b in zip(g2, data[group2]["values"][t])] 1482 elif comb == "/": 1483 g2 = [a / b for a, b in zip(g2, data[group2]["values"][t])] 1484 1485 stat, p_val = stats.ttest_ind( 1486 g1, g2, alternative="two-sided", equal_var=False 1487 ) 1488 g = sorted([group1, group2]) 1489 final_results["group1"].append(g[0]) 1490 final_results["group2"].append(g[1]) 1491 final_results["stat"].append(stat) 1492 final_results["p_val"].append(p_val) 1493 adj = p_val * len(pairs) 1494 if adj > 1: 1495 final_results["adj_p_val"].append(1) 1496 else: 1497 final_results["adj_p_val"].append(adj) 1498 1499 return p_val, final_results 1500 1501 def ks_percentiles(self, input_hist): 1502 """ 1503 Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups. 1504 1505 This method extracts the percentile levels and computes the average value for 1506 each percentile to obtain a lower-dimensional representation of the data, thereby 1507 reducing the Big Data scale problem for each group. Using these metrics, it reconstructs 1508 the underlying empirical distributions to evaluate both structural proportions and scale. 1509 1510 To further mitigate the large sample size problem ("curse of Big Data") where inflating 1511 pixel counts yields artificially significant results, a controlled downsampling (resampling) 1512 is applied to standardize the sample sizes across groups. 1513 1514 A two-sample Kolmogorov-Smirnov (KS) test is then performed for every possible pair 1515 of groups to detect differences in distribution shapes. Finally, p-values are adjusted 1516 using the Bonferroni correction method to account for multiple comparisons and control 1517 the family-wise error rate. 1518 1519 Parameters 1520 ---------- 1521 input_hist : dict 1522 A nested dictionary where keys are group names. Each group must contain 1523 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1524 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1525 1526 Returns 1527 ------- 1528 final_results : dict 1529 A dictionary containing the results of the pairwise comparisons with keys: 1530 - 'group1': list of the first group names in the pairs. 1531 - 'group2': list of the second group names in the pairs. 1532 - 'K-S': list of Kolmogorov-Smirnov test statistics. 1533 - 'p_val': list of unadjusted p-values. 1534 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1535 1536 Example 1537 ------- 1538 >>> results = self.ks_percentiles(input_hist) 1539 """ 1540 1541 ks_data = {} 1542 1543 for d in input_hist.keys(): 1544 tmp_dic = {} 1545 1546 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["avg"]): 1547 tmp_dic[f"p{n+1}"] = c 1548 1549 ks_data[d] = tmp_dic 1550 1551 df_cleaned = pd.DataFrame(ks_data).T 1552 1553 pairs = list(combinations(df_cleaned.index, 2)) 1554 1555 final_results = { 1556 "group1": [], 1557 "group2": [], 1558 "K-S": [], 1559 "p_val": [], 1560 "adj_p_val": [], 1561 } 1562 1563 for group1, group2 in pairs: 1564 1565 g = sorted([group1, group2]) 1566 1567 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1568 1569 res = stats.ks_2samp(table_pair.iloc[:, 0], table_pair.iloc[:, 1]) 1570 1571 final_results["group1"].append(g[0]) 1572 final_results["group2"].append(g[1]) 1573 final_results["K-S"].append(res.statistic) 1574 final_results["p_val"].append(res.pvalue) 1575 adj = res.pvalue * len(pairs) 1576 if adj > 1: 1577 final_results["adj_p_val"].append(1) 1578 else: 1579 final_results["adj_p_val"].append(adj) 1580 1581 return final_results 1582 1583 def fisher_percentiles(self, input_hist): 1584 """ 1585 Perform pairwise Fisher's exact tests on percentile data across all groups. 1586 1587 This method extracts the raw counts (N) for each percentile bin across all 1588 groups to construct a contingency table representation of the data. By utilizing 1589 the discrete frequency counts per bin rather than continuous average values, it 1590 evaluates both structural distribution proportions and sample size scaling 1591 differences simultaneously. 1592 1593 An exact testing approach is applied to every unique pair of groups by extracting 1594 their corresponding sub-tables. For each pair, a Fisher's exact test (or its 1595 extension for larger contingency tables) is performed to detect statistically 1596 significant deviations in distribution profiles. 1597 1598 Finally, p-values are manually adjusted using the Bonferroni correction method 1599 by multiplying the raw p-values by the total number of comparisons to control 1600 the family-wise error rate across multiple pair-wise tests. 1601 the family-wise error rate. 1602 1603 Parameters 1604 ---------- 1605 input_hist : dict 1606 A nested dictionary where keys are group names. Each group must contain 1607 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1608 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1609 1610 Returns 1611 ------- 1612 final_results : dict 1613 A dictionary containing the results of the pairwise comparisons with keys: 1614 - 'group1': list of the first group names in the pairs. 1615 - 'group2': list of the second group names in the pairs. 1616 - 'fish': list of Fisher's exact test statistics. 1617 - 'p_val': list of unadjusted p-values. 1618 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1619 1620 Example 1621 ------- 1622 >>> results = self.fisher_percentiles(input_hist) 1623 """ 1624 1625 fish_data = {} 1626 1627 for d in input_hist.keys(): 1628 tmp_dic = {} 1629 1630 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["n"]): 1631 tmp_dic[f"p{n+1}"] = c 1632 1633 fish_data[d] = tmp_dic 1634 1635 df_cleaned = pd.DataFrame(fish_data).T 1636 1637 pairs = list(combinations(df_cleaned.index, 2)) 1638 1639 final_results = { 1640 "group1": [], 1641 "group2": [], 1642 "fish": [], 1643 "p_val": [], 1644 "adj_p_val": [], 1645 } 1646 1647 for group1, group2 in pairs: 1648 1649 g = sorted([group1, group2]) 1650 1651 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1652 1653 res = stats.fisher_exact(table_pair) 1654 1655 final_results["group1"].append(g[0]) 1656 final_results["group2"].append(g[1]) 1657 final_results["fish"].append(res.statistic) 1658 final_results["p_val"].append(res.pvalue) 1659 adj = res.pvalue * len(pairs) 1660 if adj > 1: 1661 final_results["adj_p_val"].append(1) 1662 else: 1663 final_results["adj_p_val"].append(adj) 1664 1665 return final_results 1666 1667 def to_wasserstein_distance(self, data): 1668 """ 1669 Calculate scaled pairwise Wasserstein distances for grouped distributions. 1670 1671 This method computes the 1D Wasserstein distance (Earth Mover's Distance) 1672 between all possible combinations of groups in the provided dataset. 1673 Before calculating the distance, the standardized frequencies are scaled 1674 by a factor representing the average total count (sample size) of the 1675 two compared groups. 1676 1677 Parameters 1678 ---------- 1679 data : dict 1680 A nested dictionary where keys are group names. For each group, the 1681 method expects the following internal data structure: 1682 - `data[group_name]['percentiles']['mutual']['n']` : list-like 1683 Absolute counts or sample sizes for the distribution. 1684 - `data[group_name]['percentiles']['mutual']['n_standarized']` : list-like 1685 Standardized frequencies or probabilities to be compared. 1686 1687 Returns 1688 ------- 1689 final_results : dict 1690 A dictionary containing the results of the pairwise distance calculations: 1691 - 'group1' : list of str 1692 The name of the first group in the comparison. 1693 - 'group2' : list of str 1694 The name of the second group in the comparison. 1695 - 'wasserstein_distance' : list of float 1696 The computed scaled Wasserstein distance for each pair. 1697 """ 1698 1699 pairs = list(combinations(data.keys(), 2)) 1700 1701 final_results = {"group1": [], "group2": [], "wasserstein_distance": []} 1702 1703 for group1, group2 in pairs: 1704 1705 factor = ( 1706 sum(data[group1]["percentiles"]["mutual"]["n"]) 1707 + sum(data[group2]["percentiles"]["mutual"]["n"]) 1708 ) / 2 1709 1710 dist = wasserstein_distance( 1711 [ 1712 x * factor 1713 for x in data[group1]["percentiles"]["mutual"]["n_standarized"] 1714 ], 1715 [ 1716 x * factor 1717 for x in data[group2]["percentiles"]["mutual"]["n_standarized"] 1718 ], 1719 ) 1720 1721 g = sorted([group1, group2]) 1722 final_results["group1"].append(g[0]) 1723 final_results["group2"].append(g[1]) 1724 final_results["wasserstein_distance"].append(dist) 1725 1726 return final_results 1727 1728 def to_fold_change(self, data, tested_value): 1729 """ 1730 Calculate the Fold Change (FC) between all permutations of groups. 1731 1732 This method computes the ratio of the mean values of a specified feature 1733 (`tested_value`) for every directed pair of groups. Because permutations 1734 are used, the calculation is directional (i.e., both Group A / Group B 1735 and Group B / Group A are computed). 1736 1737 Parameters 1738 ---------- 1739 data : dict 1740 A nested dictionary where keys are group names. For each group, the 1741 method expects the following internal structure: 1742 - `data[group_name]['values'][tested_value]` : array-like 1743 Numeric values used to compute the mean for the group. 1744 1745 tested_value : str 1746 The specific key or column name within the 'values' dictionary 1747 indicating which feature's fold change should be calculated. 1748 1749 Returns 1750 ------- 1751 final_results : dict 1752 A dictionary containing the results of the pairwise fold change calculations: 1753 - 'group1' : list of str 1754 The name of the numerator group in the comparison. 1755 - 'group2' : list of str 1756 The name of the denominator group in the comparison. 1757 - 'FC' : list of float 1758 The calculated fold change (mean of group1 / mean of group2). 1759 """ 1760 1761 pairs = list(permutations(data.keys(), 2)) 1762 1763 final_results = {"group1": [], "group2": [], "FC": []} 1764 1765 values = [] 1766 for group1, group2 in pairs: 1767 1768 values = values + data[group1]["values"][tested_value] 1769 values = values + data[group2]["values"][tested_value] 1770 1771 values_min = min([x for x in values if x > 0]) 1772 values_min = values_min / 2 1773 1774 for group1, group2 in pairs: 1775 1776 g1 = np.mean(data[group1]["values"][tested_value]) 1777 g2 = np.mean(data[group2]["values"][tested_value]) 1778 1779 if g1 == 0: 1780 g1 = g1 + values_min 1781 1782 if g2 == 0: 1783 g2 = g2 + values_min 1784 1785 fc = g1 / g2 1786 1787 final_results["group1"].append(group1) 1788 final_results["group2"].append(group2) 1789 final_results["FC"].append(fc) 1790 1791 return final_results 1792 1793 def get_stats(self, data, tested_value): 1794 """ 1795 Calculate and aggregate statistical metrics (ANOVA, Fisher's Exact, 1796 Kolmogorov-Smirnov, Fold Change, Wasserstein distance). 1797 1798 This method computes overall statistics and pairwise comparisons for grouped data. 1799 To properly capture both structural proportions and total count variations across 1800 percentiles while avoiding the curse of Big Data, it runs two distinct tests: 1801 1. Fisher's exact test on discrete percentile counts to evaluate absolute scale 1802 and profile differences. 1803 2. Two-sample Kolmogorov-Smirnov (KS) test on reconstructed empirical 1804 distributions to evaluate discrepancies in distribution shapes. 1805 1806 Additionally, it calculates the Fold Change (FC) and evaluates Wasserstein 1807 distances. If the average number of replicates per group is at least 3, 1808 it conducts Welch's ANOVA. The input dictionary is modified in-place to 1809 include a new 'statistics' key containing all results. 1810 1811 Parameters 1812 ---------- 1813 data : dict 1814 A nested dictionary where keys are group names. Each group's dictionary 1815 must contain the structure `['values']['replication']` to verify sample sizes, 1816 along with the necessary data structures required by downstream statistical methods. 1817 1818 tested_value : str 1819 The key or column name representing the specific variable to evaluate 1820 (e.g., used for ANOVA and Fold Change calculations). 1821 1822 Returns 1823 ------- 1824 data : dict 1825 The original input dictionary, extended with a new `data['statistics']` key 1826 that houses the computed statistical results, including `percintiles_fish` 1827 and `percintiles_ks`. 1828 1829 Example 1830 ------- 1831 stats = self.get_stats( 1832 data, 1833 tested_value='n', 1834 ) 1835 """ 1836 1837 # parametric selected value 1838 sum_k = 0 1839 n = 0 1840 for k in data.keys(): 1841 if k != "statistics": 1842 n += 1 1843 sum_k += len(data[k]["values"]["replication"]) 1844 1845 sum_k = sum_k / n 1846 1847 if sum_k >= 3: 1848 pk, dfk = self.post_aov(data, testes_col=tested_value) 1849 1850 # fish 1851 fish = self.fisher_percentiles(data) 1852 1853 # K_S 1854 ks = self.ks_percentiles(data) 1855 1856 dw = self.to_wasserstein_distance(data) 1857 1858 fc = self.to_fold_change(data, tested_value) 1859 1860 data["statistics"] = {} 1861 1862 data["statistics"]["percintiles_fish"] = fish 1863 1864 data["statistics"]["percintiles_ks"] = ks 1865 1866 if sum_k >= 3: 1867 data["statistics"]["ANOVA"] = {} 1868 1869 data["statistics"]["ANOVA"]["p_value"] = pk 1870 data["statistics"]["ANOVA"]["pair-comparison"] = dfk 1871 else: 1872 import warnings 1873 1874 warnings.warn( 1875 f"Insufficient replicates for statistical analysis. " 1876 f"At least 3 replicates per group (3 vs 3) are required. " 1877 f"The average number of samples per probe in this dataset was {n}.", 1878 RuntimeWarning, 1879 ) 1880 1881 data["statistics"]["FC"] = fc 1882 1883 data["statistics"]["wasserstein_distance"] = dw 1884 1885 data["statistics"]["tested_value"] = tested_value 1886 1887 return data 1888 1889 def hist_compare_plot( 1890 self, data, queue=None, p_adj: bool = True, txt_size: int = 20 1891 ): 1892 """ 1893 Generate comparative histograms and display results of statistical tests 1894 (ANOVA / T-test [if min 3 vs. 3 comparison available], Kolmogorov-Smirnov, Fisher percentiles) 1895 and statistics (FC, Wasserstein distance). 1896 1897 1898 Parameters 1899 ---------- 1900 data : dict 1901 Dictionary where keys are group names and values are containing histogram data. 1902 Each DataFrame should include the column specified by `tested_value`. 1903 1904 queue : list of str or None 1905 Defines the order of groups to be plotted. 1906 1907 p_adj : bool, optional 1908 If True, applies Bonferroni correction for multiple comparisons (default is True). 1909 1910 txt_size : int, optional 1911 Font size for text annotations in the plot (default is 20). 1912 1913 Returns 1914 ------- 1915 fig : matplotlib.figure.Figure 1916 Matplotlib figure object containing the generated histograms and statistical test results. 1917 1918 Example 1919 ------- 1920 fig = self.hist_compare_plot( 1921 data, 1922 queue=['group1', 'group2', 'group3'], 1923 p_adj=True, 1924 txt_size=18 1925 ) 1926 plt.show() 1927 """ 1928 1929 if queue is None: 1930 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1931 1932 if sorted(queue) != [x for x in sorted(data.keys()) if x != "statistics"]: 1933 print( 1934 "\n Wrong queue provided! The queue will be sorted with default settings!" 1935 ) 1936 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1937 1938 # parametric selected value 1939 tested_value = data["statistics"]["tested_value"] 1940 1941 ############################################################################## 1942 1943 standarized_max, standarized_min, value_max, value_min = [], [], [], [] 1944 for d in queue: 1945 standarized_max.append( 1946 max(data[d]["percentiles"]["mutual"]["n_standarized"]) 1947 ) 1948 standarized_min.append( 1949 min(data[d]["percentiles"]["mutual"]["n_standarized"]) 1950 ) 1951 value_max.append(max(data[d]["percentiles"]["mutual"][tested_value])) 1952 value_min.append(min(data[d]["percentiles"]["mutual"][tested_value])) 1953 1954 num_columns = len(queue) + 1 1955 1956 fig, axs = plt.subplots( 1957 3, 1958 num_columns, 1959 figsize=(8 * num_columns, 10), 1960 gridspec_kw={"width_ratios": [1] * len(queue) + [0.5], "wspace": 0.05}, 1961 ) 1962 1963 for i, d in enumerate(queue): 1964 tmp_data = data[d]["percentiles"]["mutual"] 1965 1966 axs[0, i].bar( 1967 [str(n) for n in range(len(tmp_data["n_standarized"]))], 1968 tmp_data["n_standarized"], 1969 width=0.95, 1970 color="gold", 1971 ) 1972 1973 # line 1974 n_groups = len(data[d]["percentiles"]["replications"].keys()) 1975 colors = plt.cm.OrRd(np.linspace(0.3, 0.9, n_groups)) 1976 1977 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 1978 1979 color = colors[ix] 1980 1981 y = data[d]["percentiles"]["replications"][dn]["n_standarized"] 1982 x = np.arange(len(y)) 1983 1984 axs[0, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 1985 1986 axs[0, i].plot( 1987 x, 1988 y, 1989 color=color, 1990 linewidth=1, 1991 marker="o", 1992 ) 1993 1994 axs[0, i].set_ylim( 1995 min(standarized_min) * 0.9995, max(standarized_max) * 1.0005 1996 ) 1997 1998 if i == 0: 1999 axs[0, i].set_ylabel("Standarized\nfrequency", fontsize=txt_size) 2000 else: 2001 axs[0, i].set_yticks([]) 2002 2003 axs[0, i].set_xticks([]) 2004 axs[0, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2005 2006 axs[1, i].bar( 2007 [str(n) for n in range(len(tmp_data[tested_value]))], 2008 tmp_data[tested_value], 2009 width=0.95, 2010 color="orange", 2011 ) 2012 2013 # line 2014 2015 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2016 2017 color = colors[ix] 2018 2019 y = data[d]["percentiles"]["replications"][dn][tested_value] 2020 x = np.arange(len(y)) 2021 2022 axs[1, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2023 2024 axs[1, i].plot( 2025 x, 2026 y, 2027 color=color, 2028 linewidth=1, 2029 marker="o", 2030 ) 2031 2032 mean_value = np.mean(data[d]["values"][tested_value]) 2033 axs[1, i].axhline(y=mean_value, color="red", linestyle="--") 2034 2035 axs[1, i].set_ylim(min(value_min) * 0.9995, max(value_max) * 1.0005) 2036 2037 if i == 0: 2038 axs[1, i].set_ylabel(f"Normalized\n{tested_value}", fontsize=txt_size) 2039 else: 2040 axs[1, i].set_yticks([]) 2041 2042 axs[1, i].set_xticks([]) 2043 axs[1, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2044 2045 axs[2, i].bar( 2046 [str(n) for n in range(len(tmp_data["n_standarized"]))], 2047 [ 2048 a * b 2049 for a, b in zip(tmp_data[tested_value], tmp_data["n_standarized"]) 2050 ], 2051 width=0.95, 2052 color="goldenrod", 2053 ) 2054 2055 # line 2056 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2057 2058 color = colors[ix] 2059 2060 y = [ 2061 a * b 2062 for a, b in zip( 2063 data[d]["percentiles"]["replications"][dn][tested_value], 2064 data[d]["percentiles"]["replications"][dn]["n_standarized"], 2065 ) 2066 ] 2067 x = np.arange(len(y)) 2068 2069 axs[2, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2070 2071 axs[2, i].plot( 2072 x, 2073 y, 2074 color=color, 2075 linewidth=1, 2076 marker="o", 2077 ) 2078 2079 mean_value = np.mean( 2080 data[d]["values"][data["statistics"]["tested_value"]] 2081 ) * np.mean(tmp_data["n_standarized"]) 2082 2083 axs[2, i].axhline(y=mean_value, color="red", linestyle="--") 2084 2085 axs[2, i].set_ylim( 2086 (min(standarized_min) * min(value_min)) * 0.9995, 2087 (max(standarized_max) * max(value_max) * 1.0005), 2088 ) 2089 axs[2, i].set_xlabel(d, fontsize=txt_size) 2090 2091 if i == 0: 2092 axs[2, i].set_ylabel( 2093 f"Standarized\nnorm_{tested_value}", fontsize=txt_size 2094 ) 2095 else: 2096 axs[2, i].set_yticks([]) 2097 2098 axs[2, i].set_xticks([]) 2099 axs[2, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2100 2101 # statistics 2102 2103 # ANOVA / t-test 2104 2105 if "ANOVA" in data["statistics"].keys(): 2106 pk = data["statistics"]["ANOVA"]["p_value"] 2107 dfk = data["statistics"]["ANOVA"]["pair-comparison"] 2108 dfk = pd.DataFrame(dfk) 2109 2110 dfk = dfk.sort_values( 2111 by=["group1", "group2"], 2112 key=lambda col: [ 2113 queue.index(val) if val in queue else -1 for val in col 2114 ], 2115 ).reset_index(drop=True) 2116 2117 sign = "ns" 2118 if float(self.round_to_scientific_notation(pk)) < 0.001: 2119 sign = "***" 2120 elif float(self.round_to_scientific_notation(pk)) < 0.01: 2121 sign = "**" 2122 elif float(self.round_to_scientific_notation(pk)) < 0.05: 2123 sign = "*" 2124 2125 text = f"Test Welch's ANOVA\non '{tested_value}' values\np-value: {self.round_to_scientific_notation(pk)} - {sign}\n" 2126 2127 if p_adj == True: 2128 for i in range(len(dfk["group1"])): 2129 sign = "ns" 2130 if dfk["adj_p_val"][i] < 0.001: 2131 sign = "***" 2132 elif dfk["adj_p_val"][i] < 0.01: 2133 sign = "**" 2134 elif dfk["adj_p_val"][i] < 0.05: 2135 sign = "*" 2136 2137 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['adj_p_val'][i])} - {sign}\n" 2138 else: 2139 for i in range(len(dfk["group1"])): 2140 sign = "ns" 2141 if dfk["p_val"][i] < 0.001: 2142 sign = "***" 2143 elif dfk["p_val"][i] < 0.01: 2144 sign = "**" 2145 elif dfk["p_val"][i] < 0.05: 2146 sign = "*" 2147 2148 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['p_val'][i])} - {sign}\n" 2149 2150 axs[2, -1].text( 2151 0.5, 2152 0.5, 2153 text, 2154 ha="center", 2155 va="center", 2156 fontsize=txt_size * 0.7, 2157 wrap=True, 2158 ) 2159 axs[2, -1].set_axis_off() 2160 else: 2161 axs[2, -1].set_axis_off() 2162 2163 # FC / Distance 2164 2165 ranking_FC = pd.DataFrame(data["statistics"]["FC"]) 2166 2167 ranking_dw = pd.DataFrame(data["statistics"]["wasserstein_distance"]) 2168 2169 ranking_combined = pd.merge( 2170 ranking_FC, ranking_dw, on=["group1", "group2"], how="right" 2171 ) 2172 2173 ranking_combined = ranking_combined.sort_values( 2174 by=["group1", "group2"], 2175 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2176 ).reset_index(drop=True) 2177 2178 text = "FC / Wasserstein distance\n" 2179 for i in range(len(ranking_combined)): 2180 group1 = ranking_combined["group1"][i] 2181 group2 = ranking_combined["group2"][i] 2182 fc_val = ranking_combined["FC"][i] 2183 wasserstein_val = ranking_combined["wasserstein_distance"][i] 2184 2185 text += f"{group1} vs. {group2}\n {fc_val:.2f} | {wasserstein_val:.2f}\n" 2186 2187 axs[1, -1].text( 2188 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2189 ) 2190 axs[1, -1].set_axis_off() 2191 2192 # fish 2193 2194 fish = pd.DataFrame(data["statistics"]["percintiles_fish"]) 2195 2196 # K-S 2197 2198 ks = pd.DataFrame(data["statistics"]["percintiles_ks"]) 2199 2200 fish = fish.sort_values( 2201 by=["group1", "group2"], 2202 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2203 ).reset_index(drop=True) 2204 2205 ks = ks.sort_values( 2206 by=["group1", "group2"], 2207 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2208 ).reset_index(drop=True) 2209 2210 text = f"Tests\nFisher's exact / Kolmogorov-Smirnov\n" 2211 2212 if p_adj == True: 2213 for i in range(len(fish["group1"])): 2214 sign1 = "ns" 2215 if fish["adj_p_val"][i] < 0.001: 2216 sign1 = "***" 2217 elif fish["adj_p_val"][i] < 0.01: 2218 sign1 = "**" 2219 elif fish["adj_p_val"][i] < 0.05: 2220 sign1 = "*" 2221 2222 sign2 = "ns" 2223 if ks["adj_p_val"][i] < 0.001: 2224 sign2 = "***" 2225 elif ks["adj_p_val"][i] < 0.01: 2226 sign2 = "**" 2227 elif ks["adj_p_val"][i] < 0.05: 2228 sign2 = "*" 2229 2230 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2231 2232 else: 2233 for i in range(len(fish["group1"])): 2234 sign1 = "ns" 2235 if fish["p_val"][i] < 0.001: 2236 sign1 = "***" 2237 elif fish["p_val"][i] < 0.01: 2238 sign1 = "**" 2239 elif fish["p_val"][i] < 0.05: 2240 sign1 = "*" 2241 2242 sign2 = "ns" 2243 if ks["p_val"][i] < 0.001: 2244 sign2 = "***" 2245 elif ks["p_val"][i] < 0.01: 2246 sign2 = "**" 2247 elif ks["p_val"][i] < 0.05: 2248 sign2 = "*" 2249 2250 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2251 2252 axs[0, -1].text( 2253 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2254 ) 2255 axs[0, -1].set_axis_off() 2256 2257 plt.tight_layout() 2258 2259 if cfg._DISPLAY_MODE: 2260 plt.show() 2261 2262 return fig
Class for performing percentile-based statistical analysis on grouped data.
This class provides methods to calculate percentiles, remove outliers, aggregate data into percentile bins, perform Welch's ANOVA, Kolmogorov-Smirnov and Fisher's exact tests, evaluate Wasserstein distances, calculate Fold Change, and visualize results via comparative histograms. It is designed to handle both single-column and multi-column combinations of values for group-based analysis.
Methods
drop_up_df(data, group_col, values_col) Removes upper outliers from a DataFrame based on a grouping column.
percentiles_calculation(values, sep_perc=1) Calculates percentiles and creates loopable percentile ranges.
to_percentil(values, percentiles, percentiles_loop, values_col, replication_col) Aggregates statistics based on percentile ranges.
df_to_percentiles(data, group_col, values_col, replication_col, sep_perc=1, drop_outlires=True) Computes percentile statistics for grouped DataFrame data.
round_to_scientific_notation(num) Formats a number in scientific notation or standard format.
aov(data, testes_col, comb="*") Performs Welch's ANOVA on percentile-based group data.
post_aov(data, testes_col, comb="*") Performs Welch's ANOVA with pairwise t-tests.
ks_percentiles(input_hist) Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups.
fisher_percentiles(input_hist) Perform pairwise Fisher's exact tests on percentile data across all groups.
to_wasserstein_distance(data) Calculates scaled pairwise Wasserstein distances for grouped distributions.
to_fold_change(data, tested_value) Calculates the Fold Change (FC) between all directed permutations of groups.
get_stats(data, tested_value) Calculates and aggregates overall statistical metrics (ANOVA, Fisher's exact, Kolmogorov-Smirnov, FC, Wasserstein distance).
hist_compare_plot(data, queue=None, p_adj=True, txt_size=20) Generates comparative histograms with statistical test results and metrics.
def
drop_up_df( self, data: pandas.core.frame.DataFrame, group_col: str, values_col: str):
1007 def drop_up_df(self, data: pd.DataFrame, group_col: str, values_col: str): 1008 """ 1009 Remove upper outliers from a DataFrame based on a specified value column, grouped by a grouping column. 1010 1011 Outliers are calculated and removed separately for each group defined by `group_col`. 1012 The upper outliers are defined using the interquartile range (IQR) method: 1013 values greater than Q3 + 1.5 * IQR are considered outliers. 1014 1015 Parameters 1016 ---------- 1017 data : pd.DataFrame 1018 The input DataFrame containing the data. 1019 1020 group_col : str 1021 The name of the column used for grouping the data. 1022 1023 values_col : str 1024 The column containing the values from which upper outliers will be removed. 1025 1026 Returns 1027 ------- 1028 filtered_data : pd.DataFrame 1029 A filtered DataFrame with the upper outliers removed for each group. 1030 1031 Notes 1032 ----- 1033 - Outliers are removed separately within each group. 1034 - The original DataFrame is not modified; a new filtered DataFrame is returned. 1035 """ 1036 1037 def iqr_filter(group): 1038 q75 = np.quantile(group[values_col], 0.75) 1039 q25 = np.quantile(group[values_col], 0.25) 1040 itq = q75 - q25 1041 return group[group[values_col] <= (q75 + 1.5 * itq)] 1042 1043 filtered_data = data.groupby(group_col).apply(iqr_filter).reset_index(drop=True) 1044 1045 return filtered_data
Remove upper outliers from a DataFrame based on a specified value column, grouped by a grouping column.
Outliers are calculated and removed separately for each group defined by group_col.
The upper outliers are defined using the interquartile range (IQR) method:
values greater than Q3 + 1.5 * IQR are considered outliers.
Parameters
data : pd.DataFrame The input DataFrame containing the data.
group_col : str The name of the column used for grouping the data.
values_col : str The column containing the values from which upper outliers will be removed.
Returns
filtered_data : pd.DataFrame A filtered DataFrame with the upper outliers removed for each group.
Notes
- Outliers are removed separately within each group.
- The original DataFrame is not modified; a new filtered DataFrame is returned.
def
percentiles_calculation(self, values, sep_perc: int = 1):
1047 def percentiles_calculation(self, values, sep_perc: int = 1): 1048 """ 1049 Calculate percentiles for a set of values and generate consecutive percentile ranges. 1050 1051 This function computes percentiles from 0 to 100 at intervals defined by `sep_perc`. 1052 It also generates a list of consecutive percentile ranges that can be used for further analysis or binning. 1053 1054 Parameters 1055 ---------- 1056 values : array-like 1057 The input data values for which the percentiles are calculated. 1058 1059 sep_perc : int, optional 1060 Separation interval between percentiles (default is 1, meaning percentiles are calculated every 1%). 1061 1062 Returns 1063 ------- 1064 percentiles : np.ndarray 1065 Array of calculated percentile values. 1066 1067 percentiles_loop : list of tuple 1068 List of consecutive percentile ranges as tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1069 1070 Notes 1071 ----- 1072 - The first percentile is set to 0 to avoid issues with zero values. 1073 - `percentiles_loop` is useful for iterating through percentile ranges when aggregating statistics. 1074 """ 1075 1076 per_vector = values.copy() 1077 1078 percentiles = np.percentile(per_vector, np.arange(0, 101, sep_perc)) 1079 percentiles[0] = 0 1080 1081 percentiles_loop = [(i, i + 1) for i in range(int(100 / sep_perc))] 1082 1083 return percentiles, percentiles_loop
Calculate percentiles for a set of values and generate consecutive percentile ranges.
This function computes percentiles from 0 to 100 at intervals defined by sep_perc.
It also generates a list of consecutive percentile ranges that can be used for further analysis or binning.
Parameters
values : array-like The input data values for which the percentiles are calculated.
sep_perc : int, optional Separation interval between percentiles (default is 1, meaning percentiles are calculated every 1%).
Returns
percentiles : np.ndarray Array of calculated percentile values.
percentiles_loop : list of tuple List of consecutive percentile ranges as tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)].
Notes
- The first percentile is set to 0 to avoid issues with zero values.
percentiles_loopis useful for iterating through percentile ranges when aggregating statistics.
def
to_percentil( self, values, percentiles, percentiles_loop, values_col, replication_col):
1085 def to_percentil( 1086 self, values, percentiles, percentiles_loop, values_col, replication_col 1087 ): 1088 """ 1089 Aggregate statistics for a set of values based on percentile ranges, including replications. 1090 1091 This function calculates summary statistics (count, proportion, mean, median, 1092 standard deviation, variance) for each percentile range defined in `percentiles_loop`. 1093 It computes these statistics both for the combined data ('mutual') and separately 1094 for each individual replication. It also calculates overall metrics per replication. 1095 1096 Parameters 1097 ---------- 1098 values : pd.DataFrame 1099 Input DataFrame containing the data to be analyzed. 1100 percentiles : np.ndarray 1101 Array of percentile values used to define the boundaries of each range. 1102 percentiles_loop : list of tuple 1103 List of consecutive percentile ranges as index tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)]. 1104 values_col : str 1105 The column name in `values` containing the numeric data to aggregate. 1106 replication_col : str 1107 The column name in `values` used to identify distinct replications or samples. 1108 1109 Returns 1110 ------- 1111 full_data : dict 1112 A nested dictionary containing the calculated statistics with the following structure: 1113 - 'percentiles' : dict 1114 - 'mutual' : dict 1115 Lists of statistics ('n', 'n_standarized', 'avg', 'median', 'std', 'var') 1116 aggregated across all replications for each percentile bin. 1117 - 'replications' : dict 1118 Keys are replication names. Values are dictionaries of statistics (same as above) 1119 calculated specifically for that replication within each bin. 1120 - 'values' : dict 1121 Lists of overall statistics ('avg', 'median', 'std', 'var', 'replication') 1122 calculated for each replication as a whole (ignoring bins). 1123 1124 Notes 1125 ----- 1126 - If a percentile range contains no elements, statistics are set to 0 and count is set to 1 to avoid empty lists. 1127 """ 1128 1129 full_data = {} 1130 per_vector = values[values_col] 1131 amount = len(per_vector) 1132 1133 data_mutual = { 1134 "n": [], 1135 "n_standarized": [], 1136 "avg": [], 1137 "median": [], 1138 "std": [], 1139 "var": [], 1140 } 1141 1142 for x in percentiles_loop: 1143 subset = per_vector[ 1144 (per_vector > percentiles[x[0]]) & (per_vector <= percentiles[x[1]]) 1145 ] 1146 n_subset = len(subset) 1147 1148 if n_subset > 0: 1149 data_mutual["n"].append(n_subset) 1150 data_mutual["n_standarized"].append(n_subset / amount) 1151 data_mutual["avg"].append(np.mean(subset)) 1152 data_mutual["median"].append(np.median(subset)) 1153 data_mutual["std"].append(np.std(subset)) 1154 data_mutual["var"].append(np.var(subset)) 1155 else: 1156 data_mutual["n"].append(0) 1157 data_mutual["n_standarized"].append(0) 1158 data_mutual["avg"].append(0) 1159 data_mutual["median"].append(0) 1160 data_mutual["std"].append(0) 1161 data_mutual["var"].append(0) 1162 1163 full_data["percentiles"] = {"mutual": data_mutual, "replications": {}} 1164 1165 unique_names = set(values[replication_col]) 1166 1167 for nam in unique_names: 1168 per_vector_rep = values[values_col][values[replication_col] == nam] 1169 1170 data_rep = { 1171 "n": [], 1172 "n_standarized": [], 1173 "avg": [], 1174 "median": [], 1175 "std": [], 1176 "var": [], 1177 } 1178 1179 for x in percentiles_loop: 1180 subset = per_vector_rep[ 1181 (per_vector_rep > percentiles[x[0]]) 1182 & (per_vector_rep <= percentiles[x[1]]) 1183 ] 1184 n_subset = len(subset) 1185 1186 if n_subset > 0: 1187 data_rep["n"].append(n_subset) 1188 data_rep["n_standarized"].append(n_subset / amount) 1189 data_rep["avg"].append(np.mean(subset)) 1190 data_rep["median"].append(np.median(subset)) 1191 data_rep["std"].append(np.std(subset)) 1192 data_rep["var"].append(np.var(subset)) 1193 else: 1194 data_rep["n"].append(0) 1195 data_rep["n_standarized"].append(0) 1196 data_rep["avg"].append(0) 1197 data_rep["median"].append(0) 1198 data_rep["std"].append(0) 1199 data_rep["var"].append(0) 1200 1201 full_data["percentiles"]["replications"][nam] = data_rep 1202 1203 unique_names = set(values[replication_col]) 1204 1205 data_rep = {"avg": [], "median": [], "std": [], "var": [], "replication": []} 1206 1207 for nam in unique_names: 1208 per_vector_rep = values[values_col][values[replication_col] == nam] 1209 1210 data_rep["avg"].append(np.mean(per_vector_rep)) 1211 data_rep["median"].append(np.median(per_vector_rep)) 1212 data_rep["std"].append(np.std(per_vector_rep)) 1213 data_rep["var"].append(np.var(per_vector_rep)) 1214 data_rep["replication"].append(nam) 1215 1216 full_data["values"] = data_rep 1217 1218 return full_data
Aggregate statistics for a set of values based on percentile ranges, including replications.
This function calculates summary statistics (count, proportion, mean, median,
standard deviation, variance) for each percentile range defined in percentiles_loop.
It computes these statistics both for the combined data ('mutual') and separately
for each individual replication. It also calculates overall metrics per replication.
Parameters
values : pd.DataFrame
Input DataFrame containing the data to be analyzed.
percentiles : np.ndarray
Array of percentile values used to define the boundaries of each range.
percentiles_loop : list of tuple
List of consecutive percentile ranges as index tuples, e.g., [(0, 1), (1, 2), ..., (99, 100)].
values_col : str
The column name in values containing the numeric data to aggregate.
replication_col : str
The column name in values used to identify distinct replications or samples.
Returns
full_data : dict A nested dictionary containing the calculated statistics with the following structure: - 'percentiles' : dict - 'mutual' : dict Lists of statistics ('n', 'n_standarized', 'avg', 'median', 'std', 'var') aggregated across all replications for each percentile bin. - 'replications' : dict Keys are replication names. Values are dictionaries of statistics (same as above) calculated specifically for that replication within each bin. - 'values' : dict Lists of overall statistics ('avg', 'median', 'std', 'var', 'replication') calculated for each replication as a whole (ignoring bins).
Notes
- If a percentile range contains no elements, statistics are set to 0 and count is set to 1 to avoid empty lists.
def
df_to_percentiles( self, data: pandas.core.frame.DataFrame, group_col: str = 'individual_name', values_col: str = 'norm_intensity', replication_col: str = 'individual_number', sep_perc: int = 1, drop_outlires: bool = True):
1220 def df_to_percentiles( 1221 self, 1222 data: pd.DataFrame, 1223 group_col: str = "individual_name", 1224 values_col: str = "norm_intensity", 1225 replication_col: str = "individual_number", 1226 sep_perc: int = 1, 1227 drop_outlires: bool = True, 1228 ): 1229 """ 1230 Calculate summary statistics based on percentile ranges for each group in a DataFrame. 1231 1232 This method groups the input DataFrame by `group_col`, computes global percentile ranges 1233 based on `values_col`, and then aggregates statistics for each percentile bin. The aggregation 1234 is performed both mutually for the group and individually per replication. Optionally, 1235 upper outliers can be removed before the calculations. 1236 1237 Parameters 1238 ---------- 1239 data : pd.DataFrame 1240 Input DataFrame containing the grouped data. 1241 group_col : str, optional 1242 Column name used to define groups (default is 'individual_name'). 1243 values_col : str, optional 1244 Column name containing the numeric values for percentile calculations 1245 (default is 'norm_intensity'). 1246 replication_col : str, optional 1247 Column name used to identify separate replications within the groups 1248 (default is 'individual_number'). 1249 sep_perc : int, optional 1250 Separation interval for percentiles (default is 1, meaning 1% steps). 1251 drop_outlires : bool, optional 1252 If True, removes upper outliers from the data using the IQR method before 1253 performing calculations (default is True). 1254 1255 Returns 1256 ------- 1257 full_data : dict 1258 A dictionary where each key is a unique group name (from `group_col`). 1259 The corresponding value is the nested dictionary returned by `to_percentil()`, 1260 which includes bin-wise statistics ('mutual' and 'replications') and overall 1261 metrics ('values'). 1262 1263 Notes 1264 ----- 1265 - Outlier removal uses the IQR method within each group if `drop_outlires` is True. 1266 """ 1267 1268 full_data = {} 1269 1270 if drop_outlires == True: 1271 data = self.drop_up_df( 1272 data=data, group_col=group_col, values_col=values_col 1273 ) 1274 1275 groups = set(data[group_col]) 1276 val_dat = [x for x in data[values_col] if x > 0] 1277 1278 percentiles, percentiles_loop = self.percentiles_calculation( 1279 val_dat, sep_perc=sep_perc 1280 ) 1281 1282 for g in groups: 1283 1284 print(f"Group: {g} ...") 1285 1286 tmp_values = data[data[group_col] == g] 1287 1288 per_dat = self.to_percentil( 1289 tmp_values, percentiles, percentiles_loop, values_col, replication_col 1290 ) 1291 1292 full_data[g] = per_dat 1293 1294 return full_data
Calculate summary statistics based on percentile ranges for each group in a DataFrame.
This method groups the input DataFrame by group_col, computes global percentile ranges
based on values_col, and then aggregates statistics for each percentile bin. The aggregation
is performed both mutually for the group and individually per replication. Optionally,
upper outliers can be removed before the calculations.
Parameters
data : pd.DataFrame Input DataFrame containing the grouped data. group_col : str, optional Column name used to define groups (default is 'individual_name'). values_col : str, optional Column name containing the numeric values for percentile calculations (default is 'norm_intensity'). replication_col : str, optional Column name used to identify separate replications within the groups (default is 'individual_number'). sep_perc : int, optional Separation interval for percentiles (default is 1, meaning 1% steps). drop_outlires : bool, optional If True, removes upper outliers from the data using the IQR method before performing calculations (default is True).
Returns
full_data : dict
A dictionary where each key is a unique group name (from group_col).
The corresponding value is the nested dictionary returned by to_percentil(),
which includes bin-wise statistics ('mutual' and 'replications') and overall
metrics ('values').
Notes
- Outlier removal uses the IQR method within each group if
drop_outliresis True.
def
round_to_scientific_notation(self, num):
1296 def round_to_scientific_notation(self, num): 1297 """ 1298 Round a number to scientific notation if very small, otherwise to one decimal place. 1299 1300 Parameters 1301 ---------- 1302 num : float 1303 The number to round. 1304 1305 Returns 1306 ------- 1307 str 1308 The rounded number as a string. 1309 - If `num` is 0, returns "0.0". 1310 - If `abs(num) < 1e-4`, returns scientific notation with 1 decimal and 1-digit exponent. 1311 - Otherwise, returns the number rounded to one decimal place. 1312 """ 1313 1314 if num == 0: 1315 return "0.0" 1316 1317 if abs(num) < 0.0001: 1318 rounded_num = np.format_float_scientific(num, precision=1, exp_digits=1) 1319 return rounded_num 1320 else: 1321 return f"{num:.1f}"
Round a number to scientific notation if very small, otherwise to one decimal place.
Parameters
num : float The number to round.
Returns
str
The rounded number as a string.
- If num is 0, returns "0.0".
- If abs(num) < 1e-4, returns scientific notation with 1 decimal and 1-digit exponent.
- Otherwise, returns the number rounded to one decimal place.
def
aov(self, data, testes_col, comb: str = '*'):
1323 def aov(self, data, testes_col, comb: str = "*"): 1324 """ 1325 Perform a Welch's ANOVA analysis. 1326 1327 This function calculates group values by aggregating specified columns (testes_col) 1328 via the comb method and then conducts a Welch's ANOVA. This approach is ideal for 1329 comparing group means when data exhibits unequal variances across groups. 1330 1331 Parameters 1332 ---------- 1333 data : dict of pd.DataFrame 1334 Dictionary where keys are group names and values are DataFrames containing the data. 1335 1336 testes_col : str or list of str 1337 Column name(s) from which the group values are derived. If a list is provided, columns 1338 will be combined based on the `comb` operation. 1339 1340 comb : str, optional 1341 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1342 '*' : multiplication 1343 '+' : addition 1344 '**': exponentiation 1345 '-' : subtraction 1346 '/' : division 1347 Default is '*'. 1348 1349 Returns 1350 ------- 1351 F : float 1352 F-statistic from Welch's ANOVA. 1353 1354 p_val : float 1355 Uncorrected p-value from Welch's ANOVA, testing for significant differences between groups. 1356 1357 Notes 1358 ----- 1359 - If `testes_col` is a single string, no combination is performed, and the group values 1360 are taken directly from that column. 1361 - Welch's ANOVA is used as it accounts for unequal variances between groups. 1362 - The `df.melt()` method is used to reshape the data, allowing the ANOVA to be applied to all groups. 1363 1364 Examples 1365 -------- 1366 >>> welch_F, welch_p = self.aov(data, testes_col=['col1', 'col2'], comb='+') 1367 >>> print(f"Welch's ANOVA F-statistic: {welch_F}, p-value: {welch_p}") 1368 """ 1369 1370 groups = [] 1371 1372 for d in data.keys(): 1373 1374 if isinstance(testes_col, str): 1375 g = data[d]["values"][testes_col] 1376 elif isinstance(testes_col, list): 1377 g = [1] * len(data[d]["values"][testes_col[0]]) 1378 for t in testes_col: 1379 if comb == "*": 1380 g = [a * b for a, b in zip(g, data[d]["values"][t])] 1381 elif comb == "+": 1382 g = [a + b for a, b in zip(g, data[d]["values"][t])] 1383 elif comb == "**": 1384 g = [a**b for a, b in zip(g, data[d]["values"][t])] 1385 elif comb == "-": 1386 g = [a - b for a, b in zip(g, data[d]["values"][t])] 1387 elif comb == "/": 1388 g = [a / b for a, b in zip(g, data[d]["values"][t])] 1389 1390 groups.append(g) 1391 1392 df = pd.DataFrame({f"group_{i}": group for i, group in enumerate(groups)}) 1393 1394 df_melted = df.melt(var_name="group", value_name="value") 1395 1396 welch_results = pg.welch_anova(data=df_melted, dv="value", between="group") 1397 1398 return welch_results["F"].values[0], welch_results["p-unc"].values[0]
Perform a Welch's ANOVA analysis.
This function calculates group values by aggregating specified columns (testes_col) via the comb method and then conducts a Welch's ANOVA. This approach is ideal for comparing group means when data exhibits unequal variances across groups.
Parameters
data : dict of pd.DataFrame Dictionary where keys are group names and values are DataFrames containing the data.
testes_col : str or list of str
Column name(s) from which the group values are derived. If a list is provided, columns
will be combined based on the comb operation.
comb : str, optional
Operation used to combine multiple columns if testes_col is a list. Options include:
'' : multiplication
'+' : addition
'': exponentiation
'-' : subtraction
'/' : division
Default is ''.
Returns
F : float F-statistic from Welch's ANOVA.
p_val : float Uncorrected p-value from Welch's ANOVA, testing for significant differences between groups.
Notes
- If
testes_colis a single string, no combination is performed, and the group values are taken directly from that column. - Welch's ANOVA is used as it accounts for unequal variances between groups.
- The
df.melt()method is used to reshape the data, allowing the ANOVA to be applied to all groups.
Examples
>>> welch_F, welch_p = self.aov(data, testes_col=['col1', 'col2'], comb='+')
>>> print(f"Welch's ANOVA F-statistic: {welch_F}, p-value: {welch_p}")
def
post_aov(self, data, testes_col, comb: str = '*'):
1400 def post_aov(self, data, testes_col, comb: str = "*"): 1401 """ 1402 Perform Welch's ANOVA and pairwise Welch's t-tests on grouped data. 1403 1404 This method first conducts a Welch's ANOVA to detect significant differences 1405 in group means. It then performs pairwise Welch's t-tests across all group 1406 combinations to identify specific differences. All p-values are adjusted using 1407 the Bonferroni correction to account for multiple comparisons. 1408 1409 Parameters 1410 ---------- 1411 data : dict of pd.DataFrame 1412 Dictionary where keys are group names and values are DataFrames containing the data. 1413 1414 testes_col : str or list of str 1415 Column name(s) from which the group values are derived. If a list is provided, 1416 columns will be combined according to the `comb` operation. 1417 1418 comb : str, optional 1419 Operation used to combine multiple columns if `testes_col` is a list. Options include: 1420 '*' : multiplication 1421 '+' : addition 1422 '**': exponentiation 1423 '-' : subtraction 1424 '/' : division 1425 Default is '*'. 1426 1427 Returns 1428 ------- 1429 p_val : float 1430 Uncorrected p-value from the Welch's ANOVA. 1431 1432 final_results : dict 1433 Dictionary containing results of pairwise Welch's t-tests with keys: 1434 'group1' : list of first group names in each comparison 1435 'group2' : list of second group names in each comparison 1436 'stat' : list of t-statistics for each comparison 1437 'p_val' : list of uncorrected p-values for each comparison 1438 'adj_p_val' : list of Bonferroni-adjusted p-values for multiple comparisons 1439 """ 1440 1441 p_val = self.aov(data=data, testes_col=testes_col, comb=comb)[1] 1442 1443 pairs = list(combinations(data, 2)) 1444 final_results = { 1445 "group1": [], 1446 "group2": [], 1447 "stat": [], 1448 "p_val": [], 1449 "adj_p_val": [], 1450 } 1451 1452 for group1, group2 in pairs: 1453 if isinstance(testes_col, str): 1454 g1 = data[group1]["values"][testes_col] 1455 elif isinstance(testes_col, list): 1456 g1 = [1] * len(data[group1]["values"][testes_col[0]]) 1457 for t in testes_col: 1458 if comb == "*": 1459 g1 = [a * b for a, b in zip(g1, data[group1]["values"][t])] 1460 elif comb == "+": 1461 g1 = [a + b for a, b in zip(g1, data[group1]["values"][t])] 1462 elif comb == "**": 1463 g1 = [a**b for a, b in zip(g1, data[group1]["values"][t])] 1464 elif comb == "-": 1465 g1 = [a - b for a, b in zip(g1, data[group1]["values"][t])] 1466 elif comb == "/": 1467 g1 = [a / b for a, b in zip(g1, data[group1]["values"][t])] 1468 1469 if isinstance(testes_col, str): 1470 g2 = data[group2]["values"][testes_col] 1471 elif isinstance(testes_col, list): 1472 g2 = [1] * len(data[group2]["values"][testes_col[0]]) 1473 for t in testes_col: 1474 if comb == "*": 1475 g2 = [a * b for a, b in zip(g2, data[group2]["values"][t])] 1476 elif comb == "+": 1477 g2 = [a + b for a, b in zip(g2, data[group2]["values"][t])] 1478 elif comb == "**": 1479 g2 = [a**b for a, b in zip(g2, data[group2]["values"][t])] 1480 elif comb == "-": 1481 g2 = [a - b for a, b in zip(g2, data[group2]["values"][t])] 1482 elif comb == "/": 1483 g2 = [a / b for a, b in zip(g2, data[group2]["values"][t])] 1484 1485 stat, p_val = stats.ttest_ind( 1486 g1, g2, alternative="two-sided", equal_var=False 1487 ) 1488 g = sorted([group1, group2]) 1489 final_results["group1"].append(g[0]) 1490 final_results["group2"].append(g[1]) 1491 final_results["stat"].append(stat) 1492 final_results["p_val"].append(p_val) 1493 adj = p_val * len(pairs) 1494 if adj > 1: 1495 final_results["adj_p_val"].append(1) 1496 else: 1497 final_results["adj_p_val"].append(adj) 1498 1499 return p_val, final_results
Perform Welch's ANOVA and pairwise Welch's t-tests on grouped data.
This method first conducts a Welch's ANOVA to detect significant differences in group means. It then performs pairwise Welch's t-tests across all group combinations to identify specific differences. All p-values are adjusted using the Bonferroni correction to account for multiple comparisons.
Parameters
data : dict of pd.DataFrame Dictionary where keys are group names and values are DataFrames containing the data.
testes_col : str or list of str
Column name(s) from which the group values are derived. If a list is provided,
columns will be combined according to the comb operation.
comb : str, optional
Operation used to combine multiple columns if testes_col is a list. Options include:
'' : multiplication
'+' : addition
'': exponentiation
'-' : subtraction
'/' : division
Default is ''.
Returns
p_val : float Uncorrected p-value from the Welch's ANOVA.
final_results : dict Dictionary containing results of pairwise Welch's t-tests with keys: 'group1' : list of first group names in each comparison 'group2' : list of second group names in each comparison 'stat' : list of t-statistics for each comparison 'p_val' : list of uncorrected p-values for each comparison 'adj_p_val' : list of Bonferroni-adjusted p-values for multiple comparisons
def
ks_percentiles(self, input_hist):
1501 def ks_percentiles(self, input_hist): 1502 """ 1503 Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups. 1504 1505 This method extracts the percentile levels and computes the average value for 1506 each percentile to obtain a lower-dimensional representation of the data, thereby 1507 reducing the Big Data scale problem for each group. Using these metrics, it reconstructs 1508 the underlying empirical distributions to evaluate both structural proportions and scale. 1509 1510 To further mitigate the large sample size problem ("curse of Big Data") where inflating 1511 pixel counts yields artificially significant results, a controlled downsampling (resampling) 1512 is applied to standardize the sample sizes across groups. 1513 1514 A two-sample Kolmogorov-Smirnov (KS) test is then performed for every possible pair 1515 of groups to detect differences in distribution shapes. Finally, p-values are adjusted 1516 using the Bonferroni correction method to account for multiple comparisons and control 1517 the family-wise error rate. 1518 1519 Parameters 1520 ---------- 1521 input_hist : dict 1522 A nested dictionary where keys are group names. Each group must contain 1523 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1524 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1525 1526 Returns 1527 ------- 1528 final_results : dict 1529 A dictionary containing the results of the pairwise comparisons with keys: 1530 - 'group1': list of the first group names in the pairs. 1531 - 'group2': list of the second group names in the pairs. 1532 - 'K-S': list of Kolmogorov-Smirnov test statistics. 1533 - 'p_val': list of unadjusted p-values. 1534 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1535 1536 Example 1537 ------- 1538 >>> results = self.ks_percentiles(input_hist) 1539 """ 1540 1541 ks_data = {} 1542 1543 for d in input_hist.keys(): 1544 tmp_dic = {} 1545 1546 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["avg"]): 1547 tmp_dic[f"p{n+1}"] = c 1548 1549 ks_data[d] = tmp_dic 1550 1551 df_cleaned = pd.DataFrame(ks_data).T 1552 1553 pairs = list(combinations(df_cleaned.index, 2)) 1554 1555 final_results = { 1556 "group1": [], 1557 "group2": [], 1558 "K-S": [], 1559 "p_val": [], 1560 "adj_p_val": [], 1561 } 1562 1563 for group1, group2 in pairs: 1564 1565 g = sorted([group1, group2]) 1566 1567 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1568 1569 res = stats.ks_2samp(table_pair.iloc[:, 0], table_pair.iloc[:, 1]) 1570 1571 final_results["group1"].append(g[0]) 1572 final_results["group2"].append(g[1]) 1573 final_results["K-S"].append(res.statistic) 1574 final_results["p_val"].append(res.pvalue) 1575 adj = res.pvalue * len(pairs) 1576 if adj > 1: 1577 final_results["adj_p_val"].append(1) 1578 else: 1579 final_results["adj_p_val"].append(adj) 1580 1581 return final_results
Perform pairwise Kolmogorov-Smirnov (KS) tests on percentile data across all groups.
This method extracts the percentile levels and computes the average value for each percentile to obtain a lower-dimensional representation of the data, thereby reducing the Big Data scale problem for each group. Using these metrics, it reconstructs the underlying empirical distributions to evaluate both structural proportions and scale.
To further mitigate the large sample size problem ("curse of Big Data") where inflating pixel counts yields artificially significant results, a controlled downsampling (resampling) is applied to standardize the sample sizes across groups.
A two-sample Kolmogorov-Smirnov (KS) test is then performed for every possible pair of groups to detect differences in distribution shapes. Finally, p-values are adjusted using the Bonferroni correction method to account for multiple comparisons and control the family-wise error rate.
Parameters
input_hist : dict A nested dictionary where keys are group names. Each group must contain the following structure: input_hist[group]["percentiles"]["mutual"]["n"], which holds an iterable (e.g., list or Series) of counts per percentile/bin.
Returns
final_results : dict A dictionary containing the results of the pairwise comparisons with keys: - 'group1': list of the first group names in the pairs. - 'group2': list of the second group names in the pairs. - 'K-S': list of Kolmogorov-Smirnov test statistics. - 'p_val': list of unadjusted p-values. - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0).
Example
>>> results = self.ks_percentiles(input_hist)
def
fisher_percentiles(self, input_hist):
1583 def fisher_percentiles(self, input_hist): 1584 """ 1585 Perform pairwise Fisher's exact tests on percentile data across all groups. 1586 1587 This method extracts the raw counts (N) for each percentile bin across all 1588 groups to construct a contingency table representation of the data. By utilizing 1589 the discrete frequency counts per bin rather than continuous average values, it 1590 evaluates both structural distribution proportions and sample size scaling 1591 differences simultaneously. 1592 1593 An exact testing approach is applied to every unique pair of groups by extracting 1594 their corresponding sub-tables. For each pair, a Fisher's exact test (or its 1595 extension for larger contingency tables) is performed to detect statistically 1596 significant deviations in distribution profiles. 1597 1598 Finally, p-values are manually adjusted using the Bonferroni correction method 1599 by multiplying the raw p-values by the total number of comparisons to control 1600 the family-wise error rate across multiple pair-wise tests. 1601 the family-wise error rate. 1602 1603 Parameters 1604 ---------- 1605 input_hist : dict 1606 A nested dictionary where keys are group names. Each group must contain 1607 the following structure: input_hist[group]["percentiles"]["mutual"]["n"], 1608 which holds an iterable (e.g., list or Series) of counts per percentile/bin. 1609 1610 Returns 1611 ------- 1612 final_results : dict 1613 A dictionary containing the results of the pairwise comparisons with keys: 1614 - 'group1': list of the first group names in the pairs. 1615 - 'group2': list of the second group names in the pairs. 1616 - 'fish': list of Fisher's exact test statistics. 1617 - 'p_val': list of unadjusted p-values. 1618 - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0). 1619 1620 Example 1621 ------- 1622 >>> results = self.fisher_percentiles(input_hist) 1623 """ 1624 1625 fish_data = {} 1626 1627 for d in input_hist.keys(): 1628 tmp_dic = {} 1629 1630 for n, c in enumerate(input_hist[d]["percentiles"]["mutual"]["n"]): 1631 tmp_dic[f"p{n+1}"] = c 1632 1633 fish_data[d] = tmp_dic 1634 1635 df_cleaned = pd.DataFrame(fish_data).T 1636 1637 pairs = list(combinations(df_cleaned.index, 2)) 1638 1639 final_results = { 1640 "group1": [], 1641 "group2": [], 1642 "fish": [], 1643 "p_val": [], 1644 "adj_p_val": [], 1645 } 1646 1647 for group1, group2 in pairs: 1648 1649 g = sorted([group1, group2]) 1650 1651 table_pair = pd.DataFrame(df_cleaned).T[[group1, group2]].copy() 1652 1653 res = stats.fisher_exact(table_pair) 1654 1655 final_results["group1"].append(g[0]) 1656 final_results["group2"].append(g[1]) 1657 final_results["fish"].append(res.statistic) 1658 final_results["p_val"].append(res.pvalue) 1659 adj = res.pvalue * len(pairs) 1660 if adj > 1: 1661 final_results["adj_p_val"].append(1) 1662 else: 1663 final_results["adj_p_val"].append(adj) 1664 1665 return final_results
Perform pairwise Fisher's exact tests on percentile data across all groups.
This method extracts the raw counts (N) for each percentile bin across all groups to construct a contingency table representation of the data. By utilizing the discrete frequency counts per bin rather than continuous average values, it evaluates both structural distribution proportions and sample size scaling differences simultaneously.
An exact testing approach is applied to every unique pair of groups by extracting their corresponding sub-tables. For each pair, a Fisher's exact test (or its extension for larger contingency tables) is performed to detect statistically significant deviations in distribution profiles.
Finally, p-values are manually adjusted using the Bonferroni correction method by multiplying the raw p-values by the total number of comparisons to control the family-wise error rate across multiple pair-wise tests. the family-wise error rate.
Parameters
input_hist : dict A nested dictionary where keys are group names. Each group must contain the following structure: input_hist[group]["percentiles"]["mutual"]["n"], which holds an iterable (e.g., list or Series) of counts per percentile/bin.
Returns
final_results : dict A dictionary containing the results of the pairwise comparisons with keys: - 'group1': list of the first group names in the pairs. - 'group2': list of the second group names in the pairs. - 'fish': list of Fisher's exact test statistics. - 'p_val': list of unadjusted p-values. - 'adj_p_val': list of Bonferroni-adjusted p-values (capped at 1.0).
Example
>>> results = self.fisher_percentiles(input_hist)
def
to_wasserstein_distance(self, data):
1667 def to_wasserstein_distance(self, data): 1668 """ 1669 Calculate scaled pairwise Wasserstein distances for grouped distributions. 1670 1671 This method computes the 1D Wasserstein distance (Earth Mover's Distance) 1672 between all possible combinations of groups in the provided dataset. 1673 Before calculating the distance, the standardized frequencies are scaled 1674 by a factor representing the average total count (sample size) of the 1675 two compared groups. 1676 1677 Parameters 1678 ---------- 1679 data : dict 1680 A nested dictionary where keys are group names. For each group, the 1681 method expects the following internal data structure: 1682 - `data[group_name]['percentiles']['mutual']['n']` : list-like 1683 Absolute counts or sample sizes for the distribution. 1684 - `data[group_name]['percentiles']['mutual']['n_standarized']` : list-like 1685 Standardized frequencies or probabilities to be compared. 1686 1687 Returns 1688 ------- 1689 final_results : dict 1690 A dictionary containing the results of the pairwise distance calculations: 1691 - 'group1' : list of str 1692 The name of the first group in the comparison. 1693 - 'group2' : list of str 1694 The name of the second group in the comparison. 1695 - 'wasserstein_distance' : list of float 1696 The computed scaled Wasserstein distance for each pair. 1697 """ 1698 1699 pairs = list(combinations(data.keys(), 2)) 1700 1701 final_results = {"group1": [], "group2": [], "wasserstein_distance": []} 1702 1703 for group1, group2 in pairs: 1704 1705 factor = ( 1706 sum(data[group1]["percentiles"]["mutual"]["n"]) 1707 + sum(data[group2]["percentiles"]["mutual"]["n"]) 1708 ) / 2 1709 1710 dist = wasserstein_distance( 1711 [ 1712 x * factor 1713 for x in data[group1]["percentiles"]["mutual"]["n_standarized"] 1714 ], 1715 [ 1716 x * factor 1717 for x in data[group2]["percentiles"]["mutual"]["n_standarized"] 1718 ], 1719 ) 1720 1721 g = sorted([group1, group2]) 1722 final_results["group1"].append(g[0]) 1723 final_results["group2"].append(g[1]) 1724 final_results["wasserstein_distance"].append(dist) 1725 1726 return final_results
Calculate scaled pairwise Wasserstein distances for grouped distributions.
This method computes the 1D Wasserstein distance (Earth Mover's Distance) between all possible combinations of groups in the provided dataset. Before calculating the distance, the standardized frequencies are scaled by a factor representing the average total count (sample size) of the two compared groups.
Parameters
data : dict
A nested dictionary where keys are group names. For each group, the
method expects the following internal data structure:
- data[group_name]['percentiles']['mutual']['n'] : list-like
Absolute counts or sample sizes for the distribution.
- data[group_name]['percentiles']['mutual']['n_standarized'] : list-like
Standardized frequencies or probabilities to be compared.
Returns
final_results : dict A dictionary containing the results of the pairwise distance calculations: - 'group1' : list of str The name of the first group in the comparison. - 'group2' : list of str The name of the second group in the comparison. - 'wasserstein_distance' : list of float The computed scaled Wasserstein distance for each pair.
def
to_fold_change(self, data, tested_value):
1728 def to_fold_change(self, data, tested_value): 1729 """ 1730 Calculate the Fold Change (FC) between all permutations of groups. 1731 1732 This method computes the ratio of the mean values of a specified feature 1733 (`tested_value`) for every directed pair of groups. Because permutations 1734 are used, the calculation is directional (i.e., both Group A / Group B 1735 and Group B / Group A are computed). 1736 1737 Parameters 1738 ---------- 1739 data : dict 1740 A nested dictionary where keys are group names. For each group, the 1741 method expects the following internal structure: 1742 - `data[group_name]['values'][tested_value]` : array-like 1743 Numeric values used to compute the mean for the group. 1744 1745 tested_value : str 1746 The specific key or column name within the 'values' dictionary 1747 indicating which feature's fold change should be calculated. 1748 1749 Returns 1750 ------- 1751 final_results : dict 1752 A dictionary containing the results of the pairwise fold change calculations: 1753 - 'group1' : list of str 1754 The name of the numerator group in the comparison. 1755 - 'group2' : list of str 1756 The name of the denominator group in the comparison. 1757 - 'FC' : list of float 1758 The calculated fold change (mean of group1 / mean of group2). 1759 """ 1760 1761 pairs = list(permutations(data.keys(), 2)) 1762 1763 final_results = {"group1": [], "group2": [], "FC": []} 1764 1765 values = [] 1766 for group1, group2 in pairs: 1767 1768 values = values + data[group1]["values"][tested_value] 1769 values = values + data[group2]["values"][tested_value] 1770 1771 values_min = min([x for x in values if x > 0]) 1772 values_min = values_min / 2 1773 1774 for group1, group2 in pairs: 1775 1776 g1 = np.mean(data[group1]["values"][tested_value]) 1777 g2 = np.mean(data[group2]["values"][tested_value]) 1778 1779 if g1 == 0: 1780 g1 = g1 + values_min 1781 1782 if g2 == 0: 1783 g2 = g2 + values_min 1784 1785 fc = g1 / g2 1786 1787 final_results["group1"].append(group1) 1788 final_results["group2"].append(group2) 1789 final_results["FC"].append(fc) 1790 1791 return final_results
Calculate the Fold Change (FC) between all permutations of groups.
This method computes the ratio of the mean values of a specified feature
(tested_value) for every directed pair of groups. Because permutations
are used, the calculation is directional (i.e., both Group A / Group B
and Group B / Group A are computed).
Parameters
data : dict
A nested dictionary where keys are group names. For each group, the
method expects the following internal structure:
- data[group_name]['values'][tested_value] : array-like
Numeric values used to compute the mean for the group.
tested_value : str The specific key or column name within the 'values' dictionary indicating which feature's fold change should be calculated.
Returns
final_results : dict A dictionary containing the results of the pairwise fold change calculations: - 'group1' : list of str The name of the numerator group in the comparison. - 'group2' : list of str The name of the denominator group in the comparison. - 'FC' : list of float The calculated fold change (mean of group1 / mean of group2).
def
get_stats(self, data, tested_value):
1793 def get_stats(self, data, tested_value): 1794 """ 1795 Calculate and aggregate statistical metrics (ANOVA, Fisher's Exact, 1796 Kolmogorov-Smirnov, Fold Change, Wasserstein distance). 1797 1798 This method computes overall statistics and pairwise comparisons for grouped data. 1799 To properly capture both structural proportions and total count variations across 1800 percentiles while avoiding the curse of Big Data, it runs two distinct tests: 1801 1. Fisher's exact test on discrete percentile counts to evaluate absolute scale 1802 and profile differences. 1803 2. Two-sample Kolmogorov-Smirnov (KS) test on reconstructed empirical 1804 distributions to evaluate discrepancies in distribution shapes. 1805 1806 Additionally, it calculates the Fold Change (FC) and evaluates Wasserstein 1807 distances. If the average number of replicates per group is at least 3, 1808 it conducts Welch's ANOVA. The input dictionary is modified in-place to 1809 include a new 'statistics' key containing all results. 1810 1811 Parameters 1812 ---------- 1813 data : dict 1814 A nested dictionary where keys are group names. Each group's dictionary 1815 must contain the structure `['values']['replication']` to verify sample sizes, 1816 along with the necessary data structures required by downstream statistical methods. 1817 1818 tested_value : str 1819 The key or column name representing the specific variable to evaluate 1820 (e.g., used for ANOVA and Fold Change calculations). 1821 1822 Returns 1823 ------- 1824 data : dict 1825 The original input dictionary, extended with a new `data['statistics']` key 1826 that houses the computed statistical results, including `percintiles_fish` 1827 and `percintiles_ks`. 1828 1829 Example 1830 ------- 1831 stats = self.get_stats( 1832 data, 1833 tested_value='n', 1834 ) 1835 """ 1836 1837 # parametric selected value 1838 sum_k = 0 1839 n = 0 1840 for k in data.keys(): 1841 if k != "statistics": 1842 n += 1 1843 sum_k += len(data[k]["values"]["replication"]) 1844 1845 sum_k = sum_k / n 1846 1847 if sum_k >= 3: 1848 pk, dfk = self.post_aov(data, testes_col=tested_value) 1849 1850 # fish 1851 fish = self.fisher_percentiles(data) 1852 1853 # K_S 1854 ks = self.ks_percentiles(data) 1855 1856 dw = self.to_wasserstein_distance(data) 1857 1858 fc = self.to_fold_change(data, tested_value) 1859 1860 data["statistics"] = {} 1861 1862 data["statistics"]["percintiles_fish"] = fish 1863 1864 data["statistics"]["percintiles_ks"] = ks 1865 1866 if sum_k >= 3: 1867 data["statistics"]["ANOVA"] = {} 1868 1869 data["statistics"]["ANOVA"]["p_value"] = pk 1870 data["statistics"]["ANOVA"]["pair-comparison"] = dfk 1871 else: 1872 import warnings 1873 1874 warnings.warn( 1875 f"Insufficient replicates for statistical analysis. " 1876 f"At least 3 replicates per group (3 vs 3) are required. " 1877 f"The average number of samples per probe in this dataset was {n}.", 1878 RuntimeWarning, 1879 ) 1880 1881 data["statistics"]["FC"] = fc 1882 1883 data["statistics"]["wasserstein_distance"] = dw 1884 1885 data["statistics"]["tested_value"] = tested_value 1886 1887 return data
Calculate and aggregate statistical metrics (ANOVA, Fisher's Exact, Kolmogorov-Smirnov, Fold Change, Wasserstein distance).
This method computes overall statistics and pairwise comparisons for grouped data. To properly capture both structural proportions and total count variations across percentiles while avoiding the curse of Big Data, it runs two distinct tests:
- Fisher's exact test on discrete percentile counts to evaluate absolute scale and profile differences.
- Two-sample Kolmogorov-Smirnov (KS) test on reconstructed empirical distributions to evaluate discrepancies in distribution shapes.
Additionally, it calculates the Fold Change (FC) and evaluates Wasserstein distances. If the average number of replicates per group is at least 3, it conducts Welch's ANOVA. The input dictionary is modified in-place to include a new 'statistics' key containing all results.
Parameters
data : dict
A nested dictionary where keys are group names. Each group's dictionary
must contain the structure ['values']['replication'] to verify sample sizes,
along with the necessary data structures required by downstream statistical methods.
tested_value : str The key or column name representing the specific variable to evaluate (e.g., used for ANOVA and Fold Change calculations).
Returns
data : dict
The original input dictionary, extended with a new data['statistics'] key
that houses the computed statistical results, including percintiles_fish
and percintiles_ks.
Example
stats = self.get_stats( data, tested_value='n', )
def
hist_compare_plot(self, data, queue=None, p_adj: bool = True, txt_size: int = 20):
1889 def hist_compare_plot( 1890 self, data, queue=None, p_adj: bool = True, txt_size: int = 20 1891 ): 1892 """ 1893 Generate comparative histograms and display results of statistical tests 1894 (ANOVA / T-test [if min 3 vs. 3 comparison available], Kolmogorov-Smirnov, Fisher percentiles) 1895 and statistics (FC, Wasserstein distance). 1896 1897 1898 Parameters 1899 ---------- 1900 data : dict 1901 Dictionary where keys are group names and values are containing histogram data. 1902 Each DataFrame should include the column specified by `tested_value`. 1903 1904 queue : list of str or None 1905 Defines the order of groups to be plotted. 1906 1907 p_adj : bool, optional 1908 If True, applies Bonferroni correction for multiple comparisons (default is True). 1909 1910 txt_size : int, optional 1911 Font size for text annotations in the plot (default is 20). 1912 1913 Returns 1914 ------- 1915 fig : matplotlib.figure.Figure 1916 Matplotlib figure object containing the generated histograms and statistical test results. 1917 1918 Example 1919 ------- 1920 fig = self.hist_compare_plot( 1921 data, 1922 queue=['group1', 'group2', 'group3'], 1923 p_adj=True, 1924 txt_size=18 1925 ) 1926 plt.show() 1927 """ 1928 1929 if queue is None: 1930 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1931 1932 if sorted(queue) != [x for x in sorted(data.keys()) if x != "statistics"]: 1933 print( 1934 "\n Wrong queue provided! The queue will be sorted with default settings!" 1935 ) 1936 queue = [x for x in sorted(data.keys()) if x != "statistics"] 1937 1938 # parametric selected value 1939 tested_value = data["statistics"]["tested_value"] 1940 1941 ############################################################################## 1942 1943 standarized_max, standarized_min, value_max, value_min = [], [], [], [] 1944 for d in queue: 1945 standarized_max.append( 1946 max(data[d]["percentiles"]["mutual"]["n_standarized"]) 1947 ) 1948 standarized_min.append( 1949 min(data[d]["percentiles"]["mutual"]["n_standarized"]) 1950 ) 1951 value_max.append(max(data[d]["percentiles"]["mutual"][tested_value])) 1952 value_min.append(min(data[d]["percentiles"]["mutual"][tested_value])) 1953 1954 num_columns = len(queue) + 1 1955 1956 fig, axs = plt.subplots( 1957 3, 1958 num_columns, 1959 figsize=(8 * num_columns, 10), 1960 gridspec_kw={"width_ratios": [1] * len(queue) + [0.5], "wspace": 0.05}, 1961 ) 1962 1963 for i, d in enumerate(queue): 1964 tmp_data = data[d]["percentiles"]["mutual"] 1965 1966 axs[0, i].bar( 1967 [str(n) for n in range(len(tmp_data["n_standarized"]))], 1968 tmp_data["n_standarized"], 1969 width=0.95, 1970 color="gold", 1971 ) 1972 1973 # line 1974 n_groups = len(data[d]["percentiles"]["replications"].keys()) 1975 colors = plt.cm.OrRd(np.linspace(0.3, 0.9, n_groups)) 1976 1977 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 1978 1979 color = colors[ix] 1980 1981 y = data[d]["percentiles"]["replications"][dn]["n_standarized"] 1982 x = np.arange(len(y)) 1983 1984 axs[0, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 1985 1986 axs[0, i].plot( 1987 x, 1988 y, 1989 color=color, 1990 linewidth=1, 1991 marker="o", 1992 ) 1993 1994 axs[0, i].set_ylim( 1995 min(standarized_min) * 0.9995, max(standarized_max) * 1.0005 1996 ) 1997 1998 if i == 0: 1999 axs[0, i].set_ylabel("Standarized\nfrequency", fontsize=txt_size) 2000 else: 2001 axs[0, i].set_yticks([]) 2002 2003 axs[0, i].set_xticks([]) 2004 axs[0, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2005 2006 axs[1, i].bar( 2007 [str(n) for n in range(len(tmp_data[tested_value]))], 2008 tmp_data[tested_value], 2009 width=0.95, 2010 color="orange", 2011 ) 2012 2013 # line 2014 2015 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2016 2017 color = colors[ix] 2018 2019 y = data[d]["percentiles"]["replications"][dn][tested_value] 2020 x = np.arange(len(y)) 2021 2022 axs[1, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2023 2024 axs[1, i].plot( 2025 x, 2026 y, 2027 color=color, 2028 linewidth=1, 2029 marker="o", 2030 ) 2031 2032 mean_value = np.mean(data[d]["values"][tested_value]) 2033 axs[1, i].axhline(y=mean_value, color="red", linestyle="--") 2034 2035 axs[1, i].set_ylim(min(value_min) * 0.9995, max(value_max) * 1.0005) 2036 2037 if i == 0: 2038 axs[1, i].set_ylabel(f"Normalized\n{tested_value}", fontsize=txt_size) 2039 else: 2040 axs[1, i].set_yticks([]) 2041 2042 axs[1, i].set_xticks([]) 2043 axs[1, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2044 2045 axs[2, i].bar( 2046 [str(n) for n in range(len(tmp_data["n_standarized"]))], 2047 [ 2048 a * b 2049 for a, b in zip(tmp_data[tested_value], tmp_data["n_standarized"]) 2050 ], 2051 width=0.95, 2052 color="goldenrod", 2053 ) 2054 2055 # line 2056 for ix, dn in enumerate(data[d]["percentiles"]["replications"].keys()): 2057 2058 color = colors[ix] 2059 2060 y = [ 2061 a * b 2062 for a, b in zip( 2063 data[d]["percentiles"]["replications"][dn][tested_value], 2064 data[d]["percentiles"]["replications"][dn]["n_standarized"], 2065 ) 2066 ] 2067 x = np.arange(len(y)) 2068 2069 axs[2, i].plot(x, y, color="black", linewidth=1.5, label="_nolegend_") 2070 2071 axs[2, i].plot( 2072 x, 2073 y, 2074 color=color, 2075 linewidth=1, 2076 marker="o", 2077 ) 2078 2079 mean_value = np.mean( 2080 data[d]["values"][data["statistics"]["tested_value"]] 2081 ) * np.mean(tmp_data["n_standarized"]) 2082 2083 axs[2, i].axhline(y=mean_value, color="red", linestyle="--") 2084 2085 axs[2, i].set_ylim( 2086 (min(standarized_min) * min(value_min)) * 0.9995, 2087 (max(standarized_max) * max(value_max) * 1.0005), 2088 ) 2089 axs[2, i].set_xlabel(d, fontsize=txt_size) 2090 2091 if i == 0: 2092 axs[2, i].set_ylabel( 2093 f"Standarized\nnorm_{tested_value}", fontsize=txt_size 2094 ) 2095 else: 2096 axs[2, i].set_yticks([]) 2097 2098 axs[2, i].set_xticks([]) 2099 axs[2, i].tick_params(axis="y", labelsize=txt_size * 0.7) 2100 2101 # statistics 2102 2103 # ANOVA / t-test 2104 2105 if "ANOVA" in data["statistics"].keys(): 2106 pk = data["statistics"]["ANOVA"]["p_value"] 2107 dfk = data["statistics"]["ANOVA"]["pair-comparison"] 2108 dfk = pd.DataFrame(dfk) 2109 2110 dfk = dfk.sort_values( 2111 by=["group1", "group2"], 2112 key=lambda col: [ 2113 queue.index(val) if val in queue else -1 for val in col 2114 ], 2115 ).reset_index(drop=True) 2116 2117 sign = "ns" 2118 if float(self.round_to_scientific_notation(pk)) < 0.001: 2119 sign = "***" 2120 elif float(self.round_to_scientific_notation(pk)) < 0.01: 2121 sign = "**" 2122 elif float(self.round_to_scientific_notation(pk)) < 0.05: 2123 sign = "*" 2124 2125 text = f"Test Welch's ANOVA\non '{tested_value}' values\np-value: {self.round_to_scientific_notation(pk)} - {sign}\n" 2126 2127 if p_adj == True: 2128 for i in range(len(dfk["group1"])): 2129 sign = "ns" 2130 if dfk["adj_p_val"][i] < 0.001: 2131 sign = "***" 2132 elif dfk["adj_p_val"][i] < 0.01: 2133 sign = "**" 2134 elif dfk["adj_p_val"][i] < 0.05: 2135 sign = "*" 2136 2137 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['adj_p_val'][i])} - {sign}\n" 2138 else: 2139 for i in range(len(dfk["group1"])): 2140 sign = "ns" 2141 if dfk["p_val"][i] < 0.001: 2142 sign = "***" 2143 elif dfk["p_val"][i] < 0.01: 2144 sign = "**" 2145 elif dfk["p_val"][i] < 0.05: 2146 sign = "*" 2147 2148 text += f"{dfk['group1'][i]} vs. {dfk['group2'][i]}\np-value: {self.round_to_scientific_notation(dfk['p_val'][i])} - {sign}\n" 2149 2150 axs[2, -1].text( 2151 0.5, 2152 0.5, 2153 text, 2154 ha="center", 2155 va="center", 2156 fontsize=txt_size * 0.7, 2157 wrap=True, 2158 ) 2159 axs[2, -1].set_axis_off() 2160 else: 2161 axs[2, -1].set_axis_off() 2162 2163 # FC / Distance 2164 2165 ranking_FC = pd.DataFrame(data["statistics"]["FC"]) 2166 2167 ranking_dw = pd.DataFrame(data["statistics"]["wasserstein_distance"]) 2168 2169 ranking_combined = pd.merge( 2170 ranking_FC, ranking_dw, on=["group1", "group2"], how="right" 2171 ) 2172 2173 ranking_combined = ranking_combined.sort_values( 2174 by=["group1", "group2"], 2175 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2176 ).reset_index(drop=True) 2177 2178 text = "FC / Wasserstein distance\n" 2179 for i in range(len(ranking_combined)): 2180 group1 = ranking_combined["group1"][i] 2181 group2 = ranking_combined["group2"][i] 2182 fc_val = ranking_combined["FC"][i] 2183 wasserstein_val = ranking_combined["wasserstein_distance"][i] 2184 2185 text += f"{group1} vs. {group2}\n {fc_val:.2f} | {wasserstein_val:.2f}\n" 2186 2187 axs[1, -1].text( 2188 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2189 ) 2190 axs[1, -1].set_axis_off() 2191 2192 # fish 2193 2194 fish = pd.DataFrame(data["statistics"]["percintiles_fish"]) 2195 2196 # K-S 2197 2198 ks = pd.DataFrame(data["statistics"]["percintiles_ks"]) 2199 2200 fish = fish.sort_values( 2201 by=["group1", "group2"], 2202 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2203 ).reset_index(drop=True) 2204 2205 ks = ks.sort_values( 2206 by=["group1", "group2"], 2207 key=lambda col: [queue.index(val) if val in queue else -1 for val in col], 2208 ).reset_index(drop=True) 2209 2210 text = f"Tests\nFisher's exact / Kolmogorov-Smirnov\n" 2211 2212 if p_adj == True: 2213 for i in range(len(fish["group1"])): 2214 sign1 = "ns" 2215 if fish["adj_p_val"][i] < 0.001: 2216 sign1 = "***" 2217 elif fish["adj_p_val"][i] < 0.01: 2218 sign1 = "**" 2219 elif fish["adj_p_val"][i] < 0.05: 2220 sign1 = "*" 2221 2222 sign2 = "ns" 2223 if ks["adj_p_val"][i] < 0.001: 2224 sign2 = "***" 2225 elif ks["adj_p_val"][i] < 0.01: 2226 sign2 = "**" 2227 elif ks["adj_p_val"][i] < 0.05: 2228 sign2 = "*" 2229 2230 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2231 2232 else: 2233 for i in range(len(fish["group1"])): 2234 sign1 = "ns" 2235 if fish["p_val"][i] < 0.001: 2236 sign1 = "***" 2237 elif fish["p_val"][i] < 0.01: 2238 sign1 = "**" 2239 elif fish["p_val"][i] < 0.05: 2240 sign1 = "*" 2241 2242 sign2 = "ns" 2243 if ks["p_val"][i] < 0.001: 2244 sign2 = "***" 2245 elif ks["p_val"][i] < 0.01: 2246 sign2 = "**" 2247 elif ks["p_val"][i] < 0.05: 2248 sign2 = "*" 2249 2250 text += f"{fish['group1'][i]} vs. {fish['group2'][i]}\np-value: {sign1} / {sign2}\n" 2251 2252 axs[0, -1].text( 2253 0.5, 0.5, text, ha="center", va="center", fontsize=txt_size * 0.7, wrap=True 2254 ) 2255 axs[0, -1].set_axis_off() 2256 2257 plt.tight_layout() 2258 2259 if cfg._DISPLAY_MODE: 2260 plt.show() 2261 2262 return fig
Generate comparative histograms and display results of statistical tests (ANOVA / T-test [if min 3 vs. 3 comparison available], Kolmogorov-Smirnov, Fisher percentiles) and statistics (FC, Wasserstein distance).
Parameters
data : dict
Dictionary where keys are group names and values are containing histogram data.
Each DataFrame should include the column specified by tested_value.
queue : list of str or None Defines the order of groups to be plotted.
p_adj : bool, optional If True, applies Bonferroni correction for multiple comparisons (default is True).
txt_size : int, optional Font size for text annotations in the plot (default is 20).
Returns
fig : matplotlib.figure.Figure Matplotlib figure object containing the generated histograms and statistical test results.
Example
fig = self.hist_compare_plot( data, queue=['group1', 'group2', 'group3'], p_adj=True, txt_size=18 ) plt.show()