Description

1. Accessing the Homepage

• Open a browser (recommended: Chrome / Firefox / Edge) and navigate to the software's Web UI homepage.
• The homepage displays entry points for primary functional modules:
  a) ATRA Data Processing
  b) Virtual Cell Generation
  c) Pharmacodynamic Analysis
  d) Mechanism Interpretation
• Users may upload their own data or use the system's sample data for demonstration.

2. ATRA Data Processing Module

Function: Preprocess and perform differential gene analysis on single-cell RNA sequencing (scRNA-seq) data

Steps:

• a) Click ATRA Data Processing in the navigation bar.
• b) Upload data files (supports CSV/TSV/Excel formats; rows represent genes, columns represent cells).
• c) Configure data processing parameters:
  • QC Threshold: Minimum gene count, mitochondrial gene proportion filtering.
  • Normalization Method: TPM / RPKM / CPM.
  • Batch Effect Correction: Enable PCA/ComBat.
• d) Click Run Analysis. The system will output:
  • Differential gene table (includes logFC, p-value, significance flag).
  • Heatmap (expression distribution of differential genes across samples).
  • PCA/UMAP plot (distribution differences of cells before and after ATRA treatment).

3. Virtual Cell Generation Module

Function: Simulates transcriptional changes in cells post-ATRA treatment using machine learning and mathematical modeling

Steps:

• a) Click Virtual Cell Generation in the navigation bar.
• b) Input treatment conditions:
  o Drug concentration (µM): e.g., 0.5, 1.0, 2.0.
  o Target gene downregulation ratio (e.g., 50%).
  o Select simulation method: Random Forest / LASSO / Logistic Regression.
• c) Click Generate Virtual Cells.
• d) System Output:
  o Position of “virtual cells” in the expression matrix.
  o PCA/UMAP projection plots showing relative positions of virtual and real cells.
  o Simulated expression curves for differentially regulated genes (dose-dependent).

4. Pharmacodynamic Analysis Module

Function: Perform pharmacodynamic modeling and response prediction for differentially expressed genes and target genes

Procedure:

• a) Click Pharmacodynamic Analysis in the navigation bar.
• b) Select the gene set for analysis (can be imported from DEG results).
• c) Set simulation parameters:
  o Dose gradient (e.g., 0.1, 0.5, 1.0, 5.0 µM).
  o Effect Type: Upregulation / Downregulation.
• d) Click Run Pharmacodynamic Analysis.
• e) System Output:
  o Dose-Response Curve.
  o Population Sensitivity Distribution (Violin Plot/Box Plot).
  o Gene Expression Heatmap Before and After Drug Treatment.

5. Mechanism Analysis Module

Function: Perform enrichment analysis on differentially expressed gene sets to reveal ATRA's mechanism of action

Procedure:

• a) Click Mechanism Analysis in the navigation bar.
• b) Import the differentially expressed gene list (from the ATRA Data Processing Module).
• c) Select enrichment analysis methods:
  o GSEA (Gene Set Enrichment Analysis).
  o KEGG pathway enrichment analysis.
  o Reactome pathway analysis.
• d) Set parameters:
  o Significance threshold (p<0.05, FDR<0.25).
  o Gene set database (KEGG, MSigDB, Reactome).
• e) Click Run Mechanism Analysis.
• f) System outputs:
  o Ranked enrichment pathway list.
  o Bar chart (pathway significance comparison).
  o Network diagram (gene-pathway relationships).

6. Final Output

• After completing the above four modules, users can:
  • Generate a complete report (PDF/Word) with one click, including:
  o Data processing workflow
  o Differentially expressed gene list
  o Virtual cell simulation results
  o Pharmacodynamic analysis curves
  o Pathway mechanism analysis
• Export interactive charts (PNG/SVG) and gene expression matrices (CSV/Excel) for further research or publication.

1) Module Function

The module is the system’s first step. It cleans, standardizes, and performs differential gene analysis (DEG) on scRNA-seq inputs, producing high-quality data for virtual cell generation, pharmacodynamic analysis, and mechanism interpretation.

