CaSSiDI: novel single-cell "Cluster Similarity Scoring and Distinction Index" reveals critical functions for PirB and context-dependent Cebpb repression.

PirB is an inhibitory cell surface receptor particularly prominent on myeloid cells. PirB curtails the phenotypes of activated macrophages during inflammation or tumorigenesis, but its functions in macrophage homeostasis are obscure. To elucidate PirB-related functions in macrophages at steady-state, we generated and compared single-cell RNA-sequencing (scRNAseq) datasets obtained from myeloid cell subsets of wild type (WT) and PirB-deficient knockout (PirB KO) mice. To facilitate this analysis, we developed a novel approach to clustering parameter optimization called "Cluster Similarity Scoring and Distinction Index" (CaSSiDI). We demonstrate that CaSSiDI is an adaptable computational framework that facilitates tandem analysis of two scRNAseq datasets by optimizing clustering parameters. We further show that CaSSiDI offers more advantages than a standard Seurat analysis because it allows direct comparison of two or more independently clustered datasets, thereby alleviating the need for batch-correction while identifying the most similar and different clusters. Using CaSSiDI, we found that PirB is a novel regulator of Cebpb expression that controls the generation of Ly6Clo patrolling monocytes and the expansion properties of peritoneal macrophages. PirB's effect on Cebpb is tissue-specific since it was not observed in splenic red pulp macrophages (RPMs). However, CaSSiDI revealed a segregation of the WT RPM population into a CD68loIrf8+ "neuronal-primed" subset and an CD68hiFtl1+ "iron-loaded" subset. Our results establish the utility of CaSSiDI for single-cell assay analyses and the determination of optimal clustering parameters. Our application of CaSSiDI in this study has revealed previously unknown roles for PirB in myeloid cell populations. In particular, we have discovered homeostatic functions for PirB that are related to Cebpb expression in distinct macrophage subsets.

Distinct and opposite effects of leukemogenic Idh and Tet2 mutations in hematopoietic stem and progenitor cells.

Mutations in IDH1, IDH2, and TET2 are recurrently observed in myeloid neoplasms. IDH1 and IDH2 encode isocitrate dehydrogenase isoforms, which normally catalyze the conversion of isocitrate to α-ketoglutarate (α-KG). Oncogenic IDH1/2 mutations confer neomorphic activity, leading to the production of D-2-hydroxyglutarate (D-2-HG), a potent inhibitor of α-KG-dependent enzymes which include the TET methylcytosine dioxygenases. Given their mutual exclusivity in myeloid neoplasms, IDH1, IDH2, and TET2 mutations may converge on a common oncogenic mechanism. Contrary to this expectation, we observed that they have distinct, and even opposite, effects on hematopoietic stem and progenitor cells in genetically engineered mice. Epigenetic and single-cell transcriptomic analyses revealed that Idh2R172K and Tet2 loss-of-function have divergent consequences on the expression and activity of key hematopoietic and leukemogenic regulators. Notably, chromatin accessibility and transcriptional deregulation in Idh2R172K cells were partially disconnected from DNA methylation alterations. These results highlight unanticipated divergent effects of IDH1/2 and TET2 mutations, providing support for the optimization of genotype-specific therapies.

Cholinergic control of Th17 cell pathogenicity in experimental autoimmune encephalomyelitis.

Experimental autoimmune encephalomyelitis (EAE) is a mouse model of multiple sclerosis (MS) in which Th17 cells have a crucial but unclear function. Here we show that choline acetyltransferase (ChAT), which synthesizes acetylcholine (ACh), is a critical driver of pathogenicity in EAE. Mice with ChAT-deficient Th17 cells resist disease progression and show reduced brain-infiltrating immune cells. ChAT expression in Th17 cells is linked to strong TCR signaling, expression of the transcription factor Bhlhe40, and increased Il2, Il17, Il22, and Il23r mRNA levels. ChAT expression in Th17 cells is independent of IL21r signaling but dampened by TGFß, implicating ChAT in controlling the dichotomous nature of Th17 cells. Our study establishes a cholinergic program in which ACh signaling primes chronic activation of Th17 cells, and thereby constitutes a pathogenic determinant of EAE. Our work may point to novel targets for therapeutic immunomodulation in MS.

Three tissue resident macrophage subsets coexist across organs with conserved origins and life cycles.

Resident macrophages orchestrate homeostatic, inflammatory, and reparative activities. It is appreciated that different tissues instruct specialized macrophage functions. However, individual tissues contain heterogeneous subpopulations, and how these subpopulations are related is unclear. We asked whether common transcriptional and functional elements could reveal an underlying framework across tissues. Using single-cell RNA sequencing and random forest modeling, we observed that four genes could predict three macrophage subsets that were present in murine heart, liver, lung, kidney, and brain. Parabiotic and genetic fate mapping studies revealed that these core markers predicted three unique life cycles across 17 tissues. TLF+ (expressing TIMD4 and/or LYVE1 and/or FOLR2) macrophages were maintained through self-renewal with minimal monocyte input; CCR2+ (TIMD4−LYVE1−FOLR2−) macrophages were almost entirely replaced by monocytes, and MHC-IIhi macrophages (TIMD4−LYVE1−FOLR2−CCR2−), while receiving modest monocyte contribution, were not continually replaced. Rather, monocyte-derived macrophages contributed to the resident macrophage population until they reached a defined upper limit after which they did not outcompete pre-existing resident macrophages. Developmentally, TLF+ macrophages were first to emerge in the yolk sac and early fetal organs. Fate mapping studies in the mouse and human single-cell RNA sequencing indicated that TLF+ macrophages originated from both yolk sac and fetal monocyte precursors. Furthermore, TLF+ macrophages were the most transcriptionally conserved subset across mouse tissues and between mice and humans, despite organ- and species-specific transcriptional differences. Here, we define the existence of three murine macrophage subpopulations based on common life cycle properties and core gene signatures and provide a common starting point to understand tissue macrophage heterogeneity.

Concerted roles of PTEN and ATM in controlling hematopoietic stem cell fitness and dormancy.

In order to sustain proficient life-long hematopoiesis, hematopoietic stem cells (HSCs) must possess robust mechanisms to preserve their quiescence and genome integrity. DNA-damaging stress can perturb HSC homeostasis by affecting their survival, self-renewal, and differentiation. Ablation of the kinase ataxia telangiectasia mutated (ATM), a master regulator of the DNA damage response, impairs HSC fitness. Paradoxically, we show here that loss of a single allele of Atm enhances HSC functionality in mice. To explain this observation, we explored a possible link between ATM and the tumor suppressor phosphatase and tensin homolog (PTEN), which also regulates HSC function. We generated and analyzed a knockin mouse line (PtenS398A/S398A), in which PTEN cannot be phosphorylated by ATM. Similar to Atm+/-, PtenS398A/S398A HSCs have enhanced hematopoietic reconstitution ability, accompanied by resistance to apoptosis induced by genotoxic stress. Single-cell transcriptomic analyses and functional assays revealed that dormant PtenS398A/S398A HSCs aberrantly tolerate elevated mitochondrial activity and the accumulation of reactive oxygen species, which are normally associated with HSC priming for self-renewal or differentiation. Our results unveil a molecular connection between ATM and PTEN, which couples the response to genotoxic stress and dormancy in HSCs.

Mutant ACVR1 Arrests Glial Cell Differentiation to Drive Tumorigenesis in Pediatric Gliomas.

Diffuse intrinsic pontine gliomas (DIPGs) are aggressive pediatric brain tumors for which there is currently no effective treatment. Some of these tumors combine gain-of-function mutations in ACVR1, PIK3CA, and histone H3-encoding genes. The oncogenic mechanisms of action of ACVR1 mutations are currently unknown. Using mouse models, we demonstrate that Acvr1G328V arrests the differentiation of oligodendroglial lineage cells, and cooperates with Hist1h3bK27M and Pik3caH1047R to generate high-grade diffuse gliomas. Mechanistically, Acvr1G328V upregulates transcription factors which control differentiation and DIPG cell fitness. Furthermore, we characterize E6201 as a dual inhibitor of ACVR1 and MEK1/2, and demonstrate its efficacy toward tumor cells in vivo. Collectively, our results describe an oncogenic mechanism of action for ACVR1 mutations, and suggest therapeutic strategies for DIPGs.

Author Correction: Fcmr regulates mononuclear phagocyte control of anti-tumor immunity.

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Fcmr regulates mononuclear phagocyte control of anti-tumor immunity.

Myeloid cells contribute to tumor progression, but how the constellation of receptors they express regulates their functions within the tumor microenvironment (TME) is unclear. We demonstrate that Fcmr (Toso), the putative receptor for soluble IgM, modulates myeloid cell responses to cancer. In a syngeneic melanoma model, Fcmr ablation in myeloid cells suppressed tumor growth and extended mouse survival. Fcmr deficiency increased myeloid cell population density in this malignancy and enhanced anti-tumor immunity. Single-cell RNA sequencing of Fcmr-deficient tumor-associated mononuclear phagocytes revealed a unique subset with enhanced antigen processing/presenting properties. Conversely, Fcmr activity negatively regulated the activation and migratory capacity of myeloid cells in vivo, and T cell activation by bone marrow-derived dendritic cells in vitro. Therapeutic targeting of Fcmr during oncogenesis decreased tumor growth when used as a single agent or in combination with anti-PD-1. Thus, Fcmr regulates myeloid cell activation within the TME and may be a potential therapeutic target.

AhR controls redox homeostasis and shapes the tumor microenvironment in BRCA1-associated breast cancer.

Cancer cells have higher reactive oxygen species (ROS) than normal cells, due to genetic and metabolic alterations. An emerging scenario is that cancer cells increase ROS to activate protumorigenic signaling while activating antioxidant pathways to maintain redox homeostasis. Here we show that, in basal-like and BRCA1-related breast cancer (BC), ROS levels correlate with the expression and activity of the transcription factor aryl hydrocarbon receptor (AhR). Mechanistically, ROS triggers AhR nuclear accumulation and activation to promote the transcription of both antioxidant enzymes and the epidermal growth factor receptor (EGFR) ligand, amphiregulin (AREG). In a mouse model of BRCA1-related BC, cancer-associated AhR and AREG control tumor growth and production of chemokines to attract monocytes and activate proangiogenic function of macrophages in the tumor microenvironment. Interestingly, the expression of these chemokines as well as infiltration of monocyte-lineage cells (monocyte and macrophages) positively correlated with ROS levels in basal-like BC. These data support the existence of a coordinated link between cancer-intrinsic ROS regulation and the features of tumor microenvironment. Therapeutically, chemical inhibition of AhR activity sensitizes human BC models to Erlotinib, a selective EGFR tyrosine kinase inhibitor, suggesting a promising combinatorial anticancer effect of AhR and EGFR pathway inhibition. Thus, AhR represents an attractive target to inhibit redox homeostasis and modulate the tumor promoting microenvironment of basal-like and BRCA1-associated BC.

Cis-regulatory determinants of MyoD function.

Muscle-specific transcription factor MyoD orchestrates the myogenic gene expression program by binding to short DNA motifs called E-boxes within myogenic cis-regulatory elements (CREs). Genome-wide analyses of MyoD cistrome by chromatin immnunoprecipitation sequencing shows that MyoD-bound CREs contain multiple E-boxes of various sequences. However, how E-box numbers, sequences and their spatial arrangement within CREs collectively regulate the binding affinity and transcriptional activity of MyoD remain largely unknown. Here, by an integrative analysis of MyoD cistrome combined with genome-wide analysis of key regulatory histones and gene expression data we show that the affinity landscape of MyoD is driven by multiple E-boxes, and that the overall binding affinity-and associated nucleosome positioning and epigenetic features of the CREs-crucially depend on the variant sequences and positioning of the E-boxes within the CREs. By comparative genomic analysis of single nucleotide polymorphism (SNPs) across publicly available data from 17 strains of laboratory mice, we show that variant sequences within the MyoD-bound motifs, but not their genome-wide counterparts, are under selection. At last, we show that the quantitative regulatory effect of MyoD binding on the nearby genes can, in part, be predicted by the motif composition of the CREs to which it binds. Taken together, our data suggest that motif numbers, sequences and their spatial arrangement within the myogenic CREs are important determinants of the cis-regulatory code of myogenic CREs.

Uncovering robust patterns of microRNA co-expression across cancers using Bayesian Relevance Networks.

Co-expression networks have long been used as a tool for investigating the molecular circuitry governing biological systems. However, most algorithms for constructing co-expression networks were developed in the microarray era, before high-throughput sequencing-with its unique statistical properties-became the norm for expression measurement. Here we develop Bayesian Relevance Networks, an algorithm that uses Bayesian reasoning about expression levels to account for the differing levels of uncertainty in expression measurements between highly- and lowly-expressed entities, and between samples with different sequencing depths. It combines data from groups of samples (e.g., replicates) to estimate group expression levels and confidence ranges. It then computes uncertainty-moderated estimates of cross-group correlations between entities, and uses permutation testing to assess their statistical significance. Using large scale miRNA data from The Cancer Genome Atlas, we show that our Bayesian update of the classical Relevance Networks algorithm provides improved reproducibility in co-expression estimates and lower false discovery rates in the resulting co-expression networks. Software is available at www.perkinslab.ca.

E3 ubiquitin ligase Mule targets ß-catenin under conditions of hyperactive Wnt signaling.

Wnt signaling, named after the secreted proteins that bind to cell surface receptors to activate the pathway, plays critical roles both in embryonic development and the maintenance of homeostasis in many adult tissues. Two particularly important cellular programs orchestrated by Wnt signaling are proliferation and stem cell self-renewal. Constitutive activation of the Wnt pathway resulting from mutation or improper modulation of pathway components contributes to cancer development in various tissues. Colon cancers frequently bear inactivating mutations of the adenomatous polyposis coli (APC) gene, whose product is an important component of the destruction complex that regulates ß-catenin levels. Stabilization and nuclear localization of ß-catenin result in the expression of a panel of Wnt target genes. We previously showed that Mule/Huwe1/Arf-BP1 (Mule) controls murine intestinal stem and progenitor cell proliferation by modulating the Wnt pathway via c-Myc. Here we extend our investigation of Mule's influence on oncogenesis by showing that Mule interacts directly with ß-catenin and targets it for degradation under conditions of hyperactive Wnt signaling. Our findings suggest that Mule uses various mechanisms to fine-tune the Wnt pathway and provides multiple safeguards against tumorigenesis.

Bayesian Correlation Analysis for Sequence Count Data.

Evaluating the similarity of different measured variables is a fundamental task of statistics, and a key part of many bioinformatics algorithms. Here we propose a Bayesian scheme for estimating the correlation between different entities' measurements based on high-throughput sequencing data. These entities could be different genes or miRNAs whose expression is measured by RNA-seq, different transcription factors or histone marks whose expression is measured by ChIP-seq, or even combinations of different types of entities. Our Bayesian formulation accounts for both measured signal levels and uncertainty in those levels, due to varying sequencing depth in different experiments and to varying absolute levels of individual entities, both of which affect the precision of the measurements. In comparison with a traditional Pearson correlation analysis, we show that our Bayesian correlation analysis retains high correlations when measurement confidence is high, but suppresses correlations when measurement confidence is low-especially for entities with low signal levels. In addition, we consider the influence of priors on the Bayesian correlation estimate. Perhaps surprisingly, we show that naive, uniform priors on entities' signal levels can lead to highly biased correlation estimates, particularly when different experiments have widely varying sequencing depths. However, we propose two alternative priors that provably mitigate this problem. We also prove that, like traditional Pearson correlation, our Bayesian correlation calculation constitutes a kernel in the machine learning sense, and thus can be used as a similarity measure in any kernel-based machine learning algorithm. We demonstrate our approach on two RNA-seq datasets and one miRNA-seq dataset.

BIDCHIPS: bias decomposition and removal from ChIP-seq data clarifies true binding signal and its functional correlates.

BACKGROUND: Unraveling transcriptional regulatory networks is a central problem in molecular biology and, in this quest, chromatin immunoprecipitation and sequencing (ChIP-seq) technology has given us the unprecedented ability to identify sites of protein-DNA binding and histone modification genome wide. However, multiple systemic and procedural biases hinder harnessing the full potential of this technology. Previous studies have addressed this problem, but a thorough characterization of different, interacting biases on ChIP-seq signals is still lacking. RESULTS: Here, we present a novel framework where the genome-wide ChIP-seq signal is viewed as being quantifiably influenced by different, measurable sources of bias, which can then be computationally subtracted away. We use a compendium of 123 human ENCODE ChIP-seq datasets to build regression models that tell us how much of a ChIP-seq signal can be attributed to mappability, GC-content, chromatin accessibility, and factors represented in input DNA and IgG controls. When we use the model to separate out these non-binding influences from the ChIP-seq signal, we obtain a purified signal that associates better to TF-DNA-binding motifs than do other measures of peak significance. We also carry out a multiscale analysis that reveals how ChIP-seq signal biases differ across different scales. Finally, we investigate previously reported associations between gene expression and ChIP-seq signals at transcription start sites. We show that our model can be used to discriminate ChIP-seq signals that are truly related to gene expression from those that are merely correlated by virtue of bias-in particular, chromatin accessibility bias, which shows up in ChIP-seq signals and also relates to gene expression. CONCLUSIONS: Our study provides new insights into the behavior of ChIP-seq signal biases and proposes a novel mitigation framework that improves results compared to existing techniques. With ChIP-seq now being the central technology for studying transcriptional regulation, it is most crucial to accurately characterize, quantify, and adjust for the genome-wide effects of biases affecting ChIP-seq. Our study also emphasizes that properly accounting for confounders in ChIP-seq data is of paramount importance for obtaining biologically accurate insights into the workings of the complex regulatory mechanisms in living organisms. R and MATLAB packages implementing the framework can be obtained from http://www.perkinslab.ca/Software.html.

Chromatin tandem affinity purification sequencing.

Chromatin immunoprecipitation coupled with ultra-high-throughput sequencing (ChIP-seq) is a widely used method for mapping the interactions of proteins with DNA. However, the requirements for ChIP-grade antibodies impede wider application of this method, and variations in results can be high owing to differences in affinity and cross-reactivity of antibodies. Therefore, we developed chromatin tandem affinity purification (ChTAP) as an effective alternative to ChIP. Through the use of affinity tags and reagents that are identical for all proteins investigated, ChTAP enables one to directly compare the binding between different transcription factors and to directly assess the background in control experiments. Thus, ChTAP-seq can be used to rapidly map the genome-wide binding of multiple DNA-binding proteins in a wide range of cell types. ChTAP can be completed in 3-4 d, starting from cross-linking of chromatin to purification of ChIP DNA.

MaSC: mappability-sensitive cross-correlation for estimating mean fragment length of single-end short-read sequencing data.

