### Install edgePython Package

Source: https://github.com/pachterlab/edgepython/blob/main/README.md

Installs the edgePython package from PyPI. The optional '[all]' extra installs all dependencies.

```bash
pip install edgepython
pip install "edgepython[all]"
```

--------------------------------

### Install edgePython from Source

Source: https://github.com/pachterlab/edgepython/blob/main/README.md

Installs the edgePython package from local source code. The optional '.[all]' extra installs all dependencies.

```bash
pip install .
pip install .[all]
```

--------------------------------

### Install and Import edgePython and Libraries

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

Installs edgePython if not already present and imports necessary libraries including pandas, numpy, and matplotlib. It also defines a base URL for data retrieval.

```python
import sys, subprocess

IN_COLAB = 'google.colab' in sys.modules

try:
    import edgepython as ep
except ImportError:
    subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', 'edgepython'])
    import edgepython as ep

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

%matplotlib inline

BASE_URL = 'https://raw.githubusercontent.com/pachterlab/edgePython/main/examples/hoxa1/data'
```

--------------------------------

### Install and Import edgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/clytia/clytia_tutorial.ipynb

Installs the edgePython library if not already present and imports it along with other necessary libraries for data handling and analysis. This snippet ensures the environment is set up for the subsequent analysis.

```python
import sys, subprocess

IN_COLAB = 'google.colab' in sys.modules

try:
    import edgepython as ep
except ImportError:
    subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', 'edgepython'])
    import edgepython as ep

import os, tempfile, gzip, urllib.request
import numpy as np
import pandas as pd
import h5py
import scipy.sparse as sp
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D

%matplotlib inline
```

--------------------------------

### Quick Start: Differential Gene Expression Analysis

Source: https://github.com/pachterlab/edgepython/blob/main/README.md

Demonstrates a basic differential gene expression analysis workflow using edgePython. It includes creating a DGEList object, normalization, dispersion estimation, GLM fitting, and testing.

```python
import numpy as np
import edgepython as ep

# genes x samples count matrix
counts = np.random.poisson(lam=10, size=(1000, 6))
group = np.array(["A", "A", "A", "B", "B", "B"])

y = ep.make_dgelist(counts=counts, group=group)
y = ep.calc_norm_factors(y)
y = ep.estimate_disp(y)

design = np.column_stack([np.ones(6), (group == "B").astype(float)])
fit = ep.glm_ql_fit(y, design)
res = ep.glm_ql_ftest(fit, coef=1)
top = ep.top_tags(res, n=10)
print(top["table"].head())
```

--------------------------------

### Full RNA-seq Analysis Workflow with EdgePython

Source: https://context7.com/pachterlab/edgepython/llms.txt

Demonstrates a complete differential expression analysis workflow using edgePython, starting from simulating RNA-seq data, creating a DGEList, filtering genes, normalizing counts, estimating dispersions, fitting a quasi-likelihood generalized linear model, and testing for differential expression. It concludes with summarizing the results.

```python
import numpy as np
import pandas as pd
import edgepython as ep

# 1. Simulate RNA-seq data
np.random.seed(42)
n_genes = 2000
n_samples = 6

counts = np.random.negative_binomial(n=10, p=0.1, size=(n_genes, n_samples))
# Add 100 truly DE genes (2-fold change)
counts[:100, 3:] = (counts[:100, 3:] * 2).astype(int)

group = np.array(["Control", "Control", "Control", "Treatment", "Treatment", "Treatment"])

# Gene names
gene_names = pd.DataFrame({'GeneID': [f'Gene_{i}' for i in range(n_genes)]})

# 2. Create DGEList
y = ep.make_dgelist(counts=counts, group=group, genes=gene_names)

# 3. Filter lowly expressed genes
keep = ep.filter_by_expr(y, group=group)
y['counts'] = y['counts'][keep]
y['genes'] = y['genes'].iloc[keep].reset_index(drop=True)
print(f"Genes after filtering: {y['counts'].shape[0]}")

# 4. Normalize
y = ep.calc_norm_factors(y, method='TMM')

# 5. Estimate dispersions
design = np.column_stack([
    np.ones(n_samples),
    (group == "Treatment").astype(float)
])
y = ep.estimate_disp(y, design=design)
print(f"Common dispersion: {y['common.dispersion']:.4f}")

# 6. Fit QL model and test
fit = ep.glm_ql_fit(y, design=design)
result = ep.glm_ql_ftest(fit, coef=1)

# 7. Get results
top = ep.top_tags(result, n=20, adjust_method='BH')
print("\nTop 20 DE genes:")
print(top['table'][['logFC', 'logCPM', 'F', 'PValue', 'FDR']])

# 8. Summary
de_status = ep.decide_tests(result, p_value=0.05)
print(f"\nSummary at FDR < 0.05:")
print(f"  Up-regulated: {sum(de_status == 1)}")
print(f"  Down-regulated: {sum(de_status == -1)}")
print(f"  Not significant: {sum(de_status == 0)}")
```

--------------------------------

### Install and Import edgePython and Libraries (Python)

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_R_vs_Python.ipynb

