Skip to content

Instantly share code, notes, and snippets.

@philippbayer
Created March 26, 2021 01:34
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save philippbayer/d84482b4ad556a7254a5d1f9a5e8e1fc to your computer and use it in GitHub Desktop.
Save philippbayer/d84482b4ad556a7254a5d1f9a5e8e1fc to your computer and use it in GitHub Desktop.
circlize

First, a tab-delimited file with genome sizes

Gm01	58711475	26	60	61
Gm03	52519505	59690052	60	61

etc.

Then, to plot the thing:

library(cowplot)
library(circlize)
library(RColorBrewer)
library(tidyverse)

circos.par("track.height" = 0.2)

#### Plot chromosomes
chrom_frame <- read_tsv('genome.list', col_names =
                          c('chr', 'end', 'bla', 'blu', 'bli'))
chrom_frame$start <- 0
chrom_frame <- chrom_frame %>% select(chr, start, end)


genes <- read_delim(
  'new.pansoy.gff',
  '\t',
  col_names = c(
    'Chrom',
    'source',
    'type',
    'start',
    'end',
    'strand',
    'strand2',
    'phase',
    'names'
  )
) %>% filter(type == 'gene')

genes <- genes %>% filter(str_detect(Chrom, 'Gm'))

circos.genomicInitialize(chrom_frame)

mycol <- brewer.pal(5, 'Dark2')

### Plot regular gene density
genes_bed <- genes %>% select('Chrom', 'start', 'end') 
names(genes_bed ) <- c('chr', 'start', 'end')

newgenes_bed <- tibble('chr' = character(),'start' =numeric(), 'end'=numeric())

# some code to flip start and end coordinates around for negative strand genes
for(i in 1:nrow(genes_bed)) {
  this_gene <- genes_bed[i,]
  if (this_gene$start > this_gene$end) {
    temp = this_gene$start
    this_gene$start = this_gene$end
    this_gene$end = temp
  }
  newgenes_bed <- bind_rows(newgenes_bed, this_gene)
}

circos.genomicDensity(as.data.frame(newgenes_bed), col=mycol[1])

At this point we have one track, the gene density. The rest of my code is various gffs etc. getting fed to circos.genomicDensity().

I've also made dotplots, for some reason those need lists of data-frames:

in_up_list <- list(newwild_modern_increase_genes_bed, newwild_modern_decrease_genes_bed)
circos.genomicRainfall(in_up_list, col=c(mycol[3], mycol[4]),pch=16, cex=0.4)

To save with fancier libraries:

p <- recordPlot() # save the current plot
library(cowplot)
save_plot(filename='circos.png', 
          plot=replayPlot(p), 
          base_height=6)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment