Metacell - Single-cell RNA Sequencing Analysis¶
The metacell package implements the metacell algorithm for single-cell RNA sequencing (scRNA-seq) data analysis.
The original metacell algorithm was implemented in R. In contrast, the Python package described here implements an improved <https://todo/> version, which uses improved internal algorithms and supports much larger data sets (millions of cells).
Naively, scRNA_seq data is a set of cell profiles, where for each one, for each gene, we get a count of the mRNA molecules that existed in the cell for that gene. This serves as an indicator of how “expressed” or “active” the gene is.
As in any real world technology, the raw data may suffer from technical artifacts (counting the molecules of two cells in one profile, counting the molecules from a ruptured cells, counting only the molecules from the cell nucleus, etc.). This requires pruning the raw data to exclude such artifacts. TODO: Create view should consider these issues.
The current technology scRNA-seq data is also very sparse (typically <10%, sometimes even <5% of the RNA molecules are counted). This introduces large sampling variance on top of the original signal, which itself contains significant inherent biological noise.
Analyzing scRNA-seq data therefore requires processing the profiles in bulk. Classically, this has been done by directly clustering the cells using various methods.
In contrast, the metacell approach groups together profiles of the “same” biological state into groups with the minimal number of profiles (~100) needed for computing robust statistics (in particular, mean gene expression). Each such group is a single “metacell”.
By summing profiles together, each metacell greatly reduces the sampling variance, and provides a more robust estimation of some transcription state. In particular, a metacell is not a cell type (multiple metacells may belong to the same type), and is not a parametric model of the cell state.
The metacells should therefore be further analyzed using additional methods to classify cell types, detect cell trajectories and/or lineage, build parametric models for cell behavior, etc. Using metacells as input for such analysis techniques should benefit both from the more robust, less noisy input; and from the (~100-fold) reduction in the number of profiles to analyze.
An obvious technique is to recursively group the metacells into larger groups, and investigate the resulting clustering hierarchy. Methods for further analysis of the data may choose to build upon this readily available hierarchy (which is also provided by this package).
Given Python version 3.6 and above, run
pip install --user metacell (or,
sudo pip install
metacell if you have
sudo privileges). This will install this package and two scripts:
The metacell package is built using the DynaMake package.
metacell_flow script is simply the
dynamake script, pre-loaded with the metacell build
metacell_call is simply the
dynacall script, pre-loaded with the metacell
To use the package, create a directory somewhere with a large amount of available space (for ~2M
cells, metacell required up to 200GB of disk space). If using a servers cluster, on some shared file
system, and provide a DynaMake configuration file with an appropriate
distributing jobs across the cluster; see the DynaMake documentation for details.
raw sub-directory and download the raw_data into it. Then, run
import the gene and cell data into the fast-access format used by the metacell package. This is
unfortunately a slow process (in particular, the delaning step). Luckily, you only need to do this
This will create the following directories:
- Describes the genes used in the data.
- A COLLECTION profiles directory containing all the imported cells.
The provided automated workflow uses the following directory structure. The specifics of each directory type are provided further below.
A directory which is typically named after the data source (
HCA, etc.). Will contain:
A directory containing the raw data from the data source. Should contain one or more:
The first processing step will convert the
rawdata to the standard input format. The result will contain:
- A directory in the standard (STD) format ready to be delaned.
Creating the standard format directory is pretty fast. It typically requires either a simple copy or applying something like
sedto fix minor formatting issues. For HCA files, it is just a matter of dumping the data into a textual format.
The next processing step will split the input into distinct lanes. The result will contain:
A directory for each of the standard format directories. This would be either a standard format directory, or will contain:
- A directory in the standard format containing just the data for a single lane.
When data contains multiple lanes, delaning is a painfully slow process as it needs to split the matrix market file, processing one line at a time (without any parallelization). It would have been much easier if the raw data was already split into lanes.
A directory containing fast access per-gene data (GENES) extracted from the
The extraction process first verifies this data is identical for all the
delaneddirectories. Then one of these directories is used to actually extract the data. This is very fast, as there are only ~30,000 genes to deal with.
A directory containing fast access per-gene-per-profile (COLLECTION) data. Will contain:
A directory for each of the
delaneddirectories. This would either be a simple profiles batch (BATCH), or a collection combining:
- A profiles batch directory containing the data for a single lane.
Conversion of the (delaned) textual data to the fast-access format isn’t very fast, but this is offset by being able to run this conversion in parallel for each lane. It is therefore much less painful than delaning. Once conversion is done, the import process is complete and does not need to be repeated.
This is where the metacells package performs its analysis, recursively grouping the cells into metacells. Should contain:
Some view profiles container (VIEW) based on the full
cellsprofiles collection, to group. Typically such a view will exclude profiles with too-low total UMIs count, or that suffer from other technical issues (e.g. doublets). Genes may also be excluded.
The main function of the metacell package is to recursively compute the groups for such a view (that is, extend it to become a GROUPED directory).
The formats which can be used as the raw input for metacell processing:
A directory containing 10X Genomics data. Should contain the files:
- Matrix market UMIs data.
- The barcode of each profile (possibly with a
- The ensembl ID and name of each gene.
A directory containing Mouse Organogenesis Cell Atlas data. Should contain the files:
- Matrix market UMIs data.
- Per-cell data (see list of expected fields in
- Per-gene data (see list of expected fields in
- A file containing Human Cell Atlas in HDF format with a
.h5` suffix. Must contain a single top-level key with the sub-keys
These formats are converted to a standard input format:
A directory containing data in a standard metacell input format. Should contain the files:
- Matrix market UMIs data.
- Per-profile data.
- Per-gene data.
- A description of the fields of the CSV files.
Fast Access Data¶
Formats of data allowing fast access during the metacell processing.
A directory containing per-gene data, accessed using the
metacell.storage.genes.Genesclass. Should contain the files:
- General meta-data describing this set of genes. Should contain the keys
uuid, which must be equal to the
genes.name.txtfile. This UUID uniquely identifies this set of genes.
- Per-gene string data. These should include
genes.name.txtcontaining the (sorted, all-upper-case) gene names, and
genes.ensembl.txtcontaining the ensembl identifier of each gene.
- Per-gene numeric data, if any.
Either a BATCH, COLLECTION, or VIEW directory, providing access to per-gene-per-profile data for some genes and profiles using the
Each container directory contains at least the file:
General meta-data describing this container of profiles. Should contain at least the keys:
- The organism whose cells were sequenced (e.g.
- The UUID of the genes per profile.
- The number of profiles, and number of genes per profiles.
- The kind of profiles (
- Should be equal to the
matrix.umis.npyfile. This UUID uniquely identifies this profiles batch.
- Additional arbitrary meta-data collected for this container profile over time.
Additional keys may be specified depending on the container type.
A directory containing per-gene-per-profile data for some closely related profiles (a single batch), accessed using the
metacell.storage.profiles.ProfilesBatchclass. Barcodes can be safely assumed to be unique within each batch, but not between multiple batches. When investigating batch effects, this is the smallest group of profiles for which these effects can be assumed to be the same. Should contain the additional files:
Per-gene-per-profile data. The
mmmformat is a metacell-specific compressed sparse data format accessed using the
metacell.matrices.MemoryMappedMatrixclass. Each batch should contain
Additional per-gene-per-profile data may exist, for example:
- The log-fold-factor of the gene’s UMIs fraction (out of the profile’s total), compared to the gene’s median UMIs fraction (in all the profiles). This serves as an estimate of how strong a marker the gene is for this specific profile. This should not be used for single cells.
Per-profile string data. Common data includes:
- The unique (in this batch) per-profile barcode (without any lane suffix).
- The unique (in this batch) per-profile string name, if any.
- The name of the batch and the name of the lane (if any) per profile. These are the same for all profiles in the batch, but having these files is useful when collecting multiple batches to a single data set (see below).
Per-profile numeric data. Common data includes:
- The sum of the genes UMIs for each profile.
- The maximal fold factor of a gene in this profile.
Per-gene string and numeric data, specific for this profiles batch, if any. Common data includes:
- The median fraction of the gene’s UMIs out of the profile’s UMIs. This is used as an estimate of the global gene expression. This should not be used for single cells.
- The minimal rank of the gene’s fold factor in some profile. That is, the lower the number, the stronger the marker the gene is for distinguishing between profiles. This is used to filter genes in heatmap visualizations. This should not be used for single cells.
A directory containing a possibly restricted view of another profiles container accessed using the
profiles.yamlfile should contain the additional keys:
- The UUID of the base profiles container this is a view of.
- The path of the base profiles container this is a view of. If this is a relative path, it is relative to the view directory.
The directory may contain the additional file(s):
- If exist, contain the sorted indices of the base container profiles which are included in this view. Otherwise, the view contains all the base container profiles.
- If exist, contain the sorted indices of the base container genes which are included in this view. Otherwise, the view contains all the base container genes.
When accessing data of a view, if there exists a local file containing the data, it is used. Otherwise, if possible, the data of the base container is fetched and sliced to include only the relevant subset of profiles and/or genes. By default, data is assumed to be slicable. Non-slicable data is specified in
A directory combining a collection of other profiles container directories into a single unified profiles container, accessed using the
profiles.yamlfile should contain the additional key:
A list of mapping, one per combined profile container directories. Should contain the keys:
- The name of a sub-directory of the collection directory, containing the profiles container which is combined into the overall collection.
- The UUID identifying the combined profiles container.
- The number of profiles in the combined profiles container.
For each listed
name, a sub-directory with that name must exist, and contain a profiles container. When accessing data of a collection, if there exists a local file containing the data, it is used. Otherwise (except for per-gene data), the data is combined from the data of all the collected sub-containers.
The primary function of the metacell package is to group profiles together, in particular, to group cell profiles into metacells, metacells into clans, etc.
When groups are computed, the following files are created:
- The group index for each of the container’s profiles. A negative value indicates the profile belongs to no group (is an “outlier”).
- A profiles VIEW containing just the profiles with a non-negative group index, that is, only the non-outlier profiles.
- A profiles VIEW containing just the outlier profiles (which are not included in any
of the groups). Typically this would only hold
p.base_index.npyin addition to the
- A GROUPS profiles container for the computed groups.
- A hierarchy of ZONE (view) containers based on the computed groups.
Temporary files used to compute the groupings, but need not be used once it has been computed. The contents of this directory depends on whether we are grouping a few profiles or many profiles.
If grouping a small number of profiles:
- A view of the profiles container, including only the selected (“feature”) genes.
If grouping a large number of profiles, we can’t directly compute the groups, since the basic algorithm has a complexity of O(N^2). Computation therefore must works through a divide-and-conquer algorithm. This requires different temporary files:
A directory containing the second phase’s files.
- A profiles view of a random subset of a restricted number of to-be-grouped profiles.
To help in sorting file names, the indices are zero-padded so all have the same
00, .., 10, 11, ...). Each such pile is independently (directly) grouped.
This directory is a profiles view containing all profiles that were marked as an outlier in some indexed pile. It is also independently grouped. Unlike the indexed piles, this view may contain a large number of profiles so grouping it may require an intermediate step of its own.
Rare “types” do not have enough profiles in each of the piles to form groups, so will be marked as outliers. In this outliers-only view, there might be enough of these rare profiles to form a proper group.
When grouping the profiles in the final first phase outliers pile (small enough to be directly grouped), we disable outliers detection. That is, we insist that each of the profiles be assigned to some group by the end of the first phase, as these groups are used to compute the piles for the second phase and we’d like to give each profile a chance to be grouped in this second phase.
- A view of the grouped profiles container, with
p.group_index.npycontaining the combined group indices from all the phase-1 piles.
The intermediate grouping combines the groups from all the phase-1 piles, as well as the groups from the phase-1 outliers. This is the less-accurate grouping we will improve upon. To allow doing so, we compute the groups-of-groups (that is, we compute
p.group_index.npyhere), and use these to construct the second phase piles.
Similarly to grouping the outliers pile, we disable outlier detection when grouping the final pile, so that each group will be assigned to one of the second phase piles.
A directory containing the second phase’s files. This is similar to the first phase’s structure, with the following differences:
- The profiles contained in each pile are not chosen randomly. Instead, they are the profiles contained in one of the computed intermediate groups-of-groups.
- The outliers pile contains the phase-1 outliers in addition to the outliers of the phase-2 indexed piles. This allows very rare profile “types” a second chance at being grouped.
- There are no
groupsdirectories; the results are written directly to the original (parent) grouped profiles container. These final groups are recursively grouped as long as there is a large number of profiles, as described above.
A profiles container describing groups of profiles. The
gp.umis.npyper-profile-per-gene UMIs count in this container will contain the sum of the per-profile-per-gene of each of the grouped profiles contained in each (group) profile.
A groups profiles container has identical format to a normal profiles batch directory, except that its
profiles.yamlfile contains a
grouped_profileskey with the (relative) path to the grouped profiles container (typically
..), and a
grouped_uuidkey with the UUID of the grouped profiles container. For convenience,
mean_cells_countprovide an estimate of the number of cells and profiles in each group.
profiles_kindof the groups container should be based on the grouped container profiles kind. In particular, groups of
metacells. Such groups are expected to contain cells of the “same” state.
When grouping a large number of profiles, the
groupsdirectory might itself be grouped. This recursion continues as long as the total number of profiles is “large”. Groups of
clans, groups of
tribes, groups of
nations, and groups of
hordes. All higher-level groups contain profiles with a “similar”, but not necessarily the “same” state.
Finally, a groups container will hold the files
p.grouped_count.npywith the number of profiles summed into each group, and
p.cells_count.npywith the number of cells summed into each group. These two numbers are the same for metacells, but differ for higher level groups.
A tree directory is used to provide convenient views for exploring the hierarchy created by the computed groups. In each level of the tree, in addition to the
groupedcontainers, there is also a
zonesdirectory that contains views that group the overall profiles into groups-of-groups.
The directory structure for a small amount of cells would contain:
zones groups # All the clans (groups of metacells) ... grouped # All the metacells (grouped into clans). ... clan-0 # The 1st clan groups # All the metacells in the 1st clan ... grouped # All the cells in the 1st clan ... clan-1 # The 2nd clan groups # All the metacells in the 2nd clan ... grouped # All the cells in the 2nd clan ... ... ...
And a directory structure for a larger number (several hundred of thousands) of cells would contain:
zones groups # All the tribess (groups of clans) ... grouped # All the clans (grouped into tribes). ... tribe-0 # The 1st tribe groups # All the clans in the 1st tribe ... grouped # All the metacells in the 1st tribe ... clan-000 # The 1st clan in the 1st tribe groups # All the metacells in the 1st clan ... grouped # All the cells in the 1st clan ... clan-001 # The 2nd clan in the 1st tribe ... ... tribe-1 # The 2nd tribe groups # All the clans in the 2nd tribe ... grouped # All the metacells in the 2nd tribe ... clan # The 1st clan in the 2nd tribe groups # All the metacells in the 1st clan ... grouped # All the cells in the 1st clan ... clan-091 # The 2nd clan in the 2nd tribe ... ... ...
A much larger number (tens of millions) of cells would include an additional hierarchy level for each
nation. Even larges number of cells would (up to a few billion) would require an additional level for each
horde. Currently the largest data sets are at most of a few million so this scheme should be sufficient for the foreseeable future.
In all cases, the
groupeddirectory contains a
p.group_index.npywhich specifies the index of the group in the adjacent