Installs the edgePython library if not already present and imports necessary libraries such as numpy, pandas, and matplotlib. It also configures matplotlib for inline plotting and sets up plot aesthetics.

```python
import sys, subprocess

IN_COLAB = 'google.colab' in sys.modules

try:
    import edgepython as ep
except ImportError:
    subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', 'edgepython'])
    import edgepython as ep

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import copy

%matplotlib inline

plt.rcParams.update({
    'figure.dpi': 120, 'font.size': 10, 'axes.titlesize': 11,
    'axes.labelsize': 10, 'axes.spines.top': False, 'axes.spines.right': False,
})

group_colors = {
    'basalvirgin': '#E41A1C', 'basalpregnant': '#377EB8', 'basallactate': '#4DAF4A',
    'luminalvirgin': '#984EA3', 'luminalpregnant': '#FF7F00', 'luminallactate': '#A65628',
}

BASE_URL = 'https://raw.githubusercontent.com/pachterlab/edgePython/main/examples/mammary/data'
```

--------------------------------

### Install and Import edgePython and Load Data

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_tutorial.ipynb

Installs edgePython if not already present and imports necessary libraries (numpy, pandas, matplotlib). It then downloads and loads gene count data and sample information from a specified URL into pandas DataFrames. Finally, it prints the dimensions of the count matrix and the list of sample names.

```python
import sys, subprocess

IN_COLAB = 'google.colab' in sys.modules

try:
    import edgepython as ep
except ImportError:
    subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', 'edgepython'])
    import edgepython as ep

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

%matplotlib inline

# Download data files from the GitHub repository
BASE_URL = 'https://raw.githubusercontent.com/pachterlab/edgePython/main/examples/mammary/data'

counts_df = pd.read_csv(f'{BASE_URL}/countdata.tsv', sep='\t', index_col=0)
factor_df = pd.read_csv(f'{BASE_URL}/sampleinfo.tsv', sep='\t')

print(f'Count matrix: {counts_df.shape[0]} genes x {counts_df.shape[1]} samples')
print(f'Samples: {list(counts_df.columns)}')

```

--------------------------------

### Download and Read Kallisto Quantifications with Bootstraps in Python

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

This snippet downloads kallisto H5 quantifications including bootstrap samples from a URL and extracts them to a temporary directory. It then reads the data using edgePython's `read_data` function, specifying the source, format, sample labels, and group assignments. This is the first step in performing scaled analysis with bootstrap overdispersion.

```python
# Download kallisto H5 quantifications with bootstraps from CaltechDATA
import os, tempfile, zipfile, urllib.request

QUANT_URL = 'https://data.caltech.edu/records/3scww-j5644/files/GSE37704_kallisto_quantifications.zip?download=1'
quant_dir = tempfile.mkdtemp()
zip_path = os.path.join(quant_dir, 'quant.zip')
print('Downloading kallisto quantifications (411 MB)...', flush=True)
urllib.request.urlretrieve(QUANT_URL, zip_path)
with zipfile.ZipFile(zip_path, 'r') as zf:
    zf.extractall(quant_dir)
os.remove(zip_path)
quant_dir = os.path.join(quant_dir, 'quant')
print('Download complete.')

sample_names_sc = ['SRR493366', 'SRR493367', 'SRR493368',
                   'SRR493369', 'SRR493370', 'SRR493371']
paths = [f'{quant_dir}/{s}' for s in sample_names_sc]

y_sc = ep.read_data(paths, source='kallisto', format='h5',
                    labels=sample_names_sc,
                    group=np.array([0, 0, 0, 1, 1, 1]),
                    verbose=True)

overdisp = y_sc['genes']['Overdispersion'].values
print(f'\nTranscripts: {y_sc["counts"].shape[0]:,}')
print(f'Overdispersion range: {overdisp.min():.4f} — {overdisp.max():.4f}')
print(f'Overdispersion median: {np.median(overdisp):.4f}')
print(f'Transcripts with overdispersion > 1: {np.sum(overdisp > 1):,} '
      f'({100 * np.mean(overdisp > 1):.1f}%)')
```

--------------------------------

### Download and Read Kallisto H5 Quantifications in Python

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

This snippet downloads kallisto H5 quantification files from a URL, extracts them to a temporary directory, and then reads them into an EdgePython data object using ep.read_data. It specifies the source as 'kallisto' and format as 'h5', using provided sample names and group information.

```python
import os, tempfile, zipfile, urllib.request

QUANT_URL = 'https://data.caltech.edu/records/3scww-j5644/files/GSE37704_kallisto_quantifications.zip?download=1'
quant_dir = tempfile.mkdtemp()
zip_path = os.path.join(quant_dir, 'quant.zip')
print('Downloading kallisto quantifications (411 MB)...', flush=True)
urllib.request.urlretrieve(QUANT_URL, zip_path)
with zipfile.ZipFile(zip_path, 'r') as zf:
    zf.extractall(quant_dir)
os.remove(zip_path)
quant_dir = os.path.join(quant_dir, 'quant')
print('Download complete.')

sample_names = ['SRR493366', 'SRR493367', 'SRR493368',
                'SRR493369', 'SRR493370', 'SRR493371']
paths = [f'{quant_dir}/{s}' for s in sample_names]

y_h5 = ep.read_data(paths, source='kallisto', format='h5',
                     labels=sample_names, group=group, verbose=True)

overdisp = y_h5['genes']['Overdispersion'].values
print(f'Overdispersion range: {overdisp.min():.4f} - {overdisp.max():.4f}')
print(f'Overdispersion median: {np.median(overdisp):.4f}')
```

