lcolladotor
commented
2 years ago There are a few challenges / questions to resolve that I do not know the answer to.
aggregateAcrossCells vs manual pseudo-bulking
- Example manual pseudo-bulking code: https://github.com/LieberInstitute/spatialDLPFC/blob/ed33f1fc1124b754028f47bd655ece628c70405c/code/analysis/07_spatial_registration/07_spatial_registration.R#L28-L47
- Example
aggregateAcrossCellscode: https://github.com/LieberInstitute/Visium_IF_AD/blob/17943ea9a1659eac7992e27ab42dc79fb69dbb81/code/08_harmony_BayesSpace/04_heatmaps_dlpfc_markers.R#L86-L94 - When we manually pseudo-bulk, we don't set any
sizeFactorsso they have to be calculated bylogNormCounts(). MaybeaggregateAcrossCells()does keep some of the size factors info around. I don't know! Basically, I don't know if we get exactly the same normalized counts with both methods.
design variables
Assuming we use aggregateAcrossCells(), we would need to make sure that any variables we use later in the model like at https://github.com/LieberInstitute/spatialDLPFC/blob/ed33f1fc1124b754028f47bd655ece628c70405c/code/analysis/07_spatial_registration/07_spatial_registration.R#L72 are part of the DataFrame we pass to the function https://github.com/LieberInstitute/Visium_IF_AD/blob/17943ea9a1659eac7992e27ab42dc79fb69dbb81/code/08_harmony_BayesSpace/04_heatmaps_dlpfc_markers.R#L89-L92.
For that, I thought that it could be useful to look at other packages like ExploreModelMatrix which deal with design formulas. I see that they have this nice check to verify that the object passed is actually a formula: https://github.com/csoneson/ExploreModelMatrix/blob/master/R/ExploreModelMatrix.R#L67-L69.
Once we have a formula object, we can extract the terms like with:
fo <- y ~ x1*x2
attributes(terms(fo))$term.labels
# [1] "x1" "x2" "x1:x2"and then we can get the column names with:
terms <- attributes(terms(fo))$term.labels
terms[!grepl(":", terms)]
# [1] "x1" "x2"Function initial structure
foo <- function(sce, designFormula, block_var = NULL) {
terms <- attributes(terms(fo))$term.labels
pd <- colData(sce)[ , c(terms[!grepl(":", terms)], block_var) ]
sce_pseudo <- aggregateAcrossCells(sce, pd)
## then model
## then re-arrange results
}Re-arrange results
We would like to re-arrange the results so they look like spatialLIBD::fetch_data("modeling_results"). That can be done with code from https://github.com/LieberInstitute/HumanPilot/blob/879f11c7c57efd72334332b40feb3ad623e067c8/Analysis/Layer_Guesses/misc_numbers.R#L302-L318.
Check that it's compatible with downstream functions
We would then need to check that it works with http://research.libd.org/spatialLIBD/reference/layer_stat_cor.html and http://research.libd.org/spatialLIBD/reference/layer_stat_cor_plot.html. This new function would basically provide the new stats argument for layer_stat_cor(). Though it could also replace the modeling_results argument, if we want to do spatial registration against a new reference. Note that modeling_results is a list() object.
Add support for anova and pairwise
We could expand this new function to also run the anova and/or pairwise models. Right now this code only does the enrichment model https://github.com/LieberInstitute/spatialDLPFC/blob/ed33f1fc1124b754028f47bd655ece628c70405c/code/analysis/07_spatial_registration/07_spatial_registration.R#L84-L98.
abspangler13
It'll be easier to run the spatial registration process in your own data if we have functions for doing something like https://github.com/LieberInstitute/HumanPilot/blob/master/Analysis/Layer_Guesses/ad_snRNAseq_recast.R where we pseudo-bulk the data then compute t-stats for one cluster vs all others at a time, aka, the enrichment statistics. https://github.com/LieberInstitute/HumanPilot/blob/master/Analysis/Layer_Guesses/layer_specificity.R has the code for the whole set: anova and pairwise. Then https://github.com/LieberInstitute/HumanPilot/blob/master/Analysis/Layer_Guesses/misc_numbers.R re-shapes the results into what we can load with
fetch_data("modeling_results").This would help @abspangler13 @Nick-Eagles @sparthib and others.