cfi_toolkit.CellFunctionality

  1import copy
  2import os
  3import pickle
  4import sys
  5
  6import pandas as pd
  7from tqdm import tqdm
  8
  9_old_stdout = sys.stdout
 10sys.stdout = open(os.devnull, "w")
 11
 12from gedspy import Analysis, Enrichment
 13
 14sys.stdout.close()
 15sys.stdout = _old_stdout
 16
 17
 18class CellFunCon:
 19    """
 20    A class to perform cell-type functional analysis and enrichment based on a JDtI-COMPsc objects.
 21
 22    This class provides methods to calculate marker genes for cell types, perform functional enrichment
 23    (GO, KEGG, REACTOME, STRING, IntAct), and compute cell-cell interaction networks.
 24    Projects can also be saved and loaded via pickle.
 25
 26    Attributes
 27    ----------
 28    jdti : object
 29        JDtI-COMPsc object containing normalized single-cell data.
 30
 31    cells_markers : pd.DataFrame or None
 32        DataFrame containing marker genes per cell type after calculation.
 33
 34    enr_full_info : Enrichment
 35        Enrichment object containing all genes available for enrichment analysis.
 36
 37    cells_enrichment : dict or None
 38        Dictionary storing enrichment results per cell type.
 39
 40    cells_connection : pd.DataFrame or None
 41        DataFrame storing calculated cell-cell interaction information.
 42
 43    mt_genes : bool
 44        Whether mitochondrial genes are included (default False).
 45
 46    ribo_genes : bool
 47        Whether ribosomal genes are included (default False).
 48    """
 49
 50    def __init__(self, jdti_object, mt_genes=False, ribo_genes=False):
 51        """
 52        Initializes the CellFunCon object with a COMPsc/JDTI object.
 53
 54        Parameters
 55        ----------
 56        jdti_object : object
 57            A COMPsc or JDTI object with normalized single-cell data.
 58
 59        mt_genes : bool
 60            Whether mitochondrial genes are included (default False).
 61
 62        ribo_genes : bool
 63            Whether ribosomal genes are included (default False).
 64        """
 65
 66        self.jdti = jdti_object
 67        """JDtI-COMPsc object containing normalized single-cell data."""
 68
 69        self.cells_markers = None
 70        """DataFrame containing marker genes per cell type after calculation."""
 71
 72        self.cells_connection = None
 73        """DataFrame storing calculated cell-cell interaction information."""
 74
 75        self.cells_enrichment = None
 76        """Dictionary storing enrichment results per cell type."""
 77
 78        self.mt_genes = mt_genes
 79        """Whether mitochondrial genes are included (default False)."""
 80
 81        self.ribo_genes = ribo_genes
 82        """Whether ribosomal genes are included (default False)."""
 83
 84        names = self.jdti.normalized_data.loc[
 85            self.jdti.normalized_data.select_dtypes(include="number").sum(axis=1) > 0
 86        ].index.tolist()
 87        names = list(set(names))
 88
 89        if self.mt_genes is False:
 90            names = [x for x in names if "MT-" not in x.upper()]
 91        if self.ribo_genes is False:
 92            names = [x for x in names if "RPS" != x[:3].upper()]
 93            names = [x for x in names if "RPL" != x[:3].upper()]
 94
 95        enr = Enrichment()
 96        enr.select_features(names)
 97
 98        self.enr_full_info = enr
 99        """Enrichment object containing all genes available for enrichment analysis."""
100
101    def save_project(self, filename):
102        """
103        Saves the current CellFunCon project as a pickle file.
104
105        Parameters
106        ----------
107        filename : str
108            Path to save the project (e.g., 'project_name').
109
110        Example
111        -------
112        >>> self.save_project('my_project')
113        """
114
115        with open(f"{filename}.psc", "wb") as f:
116            pickle.dump(self, f)
117        print(f"Project saved as {filename}")
118
119    @classmethod
120    def load_project(cls, filename):
121        """
122        Loads a previously saved CellFunCon project from a pickle file.
123
124        Parameters
125        ----------
126        filename : str
127            Path to the saved pickle file.
128
129        Returns
130        -------
131        CellFunCon
132            Loaded CellFunCon self.
133
134        Raises
135        ------
136        TypeError
137            If the loaded object is not a CellFunCon self.
138
139        ValueError
140            If the file is not a valid CellFunCon project file.
141
142        Example
143        -------
144        >>> self = CellFunCon.load_project('my_project.psc')
145        """
146
147        if ".psc" in filename:
148            with open(filename, "rb") as f:
149                obj = pickle.load(f)
150            if not isinstance(obj, cls):
151                raise TypeError("File does not include project.psc")
152            print(f"Project loaded from {filename}")
153            return obj
154        else:
155            raise ValueError("Project not belong to CellFunCon project data.")
156
157    def calculate_cells_markers(self, min_exp=0, min_pct=0.05, n_proc=10):
158        """
159        Calculates marker genes for each cell type based on expression thresholds.
160
161        Perform differential gene expression (DEG) analysis on gene expression data.
162
163        The function compares groups of cells or samples (defined by `entities` or
164        `sets`) using the Mann–Whitney U test. It computes p-values, adjusted
165        p-values, fold changes, standardized effect sizes, and other statistics.
166
167
168        Parameters
169        ----------
170        min_exp : float, optional
171            Minimum expression level to consider a gene (default 0).
172
173        min_pct : float, optional
174            Minimum fraction of cells expressing a gene (default 0.05).
175
176        n_proc : int, optional
177            Number of parallel processes to use (default 10).
178
179        Notes
180        -----
181        The results are stored in the `cells_markers` attribute.
182        """
183
184        self.jdti.calculate_difference_markers(
185            min_exp=min_exp, min_pct=min_pct, n_proc=n_proc, force=True
186        )
187
188        self.cells_markers = self.jdti.var_data
189
190    def enrich_cells_fucntionality(
191        self, p_value=0.05, adj=True, log_fc=0.1, top_max=500
192    ):
193        """
194        Performs functional enrichment analysis for each cell type based on marker genes.
195
196        Parameters
197        ----------
198        p_value : float
199            Maximum p-value for significant genes (default 0.05).
200
201        adj : bool
202            If True, the adjusted p-values are used to determine significant genes.
203            Adjusted p-values are calculated using the Benjamini–Hochberg false
204            discovery rate (FDR) correction. If False, raw p-values are used instead.
205
206        log_fc : float
207            Minimum log fold-change threshold for marker genes (default 0.1).
208
209        top_max : int
210            Maximum number of top marker genes per cell type to consider (default 500).
211
212        Raises
213        ------
214        ValueError
215            If `cells_markers` is not defined.
216
217        Notes
218        -----
219        This method populates `cells_enrichment` with results for GO-TERM, KEGG, REACTOME,
220        STRING, IntAct, and specificity analyses.
221        """
222
223        if isinstance(self.cells_markers, pd.DataFrame):
224
225            markers = self.cells_markers
226            cells = set(markers["valid_group"])
227
228            data_dict = {}
229
230            max_c = len(cells)
231            for n, c in enumerate(cells):
232                print(f"\nAnalysis {n+1} of {max_c} cells --> {c} \n")
233
234                if adj:
235                    tmp = markers[
236                        (markers["valid_group"] == c)
237                        & (markers["adj_pval"] <= p_value)
238                        & (markers["log(FC)"] > log_fc)
239                    ]
240                    names = list(set(tmp["feature"]))
241
242                    tmp = tmp[tmp["feature"].isin(names)]
243
244                else:
245                    tmp = markers[
246                        (markers["valid_group"] == c)
247                        & (markers["p_val"] <= p_value)
248                        & (markers["log(FC)"] > log_fc)
249                    ]
250                    names = list(set(tmp["feature"]))
251
252                    tmp = tmp[tmp["feature"].isin(names)]
253
254                tmp = tmp.sort_values("esm", ascending=False).head(top_max)
255
256                if len(tmp.index) > 0:
257                    data_dict[c] = {}
258                    enr = copy.copy(self.enr_full_info)
259                    enr.genome = enr.genome[
260                        enr.genome["found_names"].isin(list(set(tmp["feature"])))
261                    ].reset_index(drop=True)
262                    enr.enriche_specificiti()
263                    enr.enriche_KEGG()
264                    enr.enriche_GOTERM()
265                    enr.enriche_REACTOME()
266                    enr.enriche_IntAct()
267                    enr.enriche_STRING()
268                    enr.enriche_specificiti()
269
270                    data = enr.get_results()
271                    del enr
272
273                    ans = Analysis(data)
274                    ans.gene_interaction()
275                    ans.features_specificity()
276                    ans.REACTOME_overrepresentation()
277                    ans.KEGG_overrepresentation()
278                    ans.GO_overrepresentation()
279                    ans.features_specificity()
280
281                    data_dict[c] = ans.get_full_results()
282                else:
283                    print(
284                        f"Cell {c} was not enriched. No specific markers were found in this dataset."
285                    )
286                    data_dict[c] = None
287
288            self.cells_enrichment = data_dict
289
290        else:
291            raise ValueError(
292                "`self.cells_markers` not defined. Use `self.cells_markers` to provide markers."
293            )
294
295    def get_enrichment_data(
296        self,
297        data_type="GO-TERM",
298        p_value=0.05,
299        test="FISH",
300        adj="BH",
301        parent_inc=False,
302        top_n=50,
303    ):
304        """
305        Retrieves enrichment results for all cells in a unified DataFrame.
306
307        Parameters
308        ----------
309        data_type : str
310            Type of enrichment to retrieve ('GO-TERM', 'KEGG', 'REACTOME', 'specificity').
311
312        p_value : float, optional
313            Maximum p-value threshold (default 0.05).
314
315        test : str, optional
316            Name of the statistical test column to use (default 'FISH').
317
318        adj : str, optional
319            P-value adjustment method (default 'BH').
320
321        parent_inc : bool, optional
322            Whether to include parent terms in the results (default False).
323
324        top_n : int, optional
325            Maximum number of terms per cell type to include (default 50).
326
327        Returns
328        -------
329        pd.DataFrame
330            DataFrame containing filtered enrichment results with a 'cell' column indicating cell type.
331
332        Raises
333        ------
334        ValueError
335            If `data_type` is not one of the expected values.
336        """
337
338        if not any(
339            x in data_type for x in ("GO-TERM", "KEGG", "REACTOME", "specificity")
340        ):
341            raise ValueError(
342                "Invalid value for 'data_type'. Expected: 'GO-TERM', 'KEGG', 'REACTOME' or 'specificity'."
343            )
344
345        if data_type == "GO-TERM":
346            parent_col = "parent"
347
348        elif data_type == "KEGG":
349            parent_col = "2nd"
350
351        elif data_type == "REACTOME":
352            parent_col = "top_level"
353
354        elif data_type == "specificity":
355            parent_col = "None"
356
357        pdl = []
358        for i in self.cells_enrichment.keys():
359            if self.cells_enrichment[i] is None:
360                continue
361
362            print(i)
363            if data_type == "specificity":
364                tmp_dict = self.cells_enrichment[i]["statistics"][data_type]
365                tmp = []
366                for k in tmp_dict.keys():
367                    if k != "HPA_subcellular_location":
368                        tmp.append(pd.DataFrame(tmp_dict[k]))
369
370                tmp = pd.concat(tmp)
371
372            else:
373                tmp = pd.DataFrame(self.cells_enrichment[i]["statistics"][data_type])
374
375            cols = [x for x in tmp.columns if test in x and adj in x]
376            cols = sorted(cols, reverse=True)
377            if parent_inc is False:
378                cols = [x for x in cols if parent_col not in x.lower()]
379
380            mask = (tmp[cols] <= p_value).all(axis=1)
381            tmp = tmp.loc[mask]
382            tmp["cell"] = i
383            tmp = tmp.sort_values(by=["cell"] + cols, ascending=True)
384
385            pdl.append(tmp.head(top_n))
386
387        df = pd.concat(pdl)
388        df["source"] = data_type
389        df = df.reset_index(drop=True)
390
391        return df
392
393    def get_included_cells(self):
394        """
395        Returns the list of cell types included in the enrichment analysis.
396
397        Returns
398        -------
399        list
400            List of cell type names.
401
402        Example
403        -------
404        >>> self.get_included_cells()
405        ['CellType1', 'CellType2', ...]
406        """
407
408        cl = []
409        for i in self.cells_enrichment.keys():
410            print(i)
411            cl.append(i)
412
413        return cl
414
415    def get_gene_interactions(self, cell_name):
416        """
417        Retrieves gene or protein interaction data for a specific cell type.
418
419        Parameters
420        ----------
421        cell_name : str
422            Name of the cell type.
423
424        Returns
425        -------
426        pd.DataFrame
427            DataFrame containing interactions for the specified cell.
428
429        Example
430        -------
431        >>> self.get_gene_interactions('CellType1')
432        """
433
434        tmp = pd.DataFrame(
435            self.cells_enrichment[cell_name]["statistics"]["interactions"]
436        )
437
438        return tmp
439
440    def calculate_cell_connections(self):
441        """
442        Calculates cell-cell interaction connections based on gene/protein co-expression.
443
444        Notes
445        -----
446        Populates `cells_connection` with a DataFrame containing interactions between all pairs of cells.
447
448        Each row represents an interaction between two cells and the involved genes/proteins.
449
450        Raises
451        ------
452        ValueError
453            If `normalized_data` is not defined in the JDTI object.
454        """
455
456        if isinstance(self.jdti.normalized_data, pd.DataFrame):
457
458            cells = set(self.jdti.normalized_data.columns)
459
460            data_dict = {}
461
462            for c in tqdm(cells):
463
464                tmp = self.jdti.normalized_data.loc[:, c]
465                names = tmp.loc[
466                    tmp.select_dtypes(include="number").sum(axis=1) > 0
467                ].index.tolist()
468                names = list(set(names))
469
470                enr = copy.copy(self.enr_full_info)
471                enr.genome = enr.genome[
472                    enr.genome["found_names"].isin(names)
473                ].reset_index(drop=True)
474                enr.enriche_CellCon()
475                data = enr.get_results()
476                del enr
477
478                data_dict[c] = data["CellConnections"]
479
480            full_data = []
481            for c1 in tqdm(cells):
482                for c2 in cells:
483                    if c1 != c2:
484                        c1_d = pd.DataFrame(data_dict[c1]["interactor2"])
485                        c2_d = pd.DataFrame(data_dict[c2]["interactor1"])
486
487                        mutual_lr = c1_d["interaction"][
488                            c1_d["interaction"].isin(list(c2_d["interaction"]))
489                        ]
490
491                        to_ret = (
492                            c1_d[c1_d["interaction"].isin(list(mutual_lr))]
493                            .drop(
494                                [
495                                    "Species",
496                                    "protein_id_1",
497                                    "protein_id_2",
498                                    "found_names_2",
499                                ],
500                                axis=1,
501                            )
502                            .reset_index(drop=True)
503                        )
504
505                        to_ret = to_ret.rename(columns={"found_names_1": "interactor1"})
506                        c2_subset = c2_d[["interaction", "found_names_2"]].rename(
507                            columns={"found_names_2": "interactor2"}
508                        )
509
510                        to_ret = to_ret.merge(c2_subset, on="interaction", how="left")
511                        to_ret["cell1"] = c1
512                        to_ret["cell2"] = c2
513
514                        full_data.append(to_ret)
515
516            self.cells_connection = pd.concat(full_data)
517
518        else:
519            raise ValueError(
520                "`self.cells_markers` not defined. Use `self.cells_markers` to provide markers."
521            )
522
523    def get_cell_connections(
524        self,
525        connection_type: list = [
526            "Adhesion-Adhesion",
527            "Gap-Gap",
528            "Ligand-Ligand",
529            "Ligand-Receptor",
530            "Receptor-Receptor",
531            "Undefined",
532        ],
533    ):
534        """
535        Returns the calculated cell-cell interaction connections.
536
537        Parameters
538        ----------
539        connection_type : list of str, optional
540            List of interaction types used to filter the returned cell–cell connections
541            based on the molecular directionality of the interaction. Possible values:
542
543            - "Adhesion-Adhesion" : interaction between two adhesion molecules.
544            - "Gap-Gap" : connection through gap junction proteins.
545            - "Ligand-Ligand" : interaction between two ligand molecules.
546            - "Ligand-Receptor" : directional interaction where a ligand binds to a receptor.
547            - "Receptor-Receptor" : interaction between two receptor molecules.
548            - "Undefined" : interactions where the directionality could not be determined.
549
550            By default, all interaction types are included.
551
552        Returns
553        -------
554        pd.DataFrame
555            DataFrame containing cell-cell interactions.
556
557        Example
558        -------
559        >>> connections = self.get_cell_connections()
560        """
561
562        tmp = self.cells_connection
563
564        tmp["directionality"] = [
565            x if x is not None else "Undefined" for x in tmp["directionality"]
566        ]
567
568        tmp = tmp[tmp["directionality"].isin(connection_type)]
569
570        return tmp
571
572
573def compare_connections(
574    instances_dict: dict,
575    cells_compartment: dict | None = None,
576    connection_type: list = [
577        "Adhesion-Adhesion",
578        "Gap-Gap",
579        "Ligand-Ligand",
580        "Ligand-Receptor",
581        "Receptor-Receptor",
582        "Undefined",
583    ],
584):
585    """
586    Compare gene expression between two instances based on their cell connections.
587
588    This function compares normalized gene expression data from exactly two
589    instances stored in ``instances_dict``. Optionally, the comparison can be
590    restricted to specific cell compartments for each instance. Differential
591    expression analysis is performed using ``jdti.calc_DEG``.
592
593    Parameters
594    ----------
595    instances_dict : dict
596        Dictionary containing exactly two objects. Each object must have:
597
598        - ``jdti.normalized_data`` : pandas.DataFrame
599            Gene expression matrix with genes as rows and cells as columns.
600
601        - ``cells_connection`` : pandas.DataFrame
602            DataFrame containing at least the columns ``'interactor1'`` and
603            ``'interactor2'``.
604
605        The dictionary keys are used as group labels in the comparison.
606
607    cells_compartment : dict or None, optional
608        Dictionary mapping each key in ``instances_dict`` to a list of cell names
609        to be used for the comparison. If ``None``, all cells are used and genes
610        are filtered based on cell–cell connections.
611
612    connection_type : list of str, optional
613        List of interaction types used to filter cell–cell connections that are
614        considered in the gene expression comparison. Only connections with the
615        specified molecular interaction types will be used to define interacting
616        cells between the two instances.
617
618        Possible values:
619
620        - "Adhesion-Adhesion" : interaction between two adhesion molecules.
621        - "Gap-Gap" : connection mediated by gap junction proteins.
622        - "Ligand-Ligand" : interaction between two ligand molecules.
623        - "Ligand-Receptor" : directional interaction where a ligand binds to a receptor.
624        - "Receptor-Receptor" : interaction between two receptor molecules.
625        - "Undefined" : interactions where the directionality could not be determined.
626
627        By default, all interaction types are included in the comparison.
628
629    Returns
630    -------
631    pandas.DataFrame
632        Differential expression results returned by ``calc_DEG``, filtered to
633        include only rows where ``valid_group`` matches the first key in
634        ``instances_dict``.
635
636    Raises
637    ------
638    ValueError
639        If any cell specified in ``cells_compartment`` is not present in the
640        corresponding ``normalized_data`` columns.
641
642    Notes
643    -----
644    - Only genes common to both instances are considered.
645
646    - When ``cells_compartment`` is ``None``, genes are further restricted to
647      those appearing in the cell–cell interaction networks of either instance.
648
649    - The function assumes exactly two entries in ``instances_dict``.
650
651    - Differential expression is computed with ``min_exp=0`` and ``min_pct=0.1``.
652
653    See Also
654    --------
655    jdti.calc_DEG : Function used to compute differential expression.
656    """
657
658    import pandas as pd
659    from jdti import calc_DEG
660
661    if isinstance(cells_compartment, dict):
662
663        keys_list = list(instances_dict.keys())
664        tmp1 = instances_dict[keys_list[0]].jdti.normalized_data.copy()
665        cells = cells_compartment[keys_list[0]]
666        if any(cell not in tmp1.columns for cell in cells):
667            raise ValueError(
668                'Any of {keys_list[0]} cells in dictionary "cells_compartment" do not occur!'
669            )
670        tmp1 = tmp1.loc[:, cells]
671        tmp1.columns = [keys_list[0]] * len(tmp1.columns)
672
673        tmp2 = instances_dict[keys_list[1]].jdti.normalized_data.copy()
674        cells = cells_compartment[keys_list[1]]
675        if any(cell not in tmp2.columns for cell in cells):
676            raise ValueError(
677                'Any of {keys_list[1]} cells in dictionary "cells_compartment" do not occur!'
678            )
679        tmp2 = tmp2.loc[:, cells]
680        tmp2.columns = [keys_list[1]] * len(tmp2.columns)
681
682        common_idx = tmp1.index.intersection(tmp2.index)
683
684        tmp1 = tmp1.loc[common_idx]
685        tmp2 = tmp2.loc[common_idx]
686
687        concat_df = pd.concat([tmp1, tmp2], axis=1)
688
689    else:
690
691        keys_list = list(instances_dict.keys())
692        tmp1 = instances_dict[keys_list[0]].jdti.normalized_data.copy()
693        tmp1.columns = [keys_list[0]] * len(tmp1.columns)
694
695        tmp2 = instances_dict[keys_list[1]].jdti.normalized_data.copy()
696        tmp2.columns = [keys_list[1]] * len(tmp2.columns)
697
698        common_idx = tmp1.index.intersection(tmp2.index)
699
700        tmp1 = tmp1.loc[common_idx]
701        tmp2 = tmp2.loc[common_idx]
702
703        concat_df = pd.concat([tmp1, tmp2], axis=1)
704
705    tmp_df_1 = instances_dict[keys_list[0]].cells_connection
706    tmp_df_2 = instances_dict[keys_list[1]].cells_connection
707
708    tmp_df_1["directionality"] = [
709        x if x is not None else "Undefined" for x in tmp_df_1["directionality"]
710    ]
711    tmp_df_2["directionality"] = [
712        x if x is not None else "Undefined" for x in tmp_df_2["directionality"]
713    ]
714
715    tmp_df_1 = tmp_df_1[tmp_df_1["directionality"].isin(connection_type)]
716    tmp_df_2 = tmp_df_2[tmp_df_2["directionality"].isin(connection_type)]
717
718    tmp_con1 = list(set(list(tmp_df_1["interactor1"]) + list(tmp_df_1["interactor2"])))
719
720    tmp_con2 = list(set(list(tmp_df_2["interactor1"]) + list(tmp_df_2["interactor2"])))
721
722    genes = list(set(tmp_con1 + tmp_con2))
723
724    genes2 = [x for x in genes if x in common_idx]
725
726    concat_df = concat_df.loc[genes2, :]
727
728    results = calc_DEG(
729        data=concat_df,
730        metadata_list=None,
731        entities="All",
732        sets=None,
733        min_exp=0,
734        min_pct=0,
735        n_proc=10,
736    )
737
738    results = results[results["valid_group"] == keys_list[0]]
739
740    return results