Skip to content

[Feature request] ATAC count mode for metacell construction #20

@DmitriiSeverinov

Description

@DmitriiSeverinov

Description of feature

Hi all,

For snATAC-seq processing instead of ArchR, I am using Signac. ArchR uses insertions for feature quantification, while Signac uses fragments.

In the pp/metacell.py you convert insertion counts to fragment counts:

    for cluster_name in clusters:
        cell_idx = cluster_groups.get_group(cluster_name).index

        # RNA
        rna_vals = data_rna[cell_idx].X
        if sp.issparse(rna_vals):
            rna_vals = rna_vals.toarray()
        mean_rna = np.array(rna_vals.mean(axis=0)).ravel()
        mean_rna_list.append(mean_rna)

        # ATAC
        if len(set(cell_idx).intersection(data_atac.obs_names)) == 0:
            mean_atac_list.append(np.zeros(data_atac.shape[1]))
        else:
            atac_vals = data_atac[cell_idx].X
            if sp.issparse(atac_vals):
                atac_vals = atac_vals.toarray()
            # Convert insertions to fragments
            atac_bin = (atac_vals + 1) // 2
            mean_atac = np.array(atac_bin.mean(axis=0)).ravel()
            mean_atac_list.append(mean_atac)

        cluster_names.append(cluster_name)

In case of using Signac, since we already have the fragment counts, we essentially divide the number of fragments again by 2. I am not sure if that plays a big role for the downstream data processing and subsequent model training, but it would be nice to have a simple condition asking whether we have fragments or counts:

create_metacells(..., atac_count_mode = "insertions")
...
if atac_count_mode == "insertions":
    atac_bin = (atac_vals + 1) // 2
elif atac_count_mode == "fragments":
    atac_bin = atac_vals

I have noticed the same behaviour in pp/peak_selection.py

    for c_label in clusters:
        idx_cells = cluster_groups.get_group(c_label).index
        mat = data_atac[idx_cells].X
        if sp.issparse(mat):
            mat = mat.toarray()
        mat = (mat + 1) // 2  # get fragment presence
        mean_vec = mat.mean(axis=0).A1 if hasattr(mat, "A1") else mat.mean(axis=0)
        mean_list.append(mean_vec)

Best wishes,
Dmitrii

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions