Wishbone is an algorithm to identify bifurcating developmental trajectories from single cell data. Wishbone can applied to both single cell RNA-seq and mass cytometry datasets. This notebook details the usage of Wishbone for single cell RNA-seq data.
Wishbone can read single cell RNA-seq data from a csv file. The csv file contains cells in the rows and genes in the columns. First step is to import the package. The following code snipped imports the wishbone
package along with other plotting related imports
In [1]:
import wishbone
# Plotting and miscellaneous imports
import os
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
A sample RNA-seq csv data is installed at ~/.wishbone/data/sample_scseq_data.csv
. This sample data will be used to demonstrate the utilization and capabilities of the Wishbone package. This is a sample data derived from Paul et. al. and represents the differentiation of myeloid and erythroid precursors from hematopoietic stem cells in the mouse bone marrow
The data can be loaded using the wishbone.wb.SCData.from_csv
function.
In [2]:
# Load sample data
scdata = wishbone.wb.SCData.from_csv(os.path.expanduser('~/.wishbone/data/sample_scseq_data.csv'),
data_type='sc-seq', normalize=True)
This will create an object of the type wishbone.wb.SCData
which is the base class for the analysis. This can be either be single cell RNA-seq or mass cytometry and is specified by using data_type
parameter set to sc-seq
or masscyt
respectively. The normalize
parameter is used for correcting for library size among cells.
In [3]:
scdata
Out[3]:
This shows that the data matrix contains 4423
cells and 2312
genes along with the different properties of the wishbone.wb.SCData
class.
The scdata
object can also be initialized using a pandas DataFrame
. An example is shown below
scdata = wishbone.wb.SCData(data_frame, 'sc-seq')
scdata = scdata = scdata.normalize_scseq_data()
The first step in data processing for Wishbone is to determine metagenes using principal component analysis. This representation is necessary to overcome the extensive dropouts that are pervasive in single cell RNA-seq data
In [4]:
scdata.run_pca()
Each of the analysis function updates the scdata
object. As shown below, the pca
property of the scdata
object is now changed to True
compared to `False` when the object was created
In [5]:
scdata
Out[5]:
The results of PCA i.e., the fraction of variance explained by each component can be visualized using the function plot_pca_variance_explained
. Use the ylim
and n_components
parameters to set the y-axis limits and visualize variance explained by n_components
respectively. Typically, most of the variance is explained by the first few components (No more than 15).
In [6]:
fig, ax = scdata.plot_pca_variance_explained(ylim=(0, 0.1), n_components=30)
From this, choose the appropriate number of components using the elbow method. While tSNE
visualization is sensitive to the number of components chosen, downstream results are robust to this parameter.
Wishbone uses tSNE for visualization and tSNE can be run using the run_tsne
function which takes the number of principal components as the parameter. From the above plot, 5 seems an appropriate number of components to use.
In [7]:
NO_CMPNTS = 5
scdata.run_tsne(n_components=NO_CMPNTS, perplexity=30)
perplexity
by default is set 30. This will be reduced automatically to 15 is the number of cells is less than 100.
tSNE results can be visualized by the plot_tsne
and plot_tsne_by_cell_sizes
functions. The plot_tsne_by_cell_sizes
function colors the cells by their molecule counts before normalization
In [8]:
fig, ax = scdata.plot_tsne()
In [9]:
fig = plt.figure(figsize=[5, 4])
scdata.plot_tsne_by_cell_sizes(fig=fig)
Out[9]:
Gene expression can be visualized on tSNE maps using the plot_gene_expression
function. The genes
parameter is an string iterable of genes, which are a subset of the expression of column names. The below function plots the expression of HSC gene CD34
, myeloid gene MPO
and erythroid precursor genes GATA2
and GATA1
.
In [10]:
fig, ax = scdata.plot_gene_expression(genes = ['CD34', 'GATA2', 'GATA1', 'MPO'])
Wishbone uses diffusion maps, a non-linear dimensionality reduction technique to denoise the data and capture the major axes of variation. Diffusion maps can be determined by using the run_diffusion_map
function and the diffusion components visualized on tSNE maps using plot_diffusion_components
.
In [11]:
# Run diffusion maps
scdata.run_diffusion_map()
In [12]:
fig, ax = scdata.plot_diffusion_components()
Note the component 0 is the trivial component and does not encode any information of the data
Unlike mass cytometry data, the differentiation signals in single cell RNA-seq data by virtue of measuring the transcriptome are confounded by biological processes such as stress, cell cycle and metabolism. Since diffusion components, capture the major axes of variation, Wishbone uses gene set enrichment analysis among of genes along the top diffusion components to remove components that encode confounding signals. The relevant signals are usually present in the top 2-3 diffusion components
The first step in computing these enrichments is to determine the correlation of gene expression along each component to derive a ranking of genes along each component.
In [13]:
scdata.run_diffusion_map_correlations()
The function plot_gene_component_correlations
shows the distribution of correlations along each component.
In [14]:
fig, ax = scdata.plot_gene_component_correlations()
The enrichments can be determined using the run_gsea
function. This function needs the prefix for generating GSEA reports and a gmt
file representing the different gene sets. The following invocation of the function shows the supported set of gmt
files
Note: The gmt files package with Wishbone assume all the gene names to be upper case. This can be ensured using the following code snipped.
In [15]:
scdata.data.columns = scdata.data.columns.str.upper()
In [16]:
scdata.run_gsea( output_stem= os.path.expanduser('~/.wishbone/gsea/mouse_marrow'))
Since this is data from mouse, gmt_file
parameter can be set to (mouse, gofat.bp.v1.0.gmt.txt)
In [17]:
reports = scdata.run_gsea(output_stem= os.path.expanduser('~/.wishbone/gsea/mouse_marrow'),
gmt_file=('mouse', 'gofat.bp.v1.0.gmt.txt'))
The detailed reports can be found at ~/.wishbone/gsea/
In [18]:
!open ~/.wishbone/gsea/
run_gsea
function also returns the top enrichment gene sets along each component. GSEA determines enrichments that are either positively or negatively correlated with the gene component correlations. In this datasets, components 1 and 2 show relevant enrichments and are used for running Wishbone. Please see Selection of diffusion components for single cell RNA-seq
section of the Supplementary Methods for more details.
In [19]:
# Component 1 enrichments
reports[1]['neg']
Out[19]:
In [20]:
# Component 2 enrichments
reports[2]['pos']
Out[20]:
The SCData
object can be saved to a pickle file and loaded using the save
and load
functions.
scdata.save('mouse_marrow_scdata.p')
scdata = wishbone.wb.SCdata.load('mouse_marrow_scdata.p')
wishbone.wb.Wishbone
is the class for running Wishbone. After initialization, Wishbone can be run by specifying the start cell, components to use and number of waypoints to be used. The start cell for this dataset was chosen based on high expression of CD34.
In [21]:
# Wishbone class
wb = wishbone.wb.Wishbone(scdata)
wb.run_wishbone(start_cell='W30258', components_list=[1, 2], num_waypoints=150)
Wishbone
objects contain the SCData
object along with the identified trajectory, branch associations and waypoints
In [22]:
wb
Out[22]:
Wishbone trajectory and branch results can be visualized on tSNE maps using the plot_wishbone_on_tsne
function
In [23]:
fig, ax = wb.plot_wishbone_on_tsne()
Gene expression trends along the Wishbone trajectory can be visualized using the plot_marker_trajectory
function. This function also returns the smoothed trends along with the matplotlib fig, ax handler objects.
In [24]:
vals, fig, ax = wb.plot_marker_trajectory(['CD34', 'GATA1', 'GATA2', 'MPO']);
The marker trends can be visualized as heatmaps in a given trajectory range using the following functions
In [25]:
wb.plot_marker_heatmap(vals)
Out[25]:
In [26]:
wb.plot_marker_heatmap(vals, trajectory_range=[0.1, 0.6])
Out[26]:
The change in marker trends along the trajectory or derivatives can be visualized using these functions
In [27]:
wb.plot_derivatives(vals)
Out[27]:
In [28]:
wb.plot_derivatives(vals, trajectory_range=[0.3, 0.6])
Out[28]:
Similar to SCData
objects, Wishbone
objects can also be saved and loaded using save
and load
functions.