Research_project/visualization.py at main · mathilde733/Research_project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt

########### KO #########################
print("=== Load AUCELL matrix ===")
auc = pd.read_csv("auc_mtx_KO.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("SCENIC_KO_UMAP.h5ad")

# UMAP plot colored by clusters
sc.pl.umap(adata, color="leiden", save="_KO_clusters.png")

# Identify top regulons per cluster
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20, save="_KO_regulons_ranked.png")

print("🎉 DONE : Results saved as SCENIC_KO_UMAP.h5ad + figures")

top_regs = adata.uns["rank_genes_groups"]["names"]
unique_regs = set([r[0] for r in top_regs])  # top 1 par cluster
sc.pl.umap(adata, color=list(unique_regs), ncols=4, save="_KO_top1_per_cluster.png")

markers = sc.get.rank_genes_groups_df(adata, group=None)
markers.to_csv("SCENIC_KO_cluster_markers.csv")


########### WT #########################
print("=== Load AUCELL matrix ===")
auc = pd.read_csv("auc_mtx_WT.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("SCENIC_WT_UMAP.h5ad")

# UMAP plot colored by clusters
sc.pl.umap(adata, color="leiden", save="_WT_clusters.png")

# Identify top regulons per cluster
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20, save="_WT_regulons_ranked.png")

print("🎉 DONE : Results saved as SCENIC_WT_UMAP.h5ad + figures")

top_regs = adata.uns["rank_genes_groups"]["names"]
unique_regs = set([r[0] for r in top_regs])  # top 1 par cluster
sc.pl.umap(adata, color=list(unique_regs), ncols=4, save="_WT_top1_per_cluster.png")

markers = sc.get.rank_genes_groups_df(adata, group=None)
markers.to_csv("SCENIC_WT_cluster_markers.csv")