--------------------------------

### Visualize Library Sizes

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

Generates a bar plot of library sizes for each sample, colored by experimental group (scramble or HOXA1KD). This visualization helps assess potential batch effects or differences in sequencing depth.

```python
fig, ax = plt.subplots(figsize=(8, 5))
lib_sizes = d['samples']['lib.size'].values / 1e6
colors = [group_colors[g] for g in group]
ax.bar(range(len(lib_sizes)), lib_sizes, color=colors, edgecolor='white', linewidth=0.5)
ax.set_xticks(range(len(sample_names)))
ax.set_xticklabels(sample_names, rotation=45, ha='right')
ax.set_ylabel('Library size (millions)')
ax.set_title('Library Sizes')
handles = [mpatches.Patch(color=group_colors[g], label=g) for g in group_colors]
ax.legend(handles=handles, title='Group', frameon=False)
plt.tight_layout()
plt.show()
```

--------------------------------

### Normalization Methods

Source: https://github.com/pachterlab/edgepython/blob/main/README.md

Calculates normalization factors using various methods like TMM, TMMwsp, RLE, and upper-quartile. Also provides functions to get normalized expression values.

```python
# Example usage for normalization (conceptual)
# y = ep.calc_norm_factors(y, method='TMM')
# cpm_values = ep.cpm(y)
# rpkm_values = ep.rpkm(y)
# tpm_values = ep.tpm(y)
```

--------------------------------

### Perform GLM F-tests and Treatment Analysis with edgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_R_vs_Python.ipynb

This snippet demonstrates how to perform gene-wise statistical tests using edgePython's GLM functions. It includes setting up contrast vectors, running F-tests and treatment analyses, and extracting top significant tags based on p-value. Requires a fitted model object 'fit' and numpy arrays for contrasts.

```python
import numpy as np
import edgepython as ep

# Assuming 'fit' is a pre-existing fitted model object
# Example contrast vectors for basal and luminal conditions
con_basal = np.array([-1, 1, 0, 0, 0, 0], dtype=float)
con_luminal = np.array([0, 0, 0, -1, 1, 0], dtype=float)

# Perform GLM F-tests
result_basal = ep.glm_ql_ftest(fit, contrast=con_basal)
result_luminal = ep.glm_ql_ftest(fit, contrast=con_luminal)

# Extract top tags sorted by p-value
py_qlf_basal = ep.top_tags(result_basal, n=np.inf, sort_by='PValue')['table']
py_qlf_luminal = ep.top_tags(result_luminal, n=np.inf, sort_by='PValue')['table']

# Perform GLM treatment analysis with a log fold change threshold
treat_basal = ep.glm_treat(fit, contrast=con_basal, lfc=0.58)
treat_luminal = ep.glm_treat(fit, contrast=con_luminal, lfc=0.58)

# Extract top tags sorted by p-value from treatment analysis
py_treat_basal = ep.top_tags(treat_basal, n=np.inf, sort_by='PValue')['table']
py_treat_luminal = ep.top_tags(treat_luminal, n=np.inf, sort_by='PValue')['table']

print('edgePython pipeline complete.')
```

--------------------------------

### Fit Quasi-Likelihood GLM with EdgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_tutorial.ipynb

Fits a Quasi-Likelihood (QL) Generalized Linear Model (GLM) using EdgePython's `glm_ql_fit` function. It calculates and prints the prior degrees of freedom (df.prior) from the fitted model.

```python
fit = ep.glm_ql_fit(d_filt, design=design)
print(f"Prior df: {fit['df.prior']:.4f}")
```

--------------------------------

### Generate QQ Plot of Z-Statistics

Source: https://github.com/pachterlab/edgepython/blob/main/examples/clytia/clytia_tutorial.ipynb

Creates a QQ plot to visualize the distribution of Wald z-statistics against a standard normal distribution. This plot serves to verify the calibration of the NEBULA-LN model under the null hypothesis.

```python
from scipy.stats import norm

conv = fit['convergence'] == 1
z_vals = fit['coefficients'][conv, 1] / fit['se'][conv, 1]
z_vals = z_vals[np.isfinite(z_vals)]

n_z = len(z_vals)
theoretical = norm.ppf((np.arange(1, n_z + 1) - 0.5) / n_z)
observed = np.sort(z_vals)

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(theoretical, observed, c=BLUE, s=4, alpha=0.3, edgecolors='none',
           rasterized=True)
lims = [min(theoretical.min(), observed.min()) * 1.05,
        max(theoretical.max(), observed.max()) * 1.05]
ax.plot(lims, lims, 'k-', lw=0.5, alpha=0.5)
ax.set_xlim(lims); ax.set_ylim(lims)
ax.set_xlabel('Theoretical quantiles')
ax.set_ylabel('Observed z-statistics')
ax.set_title('NEBULA-LN QQ Plot')plt.tight_layout()plt.show()

print(f'Number of z-statistics: {n_z:,}')
```

--------------------------------

### Create DGEList Object and Visualize Library Sizes

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_tutorial.ipynb

Prepares data for differential expression analysis by creating a DGEList object using edgePython's `make_dgelist` function. It then generates a bar plot visualizing the library sizes for each sample, colored according to their group, to assess sample-wise sequencing depth.

```python
group = factor_df.set_index('SampleName').loc[counts_df.columns, 'CellTypeStatus'].values
gene_ids = counts_df.index.values.astype(str)
sample_names = counts_df.columns.tolist()

genes_df = pd.DataFrame({'GeneID': gene_ids})
d = ep.make_dgelist(counts_df.values.astype(float), group=group, genes=genes_df)
print(f"DGEList: {d['counts'].shape[0]} genes, {d['counts'].shape[1]} samples")
print(f"Library sizes: {d['samples']['lib.size'].values}")

fig, ax = plt.subplots(figsize=(10, 5))
lib_sizes = d['samples']['lib.size'].values / 1e6
colors = [group_colors[g] for g in group]
ax.bar(range(len(lib_sizes)), lib_sizes, color=colors, edgecolor='white', linewidth=0.5)
ax.set_xticks(range(len(sample_names)))
ax.set_xticklabels(sample_names, rotation=45, ha='right')
ax.set_ylabel('Library size (millions)')
ax.set_title('Library Sizes by Sample')
handles = [mpatches.Patch(color=group_colors[g], label=g) for g in group_colors]
ax.legend(handles=handles, title='Group', bbox_to_anchor=(1.02, 1), loc='upper left', frameon=False)
plt.tight_layout()
plt.show()

```

--------------------------------

### Visualize Log2 CPM Density Before and After Filtering with EdgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

Generates density plots for Log2 CPM values before and after transcript filtering. It uses 'scipy.stats.gaussian_kde' for density estimation and 'matplotlib' for plotting. The function visualizes the impact of filtering on the overall expression distribution.

```python
from scipy.stats import gaussian_kde

full_lib_size = d['samples']['lib.size'].values
logcpm_before = ep.cpm(d, log=True)
logcpm_after = ep.cpm(d['counts'][keep, :], lib_size=full_lib_size, log=True)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5), sharey=True, sharex=True)
x_grid = np.linspace(-5, 18, 500)

for i in range(len(sample_names)):
    c = group_colors[group[i]]
    kde_before = gaussian_kde(logcpm_before[:, i], bw_method=0.2)
    ax1.plot(x_grid, kde_before(x_grid), color=c, alpha=0.6, linewidth=1.2)
    kde_after = gaussian_kde(logcpm_after[:, i], bw_method=0.2)
    ax2.plot(x_grid, kde_after(x_grid), color=c, alpha=0.6, linewidth=1.2)

ax1.set_xlabel('Log$_2$ CPM'); ax1.set_ylabel('Density')
ax1.set_title(f'A. Raw ({logcpm_before.shape[0]:,} transcripts)')
ax1.axvline(x=np.log2(0.5), color='black', linestyle='--', linewidth=0.8,
            alpha=0.6, label='CPM = 0.5')
ax1.legend(frameon=False, fontsize=8)
ax2.set_xlabel('Log$_2$ CPM')
ax2.set_title(f'B. Filtered ({logcpm_after.shape[0]:,} transcripts)')plt.tight_layout()plt.show()
```

--------------------------------

### Create Design Matrix with Patsy and EdgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_tutorial.ipynb

Generates a design matrix for differential expression analysis using Patsy and EdgePython. The model '~ 0 + group' creates a column for each group mean. It returns the design matrix and its column names.

```python
import patsy
design_df = pd.DataFrame({'group': group})
design = ep.model_matrix('~ 0 + group', design_df)
col_names = list(patsy.dmatrix('~ 0 + group', design_df, return_type='dataframe').columns)
print('Design matrix columns:', col_names)
print(f'Shape: {design.shape}')
pd.DataFrame(design, columns=col_names, index=sample_names)
```

--------------------------------

### Generate Quasi-Likelihood Dispersion Plot using edgeR and edgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_R_vs_Python.ipynb

This snippet generates a Quasi-Likelihood Dispersion plot, comparing dispersion estimates from edgeR and edgePython. It visualizes raw, squeezed, and trended dispersion on a quarter-root scale against average log2 CPM. Requires numpy, pandas, and matplotlib.

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Assuming r_ql_disp, d_filt, fit are defined

# R QL dispersion data
r_ql = r_ql_disp.set_index('GeneID')
r_ql_genes = r_ql.index.intersection(pd.Index(d_filt['genes']['GeneID'].values))

