Conditional density-based analysis of T cell signaling in single-cell data

Smita Krishnaswamy1, Matthew H. Spitzer2, Michael Mingueneau3, Sean C. Bendall2, Oren Litvin1,Erica Stone4, Dana Pe’er1,*,†, Garry P. Nolan2,†

Author Affiliations
1Department of Biological Sciences, Department of Systems Biology, Columbia University, New York, NY, USA.
2Baxter Laboratory in Stem Cell Biology, Department of Microbiology and Immunology, Stanford University, Stanford, CA, USA.
3Division of Immunology, Department of Microbiology and Immunobiology, Harvard Medical School, Boston, MA, USA.
4Molecular Biology Section, Division of Biological Sciences, Department of Cellular and Molecular Medicine, University of California San Diego, La Jolla, CA, USA.
*Corresponding author. E-mail: dpeer@biology.columbia.edu
† These authors contributed equally to this work.

Cellular circuits sense the environment, process signals, and compute decisions using networks of interacting proteins. To model such a system, the abundance of each activated protein species can be described as a stochastic function of the abundance of other proteins. High-dimensional single-cell technologies, such as mass cytometry, offer an opportunity to characterize signaling circuit-wide. However, the challenge of developing and applying computational approaches to interpret such complex data remains. Here, we developed computational methods, based on established statistical concepts, to characterize signaling network relationships by quantifying the strengths of network edges and deriving signaling response functions. In comparing signaling between naïve and antigen-exposed CD4+ T lymphocytes, we find that although these two cell subtypes had similarly wired networks, naïve cells transmitted more information along a key signaling cascade than did antigen-exposed cells. We validated our characterization on mice lacking the extracellular-regulated mitogen-activated protein kinase (MAPK) ERK2, which showed stronger influence of pERK on pS6 (phosphorylated-ribosomal protein S6), in naïve cells as compared with antigen-exposed cells, as predicted. We demonstrate that by using cell-to-cell variation inherent in single-cell data, we can derive response functions underlying molecular circuits and drive the understanding of how cells process signals.


CellNet: Network Biology Applied to Stem Cell Engineering

Patrick Cahan, Hu Li, Samantha A. Morris, Edroaldo Lummertz da Rocha, George Q. Daley, James J. Collins5,

Refers To
Samantha A. Morris, Patrick Cahan, Hu Li, Anna M. Zhao, Adrianna K. San Roman, Ramesh A. Shivdasani, James J. Collins, George Q. Daley
Dissecting Engineered Cell Types and Enhancing Cell Fate Conversion via CellNet
Cell, Volume 158, Issue 4, 14 August 2014, Pages 889-902
PDF (4068 K) Supplementary content
Referred to by
Kee-Pyo Kim, Hans R. Schöler
CellNet—Where Your Cells Are Standing
Cell, Volume 158, Issue 4, 14 August 2014, Pages 699-701
PDF (596 K)
Samantha A. Morris, Patrick Cahan, Hu Li, Anna M. Zhao, Adrianna K. San Roman, Ramesh A. Shivdasani, James J. Collins, George Q. Daley
Dissecting Engineered Cell Types and Enhancing Cell Fate Conversion via CellNet
Cell, Volume 158, Issue 4, 14 August 2014, Pages 889-902
PDF (4068 K) Supplementary content

Somatic cell reprogramming, directed differentiation of pluripotent stem cells, and direct conversions between differentiated cell lineages represent powerful approaches to engineer cells for research and regenerative medicine. We have developed CellNet, a network biology platform that more accurately assesses the fidelity of cellular engineering than existing methodologies and generates hypotheses for improving cell derivations. Analyzing expression data from 56 published reports, we found that cells derived via directed differentiation more closely resemble their in vivo counterparts than products of direct conversion, as reflected by the establishment of target cell-type gene regulatory networks (GRNs). Furthermore, we discovered that directly converted cells fail to adequately silence expression programs of the starting population and that the establishment of unintended GRNs is common to virtually every cellular engineering paradigm. CellNet provides a platform for quantifying how closely engineered cell populations resemble their target cell type and a rational strategy to guide enhanced cellular engineering.


Dissecting Engineered Cell Types and Enhancing Cell Fate Conversion via CellNet


Conservation of trans-acting circuitry during mammalian regulatory evolution


Andrew B. Stergachis, Shane Neph, Richard Sandstrom, Eric Haugen, Alex P. Reynolds, Miaohua Zhang, Rachel Byron, Theresa Canfield, Sandra Stelhing-Sun, Kristen Lee, Robert E. Thurman, Shinny Vong, Daniel Bates, Fidencio Neri, Morgan Diegel, Erika Giste, Douglas Dunn, Jeff Vierstra, R. Scott Hansen, Audra K. Johnson, Peter J. Sabo, Matthew S. Wilken, Thomas A. Reh, Piper M. Treuting, Rajinder Kaul et al.

Nature 515, 365–370 (20 November 2014) doi:10.1038/nature13972
Received 21 February 2014 Accepted 15 October 2014 Published online 19 November 2014


The basic body plan and major physiological axes have been highly conserved during mammalian evolution, yet only a small fraction of the human genome sequence appears to be subject to evolutionary constraint. To quantify cis- versus trans-acting contributions to mammalian regulatory evolution, we performed genomic DNase I footprinting of the mouse genome across 25 cell and tissue types, collectively defining ~8.6 million transcription factor (TF) occupancy sites at nucleotide resolution. Here we show that mouse TF footprints conjointly encode a regulatory lexicon that is ~95% similar with that derived from human TF footprints. However, only ~20% of mouse TF footprints have human orthologues. Despite substantial turnover of the cis-regulatory landscape, nearly half of all pairwise regulatory interactions connecting mouse TF genes have been maintained in orthologous human cell types through evolutionary innovation of TF recognition sequences. Furthermore, the higher-level organization of mouse TF-to-TF connections into cellular network architectures is nearly identical with human. Our results indicate that evolutionary selection on mammalian gene regulation is targeted chiefly at the level of trans-regulatory circuitry, enabling and potentiating cis-regulatory plasticity.


Regression Analysis of Combined Gene Expression Regulation in Acute Myeloid Leukemia

Yue Li, Minggao Liang, Zhaolei Zhang


Gene expression is a combinatorial function of genetic/epigenetic factors such as copy number variation (CNV), DNA methylation (DM), transcription factors (TF) occupancy, and microRNA (miRNA) post-transcriptional regulation. At the maturity of microarray/sequencing technologies, large amounts of data measuring the genome-wide signals of those factors became available from Encyclopedia of DNA Elements (ENCODE) and The Cancer Genome Atlas (TCGA). However, there is a lack of an integrative model to take full advantage of these rich yet heterogeneous data. To this end, we developed RACER (Regression Analysis of Combined Expression Regulation), which fits the mRNA expression as response using as explanatory variables, the TF data from ENCODE, and CNV, DM, miRNA expression signals from TCGA. Briefly, RACER first infers the sample-specific regulatory activities by TFs and miRNAs, which are then used as inputs to infer specific TF/miRNA-gene interactions. Such a two-stage regression framework circumvents a common difficulty in integrating ENCODE data measured in generic cell-line with the sample-specific TCGA measurements. As a case study, we integrated Acute Myeloid Leukemia (AML) data from TCGA and the related TF binding data measured in K562 from ENCODE. As a proof-of-concept, we first verified our model formalism by 10-fold cross-validation on predicting gene expression. We next evaluated RACER on recovering known regulatory interactions, and demonstrated its superior statistical power over existing methods in detecting known miRNA/TF targets. Additionally, we developed a feature selection procedure, which identified 18 regulators, whose activities clustered consistently with cytogenetic risk groups. One of the selected regulators is miR-548p, whose inferred targets were significantly enriched for leukemia-related pathway, implicating its novel role in AML pathogenesis. Moreover, survival analysis using the inferred activities identified C-Fos as a potential AML prognostic marker. Together, we provided a novel framework that successfully integrated the TCGA and ENCODE data in revealing AML-specific regulatory program at global level.


Methods for time series analysis of RNA-seq data with application to human Th17 cell differentiation

  1. Tarmo Äijö1,*,
  2. Vincent Butty2,
  3. Zhi Chen3,
  4. Verna Salo3,
  5. Subhash Tripathi3,
  6. Christopher B. Burge2,
  7. Riitta Lahesmaa3 and
  8. Harri Lähdesmäki1,3,*

+Author Affiliations

  1. 1Department of Information and Computer Science, Aalto University, FI-00076 Aalto, Finland, 2Department of Biology, Massachusetts Institute of Technology, Cambridge, MA 02139, USA and 3Turku Centre for Biotechnology, University of Turku and Åbo Akademi University, FI-20520 Turku, Finland
  1. *To whom correspondence should be addressed


Motivation: Gene expression profiling using RNA-seq is a powerful technique for screening RNA species’ landscapes and their dynamics in an unbiased way. While several advanced methods exist for differential expression analysis of RNA-seq data, proper tools to anal.yze RNA-seq time-course have not been proposed.

Results: In this study, we use RNA-seq to measure gene expression during the early human T helper 17 (Th17) cell differentiation and Tcell activation (Th0). To quantify Th17specific gene expression dynamics, we present a novel statistical methodology, DyNB, for analyzing time-course RNA-seq data. We use non-parametric Gaussian processes to model temporal correlation in gene expression and combine that with negative binomial likelihood for the count data. To account for experimentspecific biases in gene expression dynamics, such as differences in cell differentiation efficiencies, we propose a method to rescale the dynamics between replicated measurements. We develop an MCMC sampling method to make inference of differential expression dynamics between conditions. DyNB identifies several known and novel genes involved in Th17 differentiation. Analysis of differentiation efficiencies revealed consistent patterns in gene expression dynamics between different cultures. We use qRT-PCR to validate differential expression and differentiation efficiencies for selected genes. Comparison of the results with those obtained via traditional timepointwise analysis shows that time-course analysis together with time rescaling between cultures identifies differentially expressed genes which would not otherwise be detected.

Availability: An implementation of the proposed computational methods will be available at http://research.ics.aalto.fi/csb/software/

Contact: tarmo.aijo@aalto.fi or harri.lahdesmaki@aalto.fi

Supplementary information: Supplementary data are available atBioinformatics online.


Detection of active transcription factor binding sites with the combination of DNase hypersensitivity and histone modifications

  1. Eduardo G. Gusmao1,*,
  2. Christoph Dieterich2,
  3. Martin Zenke3,4 and
  4. Ivan G. Costa1,5,6,*

+Author Affiliations

  1. 1IZKF Computational Biology Research Group, Institute for Biomedical Engineering, RWTH Aachen University Medical School, 52074 Aachen, 2Computational RNA Biology Lab and Bioinformatics Core, Max Planck Institute for Biology of Ageing, 50931 Cologne, 3Department of Cell Biology, Institute for Biomedical Engineering, RWTH Aachen University Medical School, 52074, 4Helmholtz Institute for Biomedical Engineering, 52074, 5Aachen Institute for Advanced Study in Computational Engineering Science (AICES), RWTH Aachen University, 52062 Aachen, Germany and 6Center of Informatics, Federal University of Pernambuco, 50740560 Recife-PE, Brazil
  1. *To whom correspondence should be addressed
  • Received October 28, 2013.
  • Revision received June 27, 2014.
  • Accepted July 25, 2014.


Motivation: The identification of active transcriptional regulatory elements is crucial to understand regulatory networks driving cellular processes such as cell development and the onset of diseases. It has recently been shown that chromatin structure information, such as DNase I hypersensitivity (DHS) or histone modifications, significantly improves cell-specific predictions of transcription factor binding sites. However, no method has so far successfully combined both DHS and histone modification data to perform active binding site prediction.

