PRAM: a novel pooling approach for discovering intergenic transcripts from large-scale RNA sequencing experiments.

Publicly available RNA-seq data is routinely used for retrospective analysis to elucidate new biology. Novel transcript discovery enabled by joint analysis of large collections of RNA-seq data sets has emerged as one such analysis. Current methods for transcript discovery rely on a '2-Step' approach where the first step encompasses building transcripts from individual data sets, followed by the second step that merges predicted transcripts across data sets. To increase the power of transcript discovery from large collections of RNA-seq data sets, we developed a novel '1-Step' approach named Pooling RNA-seq and Assembling Models (PRAM) that builds transcript models from pooled RNA-seq data sets. We demonstrate in a computational benchmark that 1-Step outperforms 2-Step approaches in predicting overall transcript structures and individual splice junctions, while performing competitively in detecting exonic nucleotides. Applying PRAM to 30 human ENCODE RNA-seq data sets identified unannotated transcripts with epigenetic and RAMPAGE signatures similar to those of recently annotated transcripts. In a case study, we discovered and experimentally validated new transcripts through the application of PRAM to mouse hematopoietic RNA-seq data sets. We uncovered new transcripts that share a differential expression pattern with a neighboring gene Pik3cg implicated in human hematopoietic phenotypes, and we provided evidence for the conservation of this relationship in human. PRAM is implemented as an R/Bioconductor package.

Hematopoietic Signaling Mechanism Revealed from a Stem/Progenitor Cell Cistrome.

Thousands of cis-elements in genomes are predicted to have vital functions. Although conservation, activity in surrogate assays, polymorphisms, and disease mutations provide functional clues, deletion from endogenous loci constitutes the gold-standard test. A GATA-2-binding, Gata2 intronic cis-element (+9.5) required for hematopoietic stem cell genesis in mice is mutated in a human immunodeficiency syndrome. Because +9.5 is the only cis-element known to mediate stem cell genesis, we devised a strategy to identify functionally comparable enhancers ("+9.5-like") genome-wide. Gene editing revealed +9.5-like activity to mediate GATA-2 occupancy, chromatin opening, and transcriptional activation. A +9.5-like element resided in Samd14, which encodes a protein of unknown function. Samd14 increased hematopoietic progenitor levels/activity and promoted signaling by a pathway vital for hematopoietic stem/progenitor cell regulation (stem cell factor/c-Kit), and c-Kit rescued Samd14 loss-of-function phenotypes. Thus, the hematopoietic stem/progenitor cell cistrome revealed a mediator of a signaling pathway that has broad importance for stem/progenitor cell biology.

RNA-regulatory exosome complex confers cellular survival to promote erythropoiesis.

Cellular differentiation requires vast remodeling of transcriptomes, and therefore machinery mediating remodeling controls differentiation. Relative to transcriptional mechanisms governing differentiation, post-transcriptional processes are less well understood. As an important post-transcriptional determinant of transcriptomes, the RNA exosome complex (EC) mediates processing and/or degradation of select RNAs. During erythropoiesis, the erythroid transcription factor GATA1 represses EC subunit genes. Depleting EC structural subunits prior to GATA1-mediated repression is deleterious to erythroid progenitor cells. To assess the importance of the EC catalytic subunits Dis3 and Exosc10 in this dynamic process, we asked if these subunits function non-redundantly to control erythropoiesis. Dis3 or Exosc10 depletion in primary murine hematopoietic progenitor cells reduced erythroid progenitors and their progeny, while sparing myeloid cells. Dis3 loss severely compromised erythroid progenitor and erythroblast survival, rendered erythroblasts hypersensitive to apoptosis-inducing stimuli and induced γ-H2AX, indicative of DNA double-stranded breaks. Dis3 loss-of-function phenotypes were more severe than those caused by Exosc10 depletion. We innovated a genetic rescue system to compare human Dis3 with multiple myeloma-associated Dis3 mutants S447R and R750K, and only wild type Dis3 was competent to rescue progenitors. Thus, Dis3 establishes a disease mutation-sensitive, cell type-specific survival mechanism to enable a differentiation program.