r_s2_raw = r_ql.loc[r_ql_genes, 's2.raw'].values
r_s2_post = r_ql.loc[r_ql_genes, 's2.post'].values
r_s2_prior = r_ql.loc[r_ql_genes, 's2.prior'].values
r_logcpm_ql = r_ql.loc[r_ql_genes, 'AveLogCPM'].values

# Python QL dispersion data
py_logcpm_ql = fit['AveLogCPM']
py_df_res = np.asarray(fit.get('df.residual.adj', fit['df.residual']), dtype=np.float64)
py_dev = np.asarray(fit.get('deviance.adj', fit['deviance']), dtype=np.float64)
py_s2_raw = py_dev / np.maximum(py_df_res, 1e-8)
py_s2_raw[py_df_res < 1e-8] = 0
py_s2_post = fit['s2.post']
py_s2_prior = fit['s2.prior']

# Shared axis bounds
all_raw = np.concatenate([r_s2_raw ** 0.25, py_s2_raw ** 0.25])
all_x_ql = np.concatenate([r_logcpm_ql, py_logcpm_ql])
x_lim_ql = (all_x_ql.min() - 0.5, all_x_ql.max() + 0.5)
y_lim_ql = (0, np.percentile(all_raw[all_raw > 0], 99.5) * 1.1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# R panel
order_r = np.argsort(r_logcpm_ql)
ax1.scatter(r_logcpm_ql, r_s2_raw ** 0.25, s=2, alpha=0.2, color='black', rasterized=True, label='Raw')
ax1.scatter(r_logcpm_ql, r_s2_post ** 0.25, s=2, alpha=0.2, color='red', rasterized=True, label='Squeezed')
ax1.plot(r_logcpm_ql[order_r], r_s2_prior[order_r] ** 0.25, color='blue', lw=2, label='Trend')
ax1.set_xlabel('Average Log$_2$ CPM'); ax1.set_ylabel('Quarter-Root Mean Deviance')
ax1.set_title('edgeR'); ax1.legend(frameon=False, fontsize=8)
ax1.set_xlim(x_lim_ql); ax1.set_ylim(y_lim_ql)

# Python panel
order_py = np.argsort(py_logcpm_ql)
ax2.scatter(py_logcpm_ql, py_s2_raw ** 0.25, s=2, alpha=0.2, color='black', rasterized=True, label='Raw')
ax2.scatter(py_logcpm_ql, py_s2_post ** 0.25, s=2, alpha=0.2, color='red', rasterized=True, label='Squeezed')
ax2.plot(py_logcpm_ql[order_py], py_s2_prior[order_py] ** 0.25, color='blue', lw=2, label='Trend')
ax2.set_xlabel('Average Log$_2$ CPM')
ax2.set_title('edgePython'); ax2.legend(frameon=False, fontsize=8)
ax2.set_xlim(x_lim_ql); ax2.set_ylim(y_lim_ql)

fig.suptitle('Quasi-Likelihood Dispersion', fontsize=13)
plt.tight_layout()
plt.show()
```

--------------------------------

### Plot Quasi-Likelihood Dispersion with EdgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_tutorial.ipynb

Visualizes the dispersion estimates from the Quasi-Likelihood (QL) GLM fit using EdgePython's `plot_ql_disp` function. This plot helps in understanding the dispersion characteristics of the data.

```python
fig, ax = ep.plot_ql_disp(fit)
ax.set_title('Quasi-Likelihood Dispersion')
plt.tight_layout()
plt.show()
```

--------------------------------

### Load Data and R Results (Python)

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

Loads count data, sample information, and pre-computed results from edgeR (normalization factors, dispersions, QL F-test, exact test, TREAT results, DTU results, scaled analysis) from CSV files. Also loads a transcript-to-gene mapping for DTU analysis.

```python
# Count data and sample info
counts_df = pd.read_csv(f'{BASE_URL}/counts.tsv', sep='\t', index_col=0)
sample_info = pd.read_csv(f'{BASE_URL}/sample_info.tsv', sep='\t')

samples = counts_df.columns.tolist()
condition = sample_info['group'].values

# R results
r_norm = pd.read_csv(f'{BASE_URL}/R_norm_factors.csv')
r_disp = pd.read_csv(f'{BASE_URL}/R_dispersions.csv')
r_qlf = pd.read_csv(f'{BASE_URL}/R_qlf_all_results.csv', index_col=0)
r_exact = pd.read_csv(f'{BASE_URL}/R_exact_all_results.csv', index_col=0)
r_ql_fit = pd.read_csv(f'{BASE_URL}/R_ql_fit_details.csv')

# TREAT results (4 thresholds)
r_treat_lfc10 = pd.read_csv(f'{BASE_URL}/R_treat_lrt_hoxa1_lfc10.csv', index_col=0)
r_treat_lfc12 = pd.read_csv(f'{BASE_URL}/R_treat_lrt_hoxa1_lfc12.csv', index_col=0)
r_treat_lfc15 = pd.read_csv(f'{BASE_URL}/R_treat_lrt_hoxa1_lfc15.csv', index_col=0)
r_treat_worst = pd.read_csv(f'{BASE_URL}/R_treat_lrt_hoxa1_worst.csv', index_col=0)

# DTU results
r_dtu_gene = pd.read_csv(f'{BASE_URL}/R_dtu_gene_results.csv')
r_dtu_exon = pd.read_csv(f'{BASE_URL}/R_dtu_exon_results.csv')

# Scaled analysis results
r_scaled_qlf = pd.read_csv(f'{BASE_URL}/R_scaled_qlf_all_results.csv', index_col=0)
r_scaled_disp = pd.read_csv(f'{BASE_URL}/R_scaled_dispersions.csv')

# Transcript-to-gene mapping (for DTU)
t2g = pd.read_csv(f'{BASE_URL}/t2g_mapping.csv')

print(f'Count matrix: {counts_df.shape[0]:,} transcripts x {counts_df.shape[1]} samples')
print(f'R dispersions: {r_disp.shape[0]:,} transcripts (after filtering)')
print(f'R DTU genes: {r_dtu_gene.shape[0]:,}')
print(f'R scaled QLF: {r_scaled_qlf.shape[0]:,} transcripts')
```

--------------------------------

### Visualize Log2 CPM Distribution Before and After TMM Normalization with EdgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

Generates box plots to compare the Log2 CPM distribution before and after TMM normalization. It uses 'edgepython' for CPM calculation and 'matplotlib' for plotting. The box plots are colored by group to visualize the effect of normalization on sample distributions.

```python
import copy

logcpm_norm = ep.cpm(d_filt, log=True)
d_unnorm = copy.deepcopy(d_filt)
d_unnorm['samples']['norm.factors'] = 1.0
logcpm_unnorm = ep.cpm(d_unnorm, log=True)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
for ax, data, title in [(ax1, logcpm_unnorm, 'Before TMM Normalization'),
                         (ax2, logcpm_norm, 'After TMM Normalization')]:
    bp = ax.boxplot([data[:, i] for i in range(data.shape[1])],
                    patch_artist=True, medianprops=dict(color='black', linewidth=1.5),
                    flierprops=dict(marker='.', markersize=1, alpha=0.3))
    for patch, g in zip(bp['boxes'], group):
        patch.set_facecolor(group_colors[g]); patch.set_alpha(0.7)
    ax.set_xticklabels(sample_names, rotation=45, ha='right')
    ax.set_ylabel('Log$_2$ CPM'); ax.set_title(title)
plt.tight_layout()
plt.show()
```

--------------------------------

### Create DGEList Object with edgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_tutorial.ipynb

Prepares transcript IDs and creates a DGEList object using edgePython's `make_dgelist` function. It converts counts to float and includes sample group information. The output shows the dimensions of the DGEList and library sizes.

```python
transcript_ids = counts_df.index.values.astype(str)
genes_df = pd.DataFrame({'TranscriptID': transcript_ids})
d = ep.make_dgelist(counts_df.values.astype(float), group=group, genes=genes_df)
print(f'DGEList: {d["counts"].shape[0]:,} transcripts, {d["counts"].shape[1]} samples')
print(f'Library sizes: {d["samples"]["lib.size"].values.astype(int)}')
```

--------------------------------

### Visualize Library Sizes: edgeR vs edgePython

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_R_vs_Python.ipynb

This Python snippet generates a side-by-side bar plot comparing library sizes (in millions) between edgeR and edgePython. It uses matplotlib for plotting and requires pre-loaded dataframes 'r_norm' and 'd', along with sample names and group information. The plot helps visually confirm consistency in library size calculations.

```python
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np

# Assuming 'r_norm' (edgeR results) and 'd' (edgePython data) are loaded
# 'group_colors' and 'sample_names' are also assumed to be defined

r_lib = r_norm['lib.size'].values / 1e6
py_lib = d['samples']['lib.size'].values / 1e6

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5), sharey=True)
for ax, lib, title in [(ax1, r_lib, 'edgeR'), (ax2, py_lib, 'edgePython')]:
    colors = [group_colors[g] for g in group]
    ax.bar(range(len(lib)), lib, color=colors, edgecolor='white', linewidth=0.5)
    ax.set_xticks(range(len(sample_names)))
    ax.set_xticklabels(sample_names, rotation=45, ha='right')
    ax.set_title(title)
