DOCSEARCHbyAlgolia
Submit Docs

cNMF (Consensus Non-negative Matrix Factorization)

repository·main·Indexed 19 days ago

https://github.com/dylkot/cnmf

A pipeline for inferring gene expression programs (GEPs) from single-cell RNA-seq data by decomposing a count matrix into GEPs and their corresponding usage per cell. It provides a Python API and CLI for a five-step workflow: prepare, factorize, combine, k_selection_plot, and consensus. The library supports Scanpy AnnData (.h5ad), tab-delimited text, and 10x-Genomics MTX formats, and includes a Preprocess class for batch correction using an adaptation of Harmony.

Tokens
12.4K
Snippets
39
Records
45
Agent score
65%

What's inside cNMF

  1. Supported Input Data Formats for cNMF

    main

    cNMF accepts three primary input formats for the count matrix:

    1. Scanpy AnnData (.h5ad): Highly recommended for large datasets. Using sparse matrices in .h5ad files significantly reduces memory usage and speeds up loading.
    2. Tab-delimited Text: A raw text file where rows are cell IDs (barcodes) and columns are gene IDs.
    3. 10x-Genomics MTX Directory: A directory containing counts.mtx (or counts.mtx.gz), barcodes.tsv, and genes.tsv. You provide the path to the .mtx file to the counts_fn parameter.

  2. Step 2: Factorize the matrix

    main

    Run the NMF factorization for the replicates allocated during the prepare step. If you specified --total-workers greater than 1, you must run the command for each worker index separately. These commands should be submitted to distinct processors or machines to run in parallel.

    Tip: For best performance when using GNU Parallel, it is recommended to use 2 workers, as scikit-learn's NMF implementation may already utilize multiple cores on a single machine.

    # Run all jobs for worker 0 (if total-workers is 3)
cnmf factorize --output-dir ./example_data --name example_cNMF --worker-index 0 --total-workers 3

# Run all jobs for worker 1
cnmf factorize --output-dir ./example_data --name example_cNMF --worker-index 1 --total-workers 3

# Run all jobs for worker 2
cnmf factorize --output-dir ./example_data --name example_cNMF --worker-index 2 --total-workers 3

  3. Install cNMF via pip

    main

    You can install the core cNMF package using pip. It requires Python 3.7 or 3.10 and depends on scikit-learn>=1.0, scanpy>=1.8, and AnnData>=0.9.

    If you intend to use the batch correction preprocessing features, you must also install harmonypy and scikit-misc.

    # Core installation
pip install cnmf

# Required for batch correction preprocessing
pip install harmonypy
pip install scikit-misc
    pip install cnmf

  4. Step 4: Select an optimal K

    main

    Use the cnmf k_selection_plot command to iterate through the tested K values and calculate the trade-off between stability and error. This generates a PNG plot in the output directory.

    There is no single definitive way to choose K, but common practice is to select the largest value that is reasonably stable or represents a local maximum in stability.

    cnmf k_selection_plot --output-dir ./example_data --name example_cNMF

  5. Step 1: Prepare the input matrix and run parameters

    main

    Use the cnmf prepare command to normalize your input matrix and set up the parameters for the factorization process. This step subsets the data to high-variance (over-dispersed) genes to improve factorization speed and signal quality. While factorization runs on a subset, the final spectra are re-fit to include all genes from the original input.

    Important Pre-processing Note: Ensure your input matrix does not include any cells or genes with 0 total counts. Filter out low-count cells and genes prior to running this command to avoid errors.

    cnmf prepare --output-dir ./example_data --name example_cNMF -c ./example_data/counts_prefiltered.txt -k 5 6 7 8 9 10 11 12 13 --n-iter 100 --seed 14 --numgenes 2000

  6. Step 5: Obtain consensus estimates for programs and usages

    main

    The final step clusters the spectra to produce consensus estimates for gene expression programs (GEP) and their usage. This step includes an optional outlier filtering process based on local density.

    Workflow Recommendation: Run the command twice. First, use --local-density-threshold 2.00 to observe the distribution of average distances in the diagnostic plot. Then, run it a second time with a smaller threshold (determined from the histogram) to filter out outliers.

    cnmf consensus --output-dir ./example_data --name example_cNMF --components 10 --local-density-threshold 0.2 --show-clustering

  7. Step 3: Combine individual spectra results

    main

    After factorization is complete, combine the individual replicate files for each K into a single merged file for each K.

    cnmf combine --output-dir ./example_data --name example_cNMF

# Optional: Clean up temporary files after combining
rm ./example_data/example_cNMF/cnmf_tmp/example_cNMF.spectra.k_*.iter_*.df.npz

  8. Integrate technical variables and batches using Preprocess

    main

    To handle batch effects before running cNMF, use the Preprocess class. This implements an adaptation of Harmony that corrects the underlying count matrix rather than principal components.

    Workflow

    1. Initialize Preprocess.
    2. Use preprocess_for_cnmf on an AnnData object, specifying harmony_vars (column names in adata.obs to correct).
    3. Pass the resulting corrected files into the cNMF.prepare() method.
    from cnmf import cNMF, Preprocess

# 1. Initialize Preprocess
p = Preprocess(random_seed=14)

# 2. Batch correct and save outputs
# This produces corrected counts (adata_c), TPM normalized data (adata_tpm), and high-variance genes (hvgs)
(adata_c, adata_tpm, hvgs) = p.preprocess_for_cnmf(
    adata, 
    harmony_vars=['Sex', 'Sample'], 
    n_top_rna_genes=2000, 
    librarysize_targetsum=1e6,
    save_output_base='./example_islets/batchcorrect_example_sex'
)

# 3. Run cNMF using the corrected files
cnmf_obj_corrected = cNMF(output_dir='./example_islets', name='BatchCorrected')
cnmf_obj_corrected.prepare(
    counts_fn='./example_islets/batchcorrect_example.Corrected.HVG.Varnorm.h5ad',
    tpm_fn='./example_islets/batchcorrect_example.TP10K.h5ad',
    genes_file='./example_islets/batchcorrect_example.Corrected.HVGs.txt',
    components=[15], 
    n_iter=20, 
    seed=14, 
    num_highvar_genes=2000
)
    from cnmf import cNMF, Preprocess

p = Preprocess(random_seed=14)

(adata_c, adata_tpm, hvgs) = p.preprocess_for_cnmf(adata, harmony_vars=['Sex', 'Sample'], n_top_rna_genes = 2000, librarysize_targetsum= 1e6,
                                                    save_output_base='./example_islets/batchcorrect_example_sex')

cnmf_obj_corrected = cNMF(output_dir='./example_islets', name='BatchCorrected')
cnmf_obj_corrected.prepare(counts_fn='./example_islets/batchcorrect_example.Corrected.HVG.Varnorm.h5ad',
                           tpm_fn='./example_islets/batchcorrect_example.TP10K.h5ad',
                           genes_file='./example_islets/batchcorrect_example.Corrected.HVGs.txt',
                           components=[15], n_iter=20, seed=14, num_highvar_genes=2000)

  9. How the cNMF workflow works

    main

    A standard cNMF analysis follows this lifecycle:

    1. prepare(...): Loads raw counts, calculates TPM, identifies high-variance genes, and saves normalized counts and NMF parameters to the output directory.
    2. factorize(worker_i, total_workers, ...): Runs multiple NMF factorizations for each $K$ in parallel. This step produces individual .df.npz files for each iteration's spectra and usages.
    3. combine(components=...): Merges the individual iteration files for a specific $K$ into a single merged_spectra file.
    4. consensus(k, ...): Computes the consensus spectra and usages, providing a stable estimate of the gene programs and their cell-wise usage.
    # Typical Workflow
cnmf = cNMF(output_dir="./results", name="exp1")
cnmf.prepare(counts_fn="counts.h5ad", components=[10, 20])

# In a parallel environment:
cnmf.factorize(worker_i=0, total_workers=2)
cnmf.factorize(worker_i=1, total_workers=2)

cnmf.combine(components=[10, 20])
cnmf.consensus(k=20)

  10. Run cNMF in parallel using UGER

    main

    To run factorizations in parallel across multiple workers using a UGER (Univa Grid Engine) scheduler, you can use a submission command that assigns tasks based on a --worker-index.

    Example command for 5 workers:

    qsub -cwd -b y -l h_vmem=2g,h_rt=3:00:00 -o ./log -e ./log -N cnmf -t 1-5 'python ../cnmf.py factorize --output-dir ./simulated_example_data --name example_cNMF --worker-index $SGE_TASK_ID'

    Note: You must first run the prepare step with the total number of workers specified to ensure tasks are correctly allocated.

  11. Use the cNMF Python Class Interface

    main

    The cNMF class provides a high-level Python interface for performing consensus Non-negative Matrix Factorization. The typical workflow involves initializing the object, preparing the data, running factorizations, combining results, and performing consensus analysis.

    Workflow Steps:

    1. Initialize: Create a cNMF object with an output_dir and a name.
    2. Prepare: Use .prepare() to normalize counts and select high-variance genes.
    3. Factorize: Use .factorize() to run the NMF iterations.
    4. Combine: Use .combine() to merge replicate spectra.
    5. Consensus: Use .consensus() to obtain the final consensus matrix and clusterings.
    6. Analyze: Use .load_results() to retrieve processed data for plotting.
    from cnmf import cNMF

# 1. Initialize
cnmf_obj = cNMF(output_dir='./simulated_example_data', name='example_cNMF')

# 2. Prepare
cnmf_obj.prepare(counts_fn='path/to/counts.txt', components=[5, 6, 7], n_iter=20, seed=14, num_highvar_genes=1500)

# 3. Factorize
cnmf_obj.factorize()

# 4. Combine
cnmf_obj.combine()

# 5. Consensus
cnmf_obj.consensus(k=7, density_threshold=0.1, show_clustering=True)

# 6. Load Results
normalized_usage, gep_scores, gep_tpm, top_genes = cnmf_obj.load_results(K=7, density_threshold=0.1)

  12. Download example simulation data

    main

    To follow the simulation guide, download the example data package using wget and extract it with tar:

    ! wget -O ./example_simulated_data.tar.gz https://storage.googleapis.com/sabeti-public/dkotliar/cNMF/example_data_20191024.tar.gz
! tar -zxvf ./example_simulated_data.tar.gz && rm ./example_simulated_data.tar.gz