Complete Analysis Pipeline for Hi-C 3D Genome

Description

Hi-C technology is a high-throughput sequencing technique used to capture three-dimensional structural information of chromatin across the entire genome. Through this technology, we can study the spatial organization of chromatin within the cell nucleus, such as chromatin compartments, topologically associated domains (TADs), and chromatin loops. These structures play important roles in key life activities such as gene expression regulation, DNA replication, and repair.

This analysis pipeline is a complete Hi-C data analysis process, starting from the preprocessing of raw sequencing data to the final identification and visualization of advanced chromatin structures. The entire process mainly relies on bioinformatics tools such as <span>BWA-MEM2</span>, <span>SAMtools</span>, <span>HiCExplorer</span>, <span>Cooler</span>, and <span>Mustache</span>.

The analysis is primarily conducted using HiCExplorer, which provides various tools, including:

Creation and correction of contact matrices: Transforming raw Hi-C sequencing data into contact matrices and correcting them to remove experimental biases.TAD detection: Identifying topologically associated domains (TADs), which are regions that are closely connected in the three-dimensional genome structure.A/B compartment classification: Dividing the genome into different types of compartments based on the open and closed states of chromatin.Chromosome merging and rearrangement: Merging or rearranging chromosomes as needed in comparative analysis across multiple cell types.Format conversion: Supporting various data formats, such as cooler format, to facilitate data conversion between different tools.Long-range contact detection: Identifying long-range interactions on chromatin fibers.

Section 1: Environment and Parameter Setup

Before starting any analysis, it is essential to configure the working environment. This includes defining the working directory, data storage paths, reference genome locations, and computational resources required for the analysis (such as the number of threads). To ensure the reproducibility and consistency of the software environment, the script uses <span>Singularity</span> containers to run core analysis software. This approach avoids discrepancies in analysis results due to different software versions or dependency libraries.

Code

#!/bin/bash

# === Basic Path Parameter Settings ===
THREADS=36
WORK_DIR="/data/HIC_3Dgenome"
DATA_DIR="$WORK_DIR/01data"
REF="$DATA_DIR/genome.fa"

# Singularity container command aliases
SIF="singularity exec -B /data:/data hicexplorer_20240327.sif"
SIFc="singularity exec -B /data:/data cooler_20240402.sif"
SIFm="singularity exec -B /data:/data mustache_20240410.sif"

# Define output directories for each stage
MAPPED_DIR="$WORK_DIR/02mapped"
MATRIX_DIR="$WORK_DIR/03matrix"
COMPARTMENTS_DIR="$WORK_DIR/04Compartments"
TAD_DIR="$WORK_DIR/05TAD"
QC_DIR="$WORK_DIR/06QC"
LOOPS="$WORK_DIR/07Loops"

# Automatically create required directories
mkdir -p $MAPPED_DIR$MATRIX_DIR$COMPARTMENTS_DIR$TAD_DIR$QC_DIR$LOOPS

Section 2: Preparation of Reference Genome and Annotation Files

The main purposes of this step are as follows:

  1. Filtering chromosomes: Retaining only the main chromosome sequences from the original genome file (which may contain scaffolds and contigs).
  2. Generating index: Creating a <span>samtools faidx</span> index (<span>.fai</span> file) for quick access to genome sequences.
  3. Generating <span>chrom.sizes</span> file: Creating a file that records the size of each chromosome, which is used in many genomics tools.
  4. Extracting and formatting annotations: Extracting gene information corresponding to the selected chromosomes from annotation files in formats like GFF3 and converting them to GTF format.
  5. Building alignment index: Using <span>bwa-mem2 index</span> to create an index for the reference genome for subsequent rapid sequence alignment.

Code

ref_genome(){
    cd $DATA_DIR

    # 1. Generate .fai index
    samtools faidx ref_all_hapall.fa

    # 2. Generate chromosome size file (genome.chrom.sizes) and list (genome.chrom.list)
    grep "Chr" ref_all_hapall.fa.fai | awk '{print $1"\t"$2}' > genome.chrom.sizes
    cut -f 1 genome.chrom.sizes > genome.chrom.list

    # 3. Extract genome sequences by chromosome list and format to 60 bases per line
    seqtk subseq ref_all_hapall.fa genome.chrom.list | seqtk seq -l 60 > genome.fa

    # 4. Extract GTF/GFF annotations by chromosome list
    grep -wFf genome.chrom.list ref_all_hapall.fa.gff3 > genome.gff
    # gffread -T -o ref_hap01.rename.gtf ref_hap01.rename.gff3 # (Example, adjust according to actual annotation file)
    # grep -wFf genome.chrom.list ref_hapall.rename.gtf > genome.gtf

    # 5. Build bwa-mem2 alignment index
    bwa-mem2 index -p genome genome.fa

    cd $WORK_DIR
}

