Development branch usage

[angeyang123]

Hi,

I'm interested in trying out features (and bugfixes :) ) on the development branch. However, whenever I try to use it, I run into errors and pretty significant changes. For example, the run_scenicplus wrapper has been removed. Are there any examples of how to use the development branch?

Thanks!

[SeppeDeWinter]

Hi @angeyang123

Thank you for your interest.

I'm indeed quite drastically changing things in that branch. The wrapper function will be removed indeed (it's a bit buggy and not super efficient). We are switching to a snakemake pipeline. This pipeline can already be used and should be functional.

In order to use it (after installing SCENC+ from the developmental branch) you should first adapt the config.yaml file so it points to your data, see https://github.com/aertslab/scenicplus/blob/development/Snakemake/config/config.yaml

Once you filled in the file you can run the pipeline by moving into the Snakemake directory and running snakemake like this:

$ ls
Dockerfile  Dockerfile_meme  README.md  Snakemake  build  cytoscape_styles  dist  docs  notebooks  pytest.ini  requirements.txt  resources  setup.py  src  tests

$ cd Snakemake

$ snakemake --cores <number of cores> # e.g. 20 cores
 

This should start the pipeline.

Once you managed to run the pipeline (which should run drastically faster compared the the old wrapper b.t.w.) you will notice that the output format is also changed. We are moving away from the old scplus_obj object and using MuData instead (https://mudata.readthedocs.io/en/latest/). I'm currently in the progress of updating the plotting function (and other downstream analysis functons) to use that datatype instead. Currently the dotplot is already working. When that is done I will make a new tutorial on how to use the new version of SCENIC+ and merge the development branch.

Feel free to already try, feedback is always appreciated.

Best,

Seppe

[angeyang123]

Thank you @SeppeDeWinter for your response!

I did try out the Snakemake function, and ran into this error related to the metaclass:

Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 20
Rules claiming more threads will be scaled down.
Job stats:
job                            count
---------------------------  -------
AUCell_direct                      1
AUCell_extended                    1
all                                1
download_genome_annotations        1
eGRN_direct                        1
eGRN_extended                      1
get_search_space                   1
prepare_GEX_ACC_multiome           1
prepare_menr                       1
region_to_gene                     1
scplus_mudata                      1
tf_to_gene                         1
total                             12

Select jobs to execute...

[Wed Aug 16 13:13:50 2023]
rule prepare_GEX_ACC_multiome:
    input: /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc/scATAC/cistopic_obj.pkl, /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc.adata.h5ad
    output: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/ACC_GEX.h5mu
    jobid: 2
    reason: Missing output files: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/ACC_GEX.h5mu
    resources: tmpdir=/tmp


[Wed Aug 16 13:13:50 2023]
rule download_genome_annotations:
    output: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/genome_annotation.tsv, /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/chromsizes.tsv
    jobid: 9
    reason: Missing output files: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/genome_annotation.tsv,/home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/chromsizes.tsv
    resources: tmpdir=/tmp

Traceback (most recent call last):
  File "/home/jupyter-ayang/.conda/envs/scenicplus/bin/scenicplus", line 33, in <module>
    sys.exit(load_entry_point('scenicplus', 'console_scripts', 'scenicplus')())
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/scenicplus.py", line 631, in main
    args.func(args)
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/scenicplus.py", line 25, in command_prepare_GEX_ACC
    from scenicplus.cli.commands import prepare_GEX_ACC
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/commands.py", line 24, in <module>
    from scenicplus.TF_to_gene import calculate_TFs_to_genes_relationships
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/TF_to_gene.py", line 21, in <module>
    from scenicplus.scenicplus_class import SCENICPLUS
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/scenicplus_class.py", line 20, in <module>
    from scanpy import AnnData
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/__init__.py", line 16, in <module>
    from . import plotting as pl
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/__init__.py", line 1, in <module>
    from ._anndata import (
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/_anndata.py", line 28, in <module>
    from . import _utils
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/_utils.py", line 35, in <module>
    class _AxesSubplot(Axes, axes.SubplotBase, ABC):
TypeError: metaclass conflict: the metaclass of a derived class must be a (non-strict) subclass of the metaclasses of all its bases
Traceback (most recent call last):
  File "/home/jupyter-ayang/.conda/envs/scenicplus/bin/scenicplus", line 33, in <module>
    sys.exit(load_entry_point('scenicplus', 'console_scripts', 'scenicplus')())
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/scenicplus.py", line 631, in main
    args.func(args)
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/scenicplus.py", line 148, in download_command
    from scenicplus.cli.commands import download_gene_annotation_chromsizes
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/cli/commands.py", line 24, in <module>
    from scenicplus.TF_to_gene import calculate_TFs_to_genes_relationships
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/TF_to_gene.py", line 21, in <module>
    from scenicplus.scenicplus_class import SCENICPLUS
  File "/home/jupyter-ayang/notebooks/scenic_plus/scenicplus/src/scenicplus/scenicplus_class.py", line 20, in <module>
    from scanpy import AnnData
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/__init__.py", line 16, in <module>
    from . import plotting as pl
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/__init__.py", line 1, in <module>
    from ._anndata import (
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/_anndata.py", line 28, in <module>
    from . import _utils
  File "/home/jupyter-ayang/.conda/envs/scenicplus/lib/python3.8/site-packages/scanpy/plotting/_utils.py", line 35, in <module>
    class _AxesSubplot(Axes, axes.SubplotBase, ABC):
TypeError: metaclass conflict: the metaclass of a derived class must be a (non-strict) subclass of the metaclasses of all its bases
[Wed Aug 16 13:13:57 2023]
Error in rule prepare_GEX_ACC_multiome:
    jobid: 2
    input: /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc/scATAC/cistopic_obj.pkl, /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc.adata.h5ad
    output: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/ACC_GEX.h5mu
    shell:

            scenicplus prepare_data prepare_GEX_ACC                 --cisTopic_obj_fname /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc/scATAC/cistopic_obj.pkl                 --GEX_anndata_fname /home/jupyter-ayang/notebooks/scenic_plus/pbmc/pbmc.adata.h5ad                 --out_file /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/ACC_GEX.h5mu                 --bc_transform_func "lambda x: f'{x}-10x_pbmc'"

        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

[Wed Aug 16 13:13:57 2023]
Error in rule download_genome_annotations:
    jobid: 9
    output: /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/genome_annotation.tsv, /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/chromsizes.tsv
    shell:

        scenicplus prepare_data download_genome_annotations             --species hsapiens             --biomart_host http://www.ensembl.org             --genome_annotation_out_fname /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/genome_annotation.tsv             --chromsizes_out_fname /home/jupyter-ayang/notebooks/scenic_plus/pbmc_development/output/chromsizes.tsv

        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
Complete log: .snakemake/log/2023-08-16T131350.857463.snakemake.log

Also, I notice that my output directory gets deleted after running into errors with snakemake.

Thank you for your help!

[angeyang123]

It looks like installing matplotlib<3.7 fixed the issue!

[SeppeDeWinter]

Great! Did the entire pipeline run through?

[angeyang123]

Yes the entire pipeline ran successfully! I'm looking forward to using the downstream analysis and plotting functions too once they are available.

[SeppeDeWinter]

@angeyang123

Aha that's great to hear. I'm asking because you are one of the first to run it!
The dotplot function is already available. Plotting umaps/tSNE/.. is implemented in scanpy, so those you can already use as well.

Stay tuned for the next release, this one will include all features of the dev branch.

Best,

Seppe

[biobai]

Hi Seppe,
I am using Snakemake version and confusing about this parameter "region_set_folder" in config.yaml

input_data: region_set_folder: ""

what is region_set here?

[biobai]

The Snakemake pipeline was completed. I got eRegulon_direct.tsv and eRegulons_extended.tsv. But there is no tutorial from here to do downstream analysis. Any guidance you could offer would be greatly appreciated!

[SeppeDeWinter]

Hi @biobai

That's great!

Indeed, there is no tutorial yet ... I'm still working on it.
You can already generate the heatmap/dotplot though, we also plot in the tutorial.

See the function here:

scenicplus/src/scenicplus/plotting/dotplot.py

Line 67 in df11ef4

def heatmap_dotplot(

Let me know if you have troubles using it.

All the best,

Seppe

[Qotov]

Hi Seppe,

I do not have bed files. Is there a way to use menr object for running the snakemake pipeline?

Thanks,
Alex

[SeppeDeWinter]

Hi Alex,

Unfortunately this is not easy to do.
However, you can simply export the region sets used for running the motif enrichment step (and thus generating this menr.pkl file) as bed files.

Let me know if you need assistance with this.

All the best,

Seppe

[Qotov]

Thank you for your response. It worked well (after installing matplotlib<3.7 and pycistarget from the dev branch).

[SeppeDeWinter]

UPDATE: To initialise the snakemake folder run:

scenicplus init_snakemake --out_dir <PATH_TO_OUTDIR>

Best,

Seppe

[SeppeDeWinter]

Hi all

Sneak peek, I've uploaded a notebook on how to use the snakemake pipeline of SCENIC+, check it out here: https://github.com/aertslab/scenicplus/blob/development/notebooks/human_cerebellum_SCENIC%2B.ipynb

I still have to annotate the notebook before I will upload it as the official tutorial. However, you can already use this as a reference.

All the best,

Seppe

[YH-Zheng]

Hi @SeppeDeWinter,

I noticed that you updated a new version of the snakemake pipeline tutorial about ten hours ago, but the cisTopic tutorial in it seems to be inaccessible, and looking at the output in the tutorial I see that the previous topic result has changed from an html to a bed file? I have a very large dataset (currently only able to reduce the amount of data by downsampling), and if the snakemake pipeline can significantly speed up the analysis, I'm very much looking forward to having the full new tutorial and trying out the new pipeline for analysis!

Thanks for your efforts!

[SeppeDeWinter]

Hi @YH-Zheng

The pycisTopic tutorial should now also be accessible.

All the best,

Seppe

[SeppeDeWinter]

Hi all

Again an update, tutorials are now online https://scenicplus.readthedocs.io/en/development/tutorials.html.

All the best,

Seppe

[SeppeDeWinter]

Hi @YH-Zheng

Very nice to heart that it is working.

By changing the default value of the .py script "meta_cell_split" I was able to split the cell types correctly.

Could you share these changes?

With "When will the tutorial on building an eGRN network be available" are you referring to visualising the network edges (for example in Cytoscape)?

I'm planning to add some code to generate the old scenicplus_obj from the new output formats so the functions that are not updated yet can still be used while we are updating the code. Would this be useful for you?

All the best,

Seppe

[YH-Zheng]

Hi @SeppeDeWinter,

Could you share these changes?

Just change the prepare_GEX_ACC function in cli/command.py to add "meta_cell_split=" to process_non_multiome_data. But this is just a direct change to the script, if you need to add arguments directly on the command line, I think you need to add the corresponding arguments to def command_prepare_GEX_ACC(arg) in cli/scenicplus.py, and add the following to the Snakefile and config.yaml at the same time

else:
        mdata = process_non_multiome_data(
            GEX_anndata=GEX_anndata,
            cisTopic_obj=cisTopic_obj,
            key_to_group_by=key_to_group_by,
            use_raw_for_GEX_anndata=use_raw_for_GEX_anndata,
            nr_metacells=nr_metacells,
            nr_cells_per_metacells=nr_cells_per_metacells,
            meta_cell_split=" ")

I'm planning to add some code to generate the old scenicplus_obj from the new output formats so the functions that are not updated yet can still be used while we are updating the code. Would this be useful for you?

I think this is very useful for me, that means I can visualize directly with the stable version of the code. This should be the quickest and cleverest way to do it before the new visualization function is adapted.

How long will this probably take? Maybe I can finish the basic analytical calculations in the meantime.

Thanks a lot,

Yuhui

[YH-Zheng]

Hi Seppe，

I have another question, I ran the scenic+ analysis twice, once using the supplied database and the other using my own constructed database, and got the following two quantities, it seems that using my own constructed database found fewer areas instead?

use custom cisTarget database

use provided cisTarget database

How should I compare these two results to choose the best database?

[SeppeDeWinter]

Hi @angeyang123

How long will this probably take? Maybe I can finish the basic analytical calculations in the meantime.

I will try to get around to it in the next two weeks.

I have another question, I ran the scenic+ analysis twice, once using the supplied database and the other using my own constructed database, and got the following two quantities, it seems that using my own constructed database found fewer areas instead?

I'm sorry but tt's very difficult to say from these numbers what is the best result.

All the best,

Seppe

[SeppeDeWinter]

Hi @angeyang123

I just pushed this function to the development branch: f216792

All the best,

Seppe