Initial osteocar reference comparison by sjspielman · Pull Request #1426 · AlexsLemonade/OpenScPCA-analysis

sjspielman · 2026-08-05T13:16:11Z

Towards #1424

Purpose/implementation Section

The goal of this PR is to begin comparing references to each other. Originally I had one pretty large notebook, but I am breaking it up into multiple notebooks for easier review so there's more to come.
As such, this PR contains the following:

An initial script to convert references to SCE objects. This script will evolve when we implement Add script to prepare references for use #1413, since it will also need to convert to AnnData. But, this is good enough for this initial exploration now
- This is now run in the bash script too
A new exploratory notebook which just runs the match references function on aggregated (they are huge) references to get an initial sense of how similar their profiles and cell types are. I do this twice: one for patient prlmary/met references and one for xenograft primary/met. I conclude that the references indeed seem to have different profiles so we should explore them further with some very cursory cell typing in the next notebook.
- The other PR I am going to file here will have a notebook that cell types 5 patient samples with different properties to see how different reference strategies come out. This will of course be a very cursory exploration, but it should provide guidance.
I also have a renv environment building up, which does not contain more packages than are used in this PR (eg singler), since I was developing a larger notebook at the same time as this. We expect everything that is in renv to be needed anyways so this seemed fine.

Here is this PR's rendered notebook:
02-compare-references.nb.html

Author checklists

Analysis module and review

This analysis module uses the analysis template and has the expected directory structure.
The analysis module README.md has been updated to reflect code changes in this pull request.
The analytical code is documented and contains comments.
Any results and/or plots this code produces have been added to your S3 bucket for review.

Reproducibility checklist

Code in this pull request has been added to the GitHub Action workflow that runs this module.
The dependencies required to run the code in this pull request have been added to the analysis module Dockerfile.
If applicable, the dependencies required to run the code in this pull request have been added to the analysis module conda environment.yml file.
If applicable, R package dependencies required to run the code in this pull request have been added to the analysis module renv.lock file.

allyhawkins

I think this is a good start, but I had a few suggestions to add to the notebook just to help add more context and make the one large heatmap a little easier to interpret.

allyhawkins · 2026-08-05T16:17:50Z

+  patient_mets_aggr$labels$full, 
+  patient_prim_aggr$labels$full 
+)
+pheatmap::pheatmap(


Can we subset this to only show shared labels in the first place? I want to know if they have the same label in both references, are they similar. If they are different labels then they are probably distinct so we they just add noise to the heatmap and make it harder to interpret.

I actually don't think you need this heatmap at all. It's fine to keep, but if the labels don't match up then the value is going to be 0 anyways so it doesn't totally make sense to me to show both?

allyhawkins · 2026-08-05T16:21:21Z

+)
+```
+
+There's a 1:1 match of proliferative here only, but the noise around other cell types suggests that there are distinct expression profiles between shared cell types across these references.


Sure there's only a 1:1 match for proliferative, but I see some fairly strong signals for 4 out of the 6 tumor classes. It would be helpful to see a table of how many of each cell type is present in the original dataset alongside this just for added context. As in, how much of the population is assigned to this "Interactive" population vs the proliferative population?

allyhawkins · 2026-08-05T16:24:48Z

+
+## Conclusions
+
+Given the differences between primary/metastasis references for both patients and xenografts, it seems worth doing some initial exploration of annotation performance with `SingleR` next.


I think I would still run these together as a single reference with distinct labels. So you would attach primary- or met- to each of the labels prior to running. Then you can see which labels are assigned and if you get predominantly primary or met or an even mix of both within a single sample. You could then compare that to running as a single reference with harmonized labels where you don't distinguish between primary and met when building the reference model.

allyhawkins · 2026-08-05T16:26:01Z

+
+
+# export reformatted object if requested
+if (!is.null(opts$output_sce_file)) {


This seems strange to have a script solely for converting objects to SCE and then only do it optionally? I might make the output file required and remove the nesting you have here.

I see that you want to expand this script to convert to anndata here too, but I think I would still remove this as being optional and do both conversions automatically.

allyhawkins · 2026-08-05T16:29:56Z

+Rscript ${script_dir}/convert-osteocar.R --input_ref_file $mets_qs2 --output_sce_file $mets_sce
+Rscript ${script_dir}/convert-osteocar.R --input_ref_file $prim_qs2 --output_sce_file $prim_sce
+Rscript ${script_dir}/convert-osteocar.R --input_ref_file $xeno_mets_qs2 --output_sce_file $xeno_mets_sce
+Rscript ${script_dir}/convert-osteocar.R --input_ref_file $xeno_prim_qs2 --output_sce_file $xeno_prim_sce


This is fine, but there's gotta be a way in bash to turn this into a loop that's not terribly ugly and minimizes the duplication. Maybe try this?

Suggested change

Rscript ${script_dir}/convert-osteocar.R --input_ref_file $xeno_prim_qs2 --output_sce_file $xeno_prim_sce

input_files=("$mets_qs2" "$prim_qs2" "$xeno_mets_qs2" "$xeno_prim_qs2")

output_files=(

"${ref_dir}/patient_mets_sce.rds"

"${ref_dir}/patient_prim_sce.rds"

"${ref_dir}/xeno_mets_sce.rds"

"${ref_dir}/xeno_prim_sce.rds"

)

for i in "${!input_files[@]}"; do

Rscript "${script_dir}/convert-osteocar.R" \

--input_ref_file "${input_files[$i]}" \

--output_sce_file "${output_files[$i]}"

done

sjspielman · 2026-08-05T23:50:40Z

Changes include:

More heatmaps, showing both all cell types and the shared cell types (which of course are sometimes just shared in name only and not quite expression profiles...)
Barplots for distributions of cell types, as well as some more explicit printing of where they match up and don't between references
Script updates: use arrays (differently than suggested but still tighter than it was) in run-analysis.sh and remove the if from the conversion script

02-compare-references.nb.html

allyhawkins

This looks good to me. My comments are all fairly minor so I don't need to see this again. I think the biggest thing is just to make sure you add in an explanation at the top about what Ann_Level1 and Ann_Level2 are and then you probably don't need two copies of the same heatmap and can just show the shared cells.

allyhawkins · 2026-08-06T13:11:43Z

+---
+
+There are several different references in `OsteoCar`, and while we can generally map them to our sample propertoes, we have some sample which don't neatly line up.
+For example, `OsteoCar` has patient primary and metastatic (lung only) references, but we have a sample from metastatic tissue in bone which differs from the metastatic reference. 


A question for the future, is it really just one sample or is it multiple samples?

Definitely a typo, it's multiple......

> sample_df |> dplyr::filter(primary_or_metastasis=="Metastasis") |> dplyr::count(tissue_location) # A tibble: 10 × 2 tissue_location n <chr> <int> 1 Back 2 2 Femur 2 3 Kidney 1 4 Left Lung 1 5 Left upper lobe lung 1 6 Lung 29 7 Paraspinal 1 8 Right Lung 1 9 Right upper lobe lung 2 10 Thigh 1

allyhawkins · 2026-08-06T13:19:01Z

+)
+```
+
+The majority of cell types, including all classes of tumor cells, appear to be shared, but this does not guarantee they have similar expression.


It's been a while since I've read this paper, but can you remind me exactly which of these cell types they are calling "tumor" cell types vs. normal cell types? I'm actually surprised that the osteoblasts and chondrocytes population aren't overlapping with tumor, so I want to double check thats the case? Or is tumor vs normal represented by a different column in the data?

allyhawkins · 2026-08-06T13:20:26Z

+```{r}
+#| fig.height: 8
+#| fig.width: 10
+make_barplot(prim_coldata, mets_coldata) |>


Ah I see now which categories are which. Thank you for adding this since it's helpful to see them grouped by broad cell type.

allyhawkins · 2026-08-06T13:22:52Z

+We compare the references to one another to see how similar the cell type distributions and expression signatures are.
+As these references are quite large, they have to be aggregated before use, but we'll also explore their raw data too.


Can you add a note here about what's in Ann_Level1 vs Ann_Level2? That might help with some of the confusion I was experiencing (but have since resolved).

allyhawkins · 2026-08-06T13:25:14Z

+  patient_mets_aggr$labels$full, 
+  patient_prim_aggr$labels$full 
+)
+pheatmap::pheatmap(


I actually don't think you need this heatmap at all. It's fine to keep, but if the labels don't match up then the value is going to be 0 anyways so it doesn't totally make sense to me to show both?

allyhawkins · 2026-08-06T13:26:08Z

+Now, only the cell types shared between the references:
+
+```{r}
+shared_matched <- matched[shared, shared]


This is much more helpful because you can see the cell types that don't have any match that should have a match, which look to be mostly tumor cells. This tells me there probably is a difference in the marker gene expression for tumor cell subpopulations between primary and mets that we should be sure to consider.

allyhawkins · 2026-08-06T13:27:44Z

+  if (facet) {
+    p_counts <- p_counts + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)
+    p_props <- p_props + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)
+  }


Suggested change

if (facet) {

p_counts <- p_counts + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)

p_props <- p_props + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)

}

p_counts <- p_counts + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)

p_props <- p_props + facet_wrap(vars(Ann_Level1), scales = "free", nrow=3)

It looks like you don't actually use the facet option, so you can probably remove it.

allyhawkins · 2026-08-06T13:27:56Z

+## Functions
+
+```{r}
+make_barplot <- function(prim_coldata, mets_coldata, facet = TRUE) {


Suggested change

make_barplot <- function(prim_coldata, mets_coldata, facet = TRUE) {

make_barplot <- function(prim_coldata, mets_coldata) {

sjspielman · 2026-08-06T15:32:40Z

I actually don't think you need this heatmap at all. It's fine to keep, but if the labels don't match up then the value is going to be 0 anyways so it doesn't totally make sense to me to show both?

I like keeping it because this function does specific things to help you identify cases when cell types are effectively equivalent between references but the labels don't match - that's really its main purpose, but here we're using to compare references more directly. So, we actually do learn something from this heatmap - that a cell type labeled X in one reference apparently does not really correspond to cell type wit different label Y in the other. I will add this additional textual context to explain why to keep ot.

sjspielman added 6 commits August 5, 2026 08:07

add script, small for now, to convert objects to SCE

314b489

add the readme

43cde03

first comparison notebook

83fad5b

run conversion in bash script

cdad9be

update renv; this contains otehr packages coming in next PRs too

bcb8a65

readme

6103219

sjspielman requested a review from jaclyn-taroni as a code owner August 5, 2026 13:16

sjspielman requested review from allyhawkins and removed request for jaclyn-taroni August 5, 2026 13:16

allyhawkins reviewed Aug 5, 2026

View reviewed changes

sjspielman added 5 commits August 5, 2026 15:38

expand notebook with more distributions and heatmaps

c51794b

snapshot

702647d

use arrays

48edd2b

no if

d092ddd

fix typo

f73b4a0

sjspielman requested a review from allyhawkins August 5, 2026 23:50

allyhawkins approved these changes Aug 6, 2026

View reviewed changes

respond to reviews

ac757d0

sjspielman merged commit b07f64c into AlexsLemonade:main Aug 7, 2026
2 of 4 checks passed

sjspielman deleted the sjspielman/compare-references branch August 7, 2026 15:05


		## Conclusions

		Given the differences between primary/metastasis references for both patients and xenografts, it seems worth doing some initial exploration of annotation performance with `SingleR` next.



		# export reformatted object if requested
		if (!is.null(opts$output_sce_file)) {

-Rscript ${script_dir}/convert-osteocar.R --input_ref_file $xeno_prim_qs2 --output_sce_file $xeno_prim_sce
+input_files=("$mets_qs2" "$prim_qs2" "$xeno_mets_qs2" "$xeno_prim_qs2")
+output_files=(
+  "${ref_dir}/patient_mets_sce.rds"
+  "${ref_dir}/patient_prim_sce.rds"
+  "${ref_dir}/xeno_mets_sce.rds"
+  "${ref_dir}/xeno_prim_sce.rds"
+)
+for i in "${!input_files[@]}"; do
+  Rscript "${script_dir}/convert-osteocar.R" \
+    --input_ref_file "${input_files[$i]}" \
+    --output_sce_file "${output_files[$i]}"
+done

		We compare the references to one another to see how similar the cell type distributions and expression signatures are.
		As these references are quite large, they have to be aggregated before use, but we'll also explore their raw data too.

	make_barplot <- function(prim_coldata, mets_coldata, facet = TRUE) {
	make_barplot <- function(prim_coldata, mets_coldata) {

Uh oh!

Conversation

sjspielman commented Aug 5, 2026

Purpose/implementation Section

Author checklists

Analysis module and review

Reproducibility checklist

Uh oh!

allyhawkins left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

sjspielman commented Aug 5, 2026

Uh oh!

allyhawkins left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

sjspielman commented Aug 6, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants