pycisTopic

https://pycistopic.readthedocs.io/en/latest/
Admin
pycisTopic is a Python module for simultaneously identifying cell states and cis-regulatory topics...
Tokens:72,390
Snippets:453
Trust Score:9.6
Update:5 months ago
Show doc for...
Tokens:
Raw
### Install pycisTopic from Source

Source: https://pycistopic.readthedocs.io/en/latest/installation.html

Clone the pycisTopic repository and install it using pip. This method is suitable for development or when the latest version is required.

```bash
git clone https://github.com/aertslab/pycisTopic.git
cd pycisTopic
pip install -e .
```

--------------------------------

### Check pycisTopic Version

Source: https://pycistopic.readthedocs.io/en/latest/installation.html

Import the pycisTopic library and access its __version__ attribute to verify the installed version.

```python
import pycisTopic
pycisTopic.__version__
```

--------------------------------

### Import pycisTopic Clustering and Visualization Functions

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Imports necessary functions for clustering, dimensionality reduction, and plotting from the pycisTopic library. Ensure pycisTopic is installed.

```python
from pycisTopic.clust_vis import (
    find_clusters,
    run_umap,
    run_tsne,
    plot_metadata,
    plot_topic,
    cell_topic_heatmap
)
```

--------------------------------

### peak_calling

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Performs pseudobulk peak calling with MACS2. Requires MACS2 installation.

```APIDOC
## peak_calling

### Description
Performs pseudobulk peak calling with MACS2. It requires to have MACS2 installed.

### Method
`pycisTopic.pseudobulk_peak_calling.peak_calling`

### Parameters
#### Path Parameters
- **macs_path** (str) - Path to MACS binary (e.g. /xxx/MACS/xxx/bin/macs2).
- **bed_paths** (dict) - A dictionary of character strings, with sample name as names indicating the path to the fragments file/s from which pseudobulk profiles have to be created.
- **outdir** (str) - Path to the output directory.
- **genome_size** (str) - Effective genome size which is defined as the genome size which can be sequenced. Possible values: ‘hs’, ‘mm’, ‘ce’ and ‘dm’.

#### Query Parameters
- **n_cpu** (int, optional) - Number of cores to use. Default: 1.
- **input_format** (str, optional) - Format of tag file can be ELAND, BED, ELANDMULTI, ELANDEXPORT, SAM, BAM, BOWTIE, BAMPE, or BEDPE. Default is AUTO which will allow MACS to decide the format automatically. Default: ‘BEDPE’.
- **shift** (int, optional) - To set an arbitrary shift in bp. For finding enriched cutting sites (such as in ATAC-seq) a shift of 73 bp is recommended. Default: 73.
- **ext_size** (int, optional) - To extend reads in 5’->3’ direction to fix-sized fragment. For ATAC-seq data, a extension of 146 bp is recommended. Default: 146.
- **keep_dup** (str, optional) - Whether to keep duplicate tags at te exact same location. Default: ‘all’.
- **q_value** (float, optional) - The q-value (minimum FDR) cutoff to call significant regions. Default: 0.05.
- **nolambda** (bool, optional) - Do not consider the local bias/lambda at peak candidate regions.
- **skip_empty_peaks** (bool, optional) - Whether to skip empty peaks.

### Request Example
```python
peak_calling(
    macs_path="/path/to/macs2/macs2",
    bed_paths={
        "sample1": "/path/to/sample1_fragments.bed",
        "sample2": "/path/to/sample2_fragments.bed"
    },
    outdir="/path/to/output",
    genome_size="hs"
)
```
```

--------------------------------

### Load and Prepare Chromosome Sizes

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Loads chromosome sizes from a TSV file and prepares it for PyRanges, setting 'Chromosome' and 'End' columns and adding a 'Start' column with zeros.

```python
chromsizes = pd.read_table(os.path.join(out_dir, "qc", "hg38.chrom_sizes_and_alias.tsv"))
chromsizes
chromsizes.rename({"# ucsc": "Chromosome", "length": "End"}, axis = 1, inplace = True)
chromsizes["Start"] = 0
chromsizes = pr.PyRanges(chromsizes[["Chromosome", "Start", "End"]])
chromsizes
```

--------------------------------

### MACS2 Peak Calling Command Examples

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

These logs show the specific MACS2 commands executed by pycisTopic for different cell types. They detail the input treatment files, output directories, and various MACS2 parameters used for peak calling.

```bash
(macs_call_peak_ray pid=1232492) 2024-02-29 18:05:51,892 cisTopic     INFO     Calling peaks for AST_CER with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/AST_CER.fragments.tsv.gz --name AST_CER  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232485) 2024-02-29 18:05:51,919 cisTopic     INFO     Calling peaks for GC with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/GC.fragments.tsv.gz --name GC  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232488) 2024-02-29 18:05:51,975 cisTopic     INFO     Calling peaks for OPC with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/OPC.fragments.tsv.gz --name OPC  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232483) 2024-02-29 18:05:51,936 cisTopic     INFO     Calling peaks for ASTP with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/ASTP.fragments.tsv.gz --name ASTP  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232484) 2024-02-29 18:05:51,888 cisTopic     INFO     Calling peaks for INH_PVALB with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/INH_PVALB.fragments.tsv.gz --name INH_PVALB  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232489) 2024-02-29 18:05:51,895 cisTopic     INFO     Calling peaks for MGL with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/MGL.fragments.tsv.gz --name MGL  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232490) 2024-02-29 18:05:51,975 cisTopic     INFO     Calling peaks for INH_VIP with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/INH_VIP.fragments.tsv.gz --name INH_VIP  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232491) 2024-02-29 18:05:51,896 cisTopic     INFO     Calling peaks for MOL_A with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/MOL_A.fragments.tsv.gz --name MOL_A  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232486) 2024-02-29 18:05:51,961 cisTopic     INFO     Calling peaks for PURK with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/PURK.fragments.tsv.gz --name PURK  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

```bash
(macs_call_peak_ray pid=1232487) 2024-02-29 18:05:51,915 cisTopic     INFO     Calling peaks for MOL_B with macs2 callpeak --treatment outs/consensus_peak_calling/pseudobulk_bed_files/MOL_B.fragments.tsv.gz --name MOL_B  --outdir outs/consensus_peak_calling/MACS --format BEDPE --gsize hs --qvalue 0.05 --nomodel --shift 73 --extsize 146 --keep-dup all --call-summits --nolambda

```

--------------------------------

### Download and Format Chromosome Sizes

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Downloads chromosome sizes for a specified genome (hg38) from UCSC and formats it into a pandas DataFrame with 'Chromosome', 'Start', and 'End' columns.

```python
chromsizes = pd.read_table(
    "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes",
    header = None,
    names = ["Chromosome", "End"]
)
chromsizes.insert(1, "Start", 0)
chromsizes.head()
```

--------------------------------

### MALLET LDA Training Log

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Example log output from the MALLET LDA training process, showing corpus formatting, dictionary building, and topic training for a specific number of topics.

```log
2024-03-06 10:52:57,071 cisTopic     INFO     Formatting input to corpus
2024-03-06 10:52:57,690 gensim.corpora.dictionary INFO     adding document #0 to Dictionary<0 unique tokens: []>
2024-03-06 10:53:23,393 gensim.corpora.dictionary INFO     built Dictionary<436234 unique tokens: ['0', '1', '2', '3', '4']...> from 2853 documents (total 46261622 corpus positions)
2024-03-06 10:53:23,395 cisTopic     INFO     Running model with 2 topics
2024-03-06 10:53:23,409 LDAMalletWrapper INFO     Serializing temporary corpus to /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt
2024-03-06 10:54:49,354 LDAMalletWrapper INFO     Converting temporary corpus to MALLET format with Mallet-202108/bin/mallet import-file --preserve-case --keep-sequence --remove-stopwords --token-regex "\S+" --input /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt --output /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.mallet
2024-03-06 10:55:16,373 LDAMalletWrapper INFO     Training MALLET LDA with Mallet-202108/bin/mallet train-topics --input /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.mallet --num-topics 2  --alpha 50 --beta 0.1 --optimize-interval 0 --num-threads 12 --output-state /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/607c9d_state.mallet.gz --output-doc-topics /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/607c9d_doctopics.txt --output-topic-keys /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/607c9d_topickeys.txt --num-iterations 500 --inferencer-filename /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/607c9d_inferencer.mallet --doc-topics-threshold 0.0  --random-seed 555
2024-03-06 11:11:23,264 LDAMalletWrapper INFO     loading assigned topics from /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/607c9d_state.mallet.gz
2024-03-06 11:12:40,862 cisTopic     INFO     Model with 2 topics done!
2024-03-06 11:12:40,863 cisTopic     INFO     Saving model with 2 topics at /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial
2024-03-06 11:12:40,954 cisTopic     INFO     Running model with 5 topics
2024-03-06 11:12:40,967 LDAMalletWrapper INFO     Serializing temporary corpus to /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt
2024-03-06 11:14:09,509 LDAMalletWrapper INFO     Converting temporary corpus to MALLET format with Mallet-202108/bin/mallet import-file --preserve-case --keep-sequence --remove-stopwords --token-regex "\S+" --input /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt --output /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.mallet
2024-03-06 11:14:37,009 LDAMalletWrapper INFO     Training MALLET LDA with Mallet-202108/bin/mallet train-topics --input /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.mallet --num-topics 5  --alpha 50 --beta 0.1 --optimize-interval 0 --num-threads 12 --output-state /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/3efb70_state.mallet.gz --output-doc-topics /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/3efb70_doctopics.txt --output-topic-keys /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/3efb70_topickeys.txt --num-iterations 500 --inferencer-filename /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/3efb70_inferencer.mallet --doc-topics-threshold 0.0  --random-seed 555
2024-03-06 11:31:43,881 LDAMalletWrapper INFO     loading assigned topics from /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/3efb70_state.mallet.gz
2024-03-06 11:33:02,089 cisTopic     INFO     Model with 5 topics done!
2024-03-06 11:33:02,090 cisTopic     INFO     Saving model with 5 topics at /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial
2024-03-06 11:33:02,210 cisTopic     INFO     Running model with 10 topics
2024-03-06 11:33:02,224 LDAMalletWrapper INFO     Serializing temporary corpus to /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt
2024-03-06 11:34:30,049 LDAMalletWrapper INFO     Converting temporary corpus to MALLET format with Mallet-202108/bin/mallet import-file --preserve-case --keep-sequence --remove-stopwords --token-regex "\S+" --input /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.txt --output /scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial/corpus.mallet
```

--------------------------------

### Import necessary pycisTopic functions

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Imports the required functions for differential feature analysis from the pycisTopic library. Ensure pycisTopic is installed.

```python
from pycisTopic.diff_features import (
    impute_accessibility,
    normalize_scores,
    find_highly_variable_features,
    find_diff_features
)
import numpy as np
```

--------------------------------

### Get Insert Size Distribution

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Calculates the distribution of fragment insert sizes. Requires fragments that have already been filtered by cell barcodes.

```python
>>> insert_size_dist_df_pl = get_insert_size_distribution(
...     fragments_df_pl=fragments_cb_filtered_df_pl,
... )
```

--------------------------------

### Get Barcodes Passing QC for Sample

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Retrieves barcodes that pass quality control thresholds for a given sample. Allows for automatic thresholding or manual specification of thresholds for unique fragments, TSS enrichment, and FRIP. Use 'use_automatic_thresholds = True' to let pycisTopic determine thresholds.

```python
from pycisTopic.qc import get_barcodes_passing_qc_for_sample
sample_id_to_barcodes_passing_filters = {}
sample_id_to_thresholds = {}
for sample_id in fragments_dict:
    (
        sample_id_to_barcodes_passing_filters[sample_id],
        sample_id_to_thresholds[sample_id]
    ) = get_barcodes_passing_qc_for_sample(
            sample_id = sample_id,
            pycistopic_qc_output_dir = "outs/qc",
            unique_fragments_threshold = None, # use automatic thresholding
            tss_enrichment_threshold = None, # use automatic thresholding
            frip_threshold = 0,
            use_automatic_thresholds = True,
    )
```

--------------------------------

### Get TSS Annotation from Ensembl BioMart

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieve Transcription Start Site (TSS) annotations from Ensembl BioMart for specified transcript types. Supports caching and custom BioMart hosts.

```python
tss_annotation_bed_df_pl = get_tss_annotation_from_ensembl(
    biomart_name="hsapiens_gene_ensembl"
)
```

```python
tss_annotation_jul2022_bed_df_pl = get_tss_annotation_from_ensembl(
    biomart_name="hsapiens_gene_ensembl",
    biomart_host="http://jul2022.archive.ensembl.org/",
)
```

--------------------------------

### Build Docker and Singularity Images

Source: https://pycistopic.readthedocs.io/en/latest/installation.html

Steps to clone repositories, build a Docker image using Podman, export it to OCI archive, and then build a Singularity image from the archive. This is for creating containerized environments.

```bash
# Clone repositories (pycisTopic and pycistarget)
git clone https://github.com/aertslab/pycisTopic.git
git clone https://github.com/aertslab/pycistarget.git

# Build image
podman build -t aertslab/pycistopic:latest . -f pycisTopic/Dockerfile

# Export to oci
podman save --format oci-archive --output pycistopic_img.tar localhost/aertslab/pycistopic

# Build to singularity
singularity build pycistopic.sif oci-archive://pycistopic_img.tar

# Add all binding paths where you would need to access
singularity exec -B /lustre1,/staging,/data,/vsc-hard-mounts,/scratch pycistopic.sif ipython3
```

--------------------------------

### Create QC Command Lines File

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Generates a text file listing all pycistopic QC command lines. This is useful for reproducibility and auditing.

```python
with open("qc_commands.txt", "w") as f:
    for sample_id, thresholds in sample_id_to_thresholds.items():
        f.write(f"pycistopic qc --sample-id {sample_id} --thresholds {thresholds}\n")
```

--------------------------------

### get_tss_matrix

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieves the Transcription Start Site (TSS) matrix.

```APIDOC
## get_tss_matrix

### Description
Get TSS matrix.

### Function Signature
`pycisTopic.utils.get_tss_matrix(_fragments_ , _flank_window_ , _tss_space_annotation_)`
```

--------------------------------

### Download Tutorial Data

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Downloads necessary data files for the tutorial, including ATAC fragments and cell metadata. Creates a 'data' directory if it doesn't exist.

```bash
!mkdir -p data
!wget -O data/fragments.tsv.gz https://cf.10xgenomics.com/samples/cell-arc/1.0.0/human_brain_3k/human_brain_3k_atac_fragments.tsv.gz
!wget -O data/fragments.tsv.gz.tbi https://cf.10xgenomics.com/samples/cell-arc/1.0.0/human_brain_3k/human_brain_3k_atac_fragments.tsv.gz.tbi
!wget -O data/cell_data.tsv https://raw.githubusercontent.com/aertslab/pycisTopic/polars/data/cell_data_human_cerebellum.tsv
```

--------------------------------

### Generate pycisTopic QC Command Lines

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Creates a text file containing pycisTopic QC command lines for each sample, which can then be executed in parallel. Ensure the necessary filenames and output directory are defined.

```python
with open(pycistopic_qc_commands_filename, "w") as fh:
    for sample, fragment_filename in fragments_dict.items():
        print(
            "pycistopic qc",
            f"--fragments {fragment_filename}",
            f"--regions {regions_bed_filename}",
            f"--tss {tss_bed_filename}",
            f"--output {os.path.join(out_dir, "qc")}/{sample}",
            sep=" ",
            file=fh,
        )
```

--------------------------------

### macs_call_peak

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Performs pseudobulk peak calling with MACS2 in a group. Requires MACS2 installation.

```APIDOC
## macs_call_peak

### Description
Performs pseudobulk peak calling with MACS2 in a group. It requires to have MACS2 installed.

### Method
`pycisTopic.pseudobulk_peak_calling.macs_call_peak`

### Parameters
#### Path Parameters
- **macs_path** (str) - Path to MACS binary (e.g. /xxx/MACS/xxx/bin/macs2).
- **bed_path** (str) - Path to fragments file bed file.
- **name** (str) - Name of string of the group.
- **outdir** (str) - Path to the output directory.
- **genome_size** (str) - Effective genome size which is defined as the genome size which can be sequenced. Possible values: ‘hs’, ‘mm’, ‘ce’ and ‘dm’.

#### Query Parameters
- **input_format** (str, optional) - Format of tag file can be ELAND, BED, ELANDMULTI, ELANDEXPORT, SAM, BAM, BOWTIE, BAMPE, or BEDPE. Default is AUTO which will allow MACS to decide the format automatically. Default: ‘BEDPE’.
- **shift** (int, optional) - To set an arbitrary shift in bp. For finding enriched cutting sites (such as in ATAC-seq) a shift of 73 bp is recommended. Default: 73.
- **ext_size** (int, optional) - To extend reads in 5’->3’ direction to fix-sized fragment. For ATAC-seq data, a extension of 146 bp is recommended. Default: 146.
- **keep_dup** (str, optional) - Whether to keep duplicate tags at te exact same location. Default: ‘all’.
- **q_value** (float, optional) - The q-value (minimum FDR) cutoff to call significant regions. Default: 0.05.
- **nolambda** (bool, optional) - Do not consider the local bias/lambda at peak candidate regions.
- **skip_empty_peaks** (bool, optional) - Whether to skip empty peaks.

### Request Example
```python
macs_call_peak(
    macs_path="/path/to/macs2/macs2",
    bed_path="/path/to/fragments.bed",
    name="group1",
    outdir="/path/to/output",
    genome_size="hs"
)
```
```

--------------------------------

### Prepare QC Command File

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Defines file paths for consensus regions and TSS annotations, and specifies the output filename for QC commands. This is a precursor to potentially running QC in parallel.

```python
regions_bed_filename = os.path.join(out_dir, "consensus_peak_calling/consensus_regions.bed")
tss_bed_filename = os.path.join(out_dir, "qc", "tss.bed")

pycistopic_qc_commands_filename = "pycistopic_qc_commands.txt"
```

--------------------------------

### Initialize pycisTopic Object and Load R Model in Python

Source: https://pycistopic.readthedocs.io/en/latest/faqs.html

Initialize a cisTopic object in Python using a fragment matrix and add a pre-trained LDA model exported from R. Ensure the matrix path and the folder containing R data are correctly specified.

```python
import pycisTopic
## 1. Initialize cisTopic object
from pycisTopic.cistopic_class import *
# Load count matrix
matrix_path=PATH_TO_FRAGMENTS_MATRIX
fragment_matrix = pd.read_feather(matrix_path)
cisTopic_obj = create_cistopic_object(fragment_matrix)
# Also add the cell annotation, cell_data should be a pandas df with cells as rows (cell names as index) and variables as columns
cisTopic_obj.add_cell_data(cell_data)
## 2. Add model
from pycisTopic.utils import *
model = load_cisTopic_model(PATH_TO_THE_FOLDER_WHERE_YOU_SAVED_DATA_FROM_R)
cistopic_obj.add_LDA_model(model)
# You can continue with the rest of the pipeline.

```

--------------------------------

### Get Ensembl TSS Annotation

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieves TSS annotation data from Ensembl BioMart. Ensure the biomart_name is correctly specified for your organism.

```python
>>> ensembl_tss_annotation_bed_df_pl = get_tss_annotation_from_ensembl(
...     biomart_name="hsapiens_gene_ensembl"
... )
```

--------------------------------

### Configure MALLET Environment and Import Function

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Set the MALLET_MEMORY environment variable to allocate sufficient memory for MALLET. Then, import the necessary function for running LDA models with MALLET and configure the path to the MALLET binary.

```python
os.environ['MALLET_MEMORY'] = '200G'
from pycisTopic.lda_models import run_cgs_models_mallet
# Configure path Mallet
mallet_path="Mallet-202108/bin/mallet"
```

--------------------------------

### Read TSS Annotation from BED File

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Reads TSS (Transcription Start Site) annotation from a BED file. Refer to pycisTopic.gene_annotation.read_tss_annotation_from_bed() for more details.

```python
tss_annotation_bed_df_pl = read_tss_annotation_from_bed(
    tss_annotation_bed_filename="hg38.tss.bed",
)
```

--------------------------------

### Create Output Directory

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Creates an output directory named 'outs' if it does not already exist. This directory will be used to store pycisTopic results.

```python
import os
out_dir = "outs"
os.makedirs(out_dir, exist_ok = True)
```

--------------------------------

### get_tss_annotation_from_ensembl

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieves Transcription Start Site (TSS) annotation from Ensembl BioMart for specified transcript types. Returns data in BED format.

```APIDOC
## get_tss_annotation_from_ensembl

### Description
Get TSS annotation for requested transcript types from Ensembl BioMart.

### Parameters
#### Path Parameters
- **biomart_name** (string) - Required - Ensembl BioMart ID of the dataset (e.g., `hsapiens_gene_ensembl`, `mmusculus_gene_ensembl`).
- **biomart_host** (string) - Optional - BioMart host URL to use. Defaults to `http://www.ensembl.org`.
- **transcript_type** (list[string] | None) - Optional - Filter by a list of specified transcript types (e.g., `["protein_coding"]`) or `None` for all.
- **use_cache** (boolean) - Optional - Whether to cache requests to Ensembl BioMart server. Defaults to `True`.

### Returns
Polars DataFrame with TSS positions in BED format.

### Request Example
```python
tss_annotation_bed_df_pl = get_tss_annotation_from_ensembl(
    biomart_name="hsapiens_gene_ensembl"
)
tss_annotation_jul2022_bed_df_pl = get_tss_annotation_from_ensembl(
    biomart_name="hsapiens_gene_ensembl",
    biomart_host="http://jul2022.archive.ensembl.org/",
)
```
```

--------------------------------

### Download and Extract MALLET Binary

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Use these shell commands to download the MALLET binary archive and extract its contents. This is a prerequisite for using MALLET for parallel LDA.

```shell
!wget https://github.com/mimno/Mallet/releases/download/v202108/Mallet-202108-bin.tar.gz
!tar -xf Mallet-202108-bin.tar.gz
```

--------------------------------

### Create Output Directories

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Creates necessary subdirectories within the output directory for storing region sets. Ensures that directories exist before attempting to save files.

```python
os.makedirs(os.path.join(out_dir, "region_sets"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "region_sets", "Topics_otsu"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "region_sets", "Topics_top_3k"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "region_sets", "DARs_cell_type"), exist_ok = True)
```

--------------------------------

### CistopicObject Initialization

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Initializes a CistopicObject with fragment matrices, cell/region names, metadata, and path information.

```APIDOC
## CistopicObject

### Description
`CistopicObject` contains the cell by fragment matrices (stored as counts `fragment_matrix` and as binary accessibility `binary_matrix`), cell metadata `cell_data`, region metadata `region_data` and path/s to the fragments file/s `path_to_fragments`. LDA models from `CisTopicLDAModel` can be stored `selected_model` as well as cell/region projections `projections` as a dictionary.

### Parameters
- **_fragment_matrix**: csr_matrix - The cell by fragment counts matrix.
- **_binary_matrix**: csr_matrix - The cell by fragment binary accessibility matrix.
- **_cell_names**: list[str] - A list of cell names.
- **_region_names**: list[str] - A list of region names.
- **_cell_data**: DataFrame - A DataFrame containing cell metadata.
- **_region_data**: DataFrame - A DataFrame containing region metadata.
- **_path_to_fragments**: str or dict - Path(s) to the fragments file(s).
- **_project**: str, optional - Name of the cisTopic project. Defaults to 'cisTopic'.
```

--------------------------------

### Get Gene Annotation from Ensembl BioMart

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Fetches gene annotation data from Ensembl BioMart for a given organism. Ensure the biomart_name is correctly specified.

```python
>>> hg38_tss_annotation_bed_df_pl = get_tss_annotation_from_ensembl(
...     biomart_name="hsapiens_gene_ensembl",
... )
>>> hg38_tss_annotation_bed_df_pl
```

--------------------------------

### pycisTopic.clust_vis.run_tsne

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Runs tSNE and adds it to the dimensionality reduction dictionary. This function can utilize FItSNE if installed, otherwise it falls back to scikit-learn's t-SNE implementation.

```APIDOC
## pycisTopic.clust_vis.run_tsne

### Description
Run tSNE and add it to the dimensionality reduction dictionary. If FItSNE is installed it will be used, otherwise sklearn TSNE implementation will be used.

### Parameters
#### Path Parameters
- **cistopic_obj** (`class::CistopicObject`) - Required - A cisTopic object with a model in class::CistopicObject.selected_model.

#### Query Parameters
- **target** (str) - Optional - Whether cells (‘cell’) or regions (‘region’) should be used. Default: ‘cell’
- **scale** (bool) - Optional - Whether to scale the cell-topic or topic-regions contributions prior to the dimensionality reduction. Default: False
- **reduction_name** (str) - Optional - Reduction name to use as key in the dimensionality reduction dictionary. Default: ‘tSNE’
- **random_state** (int) - Optional - Seed parameter for running tSNE. Default: 555
- **perplexity** (int) - Optional - Perplexity parameter for FitSNE. Default: 30
- **selected_topics** (list[int], optional) - Optional - A list with selected topics to be used for clustering. Default: None (use all topics)
- **selected_features** (list[str], optional) - Optional - A list with selected features (cells or regions) to cluster. This is recommended when working with regions (e.g. selecting regions in binarized topics), as working with all regions can be time consuming. Default: None (use all features)
- **harmony** (bool) - Optional -
- **rna_components** (pd.DataFrame | None) - Optional -
- **rna_weight** (float) - Optional - Default: 0.5
```

--------------------------------

### Create output directory for loom files

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Ensures that the directory for saving loom files exists, creating it if necessary.

```python
os.makedirs(os.path.join(out_dir, "loom"), exist_ok=True)
```

--------------------------------

### Get TSS Annotation from Ensembl

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Fetches TSS annotations from Ensembl using the specified biomart name. This is the first step before writing to a BED file or further processing.

```python
tss_annotation_bed_df_pl = get_tss_annotation_from_ensembl(
    biomart_name="hsapiens_gene_ensembl"
)
```

--------------------------------

### Binarize topic-region distributions using Otsu's method

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Binarizes topic-region distributions using Otsu's method, which is effective for many topic-region distributions. The results are plotted for visualization.

```python
region_bin_topics_otsu = binarize_topics(
    cistopic_obj, method='otsu',
    plot=True, num_columns=5
)
```

--------------------------------

### Run cisTopic Models with MALLET

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Initiates cisTopic models using the run_cgs_models_mallet function. Specify the cisTopic object, a list of topic numbers to test, CPU count, iterations, random state, and alpha/eta parameters. Temporary and save paths for MALLET outputs are also configured.

```python
models=run_cgs_models_mallet(
    cistopic_obj,
    n_topics=[2, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50],
    n_cpu=12,
    n_iter=500,
    random_state=555,
    alpha=50,
    alpha_by_topic=True,
    eta=0.1,
    eta_by_topic=False,
    tmp_path="/scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial",
    save_path="/scratch/leuven/330/vsc33053/ray_spill/mallet/tutorial",
    mallet_path=mallet_path,
)
```

--------------------------------

### Get NCBI Assembly Accessions for Species

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Fetch NCBI assembly accession numbers and names for a given species. Useful for mapping species to their corresponding assembly identifiers.

```python
print(get_ncbi_assembly_accessions_for_species("homo sapiens"))
```

```python
print(get_ncbi_assembly_accessions_for_species("drosophila melanogaster"))
```

--------------------------------

### Get Chromosome Sizes and Alias Mapping from UCSC

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieves chromosome sizes and alias mappings for a specified UCSC assembly. This is often a prerequisite for other gene annotation tasks.

```python
>>> chrom_sizes_and_alias_hg38_df_pl = get_chrom_sizes_and_alias_mapping_from_ucsc(ucsc_assembly="hg38")
```

--------------------------------

### Import topic QC and plotting functions

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Imports functions for computing topic quality control metrics, plotting these metrics, and performing topic annotation, along with plotting utilities.

```python
from pycisTopic.topic_qc import compute_topic_metrics, plot_topic_qc, topic_annotation
import matplotlib.pyplot as plt
from pycisTopic.utils import fig2img

```

--------------------------------

### Get Fragments in Peaks

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Calculates the total and unique fragment counts within specified genomic regions (peaks). Requires pre-filtered fragments and peak region data.

```python
>>> fragments_in_peaks_df_pl = get_fragments_in_peaks(
...     fragments_df_pl=fragments_cb_filtered_df_pl,
...     regions_df_pl=regions_df_pl,
... )
```

--------------------------------

### Define Fragment Files Dictionary

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Defines a dictionary mapping sample IDs to their corresponding fragment file paths. pycisTopic automatically handles barcode collisions between samples.

```python
fragments_dict = {
    "10x_multiome_brain": "data/fragments.tsv.gz"
}
```

--------------------------------

### Get TSS Annotation

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Downloads TSS annotations from Ensembl BioMart for a specified species ('hsapiens_gene_ensembl') and UCSC genome build ('hg38'). Outputs the annotation to a BED file.

```bash
!mkdir -p outs/qc
!pycistopic tss get_tss \
    --output outs/qc/tss.bed \
    --name "hsapiens_gene_ensembl" \
    --to-chrom-source ucsc \
    --ucsc hg38
```

--------------------------------

### Intersect Genomic Regions (Default)

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Calculates the intersection of two sets of genomic regions. This example shows the default output format including information from both region sets.

```python
>>> intersection(
...     regions1_df_pl,
...     regions2_df_pl,
...     regions1_coord=True,
...     regions2_coord=True,
... )
shape: (3, 10)
┌────────────┬───────┬─────┬──────────────┬─────────┬───────┬──────────────┬─────────┬───────┬─────┐
│ Chromosome ┆ Start ┆ End ┆ Chromosome@1 ┆ Start@1 ┆ End@1 ┆ Chromosome@2 ┆ Start@2 ┆ End@2 ┆ ID  │
│ ---        ┆ ---   ┆ --- ┆ ---          ┆ ---     ┆ ---   ┆ ---          ┆ ---     ┆ ---   ┆ --- │
│ str        ┆ i64   ┆ i64 ┆ str          ┆ i64     ┆ i64   ┆ str          ┆ i64     ┆ i64   ┆ str │
╞════════════╪═══════╪═════╪══════════════╪═════════╪═══════╪══════════════╪═════════╪═══════╪═════╡
│ chr1       ┆ 2     ┆ 3   ┆ chr1         ┆ 1       ┆ 3     ┆ chr1         ┆ 2       ┆ 9     ┆ a   │
├╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌┤
│ chr1       ┆ 2     ┆ 3   ┆ chr1         ┆ 1       ┆ 3     ┆ chr1         ┆ 2       ┆ 3     ┆ a   │
├╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌┼╌╌╌╌╌┤
│ chr1       ┆ 4     ┆ 9   ┆ chr1         ┆ 4       ┆ 9     ┆ chr1         ┆ 2       ┆ 9     ┆ b   │
└────────────┴───────┴─────┴──────────────┴─────────┴───────┴──────────────┴─────────┴───────┴─────┘
```

--------------------------------

### Import binarize_topics function

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Imports the necessary function for topic binarization from the pycisTopic library.

```python
from pycisTopic.topic_binarization import binarize_topics
```

--------------------------------

### Get Fragments Per Cell Barcode

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Computes the number of fragments and duplication ratio for each cell barcode. Optionally collapses duplicate fragments and filters by minimum fragment count.

```python
>>> fragments_stats_per_cb_df_pl = get_fragments_per_cb(
...     fragments_df_pl=fragments_df_pl,
...     min_fragments_per_cb=10,
...     collapse_duplicates=True,
... )
```

--------------------------------

### Get Consensus Peaks

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Retrieves consensus peaks from narrow peak dictionaries, applying a specified half-width and blacklisting regions. Requires narrow_peak_dict, chromsizes, and a path to a blacklist file.

```python
peak_half_width=250
path_to_blacklist="pycisTopic/blacklist/hg38-blacklist.v2.bed"
consensus_peaks = get_consensus_peaks(
    narrow_peaks_dict = narrow_peak_dict,
    peak_half_width = peak_half_width,
    chromsizes = chromsizes,
    path_to_blacklist = path_to_blacklist)
```

--------------------------------

### Import pycisTopic Gene Activity Module

Source: https://pycistopic.readthedocs.io/en/latest/notebooks/human_cerebellum.html

Imports the necessary pyranges and get_gene_activity functions for gene activity inference.

```python
import pyranges as pr
from pycisTopic.gene_activity import get_gene_activity
```

--------------------------------

### Get Chromosome Sizes and Alias Mapping from UCSC

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Retrieve chromosome sizes and alias mappings from the UCSC Genome Browser for specified assemblies. Can also write the mapping to a TSV file.

```python
chrom_sizes_and_alias_hg38_df_pl = get_chrom_sizes_and_alias_mapping_from_ucsc(
    ucsc_assembly="hg38"
)
chrom_sizes_and_alias_mm10_df_pl = get_chrom_sizes_and_alias_mapping_from_ucsc(
    ucsc_assembly="mm10"
)
chrom_sizes_and_alias_dm6_df_pl = get_chrom_sizes_and_alias_mapping_from_ucsc(
    ucsc_assembly="dm6"
)
```

```python
chrom_sizes_and_alias_hg38_df_pl = get_chrom_sizes_and_alias_mapping_from_ucsc(
    ucsc_assembly="hg38",
    chrom_sizes_and_alias_tsv_filename="hg38.chrom_sizes_and_alias.tsv",
)
```

--------------------------------

### Get Ensembl BioMart Dataset Names

Source: https://pycistopic.readthedocs.io/en/latest/api.html

Fetch all available gene annotation dataset names from Ensembl BioMart. Supports specifying a custom BioMart host for the latest or archived versions.

```python
>>> biomart_latest_datasets = get_all_biomart_ensembl_dataset_names(
...    biomart_host="http://www.ensembl.org",
... )
>>> biomart_jul2022_datasets = get_all_biomart_ensembl_dataset_names(
...     biomart_host="http://jul2022.archive.ensembl.org/",
... )
```
Tokens:
Raw