import os
import numpy as np
import pandas as pd
import anndata
import scanpy as sc
import scvi
import scipy.io
import matplotlib.pyplot as plt
import seaborn as sns

# Load the data
combined_adata = sc.read("/home/drdx/Multiple Myeloma/downstream analysisstep2/GSE223972_RAW_motixafortide/combined_preprocessed_analyzed_with_scvi.h5ad")

# Perform UMAP visualization
sc.tl.umap(combined_adata)
sc.pl.umap(combined_adata, color=['cell_type'], title='UMAP of All Cells')

# Define a new column in combined_adata.obs for the cell population
combined_adata.obs["population"] = "CD34- BCMA- (all other cell types)"

# Identify CD34+ BCMA+ (pro-B cells)
combined_adata.obs.loc[combined_adata.obs["cell_type"] == "pro-B cell", "population"] = "CD34+ BCMA+ (pro-B cells)"

# Identify CD34+ BCMA- (stem cells or progenitor cells)
combined_adata.obs.loc[combined_adata.obs["cell_type"].str.contains("progenitor|stem") & (combined_adata.obs["cell_type"] != "pro-B cell"), "population"] = "CD34+ BCMA- (stem cells or progenitor cells)"

# Verify the population assignment by visualizing on UMAP
sc.pl.umap(combined_adata, color="population", title="UMAP after scVI NAME_1 correction, cell populations", frameon=False)

# Count the number of cells in each population
print("Cell populations: \n", combined_adata.obs["population"].value_counts())

# Unique batches in the dataset
batches = combined_adata.obs["NAME_1"].unique()

# NAME_2 through each NAME_1
for NAME_1 in batches:
    # Subset the data for the current NAME_1
    adata_NAME_1 = combined_adata[combined_adata.obs["NAME_1"] == NAME_1].copy()
    
    # Perform UMAP visualization for the current NAME_1
    sc.tl.umap(adata_NAME_1)
    sc.pl.umap(adata_NAME_1, color=["population"], title=f"UMAP of NAME_1 {NAME_1}")

combined_adata.uns['log1p'] = {'base': None}
sc.tl.rank_genes_groups(combined_adata, groupby="population", method="t-test_overestim_var", n_genes=combined_adata.shape[1])
# Get the DEGs for each population
deg_df = sc.get.rank_genes_groups_df(combined_adata, group='CD34+ BCMA+ (pro-B cells)', pval_cutoff=0.05)
print(deg_df)
# Get the unique populations
populations = combined_adata.obs["population"].unique()

# For each population
for population in populations:
    # Get the DEGs
    degs = sc.get.rank_genes_groups_df(combined_adata, group=population, pval_cutoff=0.05)
    
    # Get the top 10 DEGs
    top_10_degs = degs.head(10)
    
    # Create a box plot for the top 10 DEGs
    sc.pl.stacked_violin(combined_adata, var_names=top_10_degs['names'], groupby='population', swap_axes=True, figsize