Using the ENCODE ChIA-PET dataset

Updated: 2014 March 14th

From the Wikipedia article:

Chromatin Interaction Analysis by Paired-End Tag Sequencing (ChIA-PET) is a technique that incorporates chromatin immunoprecipitation (ChIP)-based enrichment, chromatin proximity ligation, Paired-End Tags, and High-throughput sequencing to determine de novo long-range chromatin interactions genome-wide.

Let's get started on using the ENCODE ChIA-PET dataset by downloading the bed files, which has the interactions:

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetK562Pol2InteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetK562Pol2InteractionsRep2.bed.gz

#get others if you want
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetHct116Pol2InteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetHelas3Pol2InteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetK562CtcfInteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7CtcfInteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7CtcfInteractionsRep2.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7EraaInteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7EraaInteractionsRep2.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7EraaInteractionsRep3.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7Pol2InteractionsRep1.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7Pol2InteractionsRep2.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7Pol2InteractionsRep3.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetMcf7Pol2InteractionsRep4.bed.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeGisChiaPet/wgEncodeGisChiaPetNb4Pol2InteractionsRep1.bed.gz

For a preview and description of these interaction bed files, have a look at the table schema, which has the definition:

ChIA-PET Chromatin Interaction PET clusters: Two different genomic regions in the chromatin are genomically far from each other or in different chromosomes, but are spatially close to each other in the nucleus and interact with each other for regulatory functions. BED12 format is used to represent the data.

One cool way of visualising these chromtain interactions would be by using Circos (instructions for installing circos). Here's a simple Perl script to parse the bed12 files and prepare them in a format readable by Circos:

#!/bin/env perl

use strict;
use warnings;

my $usage = "Usage: $0 <infile.bed>\n";
my $infile = shift or die $usage;

my %check = ();

open(IN,'<', $infile) || die "Could not open $infile: $!\n";
while(<IN>){
   chomp;
   #chr1    875400  877590  chr1:875400..877590-chr8:126185955..126188744,2 200     .       875400  877590  255,0,0 1       2190    0
   my($chr, $start, $end, $name, @rest) = split();
   if ($name =~ /^(chr[0-9xym]+):(\d+)\.\.(\d+)-(chr[0-9xym]+):(\d+)\.\.(\d+),\d+/i){

      #inter-chromosomal interactions are on two lines
      #skip the duplicated line
      if (exists $check{$name}){
         next;
      } else {
         $check{$name} = 1;
      }

      my $chr_first = $1;
      my $start_first = $2;
      my $end_first = $3;
      my $chr_second = $4;
      my $start_second = $5;
      my $end_second = $6;
      $chr_second =~ s/chr/hs/;
      $chr_first =~ s/chr/hs/;
      $chr_first = lc($chr_first);
      $chr_second = lc($chr_second);
      print join (" ", $chr_first, $start_first, $end_first, $chr_second, $start_second, $end_second),"\n";
   } else {
      die "Could not parse $name\n";
   }
}
close(IN);
exit(0);

Now to execute the script:

./to_circos.pl wgEncodeGisChiaPetK562Pol2InteractionsRep1.bed > wgEncodeGisChiaPetK562Pol2InteractionsRep1.link

I have these three files prepared for Circos (ideogram.conf, ticks.conf and test.conf). See my getting started with Circos post to get more information on these configuration files.

#ideogram.conf
cat ideogram.conf
<ideogram>

<spacing>
default = 0.005r
</spacing>

# Ideogram position, fill and outline
radius           = 0.90r
thickness        = 20p
fill             = yes
stroke_color     = dgrey
stroke_thickness = 2p

# Minimum definition for ideogram labels.

show_label       = yes
# see etc/fonts.conf for list of font names
label_font       = default
label_radius     = dims(image,radius) - 60p
label_size       = 30
label_parallel   = yes
</ideogram>

#ticks.conf
cat ticks.conf
show_ticks          = yes
show_tick_labels    = yes

<ticks>
radius           = 1r
color            = black
thickness        = 2p

# the tick label is derived by multiplying the tick position
# by 'multiplier' and casting it in 'format':
#
# sprintf(format,position*multiplier)
#

multiplier       = 1e-6

# %d   - integer
# %f   - float
# %.1f - float with one decimal
# %.2f - float with two decimals
#
# for other formats, see http://perldoc.perl.org/functions/sprintf.html

format           = %d

<tick>
spacing        = 5u
size           = 10p
</tick>

<tick>
spacing        = 25u
size           = 15p
show_label     = yes
label_size     = 20p
label_offset   = 10p
format         = %d
</tick>

</ticks>

#test.conf
cat test.conf  | grep -v "^#" | grep -v "^$"
karyotype = data/karyotype/karyotype.human.txt
chromosomes_units = 1000000
<links>
<link>
file          = wgEncodeGisChiaPetK562Pol2InteractionsRep1.link
radius        = 0.8r
bezier_radius = 0r
color         = black_a4
thickness     = 2
<rules>
<rule>
condition     = var(intrachr)
show          = no
</rule>
<rule>
condition     = 1
color         = eval(var(chr2))
flow          = continue
</rule>
<rule>
condition     = to(hs1)
radius2       = 0.99r
</rule>
</rules>
</link>
</links>
<<include ideogram.conf>>
<<include ticks.conf>>
<image>
<<include etc/image.conf>>
</image>
<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

#now run Circos assuming you have the link file in the same directory as the conf files
bin/circos -conf test.conf

If everything worked perfectly, you should get this image:

There's a whole new level of complexity when we take the spatial organisation of chromosomes into account as well.

Intra-chromosomal

If we want to focus on chromosome one and show long range interactions (over 1 mb):

#how many intra-chromosomal interactions on chromosome 1
cat wgEncodeGisChiaPetK562Pol2InteractionsRep1.link | awk '$1=="hs1" && $4=="hs1" && $5-$2>1000000 {print}' | wc -l
164
cat wgEncodeGisChiaPetK562Pol2InteractionsRep1.link | awk '$1=="hs1" && $4=="hs1" && $5-$2>1000000 {print}' > chr1_to_chr1.link

The ticks.conf and ideogram.conf are the same. Here's what the test.conf file looks like:

cat test.conf 
karyotype = data/karyotype/karyotype.human.txt
chromosomes_units = 1000000

chromosomes_display_default = no
chromosomes                 = hs1

<links>
<link>
file          = chr1_to_chr1.link
radius        = 0.8r
bezier_radius = 0r
color         = black_a4
thickness     = 2
<rules>
<rule>
condition     = var(intrachr)
show          = yes
</rule>
<rule>
condition     = 1
color         = eval(var(chr2))
flow          = continue
</rule>
<rule>
condition     = to(hs1)
radius2       = 0.99r
</rule>
</rules>
</link>
</links>
<<include ideogram.conf>>
<<include ticks.conf>>
<image>
<<include etc/image.conf>>
</image>
<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

#run circos
circos -conf test.conf

Long range intra-chromosomal interactions on chromosome one.

Conclusions

I've showed a way of visualisation the ChIA-PET dataset but not on using the dataset. One way I intend to use this dataset is to modify the Perl script above to produce a bed file containing the genomic loci that interact with another loci. Then just to get an idea of what these regions encompass, I would intersect them with some genome annotation file.