Installation

conda env create -n cigma -f envs/environment.yaml
conda activate cigma

pip install cigma

install.packages(c("optparse", "numDeriv", "Matrix"))

Testing

To verify the installation:

python3 -m cigma.tests.test

Quick start

CIGMA provides helper functions to generate model input from scRNA-seq data that have undergone quality control and normalization:

import scanpy as sc
import pandas as pd
from cigma import preprocess, fit

# Read scRNA-seq data in h5ad format.
# The data contains: a gene expression matrix (cells x genes) and cell metadata
# with columns: 'cell' (cell IDs), 'ind' (individual IDs), 'ct' (cell types).
ann = sc.read_h5ad(scRNA_h5ad_file)

# Compute pseudobulk quantities:
#   ctp  - cell type-pseudobulk matrix (individual-cell type pairs x genes)
#   ctnu - cell-to-cell variation matrix (individual-cell type pairs x genes)
#   P    - cell type proportion matrix (individuals x cell types)
#   n    - cell count matrix (individuals x cell types)
ctp, ctnu, P, n = preprocess.pseudobulk(
    ann=ann, ind_col='ind', ct_col='ct', cell_col='cell'
)
# or
# ctp, ctnu, P, n = preprocess.pseudobulk(
#    X=ann.X, obs=ann.obs, var=ann.var, ind_col='ind', ct_col='ct', cell_col='cell'
# )

# Remove genes or cell types to match the requirement of CIGMA that genes are expressed in all cell types.
# For example, keep cell types CT1, CT2, and CT3 that all genes are expressed in, and remove other cell types.
ctp = ctp.loc[ctp.index.get_level_values('ct').isin(['CT1', 'CT2', 'CT3'])]
ctnu = ctnu.loc[ctnu.index.get_level_values('ct').isin(['CT1', 'CT2', 'CT3'])]
P = P[['CT1', 'CT2', 'CT3']]
# n = n.loc[:, ['CT1', 'CT2', 'CT3']]

# Optional: scale so overall-pseudobulk has mean 0 and variance 1 across individuals.
_, _, ctp, ctnu = preprocess.std(ctp, ctnu, P)

# Fit a single gene with CIGMA-HE
gene = 'GENE1'
ctp_gene = ctp[gene].unstack().to_numpy()
ctnu_gene = ctnu[gene].unstack().to_numpy()

# Kinship matrix (individuals in same order as ctp/ctnu).
# Can be computed from genotype data using GCTA or PLINK.
K = pd.read_csv(K_file, index_col=0).to_numpy()

out, pvalues = fit.free_HE(ctp_gene, K, ctnu_gene, P)

Input data

See the test script for complete input data examples and usage.

Note

This project has been set up using PyScaffold 4.4. For details and usage information on PyScaffold see https://pyscaffold.org/.