• Remove noise and low-quality cells to ensure reliability.
• Normalize cells and batches onto a common scale for comparability.
• Identify DEGs before/after ATRA and between HCC vs. normal liver cells.

2) Module Assumptions

• Cellular heterogeneity: HCC vs. normal hepatocytes show meaningful transcriptional differences detectable by scRNA-seq.
• Drug action: ATRA modulates genes by up/down regulation with statistically significant shifts.
• Independence: Cells are treated as independent samples for DEG testing.
• Noise & batch effects: Present but mitigated by normalization, dimensionality reduction, and batch correction.

3) Mathematical Model

Data matrix. Let \(X \in \mathbb{R}^{m\times n}\) be raw expression (genes × cells):

Rows are genes \((g_i)\); columns are cells \((c_j)\).

Standardization. Log-CPM (library-size normalized then log-transform):

This addresses sequencing-depth variation and non-linear scaling.

DE testing. Group means, effect size, and hypothesis test with FDR control:

Screening rule. A gene is called a DEG if both hold:

4) Data Sources

• Experimental: scRNA-seq from HCC lines (HepG2/Huh7) and normal hepatocytes; ATRA at varying concentrations with post-treatment profiles.
• Public: GEO, TCGA, Single Cell Portal / Human Cell Atlas for controls or model references.

5) Implementation

Input formats: CSV/TSV/Excel (genes × cells), and AnnData (.h5ad) compatible with Scanpy/Seurat.

Required metadata: expression matrix + grouping labels (ATRA vs. control; HCC vs. normal).

Workflow — QC. Remove low-quality cells/genes; limit mitochondrial proportion.

import scanpy as sc
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs['percent_mito'] < 0.2]

Workflow — normalization. CPM/TPM; log; scaling.

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.scale(adata, max_value=10)

Workflow — batch correction. PCA + ComBat; or Harmony/Scanorama.

sc.tl.pca(adata, svd_solver='arpack')
# then apply Harmony / Scanorama if needed

Workflow — DEG testing. T vs C; HCC vs normal; Wilcoxon/t-test/logistic regression.

sc.tl.rank_genes_groups(adata, groupby='treatment', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)

Screening criteria. p < 0.05; log2FC > 1 or < −1.

6) Output Overview

1. DEG table: gene, log2FC, p-value, FDR, and significance marks.
2. Visualizations: PCA/UMAP (pre/post ATRA), heatmaps, violins for key targets.
3. Standardized matrix: QC’d and batch-corrected for downstream modules.

DEG Significance Illustration

Annotates significance with * / ** / *** right after p/FDR thresholds.

1) Purpose

The Virtual Cell Generation Module simulates cellular transcriptomic responses to varying ATRA treatment conditions (dose/time/downregulation ratio) based on processed single-cell RNA-seq data, generating “virtual cells” or “virtual populations.” Its primary applications include:

• Conducting in-silico drug intervention simulations when experimental coverage is insufficient;
• Expanding training datasets to enhance the robustness of downstream models (DEG, enrichment, classification);
• Comparing effects of different dosages or downregulation strategies on cell populations/single cells;
• Providing visual (PCA/UMAP) and statistical (dose-response curves, population sensitivity distributions) outputs for decision-making.

The module interfaces directly with the ATRA data processing module: input is QC/normalized/batch-corrected AnnData (or equivalent matrices).

2) Module Assumptions

• Response modelability: ATRA’s impact on gene expression can be approximated by statistical/mechanistic models (multiplicative scaling, Hill curves, or learned mappings).
• Local perturbation: Drug effects primarily alter a subset of key target genes, expandable to multi-gene coupling.
• Stability: Under identical treatments, populations are reproducible and can be Monte-Carlo simulated for uncertainty.
• Data availability: There exists at least one dataset or credible target/effect prior for calibration or supervised learning.

3) Mathematical Model

A. Per-Gene Multiplicative Perturbation

For gene \(g\) and cell \(c\), let the original expression be \(X^{\text{pre}}_{g,c}\). The virtual expression:

The perturbation coefficient can be constant or stochastic:

To preserve count statistics, optional resampling:

B. Dose-Response (Hill/Logistic) Model (for multi-dose data)

For each gene \(g\), the dose-response:

Virtual expression as additive or multiplicative modulation about a baseline \(\bar X_g\):

where \(d\) is drug concentration, \(\kappa_g\) scales gene-specific sensitivity, and \(\varepsilon_{g,c}\) captures individual variation.

C. Latent-Space Perturbation

Map each cell to latent space \(z=f(x)\) (PCA/scVI/autoencoder). Define population shift:

Shift and sample a virtual cell, then decode back:

Latent perturbation better preserves covariance structures; multiplicative is simple/transparent; Hill fits enable pharmacologic interpretation (EC50/Emax/n).

4) Data Sources

• Required: AnnData preprocessed by the ATRA module (with adata.X or adata.layers['normalized']), and adata.obs including treatment, cell_type.
• Optional: Multi-dose datasets for Hill fitting; similar public scRNA-seq (e.g., GEO) for calibration; predefined target lists or priors.

5) Implementation

import scanpy as sc
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs['percent_mito'] < 0.2]

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.scale(adata, max_value=10)

sc.tl.pca(adata, svd_solver='arpack')
# then apply Harmony / Scanorama if needed

sc.tl.rank_genes_groups(adata, groupby='treatment', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)

5.1 Input

• adata: AnnData (QC/normalized/batch-corrected); raw counts require normalization flags.
• targets: Candidate target gene list or internally screened (RF/LASSO).
• mode: 'multiplicative' / 'dose_response' / 'latent'.
• params: alpha (down-regulation ratio), dose, Emax/EC50/n, etc. (per-gene or shared).
• n_replicates: Monte-Carlo replicates per real cell.
• random_state: Random seed.

5.2 Processing Pipeline

1. Input validation (targets ∈ adata.var_names; layer availability).
2. Select population (e.g., cell_type == 'HepG2').
3. Construct coefficients:
   • multiplicative: \(s=1-\alpha\) or \(s\sim\mathrm{LogNormal}\).
   • dose_response: fit per-gene Hill or extrapolate shared parameters.
   • latent: compute \(\Delta\); choose displacement \(\lambda(d)\).
4. Apply perturbations → write to adata.layers['perturbed'] (or replicas).
5. Post-process: optional re-normalize/log1p; recompute PCA/UMAP/DEG.
6. Repeat simulations to obtain distributions/CI.
7. Save: adata, plots, DEG tables, dose-response params, adata.uns['perturbation_meta'].

def simulate_virtual_cells(adata, targets, mode="multiplicative", params=None,
                           n_replicates=1, random_state=42, layer_in="normalized",
                           layer_out="perturbed"):
    import numpy as np, datetime
    rng = np.random.default_rng(random_state)
    X = adata.layers[layer_in] if layer_in in adata.layers else adata.X.copy()

    if mode == "multiplicative":
        alpha = float(params.get("alpha", 0.5))
        sigma = float(params.get("sigma", 0.0))
        s = 1.0 - alpha
        if sigma > 0:
            s = rng.lognormal(mean=np.log(max(1e-6, 1.0 - alpha)), sigma=sigma, size=X.shape[1])
        gene_mask = adata.var_names.isin(targets)
        X[gene_mask, :] = (X[gene_mask, :].T * s).T

    elif mode == "dose_response":
        dose = float(params.get("dose", 1.0))
        Emax = float(params.get("Emax", 1.0))
        EC50 = float(params.get("EC50", 1.0))
        n = float(params.get("n", 1.0))
        kappa = float(params.get("kappa", 1.0))
        E = Emax - (Emax - 0.0) / (1 + (EC50 / dose)**n)
        gene_mask = adata.var_names.isin(targets)
        X[gene_mask, :] = X[gene_mask, :] * (1 + kappa * E)

    elif mode == "latent":
        import scanpy as sc
        if "X_pca" not in adata.obsm:
            sc.tl.pca(adata)
        lam = float(params.get("lambda", 1.0))
        tr = (adata.obs["treatment"] == "treated").values
        ct = (adata.obs["treatment"] == "control").values
        Delta = adata.obsm["X_pca"][tr].mean(0) - adata.obsm["X_pca"][ct].mean(0)
        adata.obsm["X_pca_perturbed"] = adata.obsm["X_pca"] + lam * Delta

    adata.layers[layer_out] = X
    adata.obs["is_virtual"] = True
    adata.obs["perturb_id"] = params.get("tag", "sim_001")
    adata.uns["perturbation_meta"] = {
        "alpha": params.get("alpha"),
        "dose": params.get("dose"),
        "mode": mode,
        "random_state": random_state,
        "timestamp": datetime.datetime.now().isoformat()
    }
    return adata

