local baselines script added

scripts/run_benchmark_all.sh

@@ -12,7 +12,7 @@ layer='scgen_pearson'
 metric_ids="[regression_1, regression_2]"
 cell_type_specific=false #for controls
 normalize=false
-method_ids="[pearson_corr, pearson_causal, positive_control]"
+method_ids="[pearson_corr, positive_control]"
 
 param_file="./params/${RUN_ID}.yaml"
 

scripts/sbatch/calculate_scores.sh

@@ -1,5 +1,5 @@
 #!/bin/bash
-#SBATCH --job-name=robustness
+#SBATCH --job-name=scores
 #SBATCH --time=48:00:00
 #SBATCH --output=logs/%j.out
 #SBATCH --error=logs/%j.err
@@ -8,5 +8,4 @@
 #SBATCH --mem=64G 
 #SBATCH --cpus-per-task=20  
 
-# python src/metrics/script_all.py
-python src/robustness_analysis/script_all.py
+python src/metrics/script_all.py

scripts/sbatch/robustness_analysis.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+#SBATCH --job-name=robustness
+#SBATCH --time=48:00:00
+#SBATCH --output=logs/%j.out
+#SBATCH --error=logs/%j.err
+#SBATCH --mail-type=END
+#SBATCH [email protected]
+#SBATCH --mem=64G 
+#SBATCH --cpus-per-task=20  
+
+python src/robustness_analysis/script_all.py

src/control_methods/negative_control/main.py

@@ -0,0 +1,35 @@
+
+import pandas as pd
+import anndata as ad
+import sys
+import numpy as np
+
+def process_links(net, par):
+    net = net[net.source!=net.target]
+    net = net.sample(par['max_n_links'])
+    return net
+def main(par):
+  print('Reading input data')
+  perturbation_data = ad.read_h5ad(par["perturbation_data"])
+  gene_names = perturbation_data.var_names.to_numpy()
+  tf_all = np.loadtxt(par['tf_all'], dtype=str)
+
+  n_tf = 500
+  tfs = tf_all[:n_tf]
+
+  def create_negative_control(gene_names) -> np.ndarray:
+      ratio = [.98, .01, 0.01]
+
+      net = np.random.choice([0, -1, 1], size=((len(gene_names), n_tf)),p=ratio)
+      net = pd.DataFrame(net, index=gene_names, columns=tfs)
+      return net
+  print('Inferring GRN')
+  net = create_negative_control(gene_names)
+
+  pivoted_net = net.reset_index().melt(id_vars='index', var_name='source', value_name='weight')
+
+  pivoted_net = pivoted_net.rename(columns={'index': 'target'})
+  pivoted_net = pivoted_net[pivoted_net['weight'] != 0]
+
+  pivoted_net = process_links(pivoted_net, par)
+  return pivoted_net

src/control_methods/negative_control/script.py

@@ -2,6 +2,7 @@
 import anndata as ad
 import sys
 import numpy as np
+import sys 
 
 ## VIASH START
 par = {
@@ -12,37 +13,11 @@
 }
 ## VIASH END
 print(par)
-
-def process_links(net, par):
-    net = net[net.source!=net.target]
-    net = net.sample(par['max_n_links'])
-    print(net)
-    return net
-
-print('Reading input data')
-perturbation_data = ad.read_h5ad(par["perturbation_data"])
-gene_names = perturbation_data.var_names.to_numpy()
-tf_all = np.loadtxt(par['tf_all'], dtype=str)
-
-n_tf = 500
-tfs = tf_all[:n_tf]
-
-def create_negative_control(gene_names) -> np.ndarray:
-    ratio = [.98, .01, 0.01]
-
-    net = np.random.choice([0, -1, 1], size=((len(gene_names), n_tf)),p=ratio)
-    net = pd.DataFrame(net, index=gene_names, columns=tfs)
-    return net
-print('Inferring GRN')
-net = create_negative_control(gene_names)
-
-pivoted_net = net.reset_index().melt(id_vars='index', var_name='source', value_name='weight')
-
-pivoted_net = pivoted_net.rename(columns={'index': 'target'})
-pivoted_net = pivoted_net[pivoted_net['weight'] != 0]
-
-pivoted_net = process_links(pivoted_net, par)
-
-print('Saving')
-pivoted_net.to_csv(par["prediction"])
+meta = {
+  'resources_dir':'src/control_methods/negative_control'
+}
+sys.path.append(meta['resources_dir'])
+from main import main
+prediction = main(par)
+prediction.to_csv(par["prediction"])
 

src/control_methods/positive_control/script.py

@@ -25,8 +25,6 @@
 print('Create causal corr net')
 par['causal'] = True
 par['multiomics_rna'] = par['perturbation_data']
-par['only_hvgs'] = False
-
 net = create_corr_net(par)
 
 print('Output GRN')

src/control_methods/script_all.py

@@ -0,0 +1,47 @@
+import pandas as pd
+import anndata as ad
+import sys
+import numpy as np
+import os 
+import random
+
+par = {
+  'write_dir': "resources/grn_models/d0_hvgs",
+  "multiomics_rna": "resources/grn-benchmark/multiomics_rna_d0_hvg.h5ad",
+
+  "perturbation_data": "resources/grn-benchmark/perturbation_data.h5ad",
+  "tf_all": "resources/prior/tf_all.csv",
+  "max_n_links": 50000,
+  'layer': 'scgen_pearson',
+  'cell_type_specific': False,
+  'normalize': False,
+  'causal': True
+}
+
+meta = {
+  "resources_dir": 'src/control_methods',
+  "util": 'src/utils'
+}
+sys.path.append(meta["resources_dir"])
+sys.path.append(meta["util"])
+
+os.makedirs(par['write_dir'], exist_ok=True)
+
+from util import create_corr_net
+
+
+#---- run for pearson_corr
+par['prediction'] = f"{par['write_dir']}/pearson_corr.csv"
+par['causal'] = True
+net = create_corr_net(par)
+net.to_csv(par['prediction'])
+#---- run for negative control
+from negative_control.main import main 
+par['prediction'] = f"{par['write_dir']}/negative_control.csv"
+net = main(par)
+net.to_csv(par['prediction'])
+#---- run for positive_control
+par['multiomics_rna'] = par['perturbation_data']
+par['prediction'] = f"{par['write_dir']}/positive_control.csv"
+net = create_corr_net(par)
+net.to_csv(par['prediction'])

src/metrics/script_all.py

@@ -9,18 +9,14 @@
   'reg_type': 'ridge',
   'read_dir': "resources/grn_models/d0_hvgs",
   'write_dir': "resources/results/scores",
-  'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'pearson_causal',  'portia', 'ppcor', 'genie3', 'grnboost2', 'scenic', 'scglue', 'celloracle'],
-  # 'layers': ['lognorm', 'pearson', 'seurat_lognorm', 'seurat_pearson', 'scgen_lognorm', 'scgen_pearson'],
-  'layers': ['seurat_lognorm', 'seurat_pearson'],
-
-  # 'layers': ['scgen_pearson'],
+  'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'portia', 'ppcor', 'genie3', 'grnboost2', 'scenic', 'scglue', 'celloracle'],
+  'layers': ['lognorm', 'pearson', 'seurat_lognorm', 'seurat_pearson', 'scgen_lognorm', 'scgen_pearson'],
 
   "perturbation_data": "resources/grn-benchmark/perturbation_data.h5ad",
   "tf_all": "resources/prior/tf_all.csv",
   "max_n_links": 50000,
   "apply_tf": "true",
   'subsample': -2,
-  'num_workers': 4,
   'verbose': 1,
   'binarize': True,
   'num_workers': 20,

src/robustness_analysis/script_all.py

@@ -11,7 +11,7 @@
   'write_dir': "resources/results/robustness_analysis",
   'degrees': [0, 10, 20, 50, 100],
   'noise_types': ["net", "sign"],
-  'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'pearson_causal',  'portia', 'ppcor', 'genie3', 'grnboost2', 'scenic', 'scglue', 'celloracle'],
+  'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'portia', 'ppcor', 'genie3', 'grnboost2', 'scenic', 'scglue', 'celloracle'],
 
 
   "perturbation_data": "resources/grn-benchmark/perturbation_data.h5ad",