add example folder

Author: AnneHartebrodt

examples/1-preprocess-data.ipynb

@@ -0,0 +1,80 @@
+
+import anndata as ad
+import pandas as pd
+import seaborn as sns
+import matplotlib.pyplot as plt
+from itertools import combinations
+from matplotlib.pyplot import rc_context
+import numpy as np
+import sys
+import os.path as op
+
+import anndata
+from netmap.downstream import final_downstream
+
+
+import warnings
+
+from netmap.utils.data_utils import *
+from netmap.utils.tf_utils import *
+from netmap.utils.netmap_config import NetmapConfig
+
+from netmap.model.train_model import create_model_zoo
+from netmap.grn.inferrence import inferrence, inferrence_model_wise
+from netmap.masking.internal import *
+from netmap.masking.external import *
+
+from netmap.downstream.edge_selection import *
+from netmap.downstream.clustering import *
+from netmap.downstream.final_downstream import *
+
+import scipy.sparse as scs
+import torch
+
+
+
+def train_model(adata, output_dir, model_name):
+
+   
+    gene_names = np.array(adata.var.index)
+    data_tensor = adata.X # Log normalized, but not standardized data.
+
+
+    if scs.issparse(data_tensor):
+        data_tensor = torch.tensor(data_tensor.todense(), dtype=torch.float32)
+    else:
+        data_tensor = torch.tensor(data_tensor, dtype=torch.float32)
+
+    model_zoo = create_model_zoo(data_tensor,  n_models=10, n_epochs=10000, model_type='NBAutoencoder', latent_dim= 8, dropout_rate=0.1, hidden_dim = [64] )
+
+    grn_adata = inferrence(model_zoo, data_tensor.cuda(), gene_names, xai_method='GradientShap', background_type = 'zeros', backing_file=op.join(output_dir, f'{model_name}.parquet'), return_in_memory=False)
+    
+    #Save anndata obs to grn obs for reference
+    grn_adata.obs = adata.obs
+    grn_adata.write_h5ad( op.join(output_dir, f'{model_name}_grn.h5ad'))
+    
+    grn_adata.var.to_csv(op.join(output_dir, f'{model_name}_var.tsv'), header = '\t')
+    # save the original obs
+    adata.obs.to_csv(op.join(output_dir, f'{model_name}_obs.tsv'), header = '\t')
+
+
+
+if __name__=='__main__':
+
+    # define your output dir.
+    output_dir = "netmap/case_studies/blood"
+    os.makedirs(output_dir, exist_ok=True)
+
+    # this is the folder and filename
+    model_name = 'bd-rhap-rep1-X'
+    model_output_dir = op.join(output_dir, model_name)
+    os.makedirs(model_output_dir, exist_ok=True)
+
+    ## load your data
+    adata = sc.read_h5ad('netmap/data/blood/reprocessed/bd-rhap-rep1.h5ad')
+    # Use the correct layer!
+    # Here we use the X layer. The data needs to be sc.pp.normalized(target = 10000) and sc.pp.log1() transformed
+    # But not scaled (!sc.pp.scale())
+    train_model(adata, output_dir=model_output_dir, model_name = model_name)
+
+    

examples/2-train-model.py

@@ -348,118 +348,7 @@ def filter_signatures_by_Ucell(selected_edges, adata) -> pd.DataFrame:
     data_ucell = adata.obs.filter(like='_UCell')
     return data_ucell
 
-def select_top_edges(gene_inter_adata, adata, top_per_source=10, col_cluster='leiden_remap', min_reg_size=10, verbose=True, return_copy = False):
-    """
-    Selects top gene targets per source from a clustered gene interaction AnnData.
-
-    Parameters
-    ----------
-    gene_inter_adata : AnnData
-        Gene interaction AnnData with `var` containing 'source' and 'target'.
-    adata : AnnData
-        Expression AnnData for ranking genes.
-    top_per_source : int, default=750
-        Number of top targets to select per source.
-    col_cluster : str, default='spectral'
-        Column in obs defining clusters.grn_adata3.var
-
-    Returns
-    -------
-    gene_inter_adata_filtered : AnnData
-        Filtered AnnData containing top edges.
-    reglon_sizes : list of int
-        Sizes of regulatory regions per source.
-
-    """
-
-    min_edge_support = 0.5
-
-    if verbose: print(f"Initial shape: {gene_inter_adata.shape}")
-
-    # Rank genes per cluster
-    sc.tl.rank_genes_groups(adata, groupby=col_cluster, method='t-test')
-    clusters = list(set(np.unique(gene_inter_adata.obs[col_cluster])).intersection(adata.obs[col_cluster]))
 
-    # Merge ranking data across clusters
-    rank_dfs = []
-    for c in clusters:
-        if verbose: print(f"Ranking cluster: {c}")
-        df = sc.get.rank_genes_groups_df(adata, group=str(c))
-        df = df.sort_values('scores', ascending=False)
-        df[f"rank_{c}"] = np.arange(1, len(df) + 1)
-        rank_dfs.append(df[['names', f'rank_{c}']])
-    df_rank = reduce(lambda l, r: pd.merge(l, r, on='names', how='inner'), rank_dfs)
-
-    keep_edges_dict = {}
-    # Compute differences per cluster
-    
-    for c in clusters:
-        Keep_edges, reglon_sizes = [], []
-        if verbose: print(f"Selecting targets for cluster: {c}")
-        df_rank_c = df_rank.copy()
-        rank_cols = [col for col in df_rank.columns if col != 'names']
-        rank_cols.remove(f"rank_{c}")
-        df_rank_c['avg'] = df_rank_c[rank_cols].mean(axis=1)
-        df_rank_c['diff'] = (df_rank_c[f"rank_{c}"] - df_rank_c['avg']).abs()
-        df_rank_c = df_rank_c.sort_values('diff', ascending=False)
-
-
-        for source in gene_inter_adata.var["source"].unique():
-            if df_rank_c.loc[df_rank_c['names'] == source, 'diff'].shape[0] > 0:
-                tf_rank = df_rank_c.loc[df_rank_c['names'] == source, 'diff'].values[0]
-
-                df_targets = (
-                    gene_inter_adata.var[
-                        (gene_inter_adata.var['source'] == source) &
-                        (gene_inter_adata.var[f'{c}_nonzero'] >= min_edge_support)
-                    ]
-                    .merge(df_rank_c[['names', 'diff']],
-                        left_on='target', right_on='names', how='left')
-                )
-
-                df_targets['rank_distance'] = (df_targets['diff'] - tf_rank).abs()
-                df_targets = df_targets.sort_values('rank_distance').head(top_per_source)
-
-                reglon_sizes.append(len(df_targets))
-                if len(df_targets) >= min_reg_size:
-                    Keep_edges.extend(f"{source}_{t}" for t in df_targets['target'])
-
-            keep_edges_dict[c] = Keep_edges
-    keep_edges_dict = process_cell_edges(keep_edges)
-    return keep_edges_dict
-
-
-
-def process_cell_edges(keep_edges):
-    results = {'unique': {}, 'pairwise': {}}
-    all_cells = list(keep_edges.keys())
-
-    def get_source_summary(edge_set):
-        # Handles (source, target) tuples OR strings with a separator like '->'
-        sources = []
-        for e in edge_set:
-            sources.append(e.split('_')[0])
-        
-        source_dict = dict(Counter(sources))
-        sources = pd.DataFrame({'source' :source_dict.keys(), 'count': source_dict.values()}).sort_values('count', ascending=False)
-        return sources
-
-    # Calculate Uniques
-    for cell in all_cells:
-        others = set().union(*(set(keep_edges[c]) for c in all_cells if c != cell))
-        unique = set(keep_edges[cell]) - others
-
-        df = pd.DataFrame(
-            [e.split('_', 1) for e in unique],
-            columns=['source', 'target']
-        )
-
-        results['unique'][cell] = {
-            'edges': df,
-            'summary': get_source_summary(unique)
-        }
-        
-    return results
 
 
 def compute_signatures_UCell_scores(selected_edges, adata, key='unique') -> pd.DataFrame:

examples/3-downstream.ipynb

@@ -50,6 +50,7 @@ def rank_regulon_groups_dotplot(grn_adata_filtered, adata_regl, original_cluster
     pp.dot_size_df = fractions.loc[:, pp.dot_color_df.columns]
     pp.dot_size_df = pp.dot_size_df/(pp.dot_size_df.max())
     pp.dot_size_df = pp.dot_size_df.fillna(0)
+    
     if return_fig:
         return pp
     else:

src/netmap/downstream/final_downstream.py

@@ -184,7 +184,7 @@ def aggregate_edges(selected_edges, grn_adata, key='unique') -> pd.DataFrame:
         print(ct)
         sign = selected_edges[key][ct]['edges'].groupby('source').apply(lambda x: (x['source'] + "_" + x['target']).tolist()).to_dict()
         for g in sign:
-            regulons[f'{ct}_{g}'] = grn_adata[:, sign[g]].X.sum(axis = 1)
+            regulons[f'{ct}_{g}'] = grn_adata[:, sign[g]].X.sum(axis = 1)/len(sign[g])
     regulons = pd.DataFrame(regulons)
     return regulons