6) Outputs

1. Virtual expression matrix layer: adata.layers['perturbed'] (or perturbed_rep# for multiple runs).
2. Virtual cell annotation: adata.obs['is_virtual'] (True/False), adata.obs['perturb_id'].
3. Low-dimensional embeddings: adata.obsm['X_pca_perturbed'], adata.obsm['X_umap_perturbed'] (if recomputed).
4. Virtual vs Original DEG: CSV/TSV (logFC, p-value, FDR) via sc.tl.rank_genes_groups(..., layer='perturbed').
5. Dose-response params & plots: per-gene Emax/EC50/n and fitted curves (PNG/SVG).
6. Statistical summary: means, variances, sensitivity scores; uncertainty via replicates/CI.
7. Simulation metadata: adata.uns['perturbation_meta'] stores parameters (alpha/dose/mode/random_state/timestamp).

8) Output Result Explanation

• adata.layers['perturbed']: ready for downstream (normalize → log1p → PCA → DEG); note annotations refer to perturbed data.
• DEG table: compares perturbed vs control to verify pathway shifts; can feed into mechanism analysis.
• PCA/UMAP: visual shift of virtual cells highlights population-level drug effects.
• Dose-response curve: estimates EC50/Emax for pharmacological discussion.
• Replicates: average curves with CI reflect uncertainty.

9) Features & Highlights

• Multi-strategy support: multiplicative, Hill fitting, and latent perturbation cover sparse→rich data scenarios.
• Reproducibility: all parameters/seeds stored in adata.uns['perturbation_meta'] for exact replication.
• Fidelity vs flexibility: multiplicative is simple/transparent; latent preserves covariance; Hill is pharmacologically interpretable.
• Downstream integration: seamless with Scanpy pipelines (PCA/UMAP/DEG/GSEA).
• Parallel/batch simulation: supports multi-replicate/multi-dose grid for response surfaces.
• Uncertainty quantification: Monte-Carlo replication with CI for design decisions.

1) Purpose

The Drug Efficacy Analysis Module quantitatively describes pharmacological effects of ATRA at the single-cell/population level, assessing how concentration, duration, or intervention strategies impact genes or cell states. Objectives:

• Fit dose–response curves per gene/pathway/population to estimate EC50/IC50, Emax, and Hill coefficient.
• Compute AUC and max/min response for efficacy ranking and target prioritization.
• Identify sensitive/resistant subpopulations and quantify their proportions (sensitivity distribution).
• Run statistical tests (between-dose significance, trends, population differences) with multiple-testing correction.
• Provide visualizations: dose–response plots, heatmaps, sensitivity distributions, violin plots, and dose-gradient trajectories (PCA/UMAP). Support uncertainty quantification (bootstrap/Bayesian CI) and mixed-effects to account for replicate/batch/individual effects.

2) Module Assumptions

• Dose-dependency: ATRA effects are dose-dependent and follow fitable curves (Hill/Logistic) within observed ranges.
• Reproducibility: Populations are reproducible under identical treatments; variability estimable via replicates.
• Independence & nesting: Single-cell data are nested (cells within samples/batches); mixed-effects handle correlations.
• Distribution: Transformed data meet model assumptions (log1p/CPM/GLM); counts via Poisson/NegBin GLM.
• Measurable effect: Effects are detectable at the transcriptional level with sufficient sample size.

3) Mathematical Models

3.1 Single-Gene Dose–Response (Hill / 4PL)

For gene \(g\) vs. dose \(d\) (or log-dose), model mean expression:

Parameters \(E_{0,g}\) (baseline), \(E_{\max,g}\) (max effect; negative can indicate down-regulation), \(EC_{50,g}\) (potency), \(n_g\) (steepness) are estimated via least squares / MLE.

3.2 GLMM (Generalized Linear Mixed-Effects)

For single-cell counts \(y_{i j k g}\) (sample \(i\), cell \(j\), measurement \(k\)):

with link \(g(\cdot)\) (log for counts; identity/log for continuous), random effects \(u,v\) capturing individual/batch/cell-level variation.

3.3 AUC & Sensitivity Score

Within dose range \([d_{\min}, d_{\max}]\):

Sensitivity fraction over cells with response above threshold \(\tau\):

3.4 Single-Cell Response Distribution & Clustering

For cell \(c\), define response vector \(r_c = \{ x_{c,g_1}(d)-x_{c,g_1}(0),\ldots\}\); cluster \(r_c\) to reveal sensitive/resistant groups; compute \(S(d)\) as above.

4) Data Sources

• Required: Preprocessed AnnData (from ATRA module or virtual cell module) with:
  • Multi-dose data (if available) or control + single dose;
  • adata.obs includes dose/treatment/replicate/batch/cell_type;
  • adata.layers may include perturbed layers or raw counts for GLM.
• Optional: Public dose–response data; PK/PD priors for constraints or Bayesian priors.

5) Implementation

5.1 Data Input

import scanpy as sc
adata = sc.read_h5ad('preprocessed_atra.h5ad')
print(adata.obs.columns)  # Must include 'dose' and 'treatment'

Required arguments: adata (dose/treatment in obs), genes (target list or auto-selected), and analysis params: dose_range, fit_model ('hill'/'logistic'/'glm'/'mixed'), bootstrap_reps, alpha, random_state.

5.2 Processing Flow

1. Select granularity: per gene vs. gene set/pathway.
2. Aggregate per dose (population means) or fit GLMM on single-cell data.
3. Fit model:
   • hill: non-linear least squares per gene across doses;
   • glm: NegBin regression for counts;
   • mixed: mixed-effects for nested structure.
4. Estimate parameters & uncertainty (bootstrap/SE/CIs; Bayesian: posteriors).
5. Compute derived metrics: AUC, Emax, EC50, slope, sensitivity fraction.
6. Multiple-testing correction: Benjamini–Hochberg (FDR).
7. Pattern discovery: cluster dose–response shapes.
8. Visualization: curves, heatmaps, violins, dose-colored PCA/UMAP.
9. Persist outputs: parameter tables, plots, updated adata.uns.

import numpy as np
from scipy.optimize import curve_fit

def hill_fn(d, E0, Emax, EC50, n):
    return E0 + (Emax - E0) / (1 + (EC50 / np.maximum(d, 1e-9))**n)

def fit_hill(doses, expr):
    # doses: array shape [D]; expr: mean expr per dose for one gene [D]
    p0 = [np.median(expr), np.max(expr), np.median(doses), 1.0]
    bounds = ([-np.inf, -np.inf, 1e-9, 0.1], [np.inf, np.inf, np.inf, 5.0])
    popt, pcov = curve_fit(hill_fn, doses, expr, p0=p0, bounds=bounds, maxfev=10000)
    return popt, pcov  # E0, Emax, EC50, n

def auc_from_fit(d_grid, params):
    y = hill_fn(d_grid, *params)
    auc = np.trapz(y, d_grid)
    nAUC = auc / (d_grid.max() - d_grid.min())
    return auc, nAUC

import statsmodels.api as sm
import pandas as pd
import numpy as np

def fit_glm_nb(counts, doses):
    # counts: vector; doses: vector (log-dose recommended)
    X = sm.add_constant(np.log(np.maximum(doses,1e-9)))
    model = sm.GLM(counts, X, family=sm.families.NegativeBinomial())
    res = model.fit()
    return res.params, res.bse, res.pvalues

def bootstrap_ci(doses, expr, B=200, q=(0.025,0.975), rng=None):
    rng = np.random.default_rng() if rng is None else rng
    pars = []
    for _ in range(B):
      idx = rng.choice(len(doses), len(doses), replace=True)
      p,_ = fit_hill(doses[idx], expr[idx])
      pars.append(p)
    pars = np.array(pars)
    lo = np.quantile(pars, q[0], axis=0)
    hi = np.quantile(pars, q[1], axis=0)
    return lo, hi  # element-wise CI for [E0,Emax,EC50,n]

import platform, datetime
def write_pharm_meta(adata, method="hill", seed=42, extra=None):
    meta = dict(method=method, random_seed=seed,
                time=datetime.datetime.now().isoformat(),
                python=platform.python_version())
    if extra: meta.update(extra)
    adata.uns['pharmacology_meta'] = meta

5.3 Outputs

1. Parameter table per gene/pathway: gene,E0,Emax,EC50,n,...,pvalue_dose,FDR.
2. Uncertainty intervals: bootstrap or MLE covariance (e.g., EC50_CI_low/high).
3. AUC/nAUC table for ranking.
4. Sensitivity distribution: fraction-sensitive by cell or cluster.
5. Visualization plots: dose–response curves, population heatmap (genes×doses), violin plots, UMAP/PCA with dose color bars.
6. DEG across doses: pairwise DEG; volcano plots.
7. Diagnostics: R², residuals, failure logs.
8. Meta-information and reproducibility record: fitting method, random seed, fitting time, software version, etc., written to adata.uns['pharmacology_meta'].

6) Output Interpretation

• EC50 / IC50: dose at half-max effect; lower means higher potency; report CIs.
• Emax: maximal effect (direction & magnitude) for endpoint.
• Hill n: curve steepness; \(n>1\) suggests cooperative-like behavior (interpret with caution).
• AUC: overall effect magnitude across dose range; useful for ranking across targets/pathways.
• p-value / FDR: dose effect significance; use for candidate screens.
• Single-cell sensitivity: percentage of sensitive cells reveals heterogeneity.
• Diagnostics: flag large residuals, failed fits, unstable EC50; recommend validation.

7) Module Features & Highlights

• Multi-level modeling: gene, pathway, cluster, single-cell mixed-effects.
• Flexible models: Hill, logistic, GLM (NegBin), mixed-effects, optional Bayesian.
• Uncertainty quantification: bootstrap CIs by default.
• Integrated output: parameter tables, AUC, visualizations, and metadata for reports.
• Virtual-cell integration: analyze simulated data to compare with experiments.
• Automation & interactivity: batch fitting across many genes with per-gene drill-down.

8) Best Practices and Considerations

1. Preprocess first (normalized, batch-corrected, clear dose info).
2. Choose granularity wisely (means per dose vs. GLMM at single-cell scale).
3. Use replicates and bootstrap/MCMC to get robust CIs.
4. Control multiple testing via FDR when screening many genes.
5. Use NegBin GLM (or DESeq2 in R) for raw counts.
6. Inspect residuals and extreme/unstable EC50; flag as unreliable.
7. Retain metadata: method, seed, failed fits, inits → adata.uns['pharmacology_meta'].

1) Purpose

The Mechanism Analysis Module reveals molecular mechanisms and regulatory networks underlying ATRA action using DEG and efficacy results. Scope:

• Identify significantly affected pathways and functions (GO/KEGG/Reactome).
• Quantify pathway activity at single-cell/subpopulation level (ssGSEA, AUCell, score_genes).
• Infer gene regulatory networks (GRNs) to find key TFs and upstream regulators.
• Build a three-tier Gene → Pathway → Drug network to discover mediators and upstream/downstream signaling.
• Provide pathway/network visualizations (bar/dot plots, heatmaps, GRN diagrams, trajectory plots).
• Output prioritized candidates (TFs, pathways, key nodes) for validation.

2) Module Assumptions

• DEGs reflect biological effects: ATRA-induced DEGs carry pathway-level signals.
• Pathway annotations available: KEGG/Reactome/MSigDB/GO provide coverage and interpretability.
• Regulatory relations are inferable: Co-expression and perturbation clues allow partial GRN inference (data/noise dependent).
• Independence of this stage: Depends on correctness of upstream DEG/efficacy results.

