diff --git a/dockerfiles/scenicplus/Dockerfile b/dockerfiles/scenicplus/Dockerfile
index 797d35362..de3179227 100644
--- a/dockerfiles/scenicplus/Dockerfile
+++ b/dockerfiles/scenicplus/Dockerfile
@@ -1,5 +1,12 @@
 # Use the base image
-FROM janursa/scenicplus:19-08-2024
+FROM python:3.11.8
+
+# Install scenicplus
+RUN git clone https://github.com/aertslab/scenicplus
+WORKDIR scenicplus/
+RUN git checkout fa55dae55203951c0c42e359b65a28010cd544f6
+RUN pip install .
+WORKDIR ..
 
 # Install HTSlib and tabix
 RUN git clone https://github.com/samtools/htslib.git
@@ -15,10 +22,3 @@ WORKDIR ..
 RUN apt update
 RUN apt install -y openjdk-17-jdk
 ENV JAVA_HOME /usr/lib/jvm/java-17-openjdk-amd64/
-
-# Re-install scenicplus
-RUN git clone https://github.com/aertslab/scenicplus
-WORKDIR scenicplus/
-RUN git checkout fa55dae55203951c0c42e359b65a28010cd544f6
-RUN pip install .
-WORKDIR ..
diff --git a/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather.aF53FD8D b/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather.aF53FD8D
deleted file mode 100644
index 95dd49838..000000000
Binary files a/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather.aF53FD8D and /dev/null differ
diff --git a/hg38_screen_v10_clust.regions_vs_motifs.scores.feather.e391A0Ff b/hg38_screen_v10_clust.regions_vs_motifs.scores.feather.e391A0Ff
deleted file mode 100644
index faa1792bd..000000000
Binary files a/hg38_screen_v10_clust.regions_vs_motifs.scores.feather.e391A0Ff and /dev/null differ
diff --git a/src/methods/multi_omics/scenicplus/config.vsh.yaml b/src/methods/multi_omics/scenicplus/config.vsh.yaml
index 6f842f150..f4fbcdc36 100644
--- a/src/methods/multi_omics/scenicplus/config.vsh.yaml
+++ b/src/methods/multi_omics/scenicplus/config.vsh.yaml
@@ -27,22 +27,16 @@ functionality:
       required: true
       direction: output
       description: "Cell-topics prob scores"
-    - name: --rankings_db
-      type: file
-      direction: input
-    - name: --scores_db
-      type: file
-      direction: input
-      
-
-      
   resources:
     - type: python_script
       path: script.py
 
 platforms:
   - type: docker
-    image: apassemi/scenicplus:latest
+    image: apassemi/scenicplus:1.0.0
+    setup:
+      - type: python
+        packages: [ flatbuffers ]
   - type: native
   - type: nextflow
     directives:
diff --git a/src/methods/multi_omics/scenicplus/script.py b/src/methods/multi_omics/scenicplus/script.py
index 019ea2abc..6d5160ba6 100644
--- a/src/methods/multi_omics/scenicplus/script.py
+++ b/src/methods/multi_omics/scenicplus/script.py
@@ -8,6 +8,7 @@
 import shutil
 import requests
 import traceback
+import zipfile
 import subprocess
 import gc
 import gzip
@@ -19,9 +20,8 @@
 import pandas as pd
 import anndata
 import pyranges as pr
-from pycistarget.utils import region_names_to_coordinates
-from scenicplus.wrappers.run_pycistarget import run_pycistarget
 import polars
+import mudata
 import scrublet as scr
 from sklearn.manifold import TSNE
 import pycisTopic.loom
@@ -37,6 +37,7 @@
 from pycisTopic.gene_activity import get_gene_activity
 from pycisTopic.loom import export_region_accessibility_to_loom, export_gene_activity_to_loom
 from pycisTopic.clust_vis import find_clusters, run_umap, run_tsne, plot_metadata, plot_topic, cell_topic_heatmap
+from scenicplus.wrappers.run_pycistarget import run_pycistarget
 
 
 ## VIASH START
@@ -62,13 +63,20 @@
 except NameError:
     pass
 
-par['temp_dir'] = os.path.join(os.path.dirname(par['prediction']), 'scenicplus')
+par['temp_dir'] = os.path.join(os.path.dirname(par['prediction']), 'scenicplus') # TODO
 
 out_dir = par['temp_dir']
 atac_dir = os.path.join(out_dir, 'atac')
 os.makedirs(atac_dir, exist_ok=True)
 
-par['cistopic_object'] = f"{par['temp_dir']}/cistopic_object.pkl"
+par['cistopic_object'] = f'{par["temp_dir"]}/cistopic_object.pkl'
+
+
+###############################################################################################
+###############################################################################################
+### Running pycistopic ########################################################################
+###############################################################################################
+###############################################################################################
 
 # Get list of samples (e.g., donors)
 print('Collect list of samples')
@@ -115,7 +123,7 @@
         sorted_filepath = filepath + '.sorted.tsv'
         os.system(f'sort -k1,1 -k2,2n {filepath} > {sorted_filepath}')
 
-        # Compression
+        # Decompression
         subprocess.run(['bgzip', sorted_filepath, '-o', compressed_filepath])
 
     fragments_dict[donor_id] = compressed_filepath
@@ -222,6 +230,16 @@
     chain=False
 )
 
+# Download Mallet
+MALLET_PATH = os.path.join(out_dir, 'Mallet-202108', 'bin', 'mallet')
+if not os.path.exists(MALLET_PATH):
+    url = 'https://github.com/mimno/Mallet/releases/download/v202108/Mallet-202108-bin.tar.gz'
+    response = requests.get(url)
+    with open(os.path.join(out_dir, 'Mallet-202108-bin.tar.gz'), 'wb') as f:
+        f.write(response.content)
+    with tarfile.open(os.path.join(out_dir, 'Mallet-202108-bin.tar.gz'), 'r:gz') as f:
+        f.extractall(path=out_dir)
+
 # Download TSS annotations
 os.makedirs(os.path.join(out_dir, 'qc'), exist_ok=True)
 if not os.path.exists(os.path.join(out_dir, 'qc', 'tss.bed')):
@@ -234,7 +252,7 @@
     ])
 
 # Create cistopic objects
-if not os.path.exists(os.path.join(out_dir, 'cistopic_obj.pkl')):
+if not os.path.exists(par['cistopic_object']):
     if par['qc']:  # Whether to perform quality control
         # Compute QC metrics
         print('Perform QC')
@@ -327,63 +345,47 @@
     else:
         cistopic_obj = merge(cistopic_obj_list, is_acc=1, copy=False, split_pattern='-')
     
+    # LDA-based topic modeling
+    print('Run LDA models', flush=True)
+    n_topics = [2, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
+    if os.path.exists(MALLET_PATH):
+        models = run_cgs_models_mallet(
+            cistopic_obj,
+            n_topics=n_topics,
+            n_cpu=par['num_workers'],
+            n_iter=500,
+            random_state=555,
+            alpha=50,
+            alpha_by_topic=True,
+            eta=0.1,
+            eta_by_topic=False,
+            mallet_path=MALLET_PATH
+        )
+    else:
+        print('Could not find Mallet. Running the sequential version of LDA instead.')
+        models = run_cgs_models(
+            cistopic_obj,
+            n_topics=n_topics,
+            n_cpu=12,
+            n_iter=500,
+            random_state=555,
+            alpha=50,
+            alpha_by_topic=True,
+            eta=0.1,
+            eta_by_topic=False
+        )
+
+    # Model selection
+    model = evaluate_models(models, select_model=40, return_model=True)
+    cistopic_obj.add_LDA_model(model)
+
     # Save cistopic objects
-    with open(os.path.join(out_dir, 'cistopic_obj.pkl'), 'wb') as f:
+    with open(par['cistopic_object'], 'wb') as f:
         pickle.dump(cistopic_obj, f)
 else:
     # Load cistopic objects
-    with open(os.path.join(out_dir, 'cistopic_obj.pkl'), 'rb') as f:
+    with open(par['cistopic_object'], 'rb') as f:
         cistopic_obj = pickle.load(f)
-
-# Download Mallet
-MALLET_PATH = os.path.join(out_dir, 'Mallet-202108', 'bin', 'mallet')
-if not os.path.exists(MALLET_PATH):
-    url = 'https://github.com/mimno/Mallet/releases/download/v202108/Mallet-202108-bin.tar.gz'
-    response = requests.get(url)
-    with open(os.path.join(out_dir, 'Mallet-202108-bin.tar.gz'), 'wb') as f:
-        f.write(response.content)
-    with tarfile.open(os.path.join(out_dir, 'Mallet-202108-bin.tar.gz'), 'r:gz') as f:
-        f.extractall(path=out_dir)
-
-# LDA-based topic modeling
-print('Run LDA models')
-
-# Topic modeling
-n_topics = [2, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
-if os.path.exists(MALLET_PATH):
-    models = run_cgs_models_mallet(
-        cistopic_obj,
-        n_topics=n_topics,
-        n_cpu=par['num_workers'],
-        n_iter=500,
-        random_state=555,
-        alpha=50,
-        alpha_by_topic=True,
-        eta=0.1,
-        eta_by_topic=False,
-        mallet_path=MALLET_PATH
-    )
-else:
-    print('Could not find Mallet. Running the sequential version of LDA instead.')
-    models = run_cgs_models(
-        cistopic_obj,
-        n_topics=n_topics,
-        n_cpu=12,
-        n_iter=500,
-        random_state=555,
-        alpha=50,
-        alpha_by_topic=True,
-        eta=0.1,
-        eta_by_topic=False
-    )
-
-# Model selection
-model = evaluate_models(models, select_model=40, return_model=True)
-cistopic_obj.add_LDA_model(model)
-
-with open(par['cistopic_object'], 'wb') as f:
-    pickle.dump(cistopic_obj, f)
-
     
 cell_topic = cistopic_obj.selected_model.cell_topic.T
 cell_topic.index = cell_topic.index.str.split('-').str[0]
@@ -450,6 +452,12 @@
     n_cpu=5,
     split_pattern='-'
 )
+
+# Remove empty markers
+for DAR in list(markers_dict.keys()):
+    if len(markers_dict[DAR]) == 0:
+        del markers_dict[DAR]
+
 # Save topics
 folder = os.path.join(out_dir, 'region_sets', 'Topics_otsu')
 os.makedirs(folder, exist_ok=True)
@@ -529,7 +537,7 @@
 )
 
 # Infer differentially accessible genes
-DAG_markers_dict= find_diff_features(
+DAG_markers_dict = find_diff_features(
     cistopic_obj,
     gene_act,
     variable='cell_type',
@@ -569,6 +577,92 @@
     split_pattern='-'
 )
 
+
+###############################################################################################
+###############################################################################################
+### Creating custom cistarget database ########################################################
+###############################################################################################
+###############################################################################################
+
+RANKINGS_DB_PATH = os.path.join(out_dir, 'cistarget-db', 'db', 'db.regions_vs_motifs.rankings.feather')
+SCORES_DB_PATH = os.path.join(out_dir, 'cistarget-db', 'db', 'db.regions_vs_motifs.scores.feather')
+
+if not (os.path.exists(RANKINGS_DB_PATH) and os.path.exists(SCORES_DB_PATH)):
+
+    # Download create_cisTarget_databases
+    os.makedirs(os.path.join(out_dir, 'cistarget-db'), exist_ok=True)
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'create_cisTarget_databases')):
+        with contextlib.chdir(os.path.join(out_dir, 'cistarget-db')):
+            subprocess.run(['git', 'clone', 'https://github.com/aertslab/create_cisTarget_databases'])
+
+    # Download cluster-buster
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'cbust')):
+        urlretrieve('https://resources.aertslab.org/cistarget/programs/cbust', os.path.join(out_dir, 'cistarget-db', 'cbust'))
+    subprocess.run(['chmod', 'a+x', os.path.join(out_dir, 'cistarget-db', 'cbust')])
+
+    # Download motif collection
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'v10nr_clust_public')):
+        urlretrieve(
+            'https://resources.aertslab.org/cistarget/motif_collections/v10nr_clust_public/v10nr_clust_public.zip',
+            os.path.join(out_dir, 'cistarget-db', 'v10nr_clust_public.zip')
+        )
+        with zipfile.ZipFile(os.path.join(out_dir, 'cistarget-db', 'v10nr_clust_public.zip'), 'r') as zip_ref:
+            zip_ref.extractall(os.path.join(out_dir))
+
+    # Download chromosome sizes
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'hg38.chrom.sizes')):
+        urlretrieve(
+            'https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes',
+            os.path.join(out_dir, 'cistarget-db', 'hg38.chrom.sizes')
+        )
+
+    # Download reference genome
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'hg38.fa')):
+        print('Downloading reference genome', flush=True)
+        if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'hg38.fa.gz')):
+            urlretrieve(
+                'https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz',
+                os.path.join(out_dir, 'cistarget-db', 'hg38.fa.gz')
+            )
+        with gzip.open(os.path.join(out_dir, 'cistarget-db', 'hg38.fa.gz'), 'rb') as f_in:
+            with open(os.path.join(out_dir, 'cistarget-db', 'hg38.fa'), 'wb') as f_out:
+                shutil.copyfileobj(f_in, f_out)
+
+    # Prepare fasta from consensus regions
+    if not os.path.exists(os.path.join(out_dir, 'cistarget-db', 'hg38.with_1kb_bg_padding.fa')):
+        subprocess.run([
+            os.path.join(out_dir, 'cistarget-db', 'create_cisTarget_databases', 'create_fasta_with_padded_bg_from_bed.sh'),
+            os.path.join(out_dir, 'cistarget-db', 'hg38.fa'),
+            os.path.join(out_dir, 'cistarget-db', 'hg38.chrom.sizes'),
+            os.path.join(out_dir, 'consensus_peak_calling', 'consensus_regions.bed'),
+            os.path.join(out_dir, 'cistarget-db', 'hg38.with_1kb_bg_padding.fa'),
+            '1000',
+            'yes'
+        ])
+
+    # Create cistarget databases
+    with open(os.path.join(out_dir, 'cistarget-db', 'motifs.txt'), 'w') as f:
+        for filename in os.listdir(os.path.join(out_dir, 'cistarget-db', 'v10nr_clust_public', 'singletons')):
+            f.write(f'{filename}\n')
+    with contextlib.chdir(os.path.join(out_dir, 'cistarget-db')):
+        subprocess.run([
+            'python',
+            os.path.join(out_dir, 'cistarget-db', 'create_cisTarget_databases', 'create_cistarget_motif_databases.py'),
+            '-f', os.path.join(out_dir, 'cistarget-db', 'hg38.with_1kb_bg_padding.fa'),
+            '-M', os.path.join(out_dir, 'cistarget-db', 'v10nr_clust_public', 'singletons'),
+            '-m', os.path.join(out_dir, 'cistarget-db', 'motifs.txt'),
+            '-c', os.path.join(out_dir, 'cistarget-db', 'cbust'),
+            '-o', 'db',
+            '--bgpadding', '1000',
+            '-t', str(par['num_workers'])
+        ])
+
+###############################################################################################
+###############################################################################################
+### Running scenicplus pipeline ###############################################################
+###############################################################################################
+###############################################################################################
+
 os.makedirs(os.path.join(out_dir, 'scRNA'), exist_ok=True)
 
 # Download databases