MOTIVATION: Reliable estimation of the mean fragment length for next-generation short-read sequencing data is an important step in next-generation sequencing analysis pipelines, most notably because of its impact on the accuracy of the enriched regions identified by peak-calling algorithms. Although many peak-calling algorithms include a fragment-length estimation subroutine, the problem has not been adequately solved, as demonstrated by the variability of the estimates returned by different algorithms. RESULTS: In this article, we investigate the use of strand cross-correlation to estimate mean fragment length of single-end data and show that traditional estimation approaches have mixed reliability. We observe that the mappability of different parts of the genome can introduce an artificial bias into cross-correlation computations, resulting in incorrect fragment-length estimates. We propose a new approach, called mappability-sensitive cross-correlation (MaSC), which removes this bias and allows for accurate and reliable fragment-length estimation. We analyze the computational complexity of this approach, and evaluate its performance on a test suite of NGS datasets, demonstrating its superiority to traditional cross-correlation analysis. AVAILABILITY: An open-source Perl implementation of our approach is available at http://www.perkinslab.ca/Software.html.

Adaptive bandwidth kernel density estimation for next-generation sequencing data.

BACKGROUND: High-throughput sequencing experiments can be viewed as measuring some sort of a "genomic signal" that may represent a biological event such as the binding of a transcription factor to the genome, locations of chromatin modifications, or even a background or control condition. Numerous algorithms have been developed to extract different kinds of information from such data. However, there has been very little focus on the reconstruction of the genomic signal itself. Such reconstructions may be useful for a variety of purposes ranging from simple visualization of the signals to sophisticated comparison of different datasets. METHODS: Here, we propose that adaptive-bandwidth kernel density estimators are well-suited for genomic signal reconstructions. This class of estimators is a natural extension of the fixed-bandwidth estimators that have been employed in several existing ChIP-Seq analysis programs. RESULTS: Using a set of ChIP-Seq datasets from the ENCODE project, we show that adaptive-bandwidth estimators have greater accuracy at signal reconstruction compared to fixed-bandwidth estimators, and that they have significant advantages in terms of visualization as well. For both fixed and adaptive-bandwidth schemes, we demonstrate that smoothing parameters can be set automatically using a held-out set of tuning data. We also carry out a computational complexity analysis of the different schemes and confirm through experimentation that the necessary computations can be readily carried out on a modern workstation without any significant issues.

Filter-based methodology for the location of hot spots in proteins and exons in DNA.

The so-called receiver operating characteristic technique is used as a tool in an optimization procedure for the improvement and assessment of a filter-based methodology for the location of hot spots in protein sequences and exons in DNA sequences. By optimizing the characteristic values of the nucleotides, high efficiency as well as improved accuracy can be achieved relative to results obtained with the electron-ion interaction potentials. On the other hand, by using the proposed filter-based methodology with binary sequences, improved accuracy can be achieved although the efficiency is somewhat compromised relative to that achieved using the optimized characteristic values. Extensive experimental results, evaluated using measures such as the g-mean, the Matthews correlation coefficient, and the chi-square statistic, show that the filter-based methodology performs much better than existing techniques using the short-time discrete Fourier transform, particularly in applications where short exons are involved.

Identification of tubulin drug binding sites and prediction of relative differences in binding affinities to tubulin isotypes using digital signal processing.

Microtubules are involved in numerous cellular processes including chromosome segregation during mitosis and, as a result, their constituent protein, tubulin, has become a successful target of several chemotherapeutic drugs. In general, these drugs bind indiscriminately to tubulin within both cancerous and healthy cells, resulting in unwanted side effects. However, differences between beta-tubulin isotypes expressed in a wide range of cell types may aid in the development of anti-tubulin drugs having increased specificity for only certain types of cells. Here, we describe a digital signal processing (DSP) method that is capable of predicting hot spots for the tubulin family of proteins as well as determining relative differences in binding affinities to these hot spots based only on the primary sequence of 10 human tubulin isotypes. Due to the fact that several drug binding sites have already been characterized within beta-tubulin, we are able to correlate hot spots with the binding sites for known chemotherapy drugs. We have also verified the accuracy of this method using the correlation between the binding affinities of characterized drugs and the tubulin isotypes. Additionally, the DSP method enables the rapid estimation of relative differences in binding affinities within the binding sites of tubulin isotypes that are yet to be experimentally determined.