3) Mathematical Models

3.1 ORA (Hypergeometric) & GSEA

Given background size \(N\), pathway gene count \(K\), DEG set size \(n\), overlap \(k\), ORA p-value:

GSEA ranks genes (e.g., by logFC/statistic) and computes an enrichment score (ES) from a weighted running-sum; significance is assessed via permutations to obtain FDR.

3.2 Single-Cell Pathway Activity

Per-cell activity \(A_{c,S}\) for set \(S\) (e.g., ssGSEA/AUCell/score_genes). A simple score_genes form (target minus reference mean):

3.3 GRN Inference

Estimate graph \(G=(V,E)\) with edge weight \(w_{ij}\). One view models gene \(i\) as a function of others:

Approaches include correlation/partial correlation, mutual information, and tree ensembles (e.g., GENIE3 / RandomForest); perturbation data aids directionality.

4) Data Sources

• Inputs: DEG table (gene, logFC, pval, FDR), preprocessed AnnData (expression + cell annotations; optional layers['perturbed']), dose–response params (EC50/Emax).
• Annotations: KEGG, Reactome, MSigDB, GO; TF-target (TRANSFAC, JASPAR, TRRUST, DoRothEA); PPIs (STRING, BioGRID).

5) Implementation

import scanpy as sc, pandas as pd
adata = sc.read_h5ad('preprocessed_atra.h5ad')   # contains normalized layer / obs annotations
deg_df = pd.read_csv('deg_results.csv')          # columns: gene, logFC, pval, FDR
background_genes = set(adata.var_names)          # N = len(background_genes)

from scipy.stats import hypergeom
def ora_pvalue(K, N, n, k):
    # pathway size K, background N, DEG size n, overlap k
    rv = hypergeom(N, K, n)
    # P(X >= k) = sf(k-1)
    return float(rv.sf(k-1))

import gseapy as gp
# ranked list: gene -> statistic (e.g., logFC)
rnk = deg_df[['gene','logFC']].sort_values('logFC', ascending=False)
pre_res = gp.prerank(rnk=rnk, gene_sets='kegg.gmt', min_size=10, max_size=500, permutation_num=1000)
gsea_results = pre_res.res2d  # NES, pval, FDR, leadingEdge

import scanpy as sc
target_genes = [...]         # from enrichment (e.g., KEGG_FATTY_ACID_MET)
ref_genes = [...]            # matched control gene set
sc.tl.score_genes(adata, gene_list=target_genes, ctrl_size=None, score_name='score_fa', use_raw=False)
# UMAP coloring:
# sc.pl.umap(adata, color=['score_fa'])

import numpy as np, pandas as pd
from sklearn.ensemble import RandomForestRegressor
def infer_grn(X_df, genes=None, n_trees=500, random_state=42):
    # X_df: genes x cells (or cells x genes, adapt as needed)
    genes = list(X_df.index) if genes is None else list(genes)
    edges = []
    for gi, tgt in enumerate(genes):
        y = X_df.loc[tgt].values
        X = X_df.drop(index=tgt).T.values
        feat_names = X_df.drop(index=tgt).index
        rf = RandomForestRegressor(n_estimators=n_trees, random_state=random_state, n_jobs=-1)
        rf.fit(X, y)
        imp = rf.feature_importances_
        for s, w in zip(feat_names, imp):
            edges.append((s, tgt, float(w)))
    grn_edges = pd.DataFrame(edges, columns=['source','target','weight'])
    return grn_edges.sort_values('weight', ascending=False)

# Map GRN key nodes to enriched pathways, compute centralities, and assemble a 3-layer view
import networkx as nx
def integrate_network(grn_edges, enriched_df, gene2path):
    G = nx.DiGraph()
    for s,t,w in grn_edges.itertuples(index=False):
        G.add_edge(s,t,weight=w,layer='gene-gene')
    # attach pathways
    for gene, paths in gene2path.items():
        for p in paths:
            G.add_edge(gene, p, weight=1.0, layer='gene-pathway')
    # drug layer (ATRA)
    G.add_node('ATRA', layer='drug')
    for p in set().union(*gene2path.values()):
        G.add_edge('ATRA', p, weight=1.0, layer='drug-pathway')
    cen = nx.pagerank(G, alpha=0.85, weight='weight')
    return G, cen