# Run this function
ref_genome

Section 3: Hi-C Reads Alignment

This step aligns the paired-end reads generated from Hi-C sequencing (<span>ref_1.fq.gz</span> and <span>ref_2.fq.gz</span>) to the reference genome prepared in the previous step. Since the Hi-C reads pairs come from spatially adjacent but linearly distant genomic regions, R1 and R2 fastq files are typically aligned separately. Here, we use <span>bwa-mem2</span>, which is a faster alignment tool than BWA-MEM, capable of efficiently completing the alignment task. The alignment results are stored in BAM format.

Code

map() {
    echo ">> Aligning samples..."
    # Align R1 reads
    bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -A1 -B4 -E50 -L0 -t $THREADS $DATA_DIR/genome \
        $DATA_DIR/ref_1.fq.gz | \
        samtools view -bSh - > $MAPPED_DIR/ref.bwa.1.bam
    
    # Align R2 reads
    bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -A1 -B4  -E50 -L0 -t $THREADS $DATA_DIR/genome \
        $DATA_DIR/ref_2.fq.gz | \
        samtools view -bSh - > $MAPPED_DIR/ref.bwa.2.bam
}

# Run alignment
map

Section 4: Constructing Interaction Matrix

Constructing the interaction matrix is the core step of Hi-C analysis. The <span>hicBuildMatrix</span> tool processes the aligned BAM files and generates the raw interaction matrix through the following sub-steps:

  1. Pairing and filtering: Pairing the alignment results of R1 and R2, and filtering out invalid read pairs (such as dangling ends, self-circles, etc.) based on alignment quality and sequence duplication.
  2. Assigning to bins: Allocating valid Hi-C read pairs to pre-defined genomic grids (bins). The size of the bin (<span>--binSize</span>) is the resolution of the analysis.
  3. Counting: Counting the number of read pairs falling between any two bins, which represents the interaction frequency between those two bins.
  4. Output: Finally, generating a sparse matrix file stored in <span>.cool</span> format, along with a QC folder containing detailed quality control information.

Code

build_matrix() {
    echo ">> Constructing contact matrix..."
    local sample="ref" # Assume sample name is ref
    $SIF hicBuildMatrix --samFiles $MAPPED_DIR/${sample}.bwa.1.bam $MAPPED_DIR/${sample}.bwa.2.bam \
                       --binSize 1000 \
                       --restrictionSequence GATC \
                       --danglingSequence GATC \
                       --threads $THREADS \
                       --inputBufferSize 500000 \
                       --outBam $MAPPED_DIR/${sample}_mapped.bam \
                       --outFileName $MATRIX_DIR/${sample}_matrix.cool \
                       --QCfolder $QC_DIR/${sample}_QC
}

# Run matrix construction
build_matrix

Note: The <span>find_rs_sites</span> function in the script is used to pre-locate restriction sites, but <span>hicBuildMatrix</span> can also run without providing the <span>--restrictionCutFile</span> parameter, directly identifying it through <span>--restrictionSequence</span> at runtime. For simplification, this is combined in the explanation.

Section 5: Matrix Normalization and Multi-resolution Conversion

The raw interaction matrix contains various experimental biases (such as GC content, restriction enzyme efficiency, PCR amplification bias, etc.), and must be normalized (or referred to as “calibrated” or “balanced”) to accurately reflect chromatin interactions.

  1. Normalization: <span>hicNormalize</span> can apply various algorithms for calibration. The script uses the <span>smallest</span> method, which adjusts the matrix to match the sum with the smallest sum among multiple matrices, suitable for inter-sample comparisons. A more commonly used method is <span>ICE</span> (Iterative Correction and Eigenvector decomposition), which effectively corrects the specificity bias of each bin.
  2. Multi-resolution conversion: To facilitate analysis and visualization at different scales (for example, from 1kb fine structure to 1Mb macro structure), we use <span>cooler zoomify</span> to aggregate a single resolution <span>.cool</span> file into a multi-resolution <span>.mcool</span> file. This is akin to creating a series of different zoom levels of “maps” for the genome.

Code

hicNormalize() {
    echo ">> Matrix normalization..."
    local sample="ref"
    $SIF hicNormalize --matrices $MATRIX_DIR/${sample}_matrix.cool \
                      --outFileName $MATRIX_DIR/${sample}_matrix_norm.cool \
                      --normalize smallest
}

zoomify(){
    echo ">> Creating multi-resolution matrix..."
    local sample="ref"
    # Create multi-resolution matrix balanced by ICE
    $SIFc cooler zoomify \
       $MATRIX_DIR/${sample}_matrix_norm.cool \
       -n $THREADS \
       --balance \
       --resolutions 1000,2000,5000,10000,25000,50000,100000 \
       -o $MATRIX_DIR/${sample}_matrix_norm.mcool

    # Create multi-resolution matrix without ICE balancing (optional)
    $SIFc cooler zoomify \
       $MATRIX_DIR/${sample}_matrix_norm.cool \
       -n $THREADS \
       --resolutions 1000,2000,5000,10000,25000,50000,100000 \
       -o $MATRIX_DIR/${sample}_matrix_norm_noblance.mcool
}

# Run sequentially
hicNormalize
zoomify

Section 6: Identification of A/B Compartments

At the megabase (Mb) scale, the genome is divided into A and B compartments. A compartments are typically regions with high gene density and transcriptional activity (open chromatin), while B compartments are the opposite, with low gene density and transcriptional inactivity (dense chromatin). The standard method for identifying A/B compartments is based on principal component analysis (PCA) of the interaction matrix. The value of PC1 usually clearly separates the genome into A and B categories. To determine the correspondence between the positive and negative values of PC1 and A/B compartments, it is often necessary to correlate with indicators reflecting transcriptional activity, such as gene density.

This script uses the <span>calder</span> tool, which integrates Hi-C data and gene density information to identify compartments.

  1. Calculating gene density: Using <span>bedtools</span> to calculate the number of genes within each genomic window (bin).
  2. Running Calder: Combining multi-resolution Hi-C matrix (<span>.mcool</span>) and gene density file to identify A/B compartments.

Code

Compartments() {
    echo ">> Identifying A/B compartments..."
    local bin_size=50000
    local gene_density_file="$DATA_DIR/ref_gene_density.${bin_size}.bedgraph"

    # Generate gene density file (this code needs to be adjusted according to your GFF file)
    awk '$3 == "gene"' $DATA_DIR/genome.gff | awk '{print $1, $4-1, $5}' OFS="\t" > $DATA_DIR/genes.bed
    bedtools makewindows -g $DATA_DIR/genome.chrom.sizes -w $bin_size > $DATA_DIR/genome_windows.${bin_size}.bed
    bedtools coverage -a $DATA_DIR/genome_windows.${bin_size}.bed -b $DATA_DIR/genes.bed -counts > $gene_density_file

    # Run Calder for compartment analysis
    Rscript scripts/calder \
      -i $MATRIX_DIR/ref_matrix_norm.mcool \
      -t mcool \
      -b $bin_size \
      -g others \
      --chr_prefix Chr \
      -f $gene_density_file \
      -p $THREADS -o $COMPARTMENTS_DIR/ref
}

# Run compartment analysis
Compartments

Section 7: Identification of Topologically Associated Domains (TADs)

Introduction to TADs

TADs are the basic folding units of the genome at the sub-megabase (sub-Mb) scale. DNA regions within a TAD interact frequently, while interactions between different TADs are suppressed, forming relatively independent domains. The boundary regions of TADs are often enriched with proteins such as CTCF and cohesin.<span>hicFindTADs</span> tool identifies TAD boundaries by calculating the “TAD insulation score.” Regions with lower insulation scores (i.e., areas with weak interactions with upstream and downstream) are considered TAD boundaries.

Code

find_TADs() {
    echo ">> Calling TADs..."
    local sample="ref"
    $SIF hicFindTADs --matrix $MATRIX_DIR/${sample}_matrix_norm.mcool::/resolutions/25000 \
                     --outPrefix $TAD_DIR/${sample}_TADs_25k \
                     --correctForMultipleTesting fdr \
                     --thresholdComparisons 0.05 \
                     --delta 0.01 \
                     --threads $THREADS
}

# Run TAD calling
find_TADs

Section 8: Identification of Chromatin Loops

Loops

Chromatin loops are stable contacts formed between two distant sites on the genome (usually enhancers and promoters) in space. This structure is crucial for the remote regulation of genes.<span>Mustache</span> is a fast and high-precision chromatin loop identification algorithm that identifies loops by analyzing the signal enrichment of local pixels in the Hi-C matrix.

Code

mustache() {
    echo ">> Calling chromatin loops..."
    $SIFm mustache --file $MATRIX_DIR/ref_matrix_norm.mcool \
      --resolution 10000 \
      --outfile $LOOPS/ref.10k.tsv \
      --processes $THREADS \
      -pt 0.1 \
      -st 0.7
}

# Run loop calling
mustache

Section 9: Differential Structure Analysis

Comparing the three-dimensional genome structure differences between different samples (e.g., treatment group vs. control group) can reveal structural dynamics during disease occurrence or cell differentiation processes.

  • Differential TAD analysis: <span>hicDifferentialTAD</span> can identify regions where TAD boundary strength or internal interaction frequency significantly changes between two samples.
  • Differential Loop analysis: <span>diff_mustache.py</span> is used to compare the strength of chromatin loops between two samples, identifying loops that are enhanced or weakened.

Code

hicDifferentialTAD() {
    echo ">> Performing differential TAD analysis..."
    # Assume there are two samples RG and RGS
    $SIF hicDifferentialTAD \
        --targetMatrix $MATRIX_DIR/RG_matrix_norm.mcool::/resolutions/10000 \
        --controlMatrix $MATRIX_DIR/RGS_matrix_norm.mcool::/resolutions/10000 \
        --tadDomains $TAD_DIR/RG_TADs_domains.bed \
        --outFileNamePrefix $TAD_DIR/RG_vs_RGS \
        --threads $THREADS
}

diff_loops() {
    echo ">> Performing differential Loop analysis..."
    $SIFm python /opt/conda/lib/python3.12/site-packages/mustache/diff_mustache.py \
      -f1 $MATRIX_DIR/RG_matrix_norm.mcool \
      -f2 $MATRIX_DIR/RGS_matrix_norm.mcool \
      --resolution 5000 \
      --outfile $LOOPS/RG_vs_RGS_5k_diff_loops \
      --processes $THREADS
}

# Run differential analysis (requires data from two samples)
# hicDifferentialTAD
# diff_loops

Section 10: Result Visualization

Visualization is a key step in understanding and presenting Hi-C analysis results.

  • Matrix heatmap: Using <span>hicPlotMatrix</span> can quickly plot Hi-C interaction heatmaps for specified chromosomes or regions, intuitively displaying structures like TADs.
  • Multi-track genome plots: <span>pyGenomeTracks</span> is a powerful visualization tool that can integrate Hi-C heatmaps, TAD boundaries, loops (represented as arcs), gene annotations, A/B compartment tracks, and other genomic data (such as ChIP-seq signals) into a single plot for comprehensive analysis.

Code

# Plotting matrix heatmap (example)
hicPlotMatrix_visualization() {
    $SIF hicPlotMatrix -m $MATRIX_DIR/ref_matrix_norm.mcool::/resolutions/50000 \
                      -o $MATRIX_DIR/ref_hic_plot_Chr01a.png \
                      --chromosomeOrder Chr01a \
                      --log1p
}

# Using pyGenomeTracks for fine visualization (requires tracks.ini configuration file)
hicPlotTADs_visualization() {
    cd $TAD_DIR # Enter the directory containing the configuration file
    # tracks.ini needs to be manually created by the user to define the tracks and data to display
    $SIF pyGenomeTracks --tracks tracks.ini --region Chr01a:21000000-22500000 -o Chr01a_region_TADs.pdf
    cd $WORK_DIR
}

# Run visualization
# hicPlotMatrix_visualization
# hicPlotTADs_visualization

The involved Singularity container files and integrated scripts can be obtained from the author.

Leave a Comment