Thursday 20 June 2019 — Day 1

09h00 – 10h30 | Lecture 1 — Introduction to Circos

11h00 – 12h30 | Lecture (practical) 2 — Visualizing gene distribution and size in Yeast—the histogram data track

14h00 – 15h30 | Lecture (practical) 3 — Conservation in Yeast—the link data track

16h00 – 18h00 | Lecture (practical) 4 — Drawing the human genome

18h15 – 19h30 | Lecture (practical) 5 — Afterhours—Perl refresher

18h15 – 19h30 | Lecture (practical) 6 — Afterhours—Visualizing an Ebola strain

Let's draw chromosomes from a couple of the genomes, SACE and CAGL.

Here the regular expression is using the alternation pipe, meaning that chromosome names can match either sace or cagl.

Once you've worked through this example, go back here and limit the display to only zyro-d and sace-o. Turn off the rule that hides 95% of the data.

Which chromosome combination in ../../data/links.conservation.10000.txt appears most frequently? Hint: parse the file with cut, sort and uniq to answer this question.

chromosomes_display_default = no
chromosomes = /sace|cagl/

Like plots, links are defined in blocks, except they go in <links> blocks. Individual link tracks have their own <link> block. Each line of a link data file defines two intervals: the start and end of the link.

# CHR1 START1 END1 CHR2 START2 END2
zyro-d 540803 555583 sace-l 349006 363738
zyro-d 540803 555583 sace-l 349006 363738
zyro-d 540803 555583 cagl-m 1143091 1157733
zyro-d 540803 555583 cagl-m 1143091 1157733
...

<links>
<link>
file          = ../../data/links.conservation.10000.txt

As you saw with histograms and heat maps, different plots expect different parameters. For links, you need to define the position of the ends of the links with the radius parameter. This is best done relative to the ideogram of the circle.

So here our links will go as far as 98% of the radius of the ideogram circle.

radius        = 0.98r

This controls how curved the links are. Specifically, this is the radius of the control point. Try setting this to a larger value like 0.25r or 0.5r.

bezier_radius = 0r

Since there are 10,000 links in the file, the image might take a while to draw. The record_limit parameter is useful for limiting the number of entries read from the file for debugging. You activate this if the image is taking a long time to draw. record_limit = 1000

color         = black
thickness     = 1

If you have a large data set and only want to show a random subset of the data, you can define a rule that turns the display of the data point off show = no for a fraction of the data. Here, the rand() function returns a uniformly distributed number in the interval [0,1). This number will be < 0.95 exactly 95% of the time. The condition is evaluated for each data point and will be true, on average, for 95% of the data points, which will not be shown.

The difference between this and record_limit is that here we're sampling a random subset of the data, whereas record_limit draws the data subset from the start of the file.

Activate the rules by commenting out use = no or setting use = yes and play with the cutoff fraction.

<rules>
use       = no
<rule>
condition = rand() < 0.95
show      = no
</rule>
</rules>

Uncomment these when you're drawing zyro-d vs sace-o. The links will be drawn as ribbons (with no twist in the middle). ribbon = yes flat = yes

</link>
</links>

▲ sessions/day.1/lecture.3/1/circos.png (zoom)

# read the values from STDIN and use the -col COLUMN, by default 0
my @values;
while(my $line = <>) {
  chomp $line;
  next if $line =~ /^\#/;
  my @tok = split(" ",$line);
  if(defined $CONF{col} && ! defined $tok[ $CONF{col} ]) {
    die "The line [$line] has no number in column [$CONF{col}]";
  }
  push @values, $tok[ $CONF{col}];
}

# get min and max of values
my ($min,$max) = minmax(@values);

# change min/max range if redefined via -min and/or -max
$min = $CONF{min} if defined $CONF{min};# && $CONF{min} > $min;
$max = $CONF{max} if defined $CONF{max};# && $CONF{max} < $max;

# *very* quick and dirty percentile if -pmin or -pmax used
if($CONF{pmin} || $CONF{pmax}) {
  my @sorted = sort {$a <=> $b} @values;
  $min = $sorted[ $CONF{pmin}/100 * (@sorted-1) ] if $CONF{pmin};
  $max = $sorted[ $CONF{pmax}/100 * (@sorted-1) ] if $CONF{pmax};
}

# determine bin size from either the -binsize command-line parameter
# or from -nbins. At least one of these is defined — we made sure of
# that in validateconfiguration()
my $bin_size;
if($CONF{binsize}) {
  $bin_size = $CONF{binsize};
  $CONF{nbins} = round(($max-$min)/$CONF{binsize});
} else {
  $bin_size = ($max-$min)/$CONF{nbins};
}

# populate bins
my $bins;
for my $i (0..$CONF{nbins}-1) {
  my $start = $min + $i * $bin_size;
  my $end   = $min + ($i+1) * $bin_size;
  my $n     = grep($_ >= $start && ($_ < $end || ( $i == $CONF{nbins}-1 && $_ <= $end)), @values);
  push @$bins, { start=>$start,
                 i=>$i,
                 end=>$end,
                 size=>$bin_size,
                 n=>$n };
}

my $histogram_width = $CONF{height};
my $maxn   = max(map { $_->{n} } @$bins);
my $sumn   = sum(map { $_->{n} } @$bins);
my $cumuln = 0;

# report number and fraction of values smaller than requested minimum
printinfo(sprintf("%10s %10.4f>%6d %6.3f","",
                  $min,
                  int(grep($_ < $min, @values)),
                  int(grep($_ < $min, @values))/@values));

# report each bin: min, max, count, fractional count, cumulative count
for my $bin (@$bins) {
  $cumuln += $bin->{n};
  my $binwidth = $bin->{n}*$histogram_width/$maxn;
  printinfo(sprintf("%10.4f %10.4f %6d %6.3f %6.3f %-${histogram_width}s",
                    $bin->{start},
                    $bin->{end},
                    $bin->{n},
                    $bin->{n}/$sumn,
                    $cumuln/$sumn,
                    "*"x$binwidth));
}

# report number and fraction of values larger than requested maximum
printinfo(sprintf("%10s %10.4f<%6d %6.3f","",
                  $max,
                  int(grep($_ > $max, @values)),
                  int(grep($_ > $max, @values))/@values));

# aggregate stats for full data set
printinfo(sprintf("%8s %12d","n",int(@values)));
printinfo(sprintf("%8s %12.5f","average",average(@values)));
printinfo(sprintf("%8s %12.5f","sd",stddev(@values)));
printinfo(sprintf("%8s %12.5f","min",min(@values)));
printinfo(sprintf("%8s %12.5f","max",max(@values)));
printinfo(sprintf("%8s %12.5f","sum",sum(@values)));

chromosomes_display_default = no
chromosomes = /sace|cagl/

<links>
<link>
file          = ../../data/links.conservation.10000.txt
radius        = 0.98r
bezier_radius = 0r
thickness     = 1

record_limit = 1000 If you drew a lot of links in the previous section, the entire image was saturated with black because all the links overlapped and filled the image.

You can set a color's transparency by using the _aN suffix on the color's name. In the underlying Circos configuration files (which are included at the bottom of this file) I've defined a parameter

auto_alpha_steps = 20

This gives you 20 levels of transparency for each color.

black     -> opaque
black_a1  -> 1/21 transparent
black_a2  -> 2/21 transparent
...
black_a20 -> 20/21 transparent

For obvious reasons, you don't need a 100% transparent version of each color, so the denominator in the fractions above is N+1. If you're interested look in the sessions/day.1/etc directory. If you follow all the include directives up the tree, you'll find that this file is included by all the parts of all the lectures for this day.

color         = black_a18

Even with _a18 (85% = 18/21 transparency) the image looks quite busy. Try _a20. In the last lecture, you saw how rules were used to change the color of histogram and heat map bins. In the same way, they can be used to change how links are drawn.

Turn the rules on below by setting use=yes.

<rules>
use = no

Let's remove some links form the image. You discovered in the last section that the minimum size of the start of the links in this top 10,000 set is about 3.4 kb.

Let's only show those that have a size of at least 5 kb.

To do this, we set a condition that triggers for links smaller than 5 kb. Any other directives or parameters in the rule will be applied to these links before they're drawn. In this case setting show=no will hide the links.

size1 is the size of the interval of the start of the link and size2 is the size of its end. Here we're testing whether either of these is smaller than 4kb.

<rule>
condition     = var(size1) < 5000 || var(size2) < 5000
show          = no
</rule>

The rule below maps the color of the link based on the size of its start given by var(size1) onto the 11-color diverging spectral Brewer palette. Two additional lines change the thickness of the link and its z-depth.

Circos defines a large number of named colors. All these definitions are in the Circos installation directory under etc/. For convenience, I have included these files in colornames/ in this lecture's section. You'll see the contents of colors.conf below.

<rule>
use           = no

This condition applies to all data points that reach this rule.

condition     = 1

The function remap_int(var,min,max,mapmin,mapmax) linearly maps the value given by var from the range [min,max] to [mapmin,mapmax]. We're using sprintf(fmt,list) to format the value as a string that corresponds to the name of a Brewer color.

The eval() is required to parse the contents as code as opposed to treating them as a string literal. Without eval(), z would be set to the literal string "var(size1") and not to the actual value of the size1 parameter.

color         = eval(sprintf("spectral-11-div-%d",remap_int(var(size1),5000,10000,1,11)))
thickness     = eval(sprintf("%d",remap_int(var(size1),5000,10000,1,3)))

The z parameter determines the order in which links are drawn. Links with lower z are drawn before (i.e. below) those with higher z.

z             = eval(var(size1))

</rule>
</rules>

</link>
</links>

▲ sessions/day.1/lecture.3/3/circos.png (zoom)

RGB color definition. Colors are refered to within configuration files by their name. In order to use a color, you must define it here.

Default triplet is rgb

red = 255,0,0
red = rgb(255,0,0)

You can use LCH and HSV also

red = hsv(0,1,1)
red = lch(54,105,40)

Common colors (red, green, blue, etc) are already defined in 7 tones, based on 7 color Brewer palettes, where available.

vvlHUE (very very light HUE, e.g. vvlred)
vlHUE (very light HUE, e.g. vlred)
lHUE  (light HUE, e.g. red)
HUE   (e.g. red)
dHUE  (dark HUE, e.g. dred)
vdHUE (very dark HUE, e.g. vdred))
vvdHUE (very very dark HUE, e.g. vvdred))

Pure colors are prefixed with "p" (e.g. pred, dpred, vdpred) and are based on 10% luminance steps from a pure hue color. Use these at your discretion.

Boutique colors are not. If you really must use 'bisque', then add

bisque = 255,228,196

See colors.unixnames.txt for a large number of named color definitions.

In addition to hues, two other color groups are defined: (1) Cytogenetic band colors (e.g. gposNNN, acen, stalk, etc.) which correspond to colors on ideogram bands. (2) UCSC chromosome color palette (e.g. chrNN, chrUn, chrNA), included from colors.ucsc.conf.

Another useful list of colors is in colors.brewer.conf, which defines Brewer Palettes, and corresponding lists.

white       = 255,255,255
vvvvlgrey   = greys-9-seq-1
vvvlgrey    = greys-9-seq-1
vvlgrey     = greys-9-seq-2
vlgrey      = greys-9-seq-3
lgrey       = greys-9-seq-4
grey        = greys-9-seq-5
dgrey       = greys-9-seq-6
vdgrey      = greys-9-seq-7
vvdgrey     = greys-9-seq-8
vvvdgrey    = greys-9-seq-9
vvvvdgrey   = greys-9-seq-9
black       = 0,0,0
transparent = 0,0,1

vvlred  = reds-7-seq-1
vlred   = reds-7-seq-2
lred    = reds-7-seq-3
red     = reds-7-seq-4
dred    = reds-7-seq-5
vdred   = reds-7-seq-6
vvdred  = reds-7-seq-7

vvlpred = 255,82,82
vlpred  = 255,63,63
lpred   = 255,41,41
pred    = 255,0,0
dpred   = 237,0,0
vdpred  = 219,0,0
vvdpred = 201,0,0

vvlgreen = greens-7-seq-1
vlgreen  = greens-7-seq-2
lgreen   = greens-7-seq-3
green    = greens-7-seq-4
dgreen   = greens-7-seq-5
vdgreen  = greens-7-seq-6
vvdgreen = greens-7-seq-7

vvlpgreen = 145,255,145
vlpgreen  = 112,255,112
lpgreen   = 75,255,75
pgreen    = 0,255,0
dpgreen   = 0,229,0
vdpgreen  = 0,204,0
vvdpgreen = 0,180,0

vvlblue = blues-7-seq-1
vlblue  = blues-7-seq-2
lblue   = blues-7-seq-3
blue    = blues-7-seq-4
dblue   = blues-7-seq-5
vdblue  = blues-7-seq-6
vvdblue = blues-7-seq-7

vvlpblue = 120,224,250
vlpblue  = 95,206,255
lpblue   = 64,188,255
pblue    = 0,170,255
dpblue   = 0,153,237
vdpblue  = 0,136,220
vvdpblue = 0,120,204

vvlpurple = purples-7-seq-1
vlpurple  = purples-7-seq-2
lpurple   = purples-7-seq-3
purple    = purples-7-seq-4
dpurple   = purples-7-seq-5
vdpurple  = purples-7-seq-6
vvdpurple = purples-7-seq-7

vvlppurple = 166,158,255
vlppurple  = 151,143,255
lppurple   = 136,128,255
ppurple    = 122,112,255
dppurple   = 107,97,243
vdppurple  = 92,80,232
vvdppurple = 77,62,224

vvlorange = oranges-7-seq-1
vlorange  = oranges-7-seq-2
lorange   = oranges-7-seq-3
orange    = oranges-7-seq-4
dorange   = oranges-7-seq-5
vdorange  = oranges-7-seq-6
vvdorange = oranges-7-seq-7

vvlporange = 255,182,106
vlporange  = 255,164,82
lporange   = 255,146,54
porange    = 255,127,0
dporange   = 234,110,0
vdporange  = 213,92,0
vvdporange = 193,75,0

Since there is no sequential yellow Brewer palette, a CMYK-based ladder is used here.

vvlyellow = 255,252,213
vlyellow  = 255,249,174 
lyellow   = 255,246,133 
yellow    = 255,244,80
dyellow   = 255,242,0
vdyellow  = 227,210,0
vvdyellow = 193,180,0

vvlpyellow = vvlyellow
vlpyellow  = vlyellow
lpyellow   = lyellow
pyellow    = yellow
dpyellow   = dyellow
vdpyellow  = vdyellow
vvdpyellow = vvdyellow

Brewer palettes

<<include colors.brewer.conf>>

UCSC chromosome colors

<<include colors.ucsc.conf>>

HSV pure colors

<<include colors.hsv.conf>>

chromosomes_display_default = no
chromosomes = /sace|cagl/

<links>
<link>
file          = ../../data/links.conservation.10000.txt
radius        = 0.70r
bezier_radius = 0r
thickness     = 1
color         = black_a15

<rules>

<rule>
condition     = var(size1) < 5500
show          = no
</rule>

Try other Brewer palettes here, such as piyg-11-div and reds-11-seq. Make sure you understand how remap_int is working.

<rule>
condition     = 1
color         = eval(sprintf("spectral-11-div-%d",remap_int(var(size1),5500,10000,1,11)))
z             = eval(var(size1))
</rule>

</rules>

</link>
</links>

These are histograms of the counts of the conservation links in 50 kb windows.

<plots>
<plot>
file        = ../../data/links.conservation.count.50kb.txt
type        = histogram
fill_color  = black
stroke_thickness  = 0
r1 = 0.79r
r0 = 0.70r 
orientation = out

You can reference the min, max, average (avg), standard deviation (sd) of a track using var(plot_STATISTIC). For example, var(plot_avg) gives the average and the rule above, when activated, fills assigns a black fill to bins whose size smaller than track average.

<rules>
use        = no
<rule>
condition  = var(value) < var(plot_avg)
fill_color = grey
</rule>
</rules>

Use the histogram script to generate a distribution of link counts from the data file used here. What is the 80% percentile count? Change the y0 and y1 values below to (a) color background light red above 80% percentile and (b) color background light blue between 50% and 80% percentile.

<backgrounds>
<background>
color = vvlgrey

y1 = 300 # change this to 80% percentile y0 = 50 # change this to 50% percentile

</background>

</backgrounds>

</plot>

</plots>

▲ sessions/day.1/lecture.3/4/circos.png (zoom)

chromosomes_display_default = no

chromosomes = /[scz]/ # why does this display all the chromosomes?

In this section, we'll extend the previous example and add intrachromosomal links, which will be oriented outwards.

<links>

<link>
file          = ../../data/links.conservation.10000.txt
radius        = 0.70r
bezier_radius = 0r
thickness     = 1
color         = black_a15

<rules>
use           = no

<rule>
condition     = var(size1) < 6000
show          = no
</rule>

<rule>
condition     = 1
color         = eval(sprintf("spectral-11-div-%d",remap_int(var(size1),5000,10000,1,11)))
z             = eval(var(size1))

</rule>
</rules>

</link>

If you set the bezier_radius to be larger than the radius, the control point for the link curve will be placed further from the center of the circle than the ends and cause the link to face away from the circle. For more details about link geometry, see Link geometry Circos tutorial.

<link>
file          = ../../data/links.duplication.10000.txt
radius        = 0.90r
bezier_radius = 1.15r
thickness     = 1
color         = black_a15

Create an image without any rules first&mdash;what a mess! Activete the rules one by one.

<rules>
use           = no

What does the rule below do?

<rule>
condition     = var(chr1) ne var(chr2) || var(size1) < 4000
show          = no
</rule>

<rule>
condition     = 1
color         = eval(sprintf("spectral-11-div-%d",remap_int(var(size1),5000,10000,1,11)))
z             = eval(var(size1))
</rule>
</rules>

</link>

</links>

These are histograms of the counts of the conservation links in 50 kb windows.

<plots>
<plot>
file        = ../../data/links.conservation.count.50kb.txt
type        = histogram
fill_color  = black
stroke_thickness  = 0
r1 = 0.77r
r0 = 0.70r 
orientation = out

<rules>
use = no

<rule>
condition  = var(value) < var(plot_avg)
fill_color = grey
</rule>
</rules>

<backgrounds>
<background>
color = vvlgrey
</background>
</backgrounds>

</plot>

<plot>
file             = ../../data/links.duplication.count.50kb.txt
type             = histogram
fill_color       = black
stroke_thickness = 0
r1 = 0.90r 
r0 = 0.83r 
orientation = in

You can reference the min, max, average (avg), standard deviation (sd) of a track using var(plot_STATISTIC). For example, var(plot_avg) gives the average and the rule above, when activated, fills assigns a black fill to bins whose size smaller than track average.

<rules>
use = no

<rule>
condition  = var(value) < var(plot_avg)
fill_color = grey
</rule>
</rules>

<backgrounds>
<background>
color = vvlgrey
</background>
</backgrounds>

</plot>

</plots>

<<include ideogram.conf>>

We're not going to show ticks in this image.

<<include ../../../etc/ticks.conf>>

▲ sessions/day.1/lecture.3/5/circos.png (zoom)