@@ -579,49 +673,10 @@ def download(url: str, filepath: str) -> None:
         return
     print(f'Download {url}...')
     urlretrieve(url, filepath)
-def download_and_checksum(url: str, filepath: str, digest: str) -> None:
-    download(url, filepath)
-    #with open(filepath, 'rb') as f:
-    #    file_hash = hashlib.file_digest(f, 'sha1')
-    #if file_hash.hexdigest() != digest:
-    #    os.remove(filepath)
-    #print(file_hash.hexdigest(), digest)
-    #assert file_hash.hexdigest() == digest
-def download_checksum(url: str, filepath: str) -> str:
-    if not os.path.exists(filepath):
-        response = requests.get(url)
-        with open(filepath, 'w') as f:
-            f.write(response.text)
-    with open(filepath, 'r') as f:
-        s = f.read()
-    return s.split()[0]
-# Define rankings, score and motif annotation databases
-if par['rankings_db']:
-    print("Read tf motif db locally")
-    rankings_db = par['rankings_db']
-    scores_db =  par['scores_db']
-else:
-    print("Downloading tf motif db")
-    url = 'https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/screen/mc_v10_clust/region_based/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather.sha1sum.txt'
-    digest = download_checksum(url, os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.rankings.feather.sha1sum.txt'))
-
-    url = 'https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/screen/mc_v10_clust/region_based/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather'
-    download_and_checksum(url, os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.rankings.feather'), digest)
-
-    url = 'https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/screen/mc_v10_clust/region_based/hg38_screen_v10_clust.regions_vs_motifs.scores.feather.sha1sum.txt'
-    digest = download_checksum(url, os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.scores.feather.sha1sum.txt'))
-
-    url = 'https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/screen/mc_v10_clust/region_based/hg38_screen_v10_clust.regions_vs_motifs.scores.feather'
-    download_and_checksum(url, os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.scores.feather'), digest)
-
-    rankings_db = os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.rankings.feather')
-    scores_db =  os.path.join(DB_PATH, 'hg38_screen_v10_clust.regions_vs_motifs.scores.feather')
-
 url = 'https://resources.aertslab.org/cistarget/motif_collections/v10nr_clust_public/snapshots/motifs-v10-nr.hgnc-m0.00001-o0.0.tbl'
 download(url, os.path.join(DB_PATH, 'motifs-v10-nr.hgnc-m0.00001-o0.0.tbl'))
 motif_annotation = os.path.join(DB_PATH, 'motifs-v10-nr.hgnc-m0.00001-o0.0.tbl')
 
-
 print("Preprocess RNA-seq", flush=True)
 # Load scRNA-seq data
 adata_rna = anndata.read_h5ad(par['multiomics_rna'])
@@ -641,47 +696,31 @@ def download_checksum(url: str, filepath: str) -> str:
 adata_rna.raw = adata_rna
 sc.pp.normalize_total(adata_rna, target_sum=1e4)
 sc.pp.log1p(adata_rna)
-adata_rna.write_h5ad(os.path.join(out_dir, 'rna.h5ad'))
 
-# Load candidate enhancer regions
-with open(os.path.join(out_dir, f'candidate_enhancers/region_bin_topics_otsu.pkl'), 'rb') as f:
-    region_bin_topics_otsu = pickle.load(f)
-with open(os.path.join(out_dir, f'candidate_enhancers/region_bin_topics_top3k.pkl'), 'rb') as f:
-    region_bin_topics_top3k = pickle.load(f)
-with open(os.path.join(out_dir, f'candidate_enhancers/markers_dict.pkl'), 'rb') as f:
-    markers_dict = pickle.load(f)
-
-
-# Convert to dictionary of pyrange objects
-region_sets = {}
-region_sets['topics_otsu'] = {}
-region_sets['topics_top_3'] = {}
-region_sets['DARs'] = {}
-for topic in region_bin_topics_otsu.keys():
-    regions = region_bin_topics_otsu[topic].index[region_bin_topics_otsu[topic].index.str.startswith('chr')] #only keep regions on known chromosomes
-    region_sets['topics_otsu'][topic] = pr.PyRanges(region_names_to_coordinates(regions))
-for topic in region_bin_topics_top3k.keys():
-    regions = region_bin_topics_top3k[topic].index[region_bin_topics_top3k[topic].index.str.startswith('chr')] #only keep regions on known chromosomes
-    region_sets['topics_top_3'][topic] = pr.PyRanges(region_names_to_coordinates(regions))
-for DAR in markers_dict.keys():
-    regions = markers_dict[DAR].index[markers_dict[DAR].index.str.startswith('chr')]
-    region_sets['DARs'][DAR] = pr.PyRanges(region_names_to_coordinates(regions))
+# Change barcodes to match the barcodes in the scATAC-seq data
+bar_codes = [f'{obs_name.replace("-", "")}-{donor_id}' for obs_name, donor_id in zip(adata_rna.obs_names, adata_rna.obs.donor_id)]
+adata_rna.obs_names = bar_codes
+
+# Save scRNA-seq data
+adata_rna.write_h5ad(os.path.join(out_dir, 'rna.h5ad'))
 
 # Init scenicplus pipeline
 os.makedirs(os.path.join(out_dir, 'scplus_pipeline'), exist_ok=True)
 os.makedirs(os.path.join(out_dir, 'scplus_pipeline', 'temp'), exist_ok=True)
 subprocess.run(['scenicplus', 'init_snakemake', '--out_dir', os.path.join(out_dir, 'scplus_pipeline')])
-print("snake make initialized", flush=True)
+print('snake make initialized', flush=True)
+
 # Load pipeline settings
 with open(os.path.join(out_dir, 'scplus_pipeline', 'Snakemake', 'config', 'config.yaml'), 'r') as f:
     settings = yaml.safe_load(f)
 print('output_data:', settings['output_data'], flush=True)
+
 # Update settings
 settings['input_data']['cisTopic_obj_fname'] = par['cistopic_object']
 settings['input_data']['GEX_anndata_fname'] = os.path.join(out_dir, 'rna.h5ad')
 settings['input_data']['region_set_folder'] = os.path.join(out_dir, 'region_sets')
-settings['input_data']['ctx_db_fname'] = rankings_db
-settings['input_data']['dem_db_fname'] = scores_db
+settings['input_data']['ctx_db_fname'] = RANKINGS_DB_PATH
+settings['input_data']['dem_db_fname'] = SCORES_DB_PATH
 settings['input_data']['path_to_motif_annotations'] = motif_annotation
 settings['params_general']['temp_dir'] = os.path.join(out_dir, 'scplus_pipeline', 'temp')
 settings['params_general']['n_cpu'] = par['num_workers']
@@ -701,34 +740,32 @@ def download_checksum(url: str, filepath: str) -> str:
 settings['params_motif_enrichment']['dem_max_bg_regions'] = 3000
 settings['params_motif_enrichment']['ctx_auc_threshold'] = 0.005
 settings['params_motif_enrichment']['ctx_nes_threshold'] = 3.0
+settings['params_data_preparation']['bc_transform_func'] = '\"lambda x: f\'{x}\'\"'  # '"lambda x: f''{x}''"'
 settings['params_data_preparation']['is_multiome'] = True
-settings['params_data_preparation']['is_multiome'] = 'key_to_group_by'
+settings['params_data_preparation']['key_to_group_by'] = ''
 settings['params_data_preparation']['nr_cells_per_metacells'] = 5
 settings['params_data_preparation']['species'] = 'hsapiens'
 settings['output_data']['scplus_mdata'] = par['scplus_mdata']
 
-print('output_data:', settings['output_data'], flush=True)
 # Save pipeline settings
+print('output_data:', settings['output_data'], flush=True)
 with open(os.path.join(out_dir, 'scplus_pipeline', 'Snakemake', 'config', 'config.yaml'), 'w') as f:
     yaml.dump(settings, f)
 
-# TODO: from this line onward, the code is untested (could not run it locally due to excessive memory requirements)
-print('run snakemake ', flush=True)
 # Run pipeline
-with contextlib.chdir(os.path.join(out_dir, 'scplus_pipeline', 'Snakemake')):
-    subprocess.run([
-        'snakemake',
-        '--cores', str(par['num_workers'])
-    ])
-import mudata
+print('run snakemake ', flush=True)
+if not os.path.exists(par['scplus_mdata']):
+    with contextlib.chdir(os.path.join(out_dir, 'scplus_pipeline', 'Snakemake')):
+        subprocess.run([
+            'snakemake',
+            '--cores', str(par['num_workers']),
+            #'--unlock'
+        ])
 scplus_mdata = mudata.read(par['scplus_mdata'])
-prediction = scplus_mdata.uns["direct_e_regulon_metadata"]
 
+# Save results
+prediction = scplus_mdata.uns['direct_e_regulon_metadata']
+prediction.insert(0, 'source', prediction['TF'])
+prediction.insert(1, 'target', prediction['Gene'])
+prediction.insert(2, 'weight', prediction['importance_x_abs_rho'])
 prediction.to_csv(par['prediction'])
-
-# # Make sure the file is properly formatted, and re-format it if needed
-# filepath = os.path.join(out_dir, 'tf_to_gene_adj.tsv')
-# shutil.copyfile(filepath, par['prediction'])
-
-
-
diff --git a/src/methods/multi_omics/scenicplus/test.sh b/src/methods/multi_omics/scenicplus/test.sh
index 312b6bb92..70b54f1f8 100644
--- a/src/methods/multi_omics/scenicplus/test.sh
+++ b/src/methods/multi_omics/scenicplus/test.sh
@@ -1,8 +1,2 @@
-viash run src/methods/multi_omics/scenicplus/config.vsh.yaml -- --multiomics_atac resources_test/grn-benchmark/multiomics_atac.h5ad \
-    --multiomics_rna resources_test/grn-benchmark/multiomics_rna.h5ad \
-    --temp_dir output/scenicplus \
-    --prediction output/scenicplus/prediction.csv \
-    --cell_topic output/scenicplus/cell_topic.csv \
-    --scplus_mdata output/scenicplus/scplus_mdata.h5mu \
-    --scores_db resources_test/supplementary/scenicplus_data/hg38_screen_v10_clust.regions_vs_motifs.scores.feather \
-    --rankings_db resources_test/supplementary/scenicplus_data/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather
+#viash run src/methods/multi_omics/scenicplus/config.vsh.yaml -p docker -- ---setup build
+viash run src/methods/multi_omics/scenicplus/config.vsh.yaml -- --multiomics_atac resources_test/grn-benchmark/multiomics_atac.h5ad --multiomics_rna resources_test/grn-benchmark/multiomics_rna.h5ad --temp_dir output/scenicplus --prediction output/scenicplus/prediction.csv --cell_topic output/scenicplus/cell_topic.csv --scplus_mdata output/scenicplus/scplus_mdata.h5mu
\ No newline at end of file