Dots is a Python package for working with microarray data. Its back-end is a standalone package for reading in, normalisation, statistical analysis and plotting of Agilent single-colour microarray data. Its front-end isn't finished yet (more on that below).
IMPORTANT Dots requires PhantomJS to render the html plots produced by Bokeh to png image files. The most recent version of PhantomJS (2.0) does not seem to work properly on OS X El Capitan, so I'd recommend using PhantomJS version 1.9.8. Download links for OS X, Windows and Linux are below.
Dots is now on PyPi! You can install it as follows:
sudo pip install dots_for_microarrays
Dots has a number of dependencies including NumPy and SciPy and the least painful way of getting these is to use the Anaconda Python distribution which includes NumPy and SciPy and a couple of the other required dependencies like Pandas, Scikit-learn and Bokeh.
Once you've downloaded and installed Anaconda, you can install the Dots package as follows:
git clone https://github.com/sandyjmacdonald/dots_for_microarrays cd dots_for_microarrays sudo python setup.py install
Setuptools should take care of the dependencies but, in testing, I've found the NumPy, SciPy and Scikit-learn installations to be problematic, hence my recommendation of using Anaconda to relieve those headaches.
Once you have Anaconda, if you'd like to install and use Dots in a fenced-off virtual environment that won't interfere with anything else, then you can do so as follows:
conda create --yes -n dots_env python=2.7 source activate dots_env git clone https://github.com/sandyjmacdonald/dots_for_microarrays cd dots_for_microarrays sudo python setup.py install
To test that the installation has worked, you can run Nosetests as follows:
sudo python setup.py nosetests
What Dots does
- Reads in a series of Agilent single-color array files.
**It's important that your array files are named correctly, in order for Dots to work out
to which group and replicate they belong e.g. for treated and untreated groups each with
three replicates name the files
treated_1.txt, treated_2.txt, treated_3.txt, untreated_1.txt, untreated_2.txt, untreated_3.txt.
- Normalises the data by log2-transforming, 75th percentile-shifting and setting the baseline to the median for each gene across all samples.
- Calculates fold changes and log fold changes for all of the pairs of groups.
- Runs either a T-test or ANOVA (determined automagically by the number of groups) with Benjamini-Hochberg p-value adjustment and a Tukey HSD post hoc test to determine signifcant pairs from the ANOVA.
- Provides a number of different visualisations of the data: box and whisker plots of the normalised data for each sample, a PCA plot of all of the samples, a hierarchically-clustered (by gene) heatmap for the significantly differentially expressed genes (> +/- 2-fold, p < 0.05), a plot of k-means clustered groups of genes with similar expression patterns across the samples, and volcano plots for each pair of samples.
- Generates tab-separated tables of the normalised data and fold change and statstical analysis.
- Did I mention the hover-y plots? It has hover-y plots.
What Dots will do in the future
- Read the array data into an SQLite3 database, signifcantly speeding the whole workflow if you re-analyse your array data at a later date.
- Assess the quality of the arrays.
- Provide a web front-end to guide you through the workflow.
- (Possibly) read in Affymetrix array data.
- Use linear models (similar to Limma) for the stats.
- Run functional analyses (Gene Ontology, pathways).
Make you a cappuccino while you wait.
Dots has a handy workflow script that takes as input a folder containing some Agilent array files (labelled correctly as explained above) and reads in the data, normalises it, and produces tables of data and all of the various volcano plots, etc.
You can run it, for example, on the sample data included here (in
python dots_workflow.py dots_sample_data -o sample_data_output
-o is an optional argument and, if you don't include it, then they'll be put in a
Getting your hands dirty
I've tried to comment the code as thoroughly as possible, so the best way to find out everything
that it can do is to dig into the code. Currently, it's organised in three modules that handle
reading in the arrays, analysing them and the plotting. The docstrings allow you to get information
about a function or class by typing, e.g.
The three modules -
dots_plotting - are all part of the
dots_backend package. As an example, you can import
dots_arrays by typing
from dots_backend import dots_arrays
or to import the
read_experiment function from
dots_arrays you can type
from dots_backend.dots_arrays import read_experiment
The dots_arrays module
This module handles reading in individual arrays or a series of arrays as an experiment. It
has classes for each of these - the
Array class and the
Experiment class - that have a
bunch of methods that you can run on them.
There are also two functions -
read_experiment - that pretty much do what
they say on the tin. Both of these return
Experiment instances that both have
Pandas data frame attributes that contain the data. Where possible, the dots modules use
Pandas data frames because... Pandas.
The Array class and read_array function
You can read in an individual array as follows:
array = read_array(filename, group, replicate)
To get a data frame with all of the array data in:
array_df = array.df
To get a list of the gene names:
genes = array.genenames
There are also attributes that store the probe names (
array.probenames), systematic names
array.systematicnames), descriptions (
array.descriptions), intensities (array.intensities),
array.group), replicate (
array.replicate) and sample id (
When arrays are read in, they aren't normalised unless you use the
norm_array = array.normalise()
And there's also a handy
get_normalised_intensities method that returns a list of the normalised
norm_intensities = array.get_normalised_intensities()
The Experiment class and read_experiment function
You can read in a whole experiment as follows:
experiment = read_experiment(array_filenames, baseline=True)
arrays_filenames should be a list of filenames and the
baseline option determines whether
the baseline is set to the median.
Experiment class is essentially a collection of
Array instances with some neat methods to,
for example, set the baseline to the median of the samples.
To set the baseline to median:
experiment_med = experiment.baseline_to_median()
To get just the normalised intensity values from your experiment:
norm_intensities = experiment_med.get_exp_values()
There's also a handy method to remove a sample from your experiment, although you can only do this if you haven't already set the baseline to median.
experiment = experiment.remove_sample('treated_1')
This method will be of more use once the quality control features are added, allowing you to remove samples that are of low quality before proceeding with the analysis and plotting.
The read_annotations function
For array files that lack gene symbols or descriptions, you can pass in an annotations file that you downloaded from the Agilent website or from the NCBI GEO database. These annotations will be added to your Array and Experiment class instances.
read_annotations function takes a filename as an argument and returns a dictionary of the
annotations with the key being the probename and the values being a dictionary of
You can read them in as part of your
read_experiment call as follows:
experiment = read_experiment(array_filenames, baseline=True, annotations_file='annotations.txt')
The dots_analysis module
This is the meat of the dots_backend.
get_fold_changes function is straightforward and just takes an experiment instance and
returns a data frame with e.g.
logFC_treated_untreated columns for
each pair of groups in the experiment. Use it as follows:
fold_changes = get_fold_changes(experiment)
run_stats function is similarly simple. It automagically decides whether to run just a
T-test (if there are two groups) or to run an ANOVA and Tukey HSD post hoc (if there are three
or more groups), and also adjusts the p values with a Benjamini-Hochberg correction. It
returns a data frame with
p_val_adj columns. The significances from the post
hoc test are in columns in the data frame named e.g.
significant_treated_untreated. Use it
stats = run_stats(experiment)
There's a simple
run_pca function that is used by the
do_pcaplot function in the
module. It returns a data frame with the x/y coordinates from the first two principal components.
Use it as follows:
pca_df = run_pca(experiment)
There are two functions for clustering your experiment,
do slightly different things. Both functions have an option to select either k-means or hierarchical
find_clusters function returns a list of cluster numbers in the same order as the rows in the
experiment data frame. If the method is hierarchical -
- then the number of
clusters is set at the square root of (number of rows divided by two), a good approximation. If the
method is k-means -how='kmeans'
- then values of k (the number of clusters) for a range of square
numbers are tested (4, 9, 16, 25)
using [silhouette analysis](https://en.wikipedia.org/wiki/Silhouette_(clustering)) and the best value
picked. An additional argument passed to the function -
k_vals=range(3,51)` allows you to increase the number of values tested to, in this example, 50. Here's
how to get a list of clusters with either hierarchical or k-means clustering:
hier_clusters = find_clusters(experiment_med.df, how='hierarchical') km_clusters = find_clusters(experiment_med.df, k_vals=range(3,11), how='kmeans')
get_clusters function includes the functionality of the
find_clustersfunction, but first
filters the data frame down to only significantly differentially expressed genes (this speeds
things up considerably, especially for the k-means clustering, and makes the heat maps more compact).
It returns the filtered data frame along with an extra column
cluster that contains the cluster
numbers. As with the
find_clusters function, it allows hierarchical or k-means clustering to be
selected. Here's how to get a filtered data frame with clusters with either hierarchical or k-means
hier_clusters_df = get_clusters(experiment_med.df, how='hierarchical') km_clusters_df = find_clusters(experiment_med.df, how='kmeans')
The last two functions in the
dots_analysis module deal with writing the fold changes and stats,
and normalised expression values to tab-separated files. You can run them as follows:
write_fcs_stats(experiment, outfile='foldchanges_stats.txt') write_normalised_expression(experiment, outfile='normalised_expression.txt')
The dots_plotting module
Dots uses (the totally awesome) Bokeh to create pretty plots of your array data. As explained above, it provides a number of different visualisations of the data:
- box and whisker plots of the normalised data for each sample
- a PCA plot of all of the samples
- a hierarchically-clustered (by gene) heatmap for the significantly differentially expressed genes
- a plot of k-means clustered groups of genes with similar expression patterns across the samples
- volcano plots for each pair of samples
You can see a couple of examples of the plots below.
All of the plot functions work similarly apart from the
do_volcanoplot function. Here are examples
of how to produce some plots from an experiment instance:
do_boxplot(experiment, show=False, image=False, html_file='boxplot.html') do_pcaplot(experiment, show=False, image=False, html_file='pcaplot.html') do_heatmap(experiment, show=False, image=False, html_file='heatmap.html') do_clusters_plot(experiment, show=True, image=False, html_file='clustersplot.html')
As you'll see, they all take an experiment instance and have a number of other optional arguments.
show=False/True argument determines whether the plot is shown in your browser after it is
generated, with the default being false. The
image=False/True argument determines whether a PNG
format image of the plot is created in addition to the HTML version. Lastly, the
allows you to specify a custom filename for your HTML plot (this is also used for the image filename).
As of version 0.2.0, the box plot outliers are limited to a total of 250,000 glyphs across all of the samples. Also, the heatmaps are limited to 2,500 rows, by tuning the fold change cutoff until there are less than 2,500 rows left. These limitations help to prevent strange behaviour when Bokeh deals with a really large number of glyphs.
All of the plots, with the exception of the clusters plot, use Bokeh's nifty hover function to show you information about the points on the plots, e.g. gene name, normalised expression values, etc.
do_volcanoplot function takes an additional argument (a tuple) that specifies the pair of groups
to plot on the volcano plot, for example:
do_volcanoplot(experiment, ('treated', 'untreated'), show=False, image=False, html_file='volcano_plot.html')