Results: We propose here a method based on hidden Markov models to integrate DHS and histone modifications occupancy for the detection of open chromatin regions and active binding sites. We have created a framework that includes treatment of genomic signals, model training and genome-wide application. In a comparative analysis, our method obtained a good trade-off between sensitivity versus specificity and superior area under the curve statistics than competing methods. Moreover, our technique does not require further training or sequence information to generate binding location predictions. Therefore, the method can be easily applied on new cell types and allow flexible downstream analysis such asde novo motif finding.

Availability and implementation: Our framework is available as part of the Regulatory Genomics Toolbox. The software information and all benchmarking data are available at http://costalab.org/wp/dh-hmm.

Contact: ivan.costa@rwth-aachen.de or eduardo.gusmao@rwth-aachen.de

Supplementary information: Supplementary data are available atBioinformatics online.


Wigwams: identifying gene modules co-regulated across multiple biological conditions

  1. Krzysztof Polanski1,,
  2. Johanna Rhodes1,,
  3. Claire Hill2,
  4. Peijun Zhang2,
  5. Dafyd J. Jenkins1,
  6. Steven J. Kiddle1,§,
  7. Aleksey Jironkin1,
  8. Jim Beynon1,2,
  9. Vicky Buchanan-Wollaston1,2,
  10. Sascha Ott1 and
  11. Katherine J. Denby1,2,*

+ Author Affiliations

  1. 1Warwick Systems Biology Centre and 2School of Life Sciences, University of Warwick, CV4 7AL, UK
  1. *To whom correspondence should be addressed.
  • Received September 17, 2013.
  • Revision received December 12, 2013.
  • Accepted December 13, 2013.


Motivation: Identification of modules of co-regulated genes is a crucial first step towards dissecting the regulatory circuitry underlying biological processes. Co-regulated genes are likely to reveal themselves by showing tight co-expression, e.g. high correlation of expression profiles across multiple time series datasets. However, numbers of up- or downregulated genes are often large, making it difficult to discriminate between dependent co-expression resulting from co-regulation and independent co-expression. Furthermore, modules of co-regulated genes may only show tight co-expression across a subset of the time series, i.e. show condition-dependent regulation.

Results: Wigwams is a simple and efficient method to identify gene modules showing evidence for co-regulation in multiple time series of gene expression data. Wigwams analyzes similarities of gene expression patterns within each time series (condition) and directly tests the dependence or independence of these across different conditions. The expression pattern of each gene in each subset of conditions is tested statistically as a potential signature of a condition-dependent regulatory mechanism regulating multiple genes. Wigwams does not require particular time points and can process datasets that are on different time scales. Differential expression relative to control conditions can be taken into account. The output is succinct and non-redundant, enabling gene network reconstruction to be focused on those gene modules and combinations of conditions that show evidence for shared regulatory mechanisms. Wigwams was run using six Arabidopsis time series expression datasets, producing a set of biologically significant modules spanning different combinations of conditions.

Availability and implementation: A Matlab implementation of Wigwams, complete with graphical user interfaces and documentation, is available at: warwick.ac.uk/wigwams.

Contact: k.j.denby@warwick.ac.uk

Supplementary Data: Supplementary data are available at Bioinformatics online.


Automatic Parameter Learning for Multiple Network Alignment

Jason Flannick1, Antal Novak1, Chuong B. Do1, Balaji S. Srinivasan2, and Serafim Batzoglou1
1Department of Computer Science, Stanford University, Stanford, CA 94305, USA
2Department of Statistics, Stanford University, Stanford, CA 94305, USA


We developed Græmlin 2.0, a new multiple network aligner with (1) a novel scoring func-
tion; (2) an algorithm that automatically learns the scoring function’s parameters; and (3) an
algorithm that uses the scoring function to globally align multiple networks. Existing alignment
tools use heuristic scoring functions, which must be hand-tuned to a given set of networks and
do not apply to multiple network alignment.
Our scoring function can use arbitrary features of a multiple network alignment, such as
protein deletions, protein duplications, protein mutations, and interaction losses. Our parameter
learning algorithm uses a training set of known network alignments to learn parameters for
our scoring function and thereby automatically adapts it to any set of networks. Our global
alignment algorithm finds approximate multiple network alignments in linear time.
We tested Græmlin 2.0’s accuracy on protein interaction networks from IntAct, DIP, and
the Stanford Network Database. We show that, on each of these datasets, Græmlin 2.0 has
higher sensitivity and specificity than existing network aligners. Græmlin 2.0 is available under
the GNU public license at http://graemlin.stanford.edu.


A Family of Algorithms for Computing Consensus about Node State from Network Data

Eleanor R. Brush, David C. Krakauer, Jessica C. Flack


Biological and social networks are composed of heterogeneous nodes that contribute differentially to network structure and function. A number of algorithms have been developed to measure this variation. These algorithms have proven useful for applications that require assigning scores to individual nodes–from ranking websites to determining critical species in ecosystems–yet the mechanistic basis for why they produce good rankings remains poorly understood. We show that a unifying property of these algorithms is that they quantify consensus in the network about a node’s state or capacity to perform a function. The algorithms capture consensus by either taking into account the number of a target node’s direct connections, and, when the edges are weighted, the uniformity of its weighted in-degree distribution (breadth), or by measuring net flow into a target node (depth). Using data from communication, social, and biological networks we find that that how an algorithm measures consensus–through breadth or depth– impacts its ability to correctly score nodes. We also observe variation in sensitivity to source biases in interaction/adjacency matrices: errors arising from systematic error at the node level or direct manipulation of network connectivity by nodes. Our results indicate that the breadth algorithms, which are derived from information theory, correctly score nodes (assessed using independent data) and are robust to errors. However, in cases where nodes “form opinions” about other nodes using indirect information, like reputation, depth algorithms, like Eigenvector Centrality, are required. One caveat is that Eigenvector Centrality is not robust to error unless the network is transitive or assortative. In these cases the network structure allows the depth algorithms to effectively capture breadth as well as depth. Finally, we discuss the algorithms’ cognitive and computational demands. This is an important consideration in systems in which individuals use the collective opinions of others to make decisions.


Inference of patient-specific pathway activities from multi-dimensional cancer genomics data using PARADIGM


Charles J. Vaske1,†, Stephen C. Benz2,†, J. Zachary Sanborn2, Dent Earl2, Christopher Szeto2, Jingchun Zhu2, David Haussler1,2 and Joshua M. Stuart2,*

+ Author Affiliations

1 Howard Hughes Medical Institute and 2 Department of Biomolecular Engineering and Center for Biomolecular Science and Engineering, UC Santa Cruz, CA, USA

* To whom correspondence should be addressed.


Motivation: High-throughput data is providing a comprehensive view of the molecular changes in cancer tissues. New technologies allow for the simultaneous genome-wide assay of the state of genome copy number variation, gene expression, DNA methylation and epigenetics of tumor samples and cancer cell lines.

Analyses of current data sets find that genetic alterations between patients can differ but often involve common pathways. It is therefore critical to identify relevant pathways involved in cancer progression and detect how they are altered in different patients.

Results: We present a novel method for inferring patient-specific genetic activities incorporating curated pathway interactions among genes. A gene is modeled by a factor graph as a set of interconnected variables encoding the expression and known activity of a gene and its products, allowing the incorporation of many types of omic data as evidence. The method predicts the degree to which a pathway’s activities (e.g. internal gene states, interactions or high-level ‘outputs’) are altered in the patient using probabilistic inference.

Compared with a competing pathway activity inference approach called SPIA, our method identifies altered activities in cancer-related pathways with fewer false-positives in both a glioblastoma multiform (GBM) and a breast cancer dataset. PARADIGM identified consistent pathway-level activities for subsets of the GBM patients that are overlooked when genes are considered in isolation. Further, grouping GBM patients based on their significant pathway perturbations divides them into clinically-relevant subgroups having significantly different survival outcomes. These findings suggest that therapeutics might be chosen that target genes at critical points in the commonly perturbed pathway(s) of a group of patients.

Availability:Source code available at http://sbenz.github.com/Paradigm

Contact: jstuart@soe.ucsc.edu

Supplementary information: Supplementary data are available at Bioinformatics online.