ax1.set_ylabel('Library size (millions)')

handles = [mpatches.Patch(color=group_colors[g], label=g) for g in group_colors]
ax2.legend(handles=handles, title='Group', bbox_to_anchor=(1.02, 1),
           loc='upper left', frameon=False, fontsize=8)
fig.suptitle('Library Sizes', fontsize=13)
plt.tight_layout()
plt.show()
```

--------------------------------

### Run edgePython LRT TREAT Pipeline (Python)

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

Initializes the LRT (Likelihood Ratio Test) fit for the TREAT (Test for REgulatory ANalysis Toolkit) pipeline using edgePython. This is a precursor to performing the TREAT analysis, which is used for testing differential transcript usage.

```python
# LRT fit for TREAT
fit_lrt = ep.glm_fit(y, design)
```

--------------------------------

### Compare R and Python QL F-test P-values

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

This snippet compares the p-values from QL F-tests in R and Python. It extracts p-values for common transcripts, filters for positive values, converts them to -log10 scale, and visualizes their correlation. It also calculates the Pearson correlation coefficient. Dependencies include pandas, numpy, and matplotlib.

```python
r_p_ql = r_qlf.loc[common_ql, 'PValue'].values
py_p_ql = py_qlf_tab.loc[common_ql, 'PValue'].values
mask_ql = (r_p_ql > 0) & (py_p_ql > 0)

