Initial preparations

First, we have to load all the required packages and prepare our working environment. I have run itno some troubles with what turned out to be missing fonts. I have found a nice solution to this problem here.

#source("https://bioconductor.org/biocLite.R")
#biocLite("karyoploteR", )
library(extrafont)
# need only do this once!
# font_import(pattern="[A/a]rial", prompt=FALSE)

library("ggbio")
library("stringr")
library("karyoploteR")
library("biovizBase")

Loading the data

Now, we can load the data and fix column names. Next, we convert the data into GRanges object required for plotting marker densities.

gdata <- read.table(file = '~/Downloads/genome_Lilly_Mendel_v4.txt',
                    comment.char = '#')
names(gdata) <- c('rsid', 'chr', 'pos', 'gtype')

regions <- GRanges(seqnames = paste0('chr', gdata$chr), 
                   ranges = IRanges(start = gdata$pos, width = 2))

Plotting markers density

Finally, there is time to plot markers density on a karyogram. A karyogram is an ideogram showing ordered chromosomes, a map of chromosomes. We plot markers density over each chromosome using kpPlotDensity. What, I think, is a nice feature in karyoploteR is that on the karyogram it shows so-called G-bands. Giemsa staining is a specific chemical reaction that stains different parts of chromosomes in different way. It is widely used in cytogenetics to visualise chromosomes and their structure. There is one important parameter here – the window size which I have empirically set to 1Mbp.

kp <- plotKaryotype(cex=.7) %>% 
  kpAddBaseNumbers(cex=.5) %>%
  kpPlotDensity(data=regions, 
                window.size = 1e6, 
                col='slateblue', 
                border='slateblue')

So what?

Well, I think the visualisation is pretty neat, although we are missing Y axis… We will fix this soon. Apart from this, we can see low markers density around centromeres (central parts of chromosomes G-stained in red). We also see a long stretch of no marker coverage on the Y chromosome. Perhaps this is the non-pseudoautosomal region of the chromosome. Similarly, the density is very low at the beginning of chromosomes 13, 14, 15, 21 and 22 and we have an exceptionally high density on the short arm of chromosome 6. I have no idea why and not much time to research right now…

Plotting the actual marker counts.

Ok, so this was the density. What if we want the actual counts per window? Easy enough!

counts <- kp$latest.plot$computed.values$density
data_counts <- GRanges(kp$latest.plot$computed.values$windows, y=counts)
max_count <- floor(max(counts) + .1 * max(counts))

# Change the default plotting parameters
# Values established in trial-and-error
pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15 # Avoid chr names to overlap with the plot
pp$topmargin <- 150 # Avoid overlap between chromosomes
pp$bottommargin <- 150

# Do the actual plotting
kp <- plotKaryotype(cex=.7, plot.params=pp) %>% 
  kpAddBaseNumbers(cex=.5) %>% 
  kpLines(data = data_counts, ymin = 0, ymax=max_count) %>% 
  kpAxis(cex=.5, ymin=0, ymax=max_count)

As you can see, I did not manage to completely avoid overlap between chromosomes and their data panels. This can probably be further refined by playing with graphical parameters though. You have got the idea, anyway.

Zooming in

We have previously spotted a region of high marker density on the short arm of chromosome 6. Let’s try to have a closer look.

pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15
pp$topmargin <- 120
to_zoom <- GRanges('chr6', IRanges(start=20e6, end=60e6))
kp <- plotKaryotype(cex=.7, plot.params=pp, zoom = to_zoom) %>% 
  kpAddBaseNumbers(cex=.5) %>% 
  kpAxis(cex=.5, ymin=0) %>%
  kpBars(data=data_counts, y1 = data_counts$y, ymax=max_count)

Well, almost perfect, only the Y axis has disappeared in a mysterious fashion… Do not know whether it is my fault or a feature of karyoploteR. Anyway, now we can see that the highest counts occur roughly between 30-35MBp. It’d be interesting to get gene names in the region together with some annotation.

Attention! There are two additional things I should have done here if it was a true research study:

• check the exact coords of high density, nt just eyeballing 30-35Mbp,
• since 23andMe is using hg19 assembly coordinates and Ensembl database is using newer hg38 assembly, (note that hg19 is the immediate predecessor of hg38!), I should have mapped the coords from one to another using, e.g. liftOver. I actually did and the coords for the region of interest are actually close enough not to bother for the purpose of this post.

Hint To find out the names of attributes and filters available in BioMart, go here and once you are done setting them, click on XML.

# Load the neccessary libraries
library("biomaRt")

# Retrieve the data
ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl")
genes <- getBM(attributes=c("ensembl_gene_id", "external_gene_name", "description", "phenotype_description"), 
               filter=c("chromosome_name","start","end"), 
               values = list(chromosome='6', start='30000000', end='35000000'), 
               mart= ensembl)

## Cache found

# Make a nice table
library("knitr")
library("kableExtra")
genes %>% kable("html") %>%
 kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
 scroll_box(height = "500px")

# Or a even nicer table...
# library("DT")
# genes %>% datatable

Conclusion

See, now we know what is sitting in the region. Any guesses why did 23andMe put a lot of markers there?