visualization_dot_plot.py

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt

def run_scenic_analysis(auc_csv, output_prefix):
    print(f"=== Load AUCELL matrix : {auc_csv} ===")
    auc = pd.read_csv(auc_csv, index_col=0)

    adata = sc.AnnData(auc)
    adata.var_names_make_unique()

    print("=== Preprocessing ===")
    sc.pp.log1p(adata)
    sc.pp.scale(adata)

    print("=== PCA / Neighbors / UMAP ===")
    sc.pp.pca(adata, n_comps=50)
    sc.pp.neighbors(adata, n_neighbors=15, use_rep="X_pca")
    sc.tl.umap(adata)

    print("=== Clustering Leiden ===")
    sc.tl.leiden(adata, resolution=0.8)

    print("=== Save results ===")
    adata.write(f"{output_prefix}_UMAP.h5ad")

    # UMAP plot
    sc.pl.umap(adata, color="leiden", save=f"_{output_prefix}_clusters.png")

    print("=== Rank regulons using Wilcoxon ===")
    sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")

    # Dotplot
    sc.pl.rank_genes_groups_dotplot(
        adata,
        groupby="leiden",
        standard_scale="var",
        n_genes=5,
        save=f"_{output_prefix}_dotplot.png"
    )

    # Volcano-style summary plot
    sc.pl.rank_genes_groups(
        adata,
        n_genes=20,
        sharey=False,
        save=f"_{output_prefix}_ranked.png"
    )

    print("\n=== Top regulons per cluster ===")
    markers = sc.get.rank_genes_groups_df(adata, group=None)
    markers.to_csv(f"{output_prefix}_cluster_markers.csv")

    # Print top 10 per cluster
    clusters = sorted(adata.obs["leiden"].unique(), key=int)

    for cl in clusters:
        top10 = sc.get.rank_genes_groups_df(adata, group=str(cl)).head(10)
        print(f"Cluster {cl} top 10 regulons:", list(top10["names"]))

    print(f"\n🎉 DONE : Output saved for {output_prefix}\n")

    return adata

adata_KO = run_scenic_analysis("auc_mtx_KO.csv", "SCENIC_KO")
adata_WT = run_scenic_analysis("auc_mtx_WT.csv", "SCENIC_WT")
adata_HET = run_scenic_analysis("auc_mtx_HET.csv", "SCENIC_HET")