Learning Circos

#!/bin/bash

# extract genes

echo "Creating gene track"
grep -v \# ebola.genes.txt | cut -f 2,5,6,9 | \
    awk '{print "ebola",$2,$3,$4,"name="$1}' > track.ebola.genes.txt

# extract SNPs

echo "Creating SNP track"
grep -v \# ebola.variation.txt | cut -f 2,3,4,10,15 | \
    perl -pe 's/\w+:\w+//' | \
    awk '{print "ebola",$1,$2,$5,"snp="$3",gene="$4}' > track.ebola.variation.txt

#!/bin/bash

DATA=.
SCRIPTS=../../../scripts

# Change this to the directory where Circos is installed 
CIRCOS=/home/martink/work/circos/svn/

# grep treats special characters as normal unless they're escaped
#
#https://stackoverflow.com/questions/6775904/grepping-using-the-alternative-operator

# About 43,000 genes that are on canonical chromosomes
echo "Creating list of 43,000 genes"
cat $DATA/genes.human.refseq.txt | grep -v "#" | grep -v -i "random\|_alt\|none\|chrun" > genes.43000.txt

# Gene entries with unique names
echo "Creating list of 19,300 genes (unique names)"
cat genes.43000.txt | sort -u -k 13,13 > genes.19300.txt

# Prefixes
echo "Counting genes for each 3-letter prefix"
cat track.genes.txt | cut -d " " -f 5 | sed 's/name=//' | cut -c 1-3 | sort | uniq -c | sort -rn > prefixes.count.txt

# Gene entries with unique first 3 letter substrings in their name
echo "Creating list of 3,500 genes (unique 3-letter prefix)"
cat genes.19300.txt | awk '{print $13,$0}' | uniq -w 3 > genes.3500.txt

# Histogram of gene size
echo "Creating histogram of gene sizes (in kb)"
cat genes.19300.txt | awk '{print ($6-$5+1)/1000}' | $SCRIPTS/histogram -min 0 -max 200 -binsize 10 > histogram.gene.size.txt

# Number of genes on each chromosome
echo "Counting genes on each chromosome"
cat genes.19300.txt | cut -d $'\t' -f 3 | sort | uniq -c | sort -nr > count.bychr.txt

# Create a gene track with fields
# CHR START END NUM_EXONS name=GENENAME
echo "Creating Circos data file of gene positions"
cat genes.19300.txt | awk '{print $3,$5,$6,$9,"name="$13}' | sed 's/chr/hs/' > track.genes.txt

# Count 
BINS="1 5 10 20"
for ws in $BINS; do
  echo "Creating Circos data file of gene counts in $ws Mb windows"
  cat track.genes.txt | $CIRCOS/tools/resample/bin/resample -bin ${ws}e6 -count > track.gene.count.${ws}mb.txt
done

#!/bin/bash

DATA=.
SCRIPTS=../../../scripts

# Change this to the directory where Circos is installed 
CIRCOS=/home/martink/work/circos/svn/

# Parse the segmental duplications into a Circos link format. Ignore all duplications on
# random or unanchored chromosomes.
echo "Creating Circos segmental duplication links"
cat $DATA/segdup.human.txt | grep -v \# | cut -d $'\t' -f 2,3,4,8,9,10 | grep -v -i 'random\|chrun\|_alt\|none' | sed 's/chr/hs/g' > track.segdup.all.txt

# Histogram of sizes
echo "Creating histogram of segmental duplication size (in kb)"
cat track.segdup.all.txt | awk '{print ($3-$2-1)/1000}' | histogram.v2 -min 1 -max 20 -binsize 1 > histogram.segdup.size.txt

# Get a list of unique chromosomes from the link track ...
# ... search for the chromosome at start of line, add link size, sort by size, add i=NR where NR
#     is the awk record number and append to file tmp.txt

echo "Creating Circos segmental duplication links with size rank for each chromosome"

for chr in `cat track.segdup.all.txt | cut -d $'\t' -f 1 | sort -u` ; do
grep -w ^$chr track.segdup.all.txt | awk 'BEGIN { OFS="\t" } {print $3-$2,$0}' | sort -nr | awk 'BEGIN { OFS="\t" } {print $0,"sizerank="NR}' >> tmp.txt
done

# Sort the indexed link by size (first field), remove the field and output to a track file

cat tmp.txt | sort -nr | cut -d $'\t' -f 2- > track.segdup.indexed.txt

\rm tmp.txt

#!/bin/bash

# Run this from day.1/data

DATA=.
SCRIPTS=../../scripts

# Change this to the directory where Circos is installed 
CIRCOS=/home/martink/work/circos/svn/

echo "Creating list of 16,000 genes"
cat conservation/coords.txt | awk '{print $2,$3,$4,$4-$3+1,"name="$1}' > genes.txt

BINS="5 10 20 50 100"

for ws in $BINS; do
  echo "Counting genes in $ws kb bins"
  bin=$((ws * 1000))
  cat genes.txt | ~/work/circos/svn/tools/resample/bin/resample -bin $bin -count > genes.count.${ws}kb.txt
done

for ws in $BINS; do
  bin=$((ws * 1000))
  echo "Counting genes in $ws kb bins"
  cat genes.txt | $CIRCOS/tools/resample/bin/resample -bin $bin -count > genes.count.${ws}kb.txt
  echo "Calculating gene average size in $ws kb bins"
  cat genes.txt | $CIRCOS/tools/resample/bin/resample -bin $bin -avg > genes.avgsize.${ws}kb.txt
done

for type in conservation duplication; do
  echo "Creating $type links"
  cat $type/C* $type/S* $type/Z* | $SCRIPTS/create.yeast.links -coord $type/coords.txt > links.$type.txt
  echo "...extracting largest 10,000 links"
  N=10000
  cat links.$type.txt | awk '{print $3-$2+1,$0}' | sort -nr | head -$N | cut -d " " -f 2- > links.$type.$N.txt

  for ws in $BINS; do
      bin=$((ws * 1000))
      echo "...counting links in $ws kb bins"
      cat links.$type.txt | $CIRCOS/tools/binlinks/bin/binlinks -bin $bin -num > links.$type.count.${ws}kb.txt
  done

done

# 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)));

#!/bin/bash

echo "parsing UCSC assembly table and making karyotype file"
tail -1 ebola.assembly.txt | cut -f 4 | awk '{print "chr - ebola ebola 0",$1,"black"}' > ebola.karyotype.txt

#!/bin/bash

# Randomize all colors except white and black. Try removing white,black and see what happens! Crazy, eh?

OPT="-param ideogram/show=no -randomcolor white,black"

# Alternatively, remap the colors to the hilarious and wonderfully useless
# watermelon-mint color palette (red-white-green), as requested in
# good fun by students from the EMBO course in Izmir. It's
# particularly fitting here because it matches the colors on the
# Italian flag.
#
# OPT="-param ideogram/show=no -wmmn"

for i in `seq 1 9` ; do
f=`echo "scale=1;(10-$i)/10" | bc`
echo "Drawing for hidden fraction $f to circos.$i.png"
(circos $OPT -outputfile circos.$i.png -param hidefraction=$f > out.$i; echo "Done tile $i fraction $f") &
done

echo "Please wait 30-60 seconds until all processes report finished..."
echo "Running..."

#!/bin/bash

for i in `seq 1 9`; do
convert circos.$i.png -gravity Center -crop 850x850+0+0 circos.$i.crop.png
done

montage -mode Concatenate -geometry 850x850 circos.*.crop.png circos.tiles.png