Matplotlib comprehensive chromosome drawing
1
1
Entering edit mode
9.6 years ago
Sakti ▴ 530

Hello Biostars,

Ryan Dale posted 4 years ago a wonderful answer to the question of how to start your own chromosome plots using Python's Matplotlib Visualize Chromosome With Python ?

I am migrating from Perl to Python and have been bashing my head for 3 days trying to add extra lines and gene annotations to Ryan's code, so that the extra gene "tracks" appear below each corresponding chromosome. I'd also like to zoom in to specific regions, and feed that as ranges in the code (i.e. from positions 1-100,000 in each chromosome). So far I have been unsuccessful :(

Anyone with matplotlib experience in here who could give me a hand?

While there are many chromosome plotters out there I would like to learn how to use Matplotlib and Python. I have seen many programmers do rather amazing publication quality figures with it.

Best,
Sakti

matplotlib python • 11k views
ADD COMMENT
11
Entering edit mode
9.6 years ago
Ryan Dale 5.0k

Upon revisiting this, I've now made it more general. You can now give it any BED file to plot along with the ideograms; here I'm just showing UCSC's knownGenes table from hg38.

To just show a particular region, you could use BEDTools outside of the script or pybedtools inside of it, to subset the input files to your region of interest. You might find pybedtools.BedTool.from_dataframe() and pybedtools.BedTool.to_dataframe() helpful for this. Also note you'll probably have to tweak the chromosome list and figure size appropriately.

I've commented the script to hopefully help you along in learning Python and matplotlib.

figure_1.png
view raw figure_1.png hosted with ❤ by GitHub
#!/usr/bin/env python
"""
Demonstrates plotting chromosome ideograms and genes (or any features, really)
using matplotlib.
1) Assumes a file from UCSC's Table Browser from the "cytoBandIdeo" table,
saved as "ideogram.txt". Lines look like this::
#chrom chromStart chromEnd name gieStain
chr1 0 2300000 p36.33 gneg
chr1 2300000 5300000 p36.32 gpos25
chr1 5300000 7100000 p36.31 gneg
2) Assumes another file, "ucsc_genes.txt", which is a BED format file
downloaded from UCSC's Table Browser. This script will work with any
BED-format file.
"""
from matplotlib import pyplot as plt
from matplotlib.collections import BrokenBarHCollection
import pandas
# Here's the function that we'll call for each dataframe (once for chromosome
# ideograms, once for genes). The rest of this script will be prepping data
# for input to this function
#
def chromosome_collections(df, y_positions, height, **kwargs):
"""
Yields BrokenBarHCollection of features that can be added to an Axes
object.
Parameters
----------
df : pandas.DataFrame
Must at least have columns ['chrom', 'start', 'end', 'color']. If no
column 'width', it will be calculated from start/end.
y_positions : dict
Keys are chromosomes, values are y-value at which to anchor the
BrokenBarHCollection
height : float
Height of each BrokenBarHCollection
Additional kwargs are passed to BrokenBarHCollection
"""
del_width = False
if 'width' not in df.columns:
del_width = True
df['width'] = df['end'] - df['start']
for chrom, group in df.groupby('chrom'):
print chrom
yrange = (y_positions[chrom], height)
xranges = group[['start', 'width']].values
yield BrokenBarHCollection(
xranges, yrange, facecolors=group['colors'], **kwargs)
if del_width:
del df['width']
# Height of each ideogram
chrom_height = 1
# Spacing between consecutive ideograms
chrom_spacing = 1
# Height of the gene track. Should be smaller than `chrom_spacing` in order to
# fit correctly
gene_height = 0.4
# Padding between the top of a gene track and its corresponding ideogram
gene_padding = 0.1
# Width, height (in inches)
figsize = (6, 8)
# Decide which chromosomes to use
chromosome_list = ['chr%s' % i for i in range(1, 23) + ['M', 'X', 'Y']]
# Keep track of the y positions for ideograms and genes for each chromosome,
# and the center of each ideogram (which is where we'll put the ytick labels)
ybase = 0
chrom_ybase = {}
gene_ybase = {}
chrom_centers = {}
# Iterate in reverse so that items in the beginning of `chromosome_list` will
# appear at the top of the plot
for chrom in chromosome_list[::-1]:
chrom_ybase[chrom] = ybase
chrom_centers[chrom] = ybase + chrom_height / 2.
gene_ybase[chrom] = ybase - gene_height - gene_padding
ybase += chrom_height + chrom_spacing
# Read in ideogram.txt, downloaded from UCSC Table Browser
ideo = pandas.read_table(
'ideogram.txt',
skiprows=1,
names=['chrom', 'start', 'end', 'name', 'gieStain']
)
# Filter out chromosomes not in our list
ideo = ideo[ideo.chrom.apply(lambda x: x in chromosome_list)]
# Add a new column for width
ideo['width'] = ideo.end - ideo.start
# Colors for different chromosome stains
color_lookup = {
'gneg': (1., 1., 1.),
'gpos25': (.6, .6, .6),
'gpos50': (.4, .4, .4),
'gpos75': (.2, .2, .2),
'gpos100': (0., 0., 0.),
'acen': (.8, .4, .4),
'gvar': (.8, .8, .8),
'stalk': (.9, .9, .9),
}
# Add a new column for colors
ideo['colors'] = ideo['gieStain'].apply(lambda x: color_lookup[x])
# Same thing for genes
genes = pandas.read_table(
'ucsc_genes.txt',
names=['chrom', 'start', 'end', 'name'],
usecols=range(4))
genes = genes[genes.chrom.apply(lambda x: x in chromosome_list)]
genes['width'] = genes.end - genes.start
genes['colors'] = '#2243a8'
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111)
# Now all we have to do is call our function for the ideogram data...
print("adding ideograms...")
for collection in chromosome_collections(ideo, chrom_ybase, chrom_height):
ax.add_collection(collection)
# ...and the gene data
print("adding genes...")
for collection in chromosome_collections(
genes, gene_ybase, gene_height, alpha=0.5, linewidths=0
):
ax.add_collection(collection)
# Axes tweaking
ax.set_yticks([chrom_centers[i] for i in chromosome_list])
ax.set_yticklabels(chromosome_list)
ax.axis('tight')
plt.show()
view raw ideograms.py hosted with ❤ by GitHub
The MIT License (MIT)
Copyright (c) 2016 Ryan Dale
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. 
view raw LICENSE hosted with ❤ by GitHub

ADD COMMENT
2
Entering edit mode

Fantastic work. I noticed a little mistake in the chromosome_collection function. line 61, facecolors should be group['colors'] and not df['colors']. thanks a lot anyway @Ryan_Dale

ADD REPLY
1
Entering edit mode

Ah, you're right. It wasn't raising an exception because it looks like matplotlib.Collection allows the colors to cycle, which means they don't have to be the same length as xrange, yrange. So the stain colors were wrong for all but the first chromosome. I updated the gist (code and plot) to fix this. Thank you, nice catch.

ADD REPLY
0
Entering edit mode

Beautiful reply, and great code to continue with my Python and matlib trainings. I really appreciate your reply, and hope it's useful for more Biostarters!!!

ADD REPLY
0
Entering edit mode

@ Ryan Dale: How can I add multiple peak sets? I guess the one here is for one peak set. Could you please highlight that part or adjust the script that it works for more than one peakset.

ADD REPLY
0
Entering edit mode

I tried to run it and it did not work the same as on the picture I would recommend a karyoploteR library in R to make these vizualizations.

ADD REPLY

Login before adding your answer.

Traffic: 3321 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6