Genome-wide effects on Escherichia coli transcription from ppGpp binding to its two sites on RNA polymerase.

The second messenger nucleotide ppGpp dramatically alters gene expression in bacteria to adjust cellular metabolism to nutrient availability. ppGpp binds to two sites on RNA polymerase (RNAP) in Escherichia coli, but it has also been reported to bind to many other proteins. To determine the role of the RNAP binding sites in the genome-wide effects of ppGpp on transcription, we used RNA-seq to analyze transcripts produced in response to elevated ppGpp levels in strains with/without the ppGpp binding sites on RNAP. We examined RNAs rapidly after ppGpp production without an accompanying nutrient starvation. This procedure enriched for direct effects of ppGpp on RNAP rather than for indirect effects on transcription resulting from starvation-induced changes in metabolism or on secondary events from the initial effects on RNAP. The transcriptional responses of all 757 genes identified after 5 minutes of ppGpp induction depended on ppGpp binding to RNAP. Most (>75%) were not reported in earlier studies. The regulated transcripts encode products involved not only in translation but also in many other cellular processes. In vitro transcription analysis of more than 100 promoters from the in vivo dataset identified a large collection of directly regulated promoters, unambiguously demonstrated that most effects of ppGpp on transcription in vivo were direct, and allowed comparison of DNA sequences from inhibited, activated, and unaffected promoter classes. Our analysis greatly expands our understanding of the breadth of the stringent response and suggests promoter sequence features that contribute to the specific effects of ppGpp.

Integrative analysis with ChIP-seq advances the limits of transcript quantification from RNA-seq.

RNA-seq is currently the technology of choice for global measurement of transcript abundances in cells. Despite its successes, isoform-level quantification remains difficult because short RNA-seq reads are often compatible with multiple alternatively spliced isoforms. Existing methods rely heavily on uniquely mapping reads, which are not available for numerous isoforms that lack regions of unique sequence. To improve quantification accuracy in such difficult cases, we developed a novel computational method, prior-enhanced RSEM (pRSEM), which uses a complementary data type in addition to RNA-seq data. We found that ChIP-seq data of RNA polymerase II and histone modifications were particularly informative in this approach. In qRT-PCR validations, pRSEM was shown to be superior than competing methods in estimating relative isoform abundances within or across conditions. Data-driven simulations suggested that pRSEM has a greatly decreased false-positive rate at the expense of a small increase in false-negative rate. In aggregate, our study demonstrates that pRSEM transforms existing capacity to precisely estimate transcript abundances, especially at the isoform level.

Zebrafish zic2 controls formation of periocular neural crest and choroid fissure morphogenesis.

The vertebrate retina develops in close proximity to the forebrain and neural crest-derived cartilages of the face and jaw. Coloboma, a congenital eye malformation, is associated with aberrant forebrain development (holoprosencephaly) and with craniofacial defects (frontonasal dysplasia) in humans, suggesting a critical role for cross-lineage interactions during retinal morphogenesis. ZIC2, a zinc-finger transcription factor, is linked to human holoprosencephaly. We have previously used morpholino assays to show zebrafish zic2 functions in the developing forebrain, retina and craniofacial cartilage. We now report that zebrafish with genetic lesions in zebrafish zic2 orthologs, zic2a and zic2b, develop with retinal coloboma and craniofacial anomalies. We demonstrate a requirement for zic2 in restricting pax2a expression and show evidence that zic2 function limits Hh signaling. RNA-seq transcriptome analysis identified an early requirement for zic2 in periocular neural crest as an activator of alx1, a transcription factor with essential roles in craniofacial and ocular morphogenesis in human and zebrafish. Collectively, these data establish zic2 mutant zebrafish as a powerful new genetic model for in-depth dissection of cell interactions and genetic controls during craniofacial complex development.

Analysis of embryonic development in the unsequenced axolotl: Waves of transcriptomic upheaval and stability.

The axolotl (Ambystoma mexicanum) has long been the subject of biological research, primarily owing to its outstanding regenerative capabilities. However, the gene expression programs governing its embryonic development are particularly underexplored, especially when compared to other amphibian model species. Therefore, we performed whole transcriptome polyA+ RNA sequencing experiments on 17 stages of embryonic development. As the axolotl genome is unsequenced and its gene annotation is incomplete, we built de novo transcriptome assemblies for each stage and garnered functional annotation by comparing expressed contigs with known genes in other organisms. In evaluating the number of differentially expressed genes over time, we identify three waves of substantial transcriptome upheaval each followed by a period of relative transcriptome stability. The first wave of upheaval is between the one and two cell stage. We show that the number of differentially expressed genes per unit time is higher between the one and two cell stage than it is across the mid-blastula transition (MBT), the period of zygotic genome activation. We use total RNA sequencing to demonstrate that the vast majority of genes with increasing polyA+ signal between the one and two cell stage result from polyadenylation rather than de novo transcription. The first stable phase begins after the two cell stage and continues until the mid-blastula transition, corresponding with the pre-MBT phase of transcriptional quiescence in amphibian development. Following this is a peak of differential gene expression corresponding with the activation of the zygotic genome and a phase of transcriptomic stability from stages 9-11. We observe a third wave of transcriptomic change between stages 11 and 14, followed by a final stable period. The last two stable phases have not been documented in amphibians previously and correspond to times of major morphogenic change in the axolotl embryo: gastrulation and neurulation. These results yield new insights into global gene expression during early stages of amphibian embryogenesis and will help to further develop the axolotl as a model species for developmental and regenerative biology.

MetaSRA: normalized human sample-specific metadata for the Sequence Read Archive.

MOTIVATION: The NCBI's Sequence Read Archive (SRA) promises great biological insight if one could analyze the data in the aggregate; however, the data remain largely underutilized, in part, due to the poor structure of the metadata associated with each sample. The rules governing submissions to the SRA do not dictate a standardized set of terms that should be used to describe the biological samples from which the sequencing data are derived. As a result, the metadata include many synonyms, spelling variants and references to outside sources of information. Furthermore, manual annotation of the data remains intractable due to the large number of samples in the archive. For these reasons, it has been difficult to perform large-scale analyses that study the relationships between biomolecular processes and phenotype across diverse diseases, tissues and cell types present in the SRA. RESULTS: We present MetaSRA, a database of normalized SRA human sample-specific metadata following a schema inspired by the metadata organization of the ENCODE project. This schema involves mapping samples to terms in biomedical ontologies, labeling each sample with a sample-type category, and extracting real-valued properties. We automated these tasks via a novel computational pipeline. AVAILABILITY AND IMPLEMENTATION: The MetaSRA is available at metasra.biostat.wisc.edu via both a searchable web interface and bulk downloads. Software implementing our computational pipeline is available at http://github.com/deweylab/metasra-pipeline. CONTACT: cdewey@biostat.wisc.edu. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Mechanism governing heme synthesis reveals a GATA factor/heme circuit that controls differentiation.

Metal ion-containing macromolecules have fundamental roles in essentially all biological processes throughout the evolutionary tree. For example, iron-containing heme is a cofactor in enzyme catalysis and electron transfer and an essential hemoglobin constituent. To meet the intense demand for hemoglobin assembly in red blood cells, the cell type-specific factor GATA-1 activates transcription of Alas2, encoding the rate-limiting enzyme in heme biosynthesis, 5-aminolevulinic acid synthase-2 (ALAS-2). Using genetic editing to unravel mechanisms governing heme biosynthesis, we discovered a GATA factor- and heme-dependent circuit that establishes the erythroid cell transcriptome. CRISPR/Cas9-mediated ablation of two Alas2 intronic cis elements strongly reduces GATA-1-induced Alas2 transcription, heme biosynthesis, and surprisingly, GATA-1 regulation of other vital constituents of the erythroid cell transcriptome. Bypassing ALAS-2 function in Alas2 cis element-mutant cells by providing its catalytic product 5-aminolevulinic acid rescues heme biosynthesis and the GATA-1-dependent genetic network. Heme amplifies GATA-1 function by downregulating the heme-sensing transcriptional repressor Bach1 and via a Bach1-insensitive mechanism. Through this dual mechanism, heme and a master regulator collaborate to orchestrate a cell type-specific transcriptional program that promotes cellular differentiation.

EBSeq-HMM: a Bayesian approach for identifying gene-expression changes in ordered RNA-seq experiments.

MOTIVATION: With improvements in next-generation sequencing technologies and reductions in price, ordered RNA-seq experiments are becoming common. Of primary interest in these experiments is identifying genes that are changing over time or space, for example, and then characterizing the specific expression changes. A number of robust statistical methods are available to identify genes showing differential expression among multiple conditions, but most assume conditions are exchangeable and thereby sacrifice power and precision when applied to ordered data. RESULTS: We propose an empirical Bayes mixture modeling approach called EBSeq-HMM. In EBSeq-HMM, an auto-regressive hidden Markov model is implemented to accommodate dependence in gene expression across ordered conditions. As demonstrated in simulation and case studies, the output proves useful in identifying differentially expressed genes and in specifying gene-specific expression paths. EBSeq-HMM may also be used for inference regarding isoform expression. AVAILABILITY AND IMPLEMENTATION: An R package containing examples and sample datasets is available at Bioconductor. CONTACT: kendzior@biostat.wisc.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Perm-seq: Mapping Protein-DNA Interactions in Segmental Duplication and Highly Repetitive Regions of Genomes with Prior-Enhanced Read Mapping.

Segmental duplications and other highly repetitive regions of genomes contribute significantly to cells' regulatory programs. Advancements in next generation sequencing enabled genome-wide profiling of protein-DNA interactions by chromatin immunoprecipitation followed by high throughput sequencing (ChIP-seq). However, interactions in highly repetitive regions of genomes have proven difficult to map since short reads of 50-100 base pairs (bps) from these regions map to multiple locations in reference genomes. Standard analytical methods discard such multi-mapping reads and the few that can accommodate them are prone to large false positive and negative rates. We developed Perm-seq, a prior-enhanced read allocation method for ChIP-seq experiments, that can allocate multi-mapping reads in highly repetitive regions of the genomes with high accuracy. We comprehensively evaluated Perm-seq, and found that our prior-enhanced approach significantly improves multi-read allocation accuracy over approaches that do not utilize additional data types. The statistical formalism underlying our approach facilitates supervising of multi-read allocation with a variety of data sources including histone ChIP-seq. We applied Perm-seq to 64 ENCODE ChIP-seq datasets from GM12878 and K562 cells and identified many novel protein-DNA interactions in segmental duplication regions. Our analysis reveals that although the protein-DNA interactions sites are evolutionarily less conserved in repetitive regions, they share the overall sequence characteristics of the protein-DNA interactions in non-repetitive regions.

Bicaudal-C spatially controls translation of vertebrate maternal mRNAs.

The Xenopus Cripto-1 protein is confined to the cells of the animal hemisphere during early embryogenesis where it regulates the formation of anterior structures. Cripto-1 protein accumulates only in animal cells because cripto-1 mRNA in cells of the vegetal hemisphere is translationally repressed. Here, we show that the RNA binding protein, Bicaudal-C (Bic-C), functioned directly in this vegetal cell-specific repression. While Bic-C protein is normally confined to vegetal cells, ectopic expression of Bic-C in animal cells repressed a cripto-1 mRNA reporter and associated with endogenous cripto-1 mRNA. Repression by Bic-C required its N-terminal domain, comprised of multiple KH motifs, for specific binding to relevant control elements within the cripto-1 mRNA and a functionally separable C-terminal translation repression domain. Bic-C-mediated repression required the 5' CAP and translation initiation factors, but not a poly(A) tail or the conserved SAM domain within Bic-C. Bic-C-directed immunoprecipitation followed by deep sequencing of associated mRNAs identified multiple Bic-C-regulated mRNA targets, including cripto-1 mRNA, providing new insights and tools for understanding the role of Bic-C in vertebrate development.

Rbm20 regulates titin alternative splicing as a splicing repressor.

Titin, a sarcomeric protein expressed primarily in striated muscles, is responsible for maintaining the structure and biomechanical properties of muscle cells. Cardiac titin undergoes developmental size reduction from 3.7 megadaltons in neonates to primarily 2.97 megadaltons in the adult. This size reduction results from gradually increased exon skipping between exons 50 and 219 of titin mRNA. Our previous study reported that Rbm20 is the splicing factor responsible for this process. In this work, we investigated its molecular mechanism. We demonstrate that Rbm20 mediates exon skipping by binding to titin pre-mRNA to repress the splicing of some regions; the exons/introns in these Rbm20-repressed regions are ultimately skipped. Rbm20 was also found to mediate intron retention and exon shuffling. The two Rbm20 speckles found in nuclei from muscle tissues were identified as aggregates of Rbm20 protein on the partially processed titin pre-mRNAs. Cooperative repression and alternative 3' splice site selection were found to be used by Rbm20 to skip different subsets of titin exons, and the splicing pathway selected depended on the ratio of Rbm20 to other splicing factors that vary with tissue type and developmental age.

Inference of alternative splicing from RNA-Seq data with probabilistic splice graphs.

MOTIVATION: Alternative splicing and other processes that allow for different transcripts to be derived from the same gene are significant forces in the eukaryotic cell. RNA-Seq is a promising technology for analyzing alternative transcripts, as it does not require prior knowledge of transcript structures or genome sequences. However, analysis of RNA-Seq data in the presence of genes with large numbers of alternative transcripts is currently challenging due to efficiency, identifiability and representation issues. RESULTS: We present RNA-Seq models and associated inference algorithms based on the concept of probabilistic splice graphs, which alleviate these issues. We prove that our models are often identifiable and demonstrate that our inference methods for quantification and differential processing detection are efficient and accurate. AVAILABILITY: Software implementing our methods is available at http://deweylab.biostat.wisc.edu/psginfer. CONTACT: cdewey@biostat.wisc.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Comparative RNA-seq analysis in the unsequenced axolotl: the oncogene burst highlights early gene expression in the blastema.

The salamander has the remarkable ability to regenerate its limb after amputation. Cells at the site of amputation form a blastema and then proliferate and differentiate to regrow the limb. To better understand this process, we performed deep RNA sequencing of the blastema over a time course in the axolotl, a species whose genome has not been sequenced. Using a novel comparative approach to analyzing RNA-seq data, we characterized the transcriptional dynamics of the regenerating axolotl limb with respect to the human gene set. This approach involved de novo assembly of axolotl transcripts, RNA-seq transcript quantification without a reference genome, and transformation of abundances from axolotl contigs to human genes. We found a prominent burst in oncogene expression during the first day and blastemal/limb bud genes peaking at 7 to 14 days. In addition, we found that limb patterning genes, SALL genes, and genes involved in angiogenesis, wound healing, defense/immunity, and bone development are enriched during blastema formation and development. Finally, we identified a category of genes with no prior literature support for limb regeneration that are candidates for further evaluation based on their expression pattern during the regenerative process.

Infection induced inflammation impairs wound healing through IL-1ß signaling.

Wound healing is impaired by infection; however, how microbe-induced inflammation modulates tissue repair remains unclear. We took advantage of the optical transparency of zebrafish and a genetically tractable microbe, Listeria monocytogenes, to probe the role of infection and inflammation in wound healing. Infection with bacteria engineered to activate the inflammasome, Lm-Pyro, induced persistent inflammation and impaired healing despite low bacterial burden. Inflammatory infections induced il1b expression and blocking IL-1R signaling partially rescued wound healing in the presence of persistent infection. We found a critical window of microbial clearance necessary to limit persistent inflammation and enable efficient wound repair. Taken together, our findings suggest that the dynamics of microbe-induced tissue inflammation impacts repair in complex tissue damage independent of bacterial load, with a critical early window for efficient tissue repair.

Positional orthology: putting genomic evolutionary relationships into context.

Orthology is a powerful refinement of homology that allows us to describe more precisely the evolution of genomes and understand the function of the genes they contain. However, because orthology is not concerned with genomic position, it is limited in its ability to describe genes that are likely to have equivalent roles in different genomes. Because of this limitation, the concept of 'positional orthology' has emerged, which describes the relation between orthologous genes that retain their ancestral genomic positions. In this review, we formally define this concept, for which we introduce the shorter term 'toporthology', with respect to the evolutionary events experienced by a gene's ancestors. Through a discussion of recent studies on the role of genomic context in gene evolution, we show that the distinction between orthology and toporthology is biologically significant. We then review a number of orthology prediction methods that take genomic context into account and thus that may be used to infer the important relation of toporthology.

Discovery of functional elements in 12 Drosophila genomes using evolutionary signatures.

Sequencing of multiple related species followed by comparative genomics analysis constitutes a powerful approach for the systematic understanding of any genome. Here, we use the genomes of 12 Drosophila species for the de novo discovery of functional elements in the fly. Each type of functional element shows characteristic patterns of change, or 'evolutionary signatures', dictated by its precise selective constraints. Such signatures enable recognition of new protein-coding genes and exons, spurious and incorrect gene annotations, and numerous unusual gene structures, including abundant stop-codon readthrough. Similarly, we predict non-protein-coding RNA genes and structures, and new microRNA (miRNA) genes. We provide evidence of miRNA processing and functionality from both hairpin arms and both DNA strands. We identify several classes of pre- and post-transcriptional regulatory motifs, and predict individual motif instances with high confidence. We also study how discovery power scales with the divergence and number of species compared, and we provide general guidelines for comparative studies.

RNA-regulatory exosome complex suppresses an apoptotic program to confer erythroid progenitor cell survival in vivo.

The RNA-regulatory exosome complex (EC) posttranscriptionally and cotranscriptionally processes and degrades RNAs in a context-dependent manner. Although the EC functions in diverse cell types, its contributions to stem and progenitor cell development are not well understood. Previously, we demonstrated that the transcriptional regulator of erythrocyte development, GATA1, represses EC subunit genes, and the EC maintains erythroid progenitors in vitro. To determine if this mechanism operates in vivo, we used the hematopoietic-specific Vav1-Cre and "conditional by inversion" mouse system to ablate Exosc3, encoding an EC structural subunit. Although Exosc3C/C Cre+ embryos developed normally until embryonic day 14.5, Exosc3 ablation was embryonic lethal and severely reduced erythromyeloid progenitor activity. RNA sequencing analysis of Exosc3-ablated burst-forming unit-erythroid revealed elevated transcripts encoding multiple proapoptotic factors, and the mutant erythroid progenitors exhibited increased apoptosis. We propose that the EC controls an ensemble of apoptosis-regulatory RNAs, thereby promoting erythroid progenitor survival and developmental erythropoiesis in vivo.

Fine-scale phylogenetic discordance across the house mouse genome.

Population genetic theory predicts discordance in the true phylogeny of different genomic regions when studying recently diverged species. Despite this expectation, genome-wide discordance in young species groups has rarely been statistically quantified. The house mouse subspecies group provides a model system for examining phylogenetic discordance. House mouse subspecies are recently derived, suggesting that even if there has been a simple tree-like population history, gene trees could disagree with the population history due to incomplete lineage sorting. Subspecies of house mice also hybridize in nature, raising the possibility that recent introgression might lead to additional phylogenetic discordance. Single-locus approaches have revealed support for conflicting topologies, resulting in a subspecies tree often summarized as a polytomy. To analyze phylogenetic histories on a genomic scale, we applied a recently developed method, Bayesian concordance analysis, to dense SNP data from three closely related subspecies of house mice: Mus musculus musculus, M. m. castaneus, and M. m. domesticus. We documented substantial variation in phylogenetic history across the genome. Although each of the three possible topologies was strongly supported by a large number of loci, there was statistical evidence for a primary phylogenetic history in which M. m. musculus and M. m. castaneus are sister subspecies. These results underscore the importance of measuring phylogenetic discordance in other recently diverged groups using methods such as Bayesian concordance analysis, which are designed for this purpose.