import datetime, platform
def write_mechanism_meta(adata, method='gsea+ora+grn', seed=42, versions=None, extra=None):
    meta = dict(method=method, random_seed=seed,
                time=datetime.datetime.now().isoformat(),
                python=platform.python_version())
    if versions: meta.update(versions)
    if extra: meta.update(extra)
    adata.uns['mechanism_meta'] = meta

6) Output Description (Outputs & Interpretation)

1. Enrichment result tables (CSV):
   • gsea_results.csv: pathway, NES, pval, FDR, leadingEdgeGenes.
   • ora_results.csv: ORA statistics and significance.
2. Single-cell pathway scores in AnnData: adata.obs['score_<pathway>'].
3. GRN edge table (CSV): source,target,weight,method,pvalue?; optional node centralities.
4. Network & pathway figures: dose–response-relevant pathways, pathway dot/bar plots, GRN diagrams, pathway scores on UMAP, gene–pathway heatmaps.
5. Candidate TF list: ranked by TF-target overlap and centrality.
6. Integrated report (PDF/HTML): enrichment plots, pathway projections, representative networks, summaries, and suggested validations.
7. Interpretation:
   • High NES / low FDR → priority pathways.
   • UMAP localization of high pathway scores → sensitive subpopulations.
   • High-centrality GRN nodes (often TFs) → upstream regulators for validation.
   • Agreement between virtual and real perturbations → stronger causal credibility.

Meta-information and reproducibility record: fitting method, random seed, fitting time, software versions, etc., are recorded in adata.uns['mechanism_meta'].

7) Module Features & Highlights

• Multi-scale analysis from genes → networks → pathways, not just lists but structured mechanisms.
• Single-cell pathway scoring exposes subpopulation heterogeneity.
• Perturbation-aware causal clues leverage virtual/real interventions to suggest directionality.
• Interactive-ready outputs for network exploration and pathway filtering.
• Integrated prioritization of TF/target/pathway candidates for downstream CRISPR/drug validation.

1) ATRA Analysis Tool — Contributions and Innovations (Summary)

Overall Positioning. This tool is an end-to-end platform for liver-targeted delivery (ATRA) single-cell RNA-seq analysis, spanning DBTL (Design–Build–Test–Learn): data cleaning, differential analysis, virtual cell generation, efficacy modeling, mechanism elucidation, and target prioritization. It integrates bioinformatics with reproducible in-silico interventions (virtual cells) to support drug delivery and target validation.

• Integrated DBTL support Links data processing → virtual experiments → pharmacodynamic modeling → mechanism analysis into an automated closed loop, enabling iterative optimization of drug design and protocols.
• Virtual Cell pipeline Generates reproducible “virtual treatment” data for multi-dose/strategy scenarios when wet coverage is limited, reducing cost while enabling sensitivity analysis and prioritization.
• Multi-scale efficacy modeling From gene → pathway → population: single-gene Hill curves, population AUC, single-cell GLMM.
• Mechanism-driven prioritization Integrates DEG, GSEA, GRN, and EC50/AUC into candidates that directly guide CRISPR/validation.
• Engineering & reproducibility Write perturbations to AnnData.layers and log adata.uns['perturbation_meta'] for traceability.
• Modular + visual output Web UI + API, interactive networks, one-click report export for collaboration.

2) Source Repositories & Resource Links

pip install scanpy anndata gseapy scvi-tools arboreto

3) Extended Applications and Future Directions

# Dependencies (example)
pip install scanpy anndata gseapy scvi-tools arboreto statsmodels pymc

# BibTeX (example skeletons to complete in paper)
@article{scanpy18, title={Scanpy}, journal={Nat Methods}, year={2018}}
@article{scvitools20, title={scvi-tools}, journal={Nat Biotechnol}, year={2021}}
@article{gseapy, title={GSEApy}, journal={Bioinformatics}, year={2023}}