kAOV: Kernel Analysis Of Variance#
[1]:
from kaov import AOV
from kaov.datasets import load_reversion
import re
[2]:
import kaov.datasets.data
[5]:
kaov.datasets
[5]:
<module 'kaov.datasets' from '/Users/polina/ktest-py3-9-venv/lib/python3.9/site-packages/kaov/datasets/__init__.py'>
[7]:
import inspect
src = inspect.getsource(kaov.datasets)
[8]:
src
[8]:
'#!/usr/bin/env python3\n# -*- coding: utf-8 -*-\n"""\nCreated on Wed Jul 29 17:40:32 2026\n\n@author: Polina Arsenteva\n"""\n\nfrom kaov.datasets._base import (\n load_reversion,\n load_rabbits_anndata,\n)\n\n'
Import data#
Various functionalities of the library will be demonstrated on a single-cell RT-qPCR dataset (Zreika et al., 2022) with a 4-level experimental condition of interest (i.e. Medium), and an 8-level batch effect. In this tutorial we will focus on different ways to specify the model and tests.
The dataset contains the expression measurements of 83 genes on 685 cells. In this context single-cell transcriptomics has been performed to investigate the cell differentiation process of chicken primary erythroid progenitor cells. To investigate the impact of the medium on cell differentiation, undifferentiated cells were initially put in a self-renewal medium (0H level of the medium effect), then put in a differentiation-inducing medium for 24h (24H). The population was then split
into a first population maintained in the same medium for an additional 24h to achieve differentiation (48HDIFF), the second population was put back in the self-renewal medium to investigate potential reversion (48HREV).
[3]:
# Importing data:
data = load_reversion()
# Specifying the formula for the dependent variables (genes):
form_exog = ' + '.join(data.columns[:-2])
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[3], line 2
1 # Importing data:
----> 2 data = load_reversion()
3 # Specifying the formula for the dependent variables (genes):
4 form_exog = ' + '.join(data.columns[:-2])
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/kaov/datasets/_base.py:24, in load_reversion()
13 """
14 Load the reversion dataset from a .csv file.
15
(...)
21
22 """
23 data_path = resources.files("kaov.datasets.data") / "reversion.csv"
---> 24 data = pd.read_csv(data_path, index_col=0)
25 return data
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/pandas/io/parsers/readers.py:1026, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)
1013 kwds_defaults = _refine_defaults_read(
1014 dialect,
1015 delimiter,
(...)
1022 dtype_backend=dtype_backend,
1023 )
1024 kwds.update(kwds_defaults)
-> 1026 return _read(filepath_or_buffer, kwds)
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/pandas/io/parsers/readers.py:620, in _read(filepath_or_buffer, kwds)
617 _validate_names(kwds.get("names", None))
619 # Create the parser.
--> 620 parser = TextFileReader(filepath_or_buffer, **kwds)
622 if chunksize or iterator:
623 return parser
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/pandas/io/parsers/readers.py:1620, in TextFileReader.__init__(self, f, engine, **kwds)
1617 self.options["has_index_names"] = kwds["has_index_names"]
1619 self.handles: IOHandles | None = None
-> 1620 self._engine = self._make_engine(f, self.engine)
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/pandas/io/parsers/readers.py:1880, in TextFileReader._make_engine(self, f, engine)
1878 if "b" not in mode:
1879 mode += "b"
-> 1880 self.handles = get_handle(
1881 f,
1882 mode,
1883 encoding=self.options.get("encoding", None),
1884 compression=self.options.get("compression", None),
1885 memory_map=self.options.get("memory_map", False),
1886 is_text=is_text,
1887 errors=self.options.get("encoding_errors", "strict"),
1888 storage_options=self.options.get("storage_options", None),
1889 )
1890 assert self.handles is not None
1891 f = self.handles.handle
File ~/ktest-py3-9-venv/lib/python3.9/site-packages/pandas/io/common.py:873, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
868 elif isinstance(handle, str):
869 # Check whether the filename is to be opened in binary mode.
870 # Binary mode does not support 'encoding' and 'newline'.
871 if ioargs.encoding and "b" not in ioargs.mode:
872 # Encoding
--> 873 handle = open(
874 handle,
875 ioargs.mode,
876 encoding=ioargs.encoding,
877 errors=errors,
878 newline="",
879 )
880 else:
881 # Binary mode
882 handle = open(handle, ioargs.mode)
FileNotFoundError: [Errno 2] No such file or directory: '/Users/polina/ktest-py3-9-venv/lib/python3.9/site-packages/kaov/datasets/data/reversion.csv'
[4]:
from importlib import resources
[5]:
data_path = resources.files("kaov.datasets.data") / "reversion.csv"
[6]:
data_path
[6]:
PosixPath('/Users/polina/ktest-py3-9-venv/lib/python3.9/site-packages/kaov/datasets/data/reversion.csv')
[ ]:
data = pd.read_csv(data_path, index_col=0)
[4]:
data[['Batch', 'Medium']].value_counts(sort=False)
[4]:
Batch Medium
REV1 0H 23
24H 21
48HDIFF 24
48HREV 21
REV2 0H 23
24H 24
48HDIFF 22
48HREV 21
REV3 0H 21
24H 24
48HDIFF 22
48HREV 22
REV4 0H 18
24H 21
48HDIFF 19
48HREV 19
REV5 0H 23
24H 21
48HDIFF 21
48HREV 23
REV6 0H 20
24H 22
48HDIFF 22
48HREV 22
REV7 0H 21
24H 20
48HDIFF 23
48HREV 20
REV8 0H 24
24H 20
48HDIFF 15
48HREV 23
Name: count, dtype: int64
An additive model#
Module AOV infers a linear model in the feature space associated with the kernel transormation of the data:
[5]:
# kAOV has a similar interface as statsmodels and uses a similar logic:
kfit = AOV.from_formula(form_exog + ' ~ C(Medium, OneHot) + C(Batch, OneHot)',
data=data)
Tests for factor effects#
Use method test in its most basic configuration in order to test for the effects of each factor in the model:
[6]:
res = kfit.test()
print(res)
Kernel Analysis of Variance (trunc. 1):
====================================
------------------------------------
Factor test | factor stat pval
------------------------------------
| Medium 67.0362 0.0000
| Batch 12.0566 0.0000
====================================
We can conclude that all the effects are significant, with the medium effect being the strongest. This can be visualized using density plots and mean embedding plots (densities and means of cellular embeddings projected onto the discriminant axes associated with the test in question):
[7]:
res.plot_density(comp=1, tests=['Medium',], colormap='jet')
[7]:
(<Figure size 800x600 with 1 Axes>,
<Axes: title={'center': 'Medium'}, xlabel='Discriminant axis', ylabel='Density'>)
[8]:
res.plot_mean_embedding_projections(colormap='jet')
[8]:
(<Figure size 1600x600 with 2 Axes>,
array([<Axes: title={'center': 'Medium'}, xlabel='Discriminant axis 1', ylabel='Discriminant axis 2'>,
<Axes: title={'center': 'Batch'}, xlabel='Discriminant axis 1', ylabel='Discriminant axis 2'>],
dtype=object))
Study of the medium effect#
We can perform pairwise comparisons of different considered media by specifying by_level=True and hypotheses='pairwise':
[9]:
res_pw_bl = kfit.test(by_level=True, hypotheses='pairwise', verbose=1)
-Computing the Gram matrix...
-Testing hypotheses:
100%|███████████████████████████████████████████| 34/34 [00:12<00:00, 2.68it/s]
[10]:
print(res_pw_bl)
Kernel Analysis of Variance (trunc. 1):
===============================================
-----------------------------------------------
Medium | factor_1_1 factor_1_2 stat pval
-----------------------------------------------
| 0H 24H 58.2115 0.0000
| 0H 48HDIFF 135.9141 0.0000
| 0H 48HREV 0.0000 0.9997
| 24H 48HDIFF 16.7631 0.0000
| 24H 48HREV 57.9055 0.0000
| 48HDIFF 48HREV 135.1675 0.0000
-----------------------------------------------
-----------------------------------------------
Batch | factor_1_1 factor_1_2 stat pval
-----------------------------------------------
| REV1 REV2 1.4273 0.2326
| REV1 REV3 21.7174 0.0000
| REV1 REV4 28.0970 0.0000
| REV1 REV5 7.9975 0.0048
| REV1 REV6 16.1508 0.0001
| REV1 REV7 0.0280 0.8671
| REV1 REV8 0.8696 0.3514
| REV2 REV3 34.4460 0.0000
| REV2 REV4 41.8119 0.0000
| REV2 REV5 16.2178 0.0001
| REV2 REV6 27.1902 0.0000
| REV2 REV7 1.0189 0.3131
| REV2 REV8 4.4279 0.0357
| REV3 REV4 0.6600 0.4168
| REV3 REV5 3.3101 0.0693
| REV3 REV6 0.3621 0.5475
| REV3 REV7 22.6567 0.0000
| REV3 REV8 13.1568 0.0003
| REV4 REV5 6.5677 0.0106
| REV4 REV6 1.9205 0.1663
| REV4 REV7 29.0581 0.0000
| REV4 REV8 18.4459 0.0000
| REV5 REV6 1.4493 0.2291
| REV5 REV7 8.7252 0.0032
| REV5 REV8 3.3787 0.0665
| REV6 REV7 17.0359 0.0000
| REV6 REV8 9.0495 0.0027
| REV7 REV8 1.1732 0.2791
===============================================
By default, pairwise by-level tests comparing levels of all factors are performed. Let us extract the results conserning the medium factor using summary:
[11]:
res_pw_bl.summary(trunc=1, factor='Medium')
[11]:
| factor_1_1 | factor_1_2 | stat | pval | |
|---|---|---|---|---|
| 1 | 0H | 24H | 58.21148 | 0.0 |
| 2 | 0H | 48HDIFF | 135.914078 | 0.0 |
| 3 | 0H | 48HREV | 0.0 | 0.999731 |
| 4 | 24H | 48HDIFF | 16.763146 | 0.000047 |
| 5 | 24H | 48HREV | 57.905546 | 0.0 |
| 6 | 48HDIFF | 48HREV | 135.167546 | 0.0 |
As expected, the most distinct medium pairs are differentiated with reverted and differentiated with non-differentiated. For the density plots, we can specify the refecrence tests of interest using the tests parameter:
[12]:
res_pw_bl.plot_density(comp=1, tests=[n for n, _ in res_pw_bl.hypotheses[:6]])
[12]:
(<Figure size 4800x600 with 6 Axes>,
array([<Axes: title={'center': 'Medium[0H] = Medium[24H]'}, xlabel='Discriminant axis', ylabel='Density'>,
<Axes: title={'center': 'Medium[0H] = Medium[48HDIFF]'}, xlabel='Discriminant axis', ylabel='Density'>,
<Axes: title={'center': 'Medium[0H] = Medium[48HREV]'}, xlabel='Discriminant axis', ylabel='Density'>,
<Axes: title={'center': 'Medium[24H] = Medium[48HDIFF]'}, xlabel='Discriminant axis', ylabel='Density'>,
<Axes: title={'center': 'Medium[24H] = Medium[48HREV]'}, xlabel='Discriminant axis', ylabel='Density'>,
<Axes: title={'center': 'Medium[48HDIFF] = Medium[48HREV]'}, xlabel='Discriminant axis', ylabel='Density'>],
dtype=object))
Study of the batch effect#
The test for the batch effect revealed that there are significant differences between batches, but are there batches that particularly stand out that could potentially introduce bias in the analysis? A first assessment can be performed by looking as the diagnostic plots, representing varios quantities of interest in the model projected onto the eigenfuctions of the residual covariance operator. By default it is the projections of the model residuals plotted against those of the predictions made by the model:
[13]:
kfit.plot_diagnostics(trunc=1, colormap='jet', figsize=(20, 5))
[13]:
(<Figure size 2000x500 with 2 Axes>,
array([<Axes: title={'center': 'Batch'}, xlabel='Predictions', ylabel='Residuals'>,
<Axes: title={'center': 'Medium'}, xlabel='Predictions', ylabel='Residuals'>],
dtype=object))
Next, we can perform a series of tests comparing each batch to the mean of all batches using the by-level option with hypotheses='one-vs-all':
[14]:
res_oa_bl = kfit.test(hypotheses='one-vs-all', by_level=True, verbose=1)
-Computing the Gram matrix...
-Testing hypotheses:
100%|███████████████████████████████████████████| 12/12 [00:04<00:00, 2.80it/s]
[15]:
res_oa_bl.summary(trunc=1, factor='Batch')
[15]:
| factor_1_1 | factor_1_2 | stat | pval | |
|---|---|---|---|---|
| 1 | REV1 | Grand Mean | 9.83492 | 0.001786 |
| 2 | REV2 | Grand Mean | 24.563171 | 0.000001 |
| 3 | REV3 | Grand Mean | 15.122884 | 0.000111 |
| 4 | REV4 | Grand Mean | 23.503089 | 0.000002 |
| 5 | REV5 | Grand Mean | 1.283709 | 0.257607 |
| 6 | REV6 | Grand Mean | 8.589473 | 0.003494 |
| 7 | REV7 | Grand Mean | 10.950288 | 0.000985 |
| 8 | REV8 | Grand Mean | 2.683988 | 0.101821 |
Batches 2 and 4 appear to be the most distinct. To detect outliers on the level of individual cells, we can examine their influences on these tests ising Cook’s distance. Method plot_influence plots these distances in function of discriminant axis projections:
[16]:
res_oa_bl.plot_influence(trunc=1, comp=1, tests=[n for n, _ in res_oa_bl.hypotheses[4:12]], colormap='jet')
[16]:
(<Figure size 6400x600 with 8 Axes>,
array([<Axes: title={'center': 'Batch[REV1] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV2] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV3] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV4] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV5] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV6] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV7] = Batch Grand Mean'}>,
<Axes: title={'center': 'Batch[REV8] = Batch Grand Mean'}>],
dtype=object))
The Cook’s distances and their p-values can be accessed using the method get_cook:
[13]:
cook_batch_bl = res_oa_bl.get_cook(factor='Batch', n_trunc=1)
[14]:
cook_batch_bl.groupby('factor_1')['cook_1'].mean()
[14]:
factor_1
REV1 0.003435
REV2 0.003386
REV3 0.002013
REV4 0.002247
REV5 0.001616
REV6 0.002871
REV7 0.00122
REV8 0.002147
Name: cook_1, dtype: object
Batches 1 and 2 have the biggest amount of significantly influential cells:
[15]:
cook_batch_bl.loc[(cook_batch_bl['cook_pval_1'] < 0.05), 'factor_1'].value_counts()
[15]:
factor_1
REV2 21
REV1 20
REV6 14
REV3 10
REV4 9
REV8 9
REV7 5
REV5 4
Name: count, dtype: int64
Model with an interaction factor#
Using the formula interface, we can specify a more complex model by adding an interaction effect:
[17]:
kfit_inter = AOV.from_formula(form_exog + ' ~ C(Medium, OneHot) * C(Batch, OneHot)',
data=data)
We can specifically test for the interaction effect by specifying the corresponding factor using hypotheses_subset:
[19]:
res_inter = kfit_inter.test(hypotheses_subset=['Medium:Batch',])
print(res_inter)
Kernel Analysis of Variance (trunc. 1):
=========================================
-----------------------------------------
Factor test | factor stat pval
-----------------------------------------
| Medium:Batch 9.8741 0.0000
=========================================
Next, we are interested in comparing different media in every batch. The corresponding set of hypotheses can be extracted from all by-level pairwise hypotheses using the method set_hypotheses:
[20]:
hyps_inter = kfit_inter.set_hypotheses(by_level=True, hypotheses='pairwise')
[27]:
hyps_to_test = []
for hyp_name, _ in hyps_inter:
lvls = re.findall(r"\[(.*?)\]", hyp_name)
if len(lvls) == 4 and lvls[1] == lvls[3]:
hyps_to_test.append(hyp_name)
[29]:
res_inter_pw_bl = kfit_inter.test(by_level=True, hypotheses='pairwise', hypotheses_subset=hyps_to_test, verbose=1)
-Computing the Gram matrix...
-Testing hypotheses:
100%|███████████████████████████████████████████| 48/48 [00:17<00:00, 2.69it/s]
[30]:
print(res_inter_pw_bl)
Kernel Analysis of Variance (trunc. 1):
=========================================================================
-------------------------------------------------------------------------
Medium:Batch | factor_1_1 factor_2_1 factor_1_2 factor_2_2 stat pval
-------------------------------------------------------------------------
| 0H REV1 24H REV1 0.0217 0.8830
| 0H REV1 48HDIFF REV1 0.3535 0.5523
| 0H REV1 48HREV REV1 1.0758 0.3000
| 24H REV1 48HDIFF REV1 0.1883 0.6645
| 24H REV1 48HREV REV1 1.3490 0.2459
| 48HDIFF REV1 48HREV REV1 2.6484 0.1041
| 0H REV2 24H REV2 0.0717 0.7890
| 0H REV2 48HDIFF REV2 0.0298 0.8630
| 0H REV2 48HREV REV2 1.1903 0.2756
| 24H REV2 48HDIFF REV2 0.1927 0.6608
| 24H REV2 48HREV REV2 1.8577 0.1733
| 48HDIFF REV2 48HREV REV2 0.8285 0.3630
| 0H REV3 24H REV3 0.0205 0.8863
| 0H REV3 48HDIFF REV3 0.0299 0.8627
| 0H REV3 48HREV REV3 0.0949 0.7581
| 24H REV3 48HDIFF REV3 0.0011 0.9731
| 24H REV3 48HREV REV3 0.2133 0.6444
| 48HDIFF REV3 48HREV REV3 0.2359 0.6274
| 0H REV4 24H REV4 1.1189 0.2905
| 0H REV4 48HDIFF REV4 0.2580 0.6117
| 0H REV4 48HREV REV4 0.9315 0.3348
| 24H REV4 48HDIFF REV4 2.5469 0.1110
| 24H REV4 48HREV REV4 0.0052 0.9423
| 48HDIFF REV4 48HREV REV4 2.2207 0.1366
| 0H REV5 24H REV5 1.6701 0.1967
| 0H REV5 48HDIFF REV5 1.2892 0.2566
| 0H REV5 48HREV REV5 0.8124 0.3677
| 24H REV5 48HDIFF REV5 0.0232 0.8790
| 24H REV5 48HREV REV5 4.7314 0.0300
| 48HDIFF REV5 48HREV REV5 4.0689 0.0441
| 0H REV6 24H REV6 0.0045 0.9467
| 0H REV6 48HDIFF REV6 0.0759 0.7830
| 0H REV6 48HREV REV6 0.7437 0.3888
| 24H REV6 48HDIFF REV6 0.1224 0.7266
| 24H REV6 48HREV REV6 0.9017 0.3427
| 48HDIFF REV6 48HREV REV6 0.3584 0.5496
| 0H REV7 24H REV7 0.7417 0.3894
| 0H REV7 48HDIFF REV7 0.9388 0.3329
| 0H REV7 48HREV REV7 1.4654 0.2265
| 24H REV7 48HDIFF REV7 0.0066 0.9355
| 24H REV7 48HREV REV7 4.2006 0.0408
| 48HDIFF REV7 48HREV REV7 4.7977 0.0288
| 0H REV8 24H REV8 2.7291 0.0990
| 0H REV8 48HDIFF REV8 2.6622 0.1032
| 0H REV8 48HREV REV8 5.0614 0.0248
| 24H REV8 48HDIFF REV8 0.0094 0.9228
| 24H REV8 48HREV REV8 0.2723 0.6020
| 48HDIFF REV8 48HREV REV8 0.1482 0.7004
=========================================================================
Bibliography#
Zreika, C. Fourneaux, E. Vallin, L. Modolo, R. Seraphin, A. Moussy, E. Ventre, M. Bouvier, A. Ozier-Lafontaine, A. Bonnaffoux, F. Picard, O. Gandrillon, and S. Gonin-Giraud. Evidence for close molecular proximity between reverting and undifferentiated cells. BMC Biology, 20(1):155, July
[ ]: