### 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/", ... ) ```