r_lp_ql = -np.log10(r_p_ql[mask_ql])
py_lp_ql = -np.log10(py_p_ql[mask_ql])

fig, ax = plt.subplots(figsize=(6, 6))
scatter_compare(ax, r_lp_ql, py_lp_ql, r'(f) QL F-test $-log_{10}(p)$')
plt.tight_layout()
plt.show()

corr_ql = np.corrcoef(r_lp_ql, py_lp_ql)[0, 1]
print(f'Pearson r: {corr_ql:.8f}')
```

--------------------------------

### Run TREAT Analysis with 4 LFC Thresholds in Python

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

This snippet demonstrates how to run the TREAT analysis using EdgePython's glm_treat function with four different log2 fold change (LFC) thresholds: 1.0, log2(1.2), log2(1.5), and a worst-case scenario. It then uses top_tags to extract the top transcripts for each threshold and maps their indices to target IDs.

```python
import math

py_treat_lfc10 = ep.glm_treat(fit_lrt, coef=1, lfc=1.0)
py_treat_lfc12 = ep.glm_treat(fit_lrt, coef=1, lfc=math.log2(1.2))
py_treat_lfc15 = ep.glm_treat(fit_lrt, coef=1, lfc=math.log2(1.5))
py_treat_worst = ep.glm_treat(fit_lrt, coef=1, lfc=math.log2(1.5), null='worst.case')

py_tt_lfc10 = ep.top_tags(py_treat_lfc10, n=ngenes, sort_by='none')['table']
py_tt_lfc12 = ep.top_tags(py_treat_lfc12, n=ngenes, sort_by='none')['table']
py_tt_lfc15 = ep.top_tags(py_treat_lfc15, n=ngenes, sort_by='none')['table']
py_tt_worst = ep.top_tags(py_treat_worst, n=ngenes, sort_by='none')['table']

# Map indices to transcript IDs
for tt in [py_tt_lfc10, py_tt_lfc12, py_tt_lfc15, py_tt_worst]:
    tt.index = target_ids_filt[tt.index.values]

print('LRT TREAT pipeline complete.')
```

--------------------------------

### Map Python QL Results to Transcript IDs and Compare Log-Fold-Changes

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

This code maps Python Quasi-Likelihood (QL) F-test results to transcript IDs and compares the log-fold-changes (logFC) with R results. It identifies common transcripts, extracts logFC values, and visualizes their correlation. It also calculates the Pearson correlation coefficient and the maximum absolute difference. Dependencies include pandas, numpy, and matplotlib.

```python
# Map Python QL results to transcript IDs
py_qlf_tab = py_qlf_all['table'].copy()
py_qlf_tab.index = target_ids_filt[py_qlf_tab.index.values]

r_qlf.index = r_qlf.index.astype(str)
py_qlf_tab.index = py_qlf_tab.index.astype(str)
common_ql = r_qlf.index.intersection(py_qlf_tab.index)

r_lfc = r_qlf.loc[common_ql, 'logFC'].values
py_lfc = py_qlf_tab.loc[common_ql, 'logFC'].values

fig, ax = plt.subplots(figsize=(6, 6))
scatter_compare(ax, r_lfc, py_lfc, '(c) QL F-test logFC')
plt.tight_layout()
plt.show()

corr_lfc = np.corrcoef(r_lfc, py_lfc)[0, 1]
print(f'Pearson r: {corr_lfc:.8f}')
print(f'Max |diff|: {np.max(np.abs(r_lfc - py_lfc)):.2e}')
```

--------------------------------

### Summary DataFrame Generation with edgePython Results

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

Creates a pandas DataFrame summarizing various statistical metrics computed by edgeR and edgePython. It calculates the difference or correlation for metrics like TMM normalization factors, Tagwise BCV, QL F-test logFC, posterior s2, exact test p-values, and others, presenting them with specific formatting.

```python
summary = pd.DataFrame({
    'Panel': ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)', '(g)', '(n)', '(o)', '(p)'],
    'Metric': [
        'TMM norm factors (max |diff|)',
        'Tagwise BCV (r)',
        'QL F-test logFC (r)',
        'QL posterior s2.post (r)',
        'Exact test -log10(p) (r)',
        'QL F-test -log10(p) (r)',
        'LRT TREAT -log10(p) pooled (r)',
        'DTU gene Simes -log10(p) (r)',
        'DTU exon -log10(p) (r)',
        'Scaled QL -log10(p) (r)',
    ],
    'Value': [
        f'{np.max(np.abs(r_nf - py_nf)):.2e}',
        f'{corr_bcv:.8f}',
        f'{corr_lfc:.8f}',
        f'{corr_s2:.8f}',
        f'{corr_et:.8f}',
        f'{corr_ql:.8f}',
        f'{corr_treat:.8f}',
        f'{corr_dtu_gene:.8f}',
        f'{corr_dtu_exon:.8f}',
        f'{corr_scaled:.8f}',
    ]
})
summary
```

--------------------------------

### Read Quantification Output with EdgePython (kallisto/Salmon)

Source: https://context7.com/pachterlab/edgepython/llms.txt

Reads transcript-level quantification data from kallisto or Salmon output directories. This function can also estimate overdispersion using bootstrap information if available in the output.

```python
import edgepython as ep

# Read kallisto output with bootstrap information
# paths = ['sample1/', 'sample2/', 'sample3/']
# result = ep.catch_kallisto(paths, verbose=True)
#
# print("Counts shape:", result['counts'].shape)
# print("Annotation columns:", result['annotation'].columns.tolist())
# # If bootstraps available: ['Length', 'EffectiveLength', 'Overdispersion']

# Read Salmon output
# result = ep.catch_salmon(paths, verbose=True)
```

--------------------------------

### Calculate Log-CPM with Effective Library Sizes in Python

Source: https://github.com/pachterlab/edgepython/blob/main/examples/mammary/mouse_mammary_R_vs_Python.ipynb

This Python code defines a helper function `logcpm_with_lib` to compute log2-Counts Per Million (Log-CPM) values. It takes raw counts and effective library sizes as input. The function is used to calculate Log-CPM before and after TMM normalization, comparing values derived from edgeR and edgePython normalization factors.

```python
import numpy as np

# Assuming 'd_filt' and 'r_norm' are loaded, containing counts and normalization factors

filt_counts = d_filt['counts']

# R effective library sizes (lib.size * norm.factors)
r_lib_raw = r_norm['lib.size'].values
r_nf_vals = r_norm['norm.factors'].values
r_eff_lib = r_lib_raw * r_nf_vals

# Python effective library sizes
py_lib_raw = d_filt['samples']['lib.size'].values
py_nf_vals = d_filt['samples']['norm.factors'].values
py_eff_lib = py_lib_raw * py_nf_vals

def logcpm_with_lib(counts, lib_sizes):
    """log2-CPM with given effective library sizes."""
    lib = lib_sizes / 1e6
    return np.log2(counts / lib[None, :] + 0.5 / lib[None, :])

# Before normalization: same for both (norm.factors = 1)
logcpm_unnorm = logcpm_with_lib(filt_counts, py_lib_raw)

# Example of calculating logCPM using R's normalization factors (if needed)
# logcpm_r_normed = logcpm_with_lib(filt_counts, r_eff_lib)

# Example of calculating logCPM using Python's normalization factors (if needed)
# logcpm_py_normed = logcpm_with_lib(filt_counts, py_eff_lib)

```

--------------------------------

### Pool and Compare LRT TREAT P-values from R and Python

Source: https://github.com/pachterlab/edgepython/blob/main/examples/hoxa1/hoxa1_R_vs_Python.ipynb

This code pools p-values from multiple LRT TREAT conditions (logFC thresholds) calculated in R and Python. It iterates through different conditions, finds common transcripts, collects p-values, filters for positive values, transforms them to -log10 scale, and visualizes their correlation. It reports the number of points and the Pearson correlation coefficient. Dependencies include numpy and matplotlib.

```python
# Pool p-values across 4 LRT TREAT conditions
r_treat_all = [r_treat_lfc10, r_treat_lfc12, r_treat_lfc15, r_treat_worst]
py_treat_all = [py_tt_lfc10, py_tt_lfc12, py_tt_lfc15, py_tt_worst]
labels = ['lfc=1.0', 'lfc=log2(1.2)', 'lfc=log2(1.5)', 'worst.case']

r_pvals_pooled, py_pvals_pooled = [], []
for r_df, py_df in zip(r_treat_all, py_treat_all):
    r_df.index = r_df.index.astype(str)
    py_df.index = py_df.index.astype(str)
    common_tr = r_df.index.intersection(py_df.index)
    r_pvals_pooled.extend(r_df.loc[common_tr, 'PValue'].values)
    py_pvals_pooled.extend(py_df.loc[common_tr, 'PValue'].values)

r_pvals_pooled = np.array(r_pvals_pooled)
py_pvals_pooled = np.array(py_pvals_pooled)
mask_tr = (r_pvals_pooled > 0) & (py_pvals_pooled > 0)

r_lp_tr = -np.log10(r_pvals_pooled[mask_tr])
py_lp_tr = -np.log10(py_pvals_pooled[mask_tr])

fig, ax = plt.subplots(figsize=(6, 6))
scatter_compare(ax, r_lp_tr, py_lp_tr,
                r'(g) LRT TREAT $-log_{10}(p)$ (pooled)',
                s=6, alpha=0.5)
plt.tight_layout()
plt.show()

corr_treat = np.corrcoef(r_lp_tr, py_lp_tr)[0, 1]
print(f'Points: {np.sum(mask_tr):,}')
print(f'Pearson r: {corr_treat:.8f}')
```