Skip to content

5. Heatmaps

Jump to bottom
Peter Chovanec edited this page Jul 1, 2021 · 4 revisions

Overview

The final step of the workflow is to transform a file of clusters into a file of contacts. A contacts file is a text file containing a simple square matrix of values representing the contact strength or contact frequency between any two points on the genome. This matrix is similar to an adjacency matrix, and can be easily plotted with MATLAB or with R.

The basic steps are as follows:

  1. Divide the genome or chromosome-of-interest into N bins of a given resolution. Finer resolutions will take longer to run and require much more memory. Initialize a contact matrix with dimensions N-by-N.

  2. For each cluster in the input clusters file, generate all pairs of reads and record the implied contacts by incrementing the corresponding matrix cells. A value of one-per-contact will be added to the cells if no downweighting is applied; otherwise the value will be some value less than one that depends on the downweighting strategy.

  3. Generate a vector of N bias factors corresponding to each row in the matrix. (The bias factor for a given row will be identical to the bias factor for a given column since the matrix is symmetric.) For each cell in the matrix, divide by (bias_factor(row_num) * bias_factor(col_num)). These bias factors are generated with Hi-Corrector. See Wenyuan Li, Ke Gong, Qingjiao Li, Frank Alber, Xianghong Jasmine Zhou; Hi-Corrector: a fast, scalable and memory-efficient package for normalizing large-scale Hi-C data. Bioinformatics 2015; 31 (6): 960-962. doi: 10.1093/bioinformatics/btu747 .

  4. Transform all values in the matrix to a value between zero and one. Currently this is done by calculating the median value of all cells one-off from the diagonal and dividing by it. Any values which were originally greater than this median value (and so would be greater than one after division) are simply set to one.

Usage

Run python get_sprite_contacts.py with the following flags:

NOTE: get_sprite_contacts.py uses HiCorrector (provided in the scripts folder), which needs to be executable before running the script e.g. chmod +x <fileName>.

  • --clusters: The input clusters file from the previous step.
  • --raw_contacts: The downweighted output (from step 2).
  • --biases: Hi-Corrector outputs a text file of bias factors. Save them to this location.
  • --iced: The ICEd output (from step 3).
  • --output: The final, transformed output (from step 4).
  • --assembly: Either "mm9" or "hg19". Used to initialize the contact matrix with the appropriate size.
  • --chromosome: For intrachromosomal heatmaps, one of "chr1", "chr2", ..., "chrX". For interchromosomal heatmaps, "genome".
  • --min_cluster_size and --max_cluster_size: Ignore clusters that fall outside these parameters, i.e, clusters that are too big or too small. Default: 2 to 1000.
  • --resolution: The binning resolution in nt. Default: 1000000 (1 Mb).
  • --downweighting: The downweighting strategy, one of "none", "n_minus_one", and "two_over_n".
    • none: Each contact contributes a value of 1 to the contact matrix in step 2. For example, a contact from a 5-cluster will contribute 1.
    • n_minus_one: Each contact contributes a value of 1/(N - 1), where N is the size of the corresponding cluster. For example, a contact from a 5-cluster will contribute 1/4.
    • n_over_two: Each contact contributes a value of 2/N, where N is the size of the corresponding cluster. For example, a contact from a 5-cluster will contribute 2/5.
  • --hicorrector: The path to the Hi-Corrector ic executable. This defaults to the correct location on the Guttman Lab workstation.
  • --iterations: The number of iterations to perform when running Hi-Corrector. Default: 100.

Visualize heatmap in R

  1. Run plot_heatmap.R.
  2. This takes a “final.txt” heatmap file and max value to define the cutoff of a heatmap in R.

Clone this wiki locally