adesomics.dge.statistics

View Source

 1import numpy as np
 2import pandas as pd
 3from pydeseq2.dds import DeseqDataSet
 4from pydeseq2.ds import DeseqStats
 5
 6from adesomics.io.metadata import Metadata
 7from adesomics.plotting.style import group_colors, annotate_study, GREY, INK, PINK
 8from adesomics.population.population_model import _save
 9
10class DifferentialExpression:
11
12    def __init__(self, count_matrix, metadata: Metadata, design=None):
13        counts = getattr(count_matrix, "counts", count_matrix)
14        self.metadata = metadata.align(counts)
15        self.counts = counts[self.metadata.samples]    
16        self.design = design or f"~{self.metadata.group_col}"
17        self.dds = None
18        self._cache = {}
19
20    @classmethod
21    def from_sra_run_table(cls, count_matrix, path, keep=None, **kwargs):
22        md = Metadata.from_sra_run_table(path, **kwargs)
23        if keep:
24            md = md.subset(keep=keep)
25        return cls(count_matrix, md)
26
27    def fit(self):
28        n = self.metadata.replicates()
29        if (n < 2).any():
30            raise ValueError(f"Need >=2 replicates per group: {n[n < 2].to_dict()}")
31
32        self.dds = DeseqDataSet(counts=self.counts.T,
33                                metadata=self.metadata.table,
34                                design=self.design)
35        self.dds.deseq2()
36        return self
37
38    def results(self, treated, reference, alpha=0.05):
39        key = (treated, reference, alpha)
40        if key not in self._cache:
41            ds = DeseqStats(
42                self.dds,
43                contrast=[self.metadata.group_col, treated, reference],
44                alpha=alpha,
45            )
46            _ = ds.summary()
47            self._cache[key] = ds.results_df
48            
49        return self._cache[key]
50
51    def volcano(self, treated, reference, alpha=0.05, lfc=1.0, ax=None):
52        import matplotlib.pyplot as plt
53
54        df = self.results(treated, reference).dropna(subset=["padj"])
55        x = df["log2FoldChange"]
56        y = -np.log10(df["padj"].clip(lower=1e-300))
57        sig = (df["padj"] < alpha) & (x.abs() >= lfc)
58        up, down = sig & (x > 0), sig & (x < 0)
59
60        lut = group_colors(self.metadata)
61        ax = ax or plt.subplots(figsize=(6, 5))[1]
62        ax.scatter(x[~sig], y[~sig], s=6, c="lightgrey", rasterized=True)
63        ax.scatter(x[up], y[up], s=8, color=lut[treated], label=treated)
64        ax.scatter(x[down], y[down], s=8, color=lut[reference], label=reference)
65        ax.axhline(-np.log10(alpha), ls="--", lw=.7, c="k")
66        ax.axvline(lfc, ls="--", lw=.7, c="k")
67        ax.axvline(-lfc, ls="--", lw=.7, c="k")
68        ax.set(xlabel="log2 fold change", ylabel="-log10 adj. p",
69               title=f"Differential Expression {treated} vs {reference}")
70        top = df[sig].nsmallest(15, "padj")
71        for gene, row in df[sig].nsmallest(15, "padj").iterrows():
72            gx = row["log2FoldChange"]
73            gy = -np.log10(max(row["padj"], 1e-300))
74            ax.annotate(str(gene), (gx, gy),
75                        xytext=(8 if gx > 0 else -8, 0),
76                        textcoords="offset points",
77                        ha="left" if gx > 0 else "right", va="center",
78                        fontsize=7, alpha=0.85,
79                        arrowprops=dict(arrowstyle="-", lw=0.4, color="grey"))
80        ax.legend(frameon=False, fontsize=7)
81        annotate_study(ax, self.metadata)
82        _save(ax.figure, f"volcano_{treated}_{reference}.png", "results")
83        return ax
84
85

class DifferentialExpression: View Source

11class DifferentialExpression:
12
13    def __init__(self, count_matrix, metadata: Metadata, design=None):
14        counts = getattr(count_matrix, "counts", count_matrix)
15        self.metadata = metadata.align(counts)
16        self.counts = counts[self.metadata.samples]    
17        self.design = design or f"~{self.metadata.group_col}"
18        self.dds = None
19        self._cache = {}
20
21    @classmethod
22    def from_sra_run_table(cls, count_matrix, path, keep=None, **kwargs):
23        md = Metadata.from_sra_run_table(path, **kwargs)
24        if keep:
25            md = md.subset(keep=keep)
26        return cls(count_matrix, md)
27
28    def fit(self):
29        n = self.metadata.replicates()
30        if (n < 2).any():
31            raise ValueError(f"Need >=2 replicates per group: {n[n < 2].to_dict()}")
32
33        self.dds = DeseqDataSet(counts=self.counts.T,
34                                metadata=self.metadata.table,
35                                design=self.design)
36        self.dds.deseq2()
37        return self
38
39    def results(self, treated, reference, alpha=0.05):
40        key = (treated, reference, alpha)
41        if key not in self._cache:
42            ds = DeseqStats(
43                self.dds,
44                contrast=[self.metadata.group_col, treated, reference],
45                alpha=alpha,
46            )
47            _ = ds.summary()
48            self._cache[key] = ds.results_df
49            
50        return self._cache[key]
51
52    def volcano(self, treated, reference, alpha=0.05, lfc=1.0, ax=None):
53        import matplotlib.pyplot as plt
54
55        df = self.results(treated, reference).dropna(subset=["padj"])
56        x = df["log2FoldChange"]
57        y = -np.log10(df["padj"].clip(lower=1e-300))
58        sig = (df["padj"] < alpha) & (x.abs() >= lfc)
59        up, down = sig & (x > 0), sig & (x < 0)
60
61        lut = group_colors(self.metadata)
62        ax = ax or plt.subplots(figsize=(6, 5))[1]
63        ax.scatter(x[~sig], y[~sig], s=6, c="lightgrey", rasterized=True)
64        ax.scatter(x[up], y[up], s=8, color=lut[treated], label=treated)
65        ax.scatter(x[down], y[down], s=8, color=lut[reference], label=reference)
66        ax.axhline(-np.log10(alpha), ls="--", lw=.7, c="k")
67        ax.axvline(lfc, ls="--", lw=.7, c="k")
68        ax.axvline(-lfc, ls="--", lw=.7, c="k")
69        ax.set(xlabel="log2 fold change", ylabel="-log10 adj. p",
70               title=f"Differential Expression {treated} vs {reference}")
71        top = df[sig].nsmallest(15, "padj")
72        for gene, row in df[sig].nsmallest(15, "padj").iterrows():
73            gx = row["log2FoldChange"]
74            gy = -np.log10(max(row["padj"], 1e-300))
75            ax.annotate(str(gene), (gx, gy),
76                        xytext=(8 if gx > 0 else -8, 0),
77                        textcoords="offset points",
78                        ha="left" if gx > 0 else "right", va="center",
79                        fontsize=7, alpha=0.85,
80                        arrowprops=dict(arrowstyle="-", lw=0.4, color="grey"))
81        ax.legend(frameon=False, fontsize=7)
82        annotate_study(ax, self.metadata)
83        _save(ax.figure, f"volcano_{treated}_{reference}.png", "results")
84        return ax

DifferentialExpression(count_matrix, metadata: adesomics.io.metadata.Metadata, design=None) View Source

13    def __init__(self, count_matrix, metadata: Metadata, design=None):
14        counts = getattr(count_matrix, "counts", count_matrix)
15        self.metadata = metadata.align(counts)
16        self.counts = counts[self.metadata.samples]    
17        self.design = design or f"~{self.metadata.group_col}"
18        self.dds = None
19        self._cache = {}

metadata

counts

design

dds

@classmethod

def from_sra_run_table(cls, count_matrix, path, keep=None, **kwargs): View Source

21    @classmethod
22    def from_sra_run_table(cls, count_matrix, path, keep=None, **kwargs):
23        md = Metadata.from_sra_run_table(path, **kwargs)
24        if keep:
25            md = md.subset(keep=keep)
26        return cls(count_matrix, md)

def fit(self): View Source

28    def fit(self):
29        n = self.metadata.replicates()
30        if (n < 2).any():
31            raise ValueError(f"Need >=2 replicates per group: {n[n < 2].to_dict()}")
32
33        self.dds = DeseqDataSet(counts=self.counts.T,
34                                metadata=self.metadata.table,
35                                design=self.design)
36        self.dds.deseq2()
37        return self

def results(self, treated, reference, alpha=0.05): View Source

39    def results(self, treated, reference, alpha=0.05):
40        key = (treated, reference, alpha)
41        if key not in self._cache:
42            ds = DeseqStats(
43                self.dds,
44                contrast=[self.metadata.group_col, treated, reference],
45                alpha=alpha,
46            )
47            _ = ds.summary()
48            self._cache[key] = ds.results_df
49            
50        return self._cache[key]

def volcano(self, treated, reference, alpha=0.05, lfc=1.0, ax=None): View Source

52    def volcano(self, treated, reference, alpha=0.05, lfc=1.0, ax=None):
53        import matplotlib.pyplot as plt
54
55        df = self.results(treated, reference).dropna(subset=["padj"])
56        x = df["log2FoldChange"]
57        y = -np.log10(df["padj"].clip(lower=1e-300))
58        sig = (df["padj"] < alpha) & (x.abs() >= lfc)
59        up, down = sig & (x > 0), sig & (x < 0)
60
61        lut = group_colors(self.metadata)
62        ax = ax or plt.subplots(figsize=(6, 5))[1]
63        ax.scatter(x[~sig], y[~sig], s=6, c="lightgrey", rasterized=True)
64        ax.scatter(x[up], y[up], s=8, color=lut[treated], label=treated)
65        ax.scatter(x[down], y[down], s=8, color=lut[reference], label=reference)
66        ax.axhline(-np.log10(alpha), ls="--", lw=.7, c="k")
67        ax.axvline(lfc, ls="--", lw=.7, c="k")
68        ax.axvline(-lfc, ls="--", lw=.7, c="k")
69        ax.set(xlabel="log2 fold change", ylabel="-log10 adj. p",
70               title=f"Differential Expression {treated} vs {reference}")
71        top = df[sig].nsmallest(15, "padj")
72        for gene, row in df[sig].nsmallest(15, "padj").iterrows():
73            gx = row["log2FoldChange"]
74            gy = -np.log10(max(row["padj"], 1e-300))
75            ax.annotate(str(gene), (gx, gy),
76                        xytext=(8 if gx > 0 else -8, 0),
77                        textcoords="offset points",
78                        ha="left" if gx > 0 else "right", va="center",
79                        fontsize=7, alpha=0.85,
80                        arrowprops=dict(arrowstyle="-", lw=0.4, color="grey"))
81        ax.legend(frameon=False, fontsize=7)
82        annotate_study(ax, self.metadata)
83        _save(ax.figure, f"volcano_{treated}_{reference}.png", "results")
84        return ax