analysis.dea

Utilities to perform differential expression analysis (DEA).

bulkdgd.analysis.dea.get_enrichment_scores(df_significant_genes: DataFrame, genes_sets: dict[str, list[str]], genes_all: list[str]) → DataFrame

Compute the enrichment scores for a set of genes of interest.

Parameters:

df_significant_genespandas.DataFrame

A data frame containing the genes that are significant at the given significance levels.

The index of the data frame is equal to the genes’ names.

genes_setsdict

A dictionary containing sets of genes of interest.

genes_alllist

A list containing all the genes in the analysis.

Returns:

df_e_scorespandas.DataFrame

A data frame containing the enrichment scores for each sample.

bulkdgd.analysis.dea.get_log2_fold_changes(obs_counts: Series, pred_means: Series, pseudocount: int = 1, scaling_factor: str = 'mean') → Series

Get the log2-fold change of the expression of a set of genes.

Parameters:

obs_countspandas.Series

The observed gene counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

pred_meanspandas.Series

The predicted means of the distributions modelling the genes’ counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

pseudocountint, 1

A pseudocount to add to both the predicted means and observed counts to avoid artifacts.

scaling_factorstr, {"mean", "median"}, "mean"

How the model computes the scaling factor of a sample - the number its predicted means are multiplied by to reach the scale of the sample’s own counts.

It must be the one the model was trained with, which is in the model’s own configuration file as "scaling_factor". The decoder is fitted against it, and undoing it here with the other one leaves every predicted mean wrong by the ratio of the two - about three, between the median and the mean - while every fold change still looks like a fold change.

Returns:

log2_fold_changespandas.Series

The log2-fold change associated with each gene in the given sample.

This is a series whose index correspond to the one of obs_counts and pred_means.

bulkdgd.analysis.dea.get_p_values(obs_counts: Series, pred_means: Series, r_values: Series | None = None, resolution: int | None = None, return_pmf_values: bool = False, pseudocount: int = 1, device: str | device = 'cpu', p_values_method: str = 'auto', max_elements: int | None = None, scaling_factor: str = 'mean') → tuple[Series, DataFrame, DataFrame]

Given the observed gene counts in a single sample, and the predicted mean gene counts in a single sample, calculate the p-value associated with the predicted mean of each distribution modeling a gene’s counts by comparing it to the actual gene count.

Parameters:

obs_countspandas.Series

The observed gene counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

pred_meanspandas.Series

The predicted means of the distributions modelling the genes’ counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

If the genes’ counts were modelled using negative binomial distributions, the predicted means are scaled by the corresponding distributions’ r-values.

r_valuespandas.Series, optional

The predicted r-values of the negative binomial distributions modelling the genes’ counts in a single sample, if the genes’ counts were modelled using negative binomial distributions.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

If r_values is not provided, it is assumed that the genes’ counts were modelled using Poisson distributions.

resolutionint, optional

How accurate the calculation of the p-values should be.

The resolution corresponds to the coarseness of the sum over the probability mass function of each distribution to compute the corresponding p-value.

The higher the resolution, the more accurate (and more computationally expensive) the calculation of the p-values will be.

If not passed, the calculation will be exact.

return_pmf_valuesbool, False

Return the points at which the log-probability mass function was evaluated and the corresponding values of the log- probability mass function, together with the p-values.

Set it to True only if you have a low resolution (for instance, 1e3 or lower) or a lot of RAM available since the arrays containing the points at which the log- probability mass function was evaluated and the corresponding values of the function will contain resolution floating-point numbers for each gene.

Setting it to True forces the p-values to be computed gene by gene on the CPU, ignoring device.

pseudocountint, 1

A pseudocount to add to both the predicted means and observed counts to avoid artifacts.

devicestr or torch.device, "cpu"

The device on which to compute the p-values.

The genes are independent of each other, so the calculation of the p-values parallelizes exactly, and running it on a GPU (for instance, by passing "cuda") speeds it up considerably.

It is ignored if return_pmf_values is True.

p_values_methodstr, "auto"

How to compute the p-values. The methods give the same p-values, but differ in how they use the machine.

"batched" computes the p-values for all the genes at once. This is what allows them to be computed on a GPU. On a CPU, torch spreads the computation over the cores by itself, so it should be used in a single process.

"per-gene" computes the p-values one gene at a time. It is meant to be parallelized over the samples - by perform_dea()’s caller, or by bulkdgd_dea’s -n option - which is how to use a machine with many cores.

The two kinds of parallelism do not compose: several processes each running a "batched" computation would fight over the CPU’s cores, and the result would be slower than either kind of parallelism on its own.

"auto" uses "batched" on a GPU, and "per-gene" on a CPU for the exact calculation (where the batch would have to be padded), and "batched" on a CPU otherwise.

max_elementsint, optional

The maximum number of points at which the log-probability mass function is evaluated in one batch, which caps the memory used by the "batched" method.

If not passed, it defaults to 2**26 on a GPU and 2**22 on a CPU, where the batch sits in RAM that several processes may be sharing.

scaling_factorstr, {"mean", "median"}, "mean"

How the model computes the scaling factor of a sample - the number its predicted means are multiplied by to reach the scale of the sample’s own counts.

It must be the one the model was trained with, which is in the model’s own configuration file as "scaling_factor". The decoder is fitted against it, and undoing it here with the other one leaves every predicted mean wrong by the ratio of the two - about three, between the median and the mean - while every p-value still looks like a p-value.

Returns:

p_valuespandas.Series

A series containing one p-value per gene.

kspandas.DataFrame

A data frame containing the count values at which the log- probability mass function was evaluated to compute the p-values.

The data frame has as many rows as the number of genes and as many columns as the number of count values.

This is an empty data frame if return_pmf_values is False.

pmfsnumpy.ndarray

A data frame containing the value of the log-probability mass function for each count value at which it was evaluated.

The data frame has as many rows as the number of genes and as many columns as the number of count values.

This is an empty data frame if return_pmf_values is False.

bulkdgd.analysis.dea.get_q_values(p_values: Series, alpha: float = 0.05, method: str = 'fdr_bh') → tuple[Series, Series]

Get the q-values associated with a set of p-values.

The q-values are the p-values adjusted for the false discovery rate.

Parameters:

p_valuespandas.Series

The p-values.

alphafloat, 0.05

The family-wise error rate for the calculation of the q-values.

methodstr, "fdr_bh"

The method used to adjust the p-values. The available methods are listed in the documentation for statsmodels.stats.multitest.multipletests.

Returns:

q_valuespandas.Series

A series containing the q-values.

The index of the series is equal to the index of the input series of p-values.

rejectedpandas.Series

A series containing booleans indicating whether a p-value in the input data frame was rejected (True) or not (False).

The index of the series is equal to the index of the input series of p-values.

bulkdgd.analysis.dea.get_significant_genes(df_stats: DataFrame, p_val: float = 0.05, q_val: float = 0.05, log2_fold_change: float = 2) → DataFrame

Get the genes that are significant at a given significance level.

Parameters:

df_statspandas.DataFrame

A data frame whose rows represent the genes on which the DEA was performed, and whose columns contain the statistics computed (p-values, q_values, log2-fold changes).

p_valfloat, 0.05

The p-value threshold to consider a gene as significant.

q_valfloat, 0.05

The q-value threshold to consider a gene as significant.

log2_fold_changefloat, 2

The log2-fold change threshold to consider a gene as significant. This value and its negative are used as the thresholds for the log2-fold change.

Returns:

df_significant_genespandas.DataFrame

A data frame containing the genes that are significant at the given significance levels.

bulkdgd.analysis.dea.get_statistics(obs_counts: Series, pred_means: Series, r_values: Series | None = None, sample_name: str | None = None, statistics: list[str] = ['p_values', 'q_values', 'log2_fold_changes'], resolution: int | None = None, alpha: float = 0.05, method: str = 'fdr_bh', pseudocount: int = 1, device: str | device = 'cpu', p_values_method: str = 'auto', max_elements: int | None = None, scaling_factor: str = 'mean') → tuple[DataFrame, str | None]

Compute p-values, q-values, and/or log2-fold changes.

Parameters:

obs_countspandas.Series

The observed gene counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

pred_meanspandas.Series

The predicted means of the distributions modelling the genes’ counts in a single sample.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

If the genes’ counts were modelled using negative binomial distributions, the predicted means are scaled by the corresponding distributions’ r-values.

r_valuespandas.Series, optional

The predicted r-values of the negative binomial distributions modelling the genes’ counts in a single sample, if the genes’ counts were modelled using negative binomial distributions.

This is a series whose index contains either the genes’ Ensembl IDs or names of fields containing additional information about the sample.

If r_values is not provided, it is assumed that the genes’ counts were modelled using Poisson distributions.

sample_namestr, optional

The name of the sample under consideration.

It is returned together with the results of the analysis to facilitate the identification of the sample when running the analysis in parallel for multiple samples (i.e., launching the function in parallel on multiple samples).

statisticslist, {["p_values", "q_values", "log2_fold_changes"]}

The statistics to be computed. By default, all of them will be computed.

resolutionint, optional

How accurate the calculation of the p-values should be.

The resolution corresponds to the coarseness of the sum over the probability mass function of each distribution to compute the corresponding p-value.

The higher the resolution, the more accurate (and more computationally expensive) the calculation of the p-values will be.

If not passed, the calculation will be exact.

alphafloat, 0.05

The family-wise error rate for the calculation of the q-values (adjusted p-values).

methodstr, "fdr_bh"

The method used to calculate the q-values (in other words, to adjust the p-values). The available methods are listed in the documentation for statsmodels.stats.multitest.multipletests.

pseudocountint, 1

A pseudocount to add to both the predicted means and observed counts to avoid artifacts.

devicestr or torch.device, "cpu"

The device on which to compute the p-values.

The genes are independent of each other, so the calculation of the p-values parallelizes exactly, and running it on a GPU (for instance, by passing "cuda") speeds it up considerably.

p_values_methodstr, "auto"

How to compute the p-values. The methods give the same p-values, but differ in how they use the machine.

"batched" computes the p-values for all the genes at once. This is what allows them to be computed on a GPU. On a CPU, torch spreads the computation over the cores by itself, so it should be used in a single process.

"per-gene" computes the p-values one gene at a time. It is meant to be parallelized over the samples - by perform_dea()’s caller, or by bulkdgd_dea’s -n option - which is how to use a machine with many cores.

The two kinds of parallelism do not compose: several processes each running a "batched" computation would fight over the CPU’s cores, and the result would be slower than either kind of parallelism on its own.

"auto" uses "batched" on a GPU, and "per-gene" on a CPU for the exact calculation (where the batch would have to be padded), and "batched" on a CPU otherwise.

max_elementsint, optional

The maximum number of points at which the log-probability mass function is evaluated in one batch, which caps the memory used by the "batched" method.

If not passed, it defaults to 2**26 on a GPU and 2**22 on a CPU, where the batch sits in RAM that several processes may be sharing.

scaling_factorstr, {"mean", "median"}, "mean"

How the model computes the scaling factor of a sample - the number its predicted means are multiplied by to reach the scale of the sample’s own counts.

It must be the one the model was trained with, which is in the model’s own configuration file as "scaling_factor". The decoder is fitted against it, and undoing it here with the other one leaves every predicted mean wrong by the ratio of the two - about three, between the median and the mean - while every p-value still looks like a p-value.

Returns:

df_statspandas.DataFrame

A data frame whose rows represent the genes on which the DEA was performed, and whose columns contain the statistics computed (p-values, q_values, log2-fold changes). If not all statistics were computed, the columns corresponding to the missing ones will be empty.

sample_namestr or None

The name of the sample under consideration.

bulkdgd.analysis.dea.perform_dea(obs_counts: DataFrame, pred_means: DataFrame, r_values: DataFrame | None = None, resolution: int | None = None, alpha: float = 0.05, method: str = 'fdr_bh', p_val: float = 0.05, q_val: float = 0.05, log2_fold_change: float = 2, genes_sets: dict[str, list[str]] | None = None, pseudocount: int = 1, device: str | device = 'cpu', p_values_method: str = 'auto', max_elements: int | None = None, scaling_factor: str = 'mean') → tuple[dict[str, DataFrame], Series, DataFrame]

Perform differential expression analysis (DEA) on multiple samples.

Parameters:

obs_countspandas.DataFrame

The observed gene counts in multiple sample.

This is a data frame whose index contains the samples’s names, and the columns contain either the genes’ Ensembl IDs or names of fields containing additional information about the samples.

pred_meanspandas.DataFrame

The predicted means of the distributions modelling the genes’ counts in each sample.

This is a data frame whose index contains the samples’ names, and the columns contain either the genes’ Ensembl IDs or names of fields containing additional information about the samples.

r_valuespandas.DataFrame, optional

The predicted r-values of the negative binomial distributions modelling the genes’ counts in each sample, if the genes’ counts were modelled using negative binomial distributions.

This is a data frame whose index contains the samples’ names, and the columns contain either the genes’ Ensembl IDs or names of fields containing additional information about the samples.

If r_values is not provided, it is assumed that the genes’ counts were modelled using Poisson distributions.

resolutionint, optional

How accurate the calculation of the p-values should be.

The resolution corresponds to the coarseness of the sum over the probability mass function of each distribution to compute the corresponding p-value.

The higher the resolution, the more accurate (and more computationally expensive) the calculation of the p-values will be.

If not passed, the calculation will be exact.

alphafloat, 0.05

The family-wise error rate for the calculation of the q-values (adjusted p-values).

methodstr, "fdr_bh"

The method used to calculate the q-values (in other words, to adjust the p-values). The available methods are listed in the documentation for statsmodels.stats.multitest.multipletests.

p_valfloat, 0.05

The p-value threshold to consider a gene as significant.

q_valfloat, 0.05

The q-value threshold to consider a gene as significant.

log2_fold_changefloat, 2

The log2-fold change threshold to consider a gene as significant. This value and its negative are used as the thresholds for the log2-fold change.

genes_setsdict, optional

A dictionary containing sets of genes of interest.

The keys are the names of the gene sets, and the values are lists of genes.

pseudocountint, 1

A pseudocount to add to both the predicted means and observed counts to avoid artifacts.

devicestr or torch.device, "cpu"

The device on which to compute the p-values.

The genes are independent of each other, so the calculation of the p-values parallelizes exactly, and running it on a GPU (for instance, by passing "cuda") speeds it up considerably.

p_values_methodstr, "auto"

How to compute the p-values. The methods give the same p-values, but differ in how they use the machine.

"batched" computes the p-values for all the genes at once. This is what allows them to be computed on a GPU. On a CPU, torch spreads the computation over the cores by itself, so it should be used in a single process.

"per-gene" computes the p-values one gene at a time. It is meant to be parallelized over the samples, which is how to use a machine with many cores.

The two kinds of parallelism do not compose: several processes each running a "batched" computation would fight over the CPU’s cores, and the result would be slower than either kind of parallelism on its own.

"auto" uses "batched" on a GPU, and "per-gene" on a CPU for the exact calculation (where the batch would have to be padded), and "batched" on a CPU otherwise.

max_elementsint, optional

The maximum number of points at which the log-probability mass function is evaluated in one batch, which caps the memory used by the "batched" method.

If not passed, it defaults to 2**26 on a GPU and 2**22 on a CPU, where the batch sits in RAM that several processes may be sharing.

scaling_factorstr, {"mean", "median"}, "mean"

How the model computes the scaling factor of a sample - the number its predicted means are multiplied by to reach the scale of the sample’s own counts.

It must be the one the model was trained with, which is in the model’s own configuration file as "scaling_factor". The decoder is fitted against it, and undoing it here with the other one leaves every predicted mean wrong by the ratio of the two - about three, between the median and the mean - while every p-value still looks like a p-value.

Returns:

dfs_statsdict

A dictionary containing the data frames with the statistics for each sample.

series_significant_genespandas.Series

A series containing the significant genes per sample.

df_e_scorespandas.DataFrame

A data frame containing the enrichment scores for each sample.

If no gene sets were passed, the data frame will be empty.