610 Add script to apply SEACells to AnnData objects #658

jashapiro · 2024-07-24T19:36:02Z

Closes #610

This PR adds a python script to apply the SEACells algorithm to AnnData objects.

The main steps and default parameters are derived from https://github.com/dpeerlab/SEACells/blob/main/notebooks/SEACell_computation.ipynb and associated SEACells source code.

One thing to note was that there are few subtleties to the expected formats and associated metadata for things like PCA and UMAP, so it was necessary to recalculate those values, so I wrote a small function to do that. I initally thought I might be able to get away with just renaming the matrices, but it seemed like some of the metadata stored when scanpy calculates the dimension reduction matrices was required downstream, so recalculation seemed the way to go. In my testing, the PCA results were nearly identical, but UMAP were subtly different, as might be expected.

I am currently writing out both the updated AnnData object and the SEACells model object (as a pickle) in order to have access to both for the later analysis.
The model object seems most useful for checking convergence, so it may not really be needed in most cases, but it is relatively small, so easy to keep.

The wrapper script currently runs only on one sample, and exists mostly just for tracking the commands used for testing.

sjspielman

This looks good to me overall! A few of my comments focus on some choices that the script makes, so more docs may be helpful there.
In addition, I do think we want to explain somewhere in the docs the values that you had to recalculate to run SEAcells. The main purpose of that is to be clear how the input data for (at least 1 script in) this module slightly differs from the processed objects.

I'm also think it makes sense to activate the module GHA before developing further, but that should probably be in a separate PR.

sjspielman · 2024-07-24T20:10:59Z

analyses/metacells/README.md


 ## Input files

 Initial analysis is expected use processed `SingleCellExperiment` objects from the ScPCA Portal, in AnnData format.

+For initial testing, files were downloaded with:


Just confirming this code is right - since you wrote testing, I wondered if you meant to include the test flag? If not, then not!

analyses/metacells/run-metacells.sh

analyses/metacells/scripts/run-seacells.py

sjspielman · 2024-07-24T20:12:57Z

analyses/metacells/scripts/run-seacells.py

+"""
+Seacells Analysis Script
+Joshua Shapiro
+2024-07-18


This seems like it will get out of date real fast, for example, now.

analyses/metacells/scripts/run-seacells.py

sjspielman · 2024-07-24T20:15:50Z

analyses/metacells/scripts/run-seacells.py

+    ----------
+    adata : anndata.AnnData
+        An annData object containing the data to run the SEACells algorithm on.
+        Should contain an X_pca field with the PCA coordinates of the cells.


We might want to formally check that X_pca if present. It doesn't seem necessary at this stage, but it might be helpful later, so can't hurt now?

I'm going to modify this to actually call the convert function from within this function.

sjspielman · 2024-07-24T20:18:29Z

analyses/metacells/scripts/run-seacells.py

+    adata : anndata.AnnData
+        An annData object containing the data to run the SEACells algorithm on.
+        Should contain an X_pca field with the PCA coordinates of the cells.
+    cell_ratio : float


Can you add a bit here about where the 75 default comes from?

sjspielman · 2024-07-24T20:35:24Z

analyses/metacells/scripts/run-seacells.py

+    # initialize and fit model
+    model.construct_kernel_matrix()
+    model.initialize_archetypes()
+    model.fit(min_iter=10, max_iter=50)


I think this is the function - any reason to not keep their max_iter=100? Is it super slow?
https://github.com/dpeerlab/SEACells/blob/3462c624ffae0df6d3930490f345f00196c3503e/SEACells/cpu.py#L653-L659

sjspielman · 2024-07-24T20:39:00Z

analyses/metacells/scripts/run-seacells.py

+    parser.add_argument(
+        "--logfile", type=pathlib.Path, help="File path for log outputs"
+    )


I'd mention that this also controls the verbosity of SEAcells

I don't think it is necessary for the argument definition. The point is that if you want logs I will put all the info there. If you don't want logs, the verbosity is just annoying.

sjspielman · 2024-07-24T20:39:34Z

analyses/metacells/README.md


 ## Output files

-TBD
+Currently, the script outputs three files per library: an updated AnnData object, a pickled SEACells model object, and a log file.


Currently it looks like it only outputs a logfile if one specified

Which it is in the run script!

yeah, fair! Since I was mostly reviewing python, "the script" meant the python script to me, but that's not what it means here! 👍

Co-authored-by: Stephanie Spielman <[email protected]>

jashapiro

Thanks for the review. I will be implementing many of your suggestions

jashapiro · 2024-07-25T15:42:01Z

analyses/metacells/README.md


 ## Output files

-TBD
+Currently, the script outputs three files per library: an updated AnnData object, a pickled SEACells model object, and a log file.


Which it is in the run script!

analyses/metacells/run-metacells.sh

jashapiro · 2024-07-25T15:44:18Z

analyses/metacells/scripts/run-seacells.py

+    )
+
+    # recompute principal components and umap (adds additional metadata)
+    sc.tl.pca(adata, n_comps=50, mask_var="highly_variable")


Intent here is to keep it the same as the values already present, but I will make that explicit and check.

jashapiro · 2024-07-25T15:45:21Z

analyses/metacells/scripts/run-seacells.py

+    ----------
+    adata : anndata.AnnData
+        An annData object containing the data to run the SEACells algorithm on.
+        Should contain an X_pca field with the PCA coordinates of the cells.


I'm going to modify this to actually call the convert function from within this function.

jashapiro · 2024-07-25T15:47:32Z

analyses/metacells/scripts/run-seacells.py

+    parser.add_argument(
+        "--logfile", type=pathlib.Path, help="File path for log outputs"
+    )


I don't think it is necessary for the argument definition. The point is that if you want logs I will put all the info there. If you don't want logs, the verbosity is just annoying.

call `convert` from within the main function other responses to review

This reverts commit 06e495f.

jashapiro · 2024-07-25T19:49:15Z

I think this is pretty much ready for another review. I briefly tried activating the action, but I ran into trouble which I expect is related to small samples sizes. I will investigate those as part of a separate PR when I activate the action for real.

sjspielman

LGTM! Left a few small comments but nothing I need to see again. Makes sense to deal with GHA separately, as well.

analyses/metacells/scripts/run-seacells.py

sjspielman · 2024-07-25T20:02:04Z

analyses/metacells/scripts/run-seacells.py

+            n_comps = min(50, adata.n_vars, adata.n_obs)
+        sc.tl.pca(adata, n_comps=n_comps, mask_var="highly_variable")
+
+    # recompute principal components and umap


left over comment?

Co-authored-by: Stephanie Spielman <[email protected]>

jashapiro added 10 commits July 18, 2024 13:41

environment update

f1128e3

start seacells script

07dec92

updates from initial testing

8269a70

add ipynb to large files

b9d4c5a

add full script

dca6ed2

Add model output

409cd2a

add bash script

ebe45a7

Merge branch 'main' into jashapiro/610-SEACells-script

e2900a5

fixes to get things running properly

6165f27

update readme

708af41

jashapiro requested a review from allyhawkins as a code owner July 24, 2024 19:36

jashapiro requested review from sjspielman and removed request for allyhawkins July 24, 2024 19:36

Merge branch 'main' into jashapiro/610-SEACells-script

1f53788

sjspielman reviewed Jul 24, 2024

View reviewed changes

Apply suggestions from code review

792d886

Co-authored-by: Stephanie Spielman <[email protected]>

jashapiro commented Jul 25, 2024

View reviewed changes

jashapiro added 4 commits July 25, 2024 15:43

only recompute PCA as needed

2ca98c7

call `convert` from within the main function other responses to review

activate action

06e495f

Fix some comments

e5487a7

Revert "activate action"

e3156b3

This reverts commit 06e495f.

Merge branch 'main' into jashapiro/610-SEACells-script

038098f

jashapiro requested a review from sjspielman July 25, 2024 19:49

sjspielman approved these changes Jul 25, 2024

View reviewed changes

jashapiro and others added 3 commits July 25, 2024 16:11

Update analyses/metacells/scripts/run-seacells.py

4646655

Co-authored-by: Stephanie Spielman <[email protected]>

remove leftover comment

9f02c79

Merge branch 'main' into jashapiro/610-SEACells-script

5bb7acf

jashapiro merged commit 3ee159b into AlexsLemonade:main Jul 25, 2024
2 checks passed

jashapiro deleted the jashapiro/610-SEACells-script branch July 25, 2024 20:58

jashapiro mentioned this pull request Jul 26, 2024

Activate GHA for metacells and fix edge case failures #669

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

610 Add script to apply SEACells to AnnData objects #658

610 Add script to apply SEACells to AnnData objects #658

jashapiro commented Jul 24, 2024

sjspielman left a comment

sjspielman Jul 24, 2024

sjspielman Jul 24, 2024

sjspielman Jul 24, 2024

jashapiro Jul 25, 2024

sjspielman Jul 24, 2024

sjspielman Jul 24, 2024

sjspielman Jul 24, 2024

jashapiro Jul 25, 2024

sjspielman Jul 24, 2024

jashapiro Jul 25, 2024

sjspielman Jul 25, 2024

jashapiro left a comment

jashapiro Jul 25, 2024

jashapiro Jul 25, 2024

jashapiro Jul 25, 2024

jashapiro Jul 25, 2024

jashapiro commented Jul 25, 2024

sjspielman left a comment

sjspielman Jul 25, 2024

610 Add script to apply SEACells to AnnData objects #658

610 Add script to apply SEACells to AnnData objects #658

Conversation

jashapiro commented Jul 24, 2024

sjspielman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jashapiro left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jashapiro commented Jul 25, 2024

sjspielman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment