R package implementation of Milo for testing for differential abundance in KNN graphs
Home Page: https://bioconductor.org/packages/release/bioc/html/miloR.html
License: GNU General Public License v3.0
At the moment this is incomplete and refers to old structure (e.g. nhoodDistances as matrix instead of list)
Matrix of cells x nhoods, something like
cell2nhoods <- function(x){
nhs <- lapply(nhoods(x), function(nh) as.vector(nh))
# nhood_mat <- matrix(nrow = ncol(x), ncol=length(nhoods(x)))
nhood_mat <- sapply(nhs, function(nh) ifelse(1:ncol(x) %in% nh, 1, 0))
return(nhood_mat)
}
Need to change the getter and setter for the graph slot to prevent name overwriting by igraph and vice versa
...to allow input from ScanPy (and Seurat?!)
It would be useful to incorporate a subsetting option in testNhood to exclude some samples in the DA analysis, basically excluding some samples from the count matrix and design matrix.
For example, a colleague is testing out Milo on a dataset where she wants to test for differences in abundance with age, but the big dataset includes a few samples that were FACS sorted differently, and DA in that subset is not of interest.
The best solution would be of course to exclude those samples at the beginning and rebuild all the graph, neighbourhoods etc, but this might be time consuming for a big dataset if you just want to do a quick check.
Would this be completely accounted for by adding a nuisance covariate in the design? Or would having a subsetting make interpretation easier?
The function complains about logcounts(x):
> x <- calcNhoodExpression(x, assay = "logcounts", subset.row = hvgs)
Error in t.default(data.set[subset.row, ]) : argument is not a matrixNot sure what's wrong, just running t(logcounts(x))[subset_rows,] works fine.
If I try to convert logcounts(x) to a non-sparse matrix the function never ends, the R session crashes, all burns 🔥
Thanks for developing this really awesome package.
I have been generating a single cell dataset, in which I have 11 time points and multiple conditions. I was wondering if it would be possible to implement an LRT test in milo to identify age:condition interactions which are DA in the graph.
Thanks!
Dylan
findNhoodMarkers can be quite slow, and the critical step seems to be the .group_nhoods_by_overlap function. I noticed that also here there is some looping in between many pairs of nhoods, similar to what was slowing down buildNhoodGraph before 42db3c4 :
https://github.com/MikeDMorgan/miloR/blob/4b192a3b90a2118055f2c4dfdeb80870f486c702/R/testDiffExp.R#L268
Maybe it could be refactored in a similar way? Or could the graph stored in the nhoodGraph slot be used here?
Try first searching for nhoods with at least some common cells using %in%, then calculate length(intersect(nh1, nh2)) just for those pairs
Any idea why the Milo object is much bigger than the SingleCellExperiment object, even before computing all the graphs and neighborhoods?
> sce
class: SingleCellExperiment
dim: 33694 50514
metadata(0):
assays(2): counts logcounts
rownames(33694): TSPAN6 TNMD ... RP11-107E5.4 RP11-299P2.2
rowData names(0):
colnames(50514): FCAImmP7179369-AAACCTGAGCCCAATT FCAImmP7179369-AAACCTGAGCCTATGT ... Human_colon_16S7985397-TTTGCGCAGAGCCCAA
Human_colon_16S7985397-TTTGGTTAGTACGCGA
colData names(14): Sample donor ... cell types umap_density_Age
reducedDimNames(1): MOFA
altExpNames(0):
> object.size(sce)
16486280 bytes
> milo <- Milo(sce)
> milo
class: Milo
dim: 33694 50514
metadata(0):
assays(1): expression
rownames(33694): TSPAN6 TNMD ... RP11-107E5.4 RP11-299P2.2
rowData names(0):
colnames(50514): FCAImmP7179369-AAACCTGAGCCCAATT FCAImmP7179369-AAACCTGAGCCTATGT ... Human_colon_16S7985397-TTTGCGCAGAGCCCAA
Human_colon_16S7985397-TTTGGTTAGTACGCGA
colData names(14): Sample donor ... cell types umap_density_Age
reducedDimNames(1): MOFA
altExpNames(0):
neighbourhoods dimensions(1): 0
neighbourhoodCounts dimensions(2): 1 1
neighbourDistances dimensions(2): 1 1
graph names(0):
> object.size(milo)
13639413808 bytesNot sure if this is a real issue, but it seems to make things like subsetting very slow and leads to frequent R crashes.
The buildGraph function is significantly slower than scran::buildKNNGraph. Here on a simulated dataset of 2000 cells:
> system.time(sim_milo <- buildGraph(sim_milo, k = 20,d = 30))
Constructing kNN graph with k:20
Retrieving distances from 20 nearest neighbours
user system elapsed
12.897 0.767 13.702
> system.time(sim2.knn <- buildKNNGraph(x=reducedDim(sim_milo, "PCA")[, c(1:30)], k=10, d=NA, transposed=TRUE))
user system elapsed
0.111 0.002 0.113
On the human thymus dataset (~ 30k cells) it takes forever. Any idea on what's the bottleneck here?
Maybe it's enough to throw an informative error message, but at the moment running with one gene gives a funny error.
> calcNhoodExpression(sim1.mylo, subset.row = c("Gene1", "Gene2"))
class: Milo
dim: 1000 1100
metadata(0):
assays(1): logcounts
rownames(1000): Gene1 Gene2 ... Gene999 Gene1000
rowData names(0):
colnames(1100): Cell1 Cell2 ... Cell1099 Cell1100
colData names(0):
reducedDimNames(1): PCA
spikeNames(0):
altExpNames(0):
nhoods dimensions(1): 100
nhoodCounts dimensions(2): 100 6
nhoodDistances dimensions(2): 1100 1100
graph names(1): graph
nhoodIndex names(1): 100
nhoodExpression dimension(2): 2 100
nhoodReducedDim names(0):
nhoodGraph names(0):> calcNhoodExpression(sim1.mylo, subset.row = c("Gene1"))
Error in t(neighbour.model) %*% t(data.set[subset.row, ]) :
non-conformable arguments
It's probably just about adding a drop=FALSE. @MikeDMorgan I can fix this in .calc_expression if that's ok with you.
Running findNhoodMarkers setting assay=counts generates this error:
Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extentTraceback:
3. | `colnames<-`(`*tmp*`, value = c("", "NGenes", "n.gene", "")) at testDiffExp.R#287
2. .perform_counts_dge(exprs, i.model, model.contrasts = i.contrast, gene.offset = gene.offset) at findNhoodMarkers.R#285
1. findNhoodMarkers(embryo_milo, da_results, assay = "counts", subset.nhoods = da_results$celltype == "Epiblast")The issue arises when setting a gene offset:
Lines 284 to 288 in 0883cd5
test.model only has one column, so colnames(test.model[, c(2:ncol(test.model))]) includes NGenes already.
The easy work-around is to set gene.offset=FALSE, but I am not super sure of the risks or of how the offset is used later in the DGE testing, since NGenes doesn't seem to be included in the design...
If setting the offset is not useful, it might be worth setting gene.offset=FALSE automatically if assay="counts"
There is an uninitialized sym1.milo object in the is_empty function called by testNeighbourhoods
which breaks testNeighbourhood
Having a look at the dataset included in the package, I wonder if we should load these as SingleCellExperiment objects to demonstrate usage the Milo constructor?
Dear authors,
First of all, the miloR package looks incredibly useful and generally well-documented.
My issue is about the lack of documentation for the output table of testNhoods(). The documentation says A data.frame of model results. In order to use and interpret miloR results correctly, I find it important to understand the table. Can you provide better documentation for the output table?
A few specific questions to get started:
logFC? What is its interpretation? How do you determine its 'direction'/'ratio' (FC=E7/E8 or the opposite)? (And does it make sense if analyzing > 2 experimental conditions (say E7, E7.5 and E8)logCPM?FDR instead of SpatialFDR? (If ever)Example output from the Mouse gastrulation example
:
da_results <- testNhoods(embryo_milo, design = ~ sequencing.batch + stage, design.df = embryo_design)
da_results
## logFC logCPM F PValue FDR Nhood SpatialFDR
## 7 -7.397680 20.20389 13.18174 2.868468e-04 0.008307757 7 0.007518813
## 15 -7.407804 19.72189 14.72980 1.263205e-04 0.008307757 15 0.007518813
## 53 7.168031 19.37642 12.80220 3.510315e-04 0.008307757 53 0.007518813
## 64 -7.164631 19.67509 14.31767 1.570634e-04 0.008307757 64 0.007518813
## 69 -7.652294 19.96889 14.84507 1.188619e-04 0.008307757 69 0.007518813
## 71 -7.357238 19.56574 15.26142 9.542589e-05 0.008307757 71 0.007518813
I have a terrible feeling that much of the time taken during graph building on large data sets, and subsequently the large memory footprint is taken up by the distance calculations.
If we can separate these two steps (we currently delegate the distance calculations to BiocNeighbours::findKNN), we could optimise these two steps independently.
To avoid getting ridiculously small p-values
We can use the ones in the main milo repo for this. I need to double-check what the format and minimal requirements are for Bioconductor here.
If labelled genes are at the top of a plot, the the legend overlaps them and they cannot be seen. Using the usual tricks outside of the function to adjust legend position does not work - I assume this is something to do with cowplot. It would be better to either not gather the legends, or to have them positioned sufficiently far from the main plot to allow gene labeling at any row of the heatmap.
Example:
A minor thing: both countCells and testNeighbourhoods have a required parameter called data, but for counting this needs to be the colData of the Milo object, while for testing it's the design matrix. I would change the name to design_df or something like that in testNeighbourhoods.
The plotNhoodGraphDA can be hard to interpret when there is a high density of neighbourhoods. Example:
It might be worth plotting 2 layers: 1) Non-DA neighbourhoods (white fill), 2) DA neighbourhoods (coloured by LFC). That way the DA neighbourhoods are never obscured by the non-DA neighbourhoods.
For gene expression signatures, is it better to have counts previously corrected by batch effect or just include them in the model?
And also, is it possible to have sort of the intermediate output of findNhoodMarkers, where you have the list of cells that were grouped together regarding the da results and the design? I just want to check there is no further batch effect in those groups
Thanks
I get this error when running testNeighbourhoods:
> milo_results <- testNeighbourhoods(sim_milo, ~ 1 + condition, design.df = design_df, fdr.weighting = "k-distance")
Error in rowSums(x.slot) :
'x' must be an array of at least two dimensionsI found this comment suggesting that it might be an underlying problem of having a matrix of class dgCMatrix in SingleCellExperiment slots, but if I convert the counts to a regular matrix I get a weirder error:
> neighbourhoodCounts(sim_milo) <- as.matrix(neighbourhoodCounts(sim_milo))
the condition has length > 1 and only the first element will be used
> milo_results <- testNeighbourhoods(sim_milo, ~ 1 + condition, design.df = design_df,
+ fdr.weighting = "k-distance")
the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedError in if (.is_empty(x, "neighbourhoodCounts")) { :
argument is of length zeroAny suggestions here
A problem with my installation?
Thanks
Mouse gastrulation example:
embryo_milo <- Milo(embryo_data)
Error in validObject(.Object) :
invalid class “Milo” object: 1: invalid object for slot "nhoodCounts" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
invalid class “Milo” object: 2: invalid object for slot "nhoodExpression" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
invalid class “Milo” object: 3: invalid object for slot "nhoodAdjacency" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
Simulated dataset example:
traj_milo <- Milo(traj_sce)
Error in validObject(.Object) :
invalid class “Milo” object: 1: invalid object for slot "nhoodCounts" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
invalid class “Milo” object: 2: invalid object for slot "nhoodExpression" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
invalid class “Milo” object: 3: invalid object for slot "nhoodAdjacency" in class "Milo": got class "ddiMatrix", should be or extend class "matrixORMatrix"
sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.1 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] MouseGastrulationData_1.4.0 patchwork_1.1.1
[3] dplyr_1.0.2 scater_1.18.3
[5] ggplot2_3.3.3 SingleCellExperiment_1.12.0
[7] SummarizedExperiment_1.20.0 Biobase_2.50.0
[9] GenomicRanges_1.42.0 GenomeInfoDb_1.26.2
[11] IRanges_2.24.1 S4Vectors_0.28.1
[13] BiocGenerics_0.36.0 MatrixGenerics_1.2.0
[15] matrixStats_0.57.0 miloR_0.1.0
[17] edgeR_3.32.0 limma_3.46.0
loaded via a namespace (and not attached):
[1] ggbeeswarm_0.6.0 colorspace_2.0-0
[3] ellipsis_0.3.1 rprojroot_2.0.2
[5] scuttle_1.0.4 XVector_0.30.0
[7] fs_1.5.0 BiocNeighbors_1.8.2
[9] farver_2.0.3 remotes_2.2.0
[11] graphlayouts_0.7.1 ggrepel_0.9.0
[13] bit64_4.0.5 RSpectra_0.16-0
[15] interactiveDisplayBase_1.28.0 AnnotationDbi_1.52.0
[17] fansi_0.4.1 codetools_0.2-18
[19] sparseMatrixStats_1.2.0 pkgload_1.1.0
[21] polyclip_1.10-0 dbplyr_2.0.0
[23] uwot_0.1.10 ggforce_0.3.2
[25] shiny_1.5.0 BiocManager_1.30.10
[27] compiler_4.0.3 httr_1.4.2
[29] assertthat_0.2.1 Matrix_1.3-0
[31] fastmap_1.0.1 cli_2.2.0
[33] later_1.1.0.1 tweenr_1.0.1
[35] BiocSingular_1.6.0 prettyunits_1.1.1
[37] htmltools_0.5.0 tools_4.0.3
[39] rsvd_1.0.3 igraph_1.2.6
[41] gtable_0.3.0 glue_1.4.2
[43] GenomeInfoDbData_1.2.4 rappdirs_0.3.1
[45] Rcpp_1.0.5 vctrs_0.3.6
[47] ExperimentHub_1.16.0 DelayedMatrixStats_1.12.1
[49] ggraph_2.0.4 stringr_1.4.0
[51] ps_1.5.0 testthat_3.0.1
[53] beachmat_2.6.4 mime_0.9
[55] lifecycle_0.2.0 irlba_2.3.3
[57] gtools_3.8.2 devtools_2.3.2
[59] AnnotationHub_2.22.0 zlibbioc_1.36.0
[61] MASS_7.3-53 scales_1.1.1
[63] tidygraph_1.2.0 promises_1.1.1
[65] yaml_2.2.1 curl_4.3
[67] memoise_1.1.0 gridExtra_2.3
[69] stringi_1.5.3 RSQLite_2.2.1
[71] BiocVersion_3.12.0 desc_1.2.0
[73] pkgbuild_1.2.0 BiocParallel_1.24.1
[75] rlang_0.4.10 pkgconfig_2.0.3
[77] bitops_1.0-6 lattice_0.20-41
[79] purrr_0.3.4 labeling_0.4.2
[81] processx_3.4.5 cowplot_1.1.1
[83] bit_4.0.4 tidyselect_1.1.0
[85] RcppAnnoy_0.0.18 magrittr_2.0.1
[87] R6_2.5.0 generics_0.1.0
[89] DelayedArray_0.16.0 DBI_1.1.0
[91] pillar_1.4.7 withr_2.3.0
[93] RCurl_1.98-1.2 tibble_3.0.4
[95] crayon_1.3.4 BiocFileCache_1.14.0
[97] usethis_2.0.0 viridis_0.5.1
[99] locfit_1.5-9.4 grid_4.0.3
[101] FNN_1.1.3 callr_3.5.1
[103] blob_1.2.1 digest_0.6.27
[105] xtable_1.8-4 tidyr_1.1.2
[107] httpuv_1.5.4 munsell_0.5.0
[109] beeswarm_0.2.3 viridisLite_0.3.0
[111] vipor_0.4.5 sessioninfo_1.1.1
For some functions we have one, for some the other.
Also it's quite long to write all the time. I would be pro using nhoods everywhere instead.
Hi!
I wanted to try Milo today, but got an error when trying to build the milo object from a graph adjacency matrix. In fact, buildFromAdjacency() throws the same error even when run with the example code in its documentation:
r <- 1000
c <- 1000
k <- 35
m <- floor(matrix(runif(r*c), r, c))
for(i in seq_along(1:r)){
m[i, sample(1:c, size=k)] <- 1
}
buildFromAdjacency(m)
output:
Casting to sparse matrix format
Inferring k from matrix
Error in validObject(x) : invalid class “Milo” object: 1:
no 'rowPairs' field in 'int_elementMetadata'
invalid class “Milo” object: 2:
no 'colPairs' field in 'int_colData'
I tried a few different matrices/sparse matrices etc. but always with the same outcome.
Thanks!
On a subset of the gastrulation data
Which step of refined sampling is slowing down the calculation?
What is the purpose of this conditional statement supposed to be? https://github.com/MikeDMorgan/miloR/blob/e86041a824c4308956ce59ce46aa7b383da9205a/R/buildGraph.R#L68
If I have a Milo object with only one slot named PCA in reducedDim(x) it wants to recompute the PCA.
If plotNhoods has been run then the between neighbourhood distances have already been calculated - use these for the neighbourhood grouping step - otherwise compute.
If you have IDs as rownames and you want to use gene names stored in rowData
I'd also add and option to subset nhoods for testing. E.g. in the liver case I want to do the test just in one lineage. I can subset manually the Milo obj and the results data.frame, but it would be a nice feature.
Running the newer makeNeighbourhoods on a data set of ~8k cells still takes 2-3 minutes to compute:
> mylo
class: Milo
dim: 27998 8419
metadata(1): log.exprs.offset
assays(2): counts logcounts
rownames(27998): ENSMUSG00000051951 ENSMUSG00000089699 ... ENSMUSG00000096730 ENSMUSG00000095742
rowData names(0):
colnames(8419): SIGAC1_AAACCTGCACAGACAG-1 SIGAC1_AAACCTGCAGCTCGCA-1 ... SIGAF7_AtohPos_TTTGTCATCAACGGCC-1
SIGAF7_AtohPos_TTTGTCATCGAATGGG-1
colData names(0):
reducedDimNames(1): PCA
spikeNames(0):
altExpNames(0):
neighbourhoods dimensions(1): 0
neighbourhoodCounts dimensions(2): 1 1
neighbourDistances dimensions(2): 8419 8419
graph names(1): graph
neighbourhoodIndex names(1): 0
> system.time(mylo <- makeNeighbourhoods(mylo, k=21, prop=0.1, refined=TRUE))
Checking valid object
Using refined sampling
user system elapsed
147.487 13.997 163.734
I don't know if this is because it only really matters for much bigger data sets. How does it behave in your hands @emdann ?
This would be useful for plotting purposes
I was running testNhoods on a small dataset and I got this error message:
Error in DGEList(counts = nhoodCounts(x)[keep.nh, ], lib.size = log(colSums(nhoodCounts(x)))) :
negative library sizes not permittedThe problem is that the cells of a one sample are not captured in any neighbourhood.
The quick (and probably the proper) fix is to remove this sample with very few cells from testing, but I think at least an error message here or a warning when running countCells could be useful. Ideas?
The reduced.dim parameter is not passed to .buildGraph, so default ("PCA") is used
Lines 118 to 119 in 7b15411
I was rerunning the liver analysis script, after reinstalling the latest miloR version, and I get an error I am not sure how to debug. Running:
nhood_markers <- findNhoodMarkers(liver_milo, milo_res, overlap=1, return.groups = TRUE,
subset.row = endo_hvgs,
subset.nhoods = str_detect(milo_res$annotation_indepth, "Endo"))
I get:
Error in UseMethod("isSymmetric") : no applicable method for 'isSymmetric' applied to an object of class "c('double', 'numeric')"
The traceback says the problem is in
.group_nhoods_from_adjacency(nhoods(x), nhood.adj = nhoodAdjacency(x), da.res = da.res, is.da = da.res$SpatialFDR < da.fdr, merge.discord = merge.discord, overlap = overlap, subset.nhoods = subset.nhoods)
Any ideas?
Minor: running buildFromAdjacency still throws a message Adding nhoodDistances to Milo object, but the nhoodDistances slot is not populated anymore.
Lines 99 to 102 in 595566b
This leads to wrong matching if the user provides subset.nhoods as indices
Lines 106 to 109 in 595566b
nhoodAdjacency(x) are in the same order as nhoods(x), but I find this is often not the case. Is nhoodAdjacency(x) added just via buildNhoodGraph ? If that's the case it might be an issue with the reordering required for plotting.Write unit tests for findNhoodMarkers and testDiffExp.
Write tests for buildNhoodGraph and plotNhoodGraph
Not urgent, but atm buildGraph looks for a dimensionality reduction named PCA, but I might want to use a precomputed dim reduction with a different name (e.g. in the mouse gastrulation atlas the corrected reduction is called pca.corrected)
testNhoods returns a data.frame with results, not a Milo object
https://github.com/MikeDMorgan/miloR/blob/7fb8c70d667536230ed009a9437b623767c54f49/R/testNhoods.R#L36
This is a place-holder, but for usability we need to make countCells handle dplyr data.frame-alikes.
Except in the case of the multiple discrete blocks, the neighbourhood grouping almost always only generates 2 groups if I use the same FDR threshold as the DA testing (relaxing this gives more granular grouping in line with expectations in the mouse thymus data). Even an FDR 100% still only returns 4 groups...
My concern is that the simple approach using weakly connected components struggles when the graph is highly connected. I propose that we use a community detection e.g. Louvain/Infomap/Walktrap or something to get more granularity in the grouping perhaps?
Hi I am very interested in trying miloR but I am a Python user. I was wondering if there is any tutorials/example using adata as input for miloR. Thank you very much!
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.