ChromBPNet

repository·master·Indexed 42 days ago

https://github.com/kundajelab/chrombpnet

Bias-factorized, base-resolution deep learning models for chromatin accessibility (ATAC-seq and DNASE-seq). It includes tools for calculating sequence contribution scores via chrombpnet_deepshap, generating bigwig tracks, calculating marginal footprints, and performing de-novo motif discovery using a MODISCO pipeline.

Tokens
15.1K
Snippets
38
Records
52
Agent score
54%

What's inside ChromBPNet

  1. How ChromBPNet handles assay bias

    master

    ChromBPNet addresses confounding enzymatic bias (from DNASE-I or Tn5) through a two-step bias factorization process:

    1. Frozen Bias Model: A simple model is first learned on the chromatin background to capture the enzyme effect.
    2. TF Model: The effect of the enzyme is regressed out from the profiles. This ensures that the subsequent ChromBPNet sequence component (the TF Model) learns to predict the biological signal without being confounded by enzymatic bias.

    This architecture uses a fully convolutional neural network with dilated convolutions and residual connections to maintain large receptive fields efficiently.

  2. Build a PWM matrix from a bigwig file

    master

    Use the build_pwm_from_bigwig.py script to build a Position Weight Matrix (PWM) centered at the non-zero entries in a bigwig file. This is used in ChromBPNet to ensure bigwigs are shifted correctly.

    Arguments

    FlagArgumentDescription
    -i[bigwig]Path to ATAC/DNASE data in bigwig format
    -g[genome]Path to reference genome fasta
    -o[output_prefix]Output prefix for the image file. If a directory path is included, it must already exist. The script appends .png to this value
    -c[chr]Chromosome name (string). Must be present in both the bigwig file and column one of the chrom_sizes file
    -cz[chrom_sizes]Path to a TSV file containing chromosome names in the first column and their sizes in the second
    -pw[pwm_width]Integer value for PWM width. Defaults to 24
    python build_pwm_from_bigwig.py -i unstranded.bw -g GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta -o /path/name_of_pwm -c chr20 -cz hg38.chrom.sizes -pw 24
  3. Reformat ChromBPNet H5 files

    master

    Use the reformat_chrombpnet_h5.py script to reformat ChromBPNet H5 files. This is typically used to process bias-factorized model outputs by providing both the no-bias ChromBPNet H5 file and the scaled bias model H5 file.

    Note: Ensure that the specified output directory exists before running the command.

    python reformat_chrombpnet_h5.py -cnb path/to/chrombpnet_nobias.h5 -bm /path/to/bias_model_scaled.h5 -o /path/to/outputdir
  4. Run MODISCO on ChromBPNet DeepSHAP outputs

    master

    The modisco.sh script automates a three-step pipeline for de-novo motif discovery and annotation from ChromBPNet DeepSHAP outputs:

    1. De-novo motif discovery: Runs run_modisco.py on the DeepSHAP output.
    2. Motif annotation: Uses TOMTOM via fetch_tomtom.py to annotate discovered motifs.
    3. Visualization: Summarizes results into an HTML format using visualize_motif_matches.py.

    Note: If you do not wish to host the results online, you can remove visualize_motif_matches.py from the run.sh script to skip the HTML generation step.

    modisco.sh [scores_prefix] [output_dir] [score_type] [seqlets] [crop] [meme_db] [meme_logos] [vier_logos] [vier_html] [html_link]
  5. Generate bigwig tracks from ChromBPNet predictions

    master

    Use the chrombpnet_predict_to_bigwig script to generate bigwig tracks from ChromBPNet predictions on a set of input regions.

    Important Requirements:

    • You must provide absolute paths to the script.
    • The input regions file (-r) must be a 10-column BED file where the number of regions matches the number of entries in the importance matrix (f['projected_shap']['seq'].shape[0]).
    • The script reads the BED file twice, so it should not be piped.
    • Region Calculation: If the BED file uses the 2nd column for start and the 10th column for summit, the input regions are calculated as [start + summit - (inputlen/2) : start + summit + (inputlen/2)].
    chrombpnet_predict_to_bigwig -bm /path/to/bias_model.h5 -cm /path/to/chrombpnet_model.h5 -cmb /path/to/chrombpnet_model_nb.h5 -r /path/to/regions.bed -g /path/to/genome.fa -c /path/to/chrom_sizes.txt -o /path/to/output_prefix
  6. Install ChromBPNet locally via Conda and Pip

    master

    For a local installation, it is recommended to use a GPU with NVIDIA drivers and CUDA already installed. Follow these steps:

    1. Create and activate a clean Conda environment with Python >= 3.8:
    conda create -n chrombpnet python=3.8
    conda activate chrombpnet
    1. Install non-Python dependencies via Conda:
    conda install -y -c conda-forge -c bioconda samtools bedtools ucsc-bedgraphtobigwig pybigwig meme
    1. Install the package using either Pip or from source:

    Via PyPI:

    pip install chrombpnet

    From Source:

    git clone https://github.com/kundajelab/chrombpnet.git
    pip install -e chrombpnet
    conda create -n chrombpnet python=3.8
    conda activate chrombpnet
    conda install -y -c conda-forge -c bioconda samtools bedtools ucsc-bedgraphtobigwig pybigwig meme
    pip install chrombpnet
  7. Install ChromBPNet via Docker

    master

    To run ChromBPNet in a pre-configured environment, use Docker. Ensure you have the NVIDIA Container Toolkit installed on your host machine to allow the container to access your system GPUs.

    Run the following command to start the container. Once inside, use cd chrombpnet to access the tutorial files.

    docker run -it --rm --memory=100g --gpus device=0  kundajelab/chrombpnet:latest
  8. Best practices for preparing input data for ChromBPNet

    master

    When preparing data for the ChromBPNet pipeline, follow these recommendations:

    • Quality Control: Ensure input files are filtered for quality metrics and duplicates are removed appropriately.
    • Duplicate Removal:
      • For Paired-End (PE) data: Duplicates can be marked and removed correctly.
      • For Single-End (SE) data: Do not remove duplicates, as this may discard actual signal information.
    • Fragment Files: Files from tools like Chromap or cellranger are suitable.
    • Single-Cell Analysis: For training cluster-specific models, consider aggregating fragments by cluster rather than by individual sample.
  9. Generate GC-matched negative regions with make_gc_matched_negatives.sh

    master

    Use the make_gc_matched_negatives.sh script to generate non-peak regions that are GC-matched to a given set of foreground (peak) regions. This is useful for training bias models or creating control sets for ChromBPNet.

    The script performs the following steps:

    1. Runs get_gc_content.py (or chrombpnet_gc_content_foreground) to find the GC content distribution of the foreground regions.
    2. Uses bedtools to filter out regions in the exclude_bed from the genome-wide candidates.
    3. Runs get_gc_matched_negatives.py (or chrombpnet_gc_matched_negatives) to sample a negative set that matches the GC distribution of the foreground.

    Note on inputlen: If a high percentage of foreground peaks are filtered out, consider reducing the inputlen value, as your genome might be too small to accommodate the requested window size.

    chmod +x make_gc_matched_negatives.sh
    ./make_gc_matched_negatives.sh [foreground_bed] [exclude_bed] [inputlen] [output_dir] [genome] [genomewide_gc] [fold] [chrom_sizes]
  10. Interpret 'counts' head using SHAP

    master

    You can calculate SHAP values for the 'counts' head of the ChromBPNet model to understand how specific bases contribute to the predicted counts.

    This process uses shap.explainers.deep.TFDeepExplainer with a specific configuration for ChromBPNet:

    1. Input: The model input and the sum of the model's second output (model.outputs[1]) across the last axis.
    2. Shuffling/Combining: Uses shap_utils.shuffle_several_times and shap_utils.combine_mult_and_diffref.

    To visualize the observed contributions (the actual contribution of the base present in the sequence), multiply the resulting counts_shap_scores by the one_hot_seqs.

    import shap
    import chrombpnet.evaluation.interpret.shap_utils as shap_utils
    
    # Setup explainer for counts head
    counts_model_input = model.input
    counts_input = one_hot_seqs
    
    counts_model_counts_explainer = shap.explainers.deep.TFDeepExplainer(
        (counts_model_input, tf.reduce_sum(model.outputs[1], axis=-1)),
        shap_utils.shuffle_several_times,
        combine_mult_and_diffref=shap_utils.combine_mult_and_diffref
    )
    
    # Generate scores
    counts_shap_scores = counts_model_counts_explainer.shap_values(counts_input, progress_message=100)
    
    # Get observed scores for visualization
    observed_counts_scores = counts_shap_scores * one_hot_seqs
  11. Interpret 'profile' head using SHAP

    master

    To understand how bases contribute to the predicted chromatin accessibility profile, use the profile head interpretation workflow.

    Unlike the counts head, the profile head explainer requires a weighted sum of mean-normed logits, which can be retrieved via shap_utils.get_weightedsum_meannormed_logits(model).

    To visualize the observed contributions, multiply the resulting profile_shap_scores by the one_hot_seqs.

    import shap
    import chrombpnet.evaluation.interpret.shap_utils as shap_utils
    
    # Setup explainer for profile head
    profile_model_input = model.input
    profile_input = one_hot_seqs
    
    weightedsum_meannormed_logits = shap_utils.get_weightedsum_meannormed_logits(model)
    
    profile_model_profile_explainer = shap.explainers.deep.TFDeepExplainer(
        (profile_model_input, weightedsum_meannormed_logits),
        shap_utils.shuffle_several_times,
        combine_mult_and_diffref=shap_utils.combine_mult_and_diffref
    )
    
    # Generate scores
    profile_shap_scores = profile_model_profile_explainer.shap_values(profile_input, progress_message=100)
    
    # Get observed scores for visualization
    observed_profile_scores = profile_shap_scores * one_hot_seqs