Benchmarking long-read RNA-sequencing analysis tools using in silico mixtures.

The lack of benchmark data sets with inbuilt ground-truth makes it challenging to compare the performance of existing long-read isoform detection and differential expression analysis workflows. Here, we present a benchmark experiment using two human lung adenocarcinoma cell lines that were each profiled in triplicate together with synthetic, spliced, spike-in RNAs (sequins). Samples were deeply sequenced on both Illumina short-read and Oxford Nanopore Technologies long-read platforms. Alongside the ground-truth available via the sequins, we created in silico mixture samples to allow performance assessment in the absence of true positives or true negatives. Our results show that StringTie2 and bambu outperformed other tools from the six isoform detection tools tested, DESeq2, edgeR and limma-voom were best among the five differential transcript expression tools tested and there was no clear front-runner for performing differential transcript usage analysis between the five tools compared, which suggests further methods development is needed for this application.

A functional genetic screen identifies aurora kinase b as an essential regulator of Sox9-positive mouse embryonic lung progenitor cells.

Development of a branching tree in the embryonic lung is crucial for the formation of a fully mature functional lung at birth. Sox9+ cells present at the tip of the primary embryonic lung endoderm are multipotent cells responsible for branch formation and elongation. We performed a genetic screen in murine primary cells and identified aurora kinase b (Aurkb) as an essential regulator of Sox9+ cells ex vivo. In vivo conditional knockout studies confirmed that Aurkb was required for lung development but was not necessary for postnatal growth and the repair of the adult lung after injury. Deletion of Aurkb in embryonic Sox9+ cells led to the formation of a stunted lung that retained the expression of Sox2 in the proximal airways, as well as Sox9 in the distal tips. Although we found no change in cell polarity, we showed that loss of Aurkb or chemical inhibition of Aurkb caused Sox9+ cells to arrest at G2/M, likely responsible for the lack of branch bifurcation. This work demonstrates the power of genetic screens in identifying novel regulators of Sox9+ progenitor cells and lung branching morphogenesis.

Toward a gold standard for benchmarking gene set enrichment analysis.

MOTIVATION: Although gene set enrichment analysis has become an integral part of high-throughput gene expression data analysis, the assessment of enrichment methods remains rudimentary and ad hoc. In the absence of suitable gold standards, evaluations are commonly restricted to selected datasets and biological reasoning on the relevance of resulting enriched gene sets. RESULTS: We develop an extensible framework for reproducible benchmarking of enrichment methods based on defined criteria for applicability, gene set prioritization and detection of relevant processes. This framework incorporates a curated compendium of 75 expression datasets investigating 42 human diseases. The compendium features microarray and RNA-seq measurements, and each dataset is associated with a precompiled GO/KEGG relevance ranking for the corresponding disease under investigation. We perform a comprehensive assessment of 10 major enrichment methods, identifying significant differences in runtime and applicability to RNA-seq data, fraction of enriched gene sets depending on the null hypothesis tested and recovery of the predefined relevance rankings. We make practical recommendations on how methods originally developed for microarray data can efficiently be applied to RNA-seq data, how to interpret results depending on the type of gene set test conducted and which methods are best suited to effectively prioritize gene sets with high phenotype relevance. AVAILABILITY: http://bioconductor.org/packages/GSEABenchmarkeR. CONTACT: ludwig.geistlinger@sph.cuny.edu.

The impact of influenza pulmonary infection and inflammation on vagal bronchopulmonary sensory neurons.

Influenza A virus (IAV) is rapidly detected in the airways by the immune system, with resident parenchymal cells and leukocytes orchestrating viral sensing and the induction of antiviral inflammatory responses. The airways are innervated by heterogeneous populations of vagal sensory neurons which also play an important role in pulmonary defense. How these neurons respond to IAV respiratory infection remains unclear. Here, we use a murine model to provide the first evidence that vagal sensory neurons undergo significant transcriptional changes following a respiratory IAV infection. RNA sequencing on vagal sensory ganglia showed that IAV infection induced the expression of many genes associated with an antiviral and pro-inflammatory response and this was accompanied by a significant increase in inflammatory cell recruitment into the vagal ganglia. Assessment of gene expression in single-vagal sensory neurons confirmed that IAV infection induced a neuronal inflammatory phenotype, which was most prominent in bronchopulmonary neurons, and also evident in some neurons innervating other organs. The altered transcriptome could be mimicked by intranasal treatment with cytokines and the lung homogenates of infected mice, in the absence of infectious virus. These data argue that IAV pulmonary infection and subsequent inflammation induces vagal sensory ganglia neuroinflammation and this may have important implications for IAV-induced morbidity.

RNA-seq mixology: designing realistic control experiments to compare protocols and analysis methods.

Carefully designed control experiments provide a gold standard for benchmarking different genomics research tools. A shortcoming of many gene expression control studies is that replication involves profiling the same reference RNA sample multiple times. This leads to low, pure technical noise that is atypical of regular studies. To achieve a more realistic noise structure, we generated a RNA-sequencing mixture experiment using two cell lines of the same cancer type. Variability was added by extracting RNA from independent cell cultures and degrading particular samples. The systematic gene expression changes induced by this design allowed benchmarking of different library preparation kits (standard poly-A versus total RNA with Ribozero depletion) and analysis pipelines. Data generated using the total RNA kit had more signal for introns and various RNA classes (ncRNA, snRNA, snoRNA) and less variability after degradation. For differential expression analysis, voom with quality weights marginally outperformed other popular methods, while for differential splicing, DEXSeq was simultaneously the most sensitive and the most inconsistent method. For sample deconvolution analysis, DeMix outperformed IsoPure convincingly. Our RNA-sequencing data set provides a valuable resource for benchmarking different protocols and data pre-processing workflows. The extra noise mimics routine lab experiments more closely, ensuring any conclusions are widely applicable.

Glimma: interactive graphics for gene expression analysis.

MOTIVATION: graphics for RNA-sequencing and microarray gene expression analyses may contain upwards of tens of thousands of points. Details about certain genes or samples of interest are easily obscured in such dense summary displays. Incorporating interactivity into summary plots would enable additional information to be displayed on demand and facilitate intuitive data exploration. RESULTS: The open-source Glimma package creates interactive graphics for exploring gene expression analysis with a few simple R commands. It extends popular plots found in the limma package, such as multi-dimensional scaling plots and mean-difference plots, to allow individual data points to be queried and additional annotation information to be displayed upon hovering or selecting particular points. It also offers links between plots so that more information can be revealed on demand. Glimma is widely applicable, supporting data analyses from a number of well-established Bioconductor workflows ( limma , edgeR and DESeq2 ) and uses D3/JavaScript to produce HTML pages with interactive displays that enable more effective data exploration by end-users. Results from Glimma can be easily shared between bioinformaticians and biologists, enhancing reporting capabilities while maintaining reproducibility. AVAILABILITY AND IMPLEMENTATION: The Glimma R package is available from http://bioconductor.org/packages/Glimma/ . CONTACT: su.s@wehi.edu.au , law@wehi.edu.au or mritchie@wehi.edu.au.

The hematopoietic oncoprotein FOXP1 promotes tumor cell survival in diffuse large B-cell lymphoma by repressing S1PR2 signaling.

Aberrant expression of the oncogenic transcription factor forkhead box protein 1 (FOXP1) is a common feature of diffuse large B-cell lymphoma (DLBCL). We have combined chromatin immunoprecipitation and gene expression profiling after FOXP1 depletion with functional screening to identify targets of FOXP1 contributing to tumor cell survival. We find that the sphingosine-1-phosphate receptor 2 (S1PR2) is repressed by FOXP1 in activated B-cell (ABC) and germinal center B-cell (GCB) DLBCL cell lines with aberrantly high FOXP1 levels; S1PR2 expression is further inversely correlated with FOXP1 expression in 3 patient cohorts. Ectopic expression of wild-type S1PR2, but not a point mutant incapable of activating downstream signaling pathways, induces apoptosis in DLBCL cells and restricts tumor growth in subcutaneous and orthotopic models of the disease. The proapoptotic effects of S1PR2 are phenocopied by ectopic expression of the small G protein Gα13 but are independent of AKT signaling. We further show that low S1PR2 expression is a strong negative prognosticator of patient survival, alone and especially in combination with high FOXP1 expression. The S1PR2 locus has previously been demonstrated to be recurrently mutated in GCB DLBCL; the transcriptional silencing of S1PR2 by FOXP1 represents an alternative mechanism leading to inactivation of this important hematopoietic tumor suppressor.

limma powers differential expression analyses for RNA-sequencing and microarray studies.

limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.

An immunohistochemical atlas of necroptotic pathway expression.

Necroptosis is a lytic form of regulated cell death reported to contribute to inflammatory diseases of the gut, skin and lung, as well as ischemic-reperfusion injuries of the kidney, heart and brain. However, precise identification of the cells and tissues that undergo necroptotic cell death in vivo has proven challenging in the absence of robust protocols for immunohistochemical detection. Here, we provide automated immunohistochemistry protocols to detect core necroptosis regulators - Caspase-8, RIPK1, RIPK3 and MLKL - in formalin-fixed mouse and human tissues. We observed surprising heterogeneity in protein expression within tissues, whereby short-lived immune barrier cells were replete with necroptotic effectors, whereas long-lived cells lacked RIPK3 or MLKL expression. Local changes in the expression of necroptotic effectors occurred in response to insults such as inflammation, dysbiosis or immune challenge, consistent with necroptosis being dysregulated in disease contexts. These methods will facilitate the precise localisation and evaluation of necroptotic signaling in vivo.

Modeling group heteroscedasticity in single-cell RNA-seq pseudo-bulk data.

Group heteroscedasticity is commonly observed in pseudo-bulk single-cell RNA-seq datasets and its presence can hamper the detection of differentially expressed genes. Since most bulk RNA-seq methods assume equal group variances, we introduce two new approaches that account for heteroscedastic groups, namely voomByGroup and voomWithQualityWeights using a blocked design (voomQWB). Compared to current gold-standard methods that do not account for group heteroscedasticity, we show results from simulations and various experiments that demonstrate the superior performance of voomByGroup and voomQWB in terms of error control and power when group variances in pseudo-bulk single-cell RNA-seq data are unequal.

Dashboard-style interactive plots for RNA-seq analysis are R Markdown ready with Glimma 2.0.

Glimma 1.0 introduced intuitive, point-and-click interactive graphics for differential gene expression analysis. Here, we present a major update to Glimma that brings improved interactivity and reproducibility using high-level visualization frameworks for R and JavaScript. Glimma 2.0 plots are now readily embeddable in R Markdown, thus allowing users to create reproducible reports containing interactive graphics. The revamped multidimensional scaling plot features dashboard-style controls allowing the user to dynamically change the colour, shape and size of sample points according to different experimental conditions. Interactivity was enhanced in the MA-style plot for comparing differences to average expression, which now supports selecting multiple genes, export options to PNG, SVG or CSV formats and includes a new volcano plot function. Feature-rich and user-friendly, Glimma makes exploring data for gene expression analysis more accessible and intuitive and is available on Bioconductor and GitHub.

Homeostatic apoptosis prevents competition-induced atrophy in follicular B cells.

While the intrinsic apoptosis pathway is thought to play a central role in shaping the B cell lineage, its precise role in mature B cell homeostasis remains elusive. Using mice in which mature B cells are unable to undergo apoptotic cell death, we show that apoptosis constrains follicular B (FoB) cell lifespan but plays no role in marginal zone B (MZB) cell homeostasis. In these mice, FoB cells accumulate abnormally. This intensifies intercellular competition for BAFF, resulting in a contraction of the MZB cell compartment, and reducing the growth, trafficking, and fitness of FoB cells. Diminished BAFF signaling dampens the non-canonical NF-κB pathway, undermining FoB cell growth despite the concurrent triggering of a protective p53 response. Thus, MZB and FoB cells exhibit a differential requirement for the intrinsic apoptosis pathway. Homeostatic apoptosis constrains the size of the FoB cell compartment, thereby preventing competition-induced FoB cell atrophy.

The long and the short of it: unlocking nanopore long-read RNA sequencing data with short-read differential expression analysis tools.

Application of Oxford Nanopore Technologies' long-read sequencing platform to transcriptomic analysis is increasing in popularity. However, such analysis can be challenging due to the high sequence error and small library sizes, which decreases quantification accuracy and reduces power for statistical testing. Here, we report the analysis of two nanopore RNA-seq datasets with the goal of obtaining gene- and isoform-level differential expression information. A dataset of synthetic, spliced, spike-in RNAs ('sequins') as well as a mouse neural stem cell dataset from samples with a null mutation of the epigenetic regulator Smchd1 was analysed using a mix of long-read specific tools for preprocessing together with established short-read RNA-seq methods for downstream analysis. We used limma-voom to perform differential gene expression analysis, and the novel FLAMES pipeline to perform isoform identification and quantification, followed by DRIMSeq and limma-diffSplice (with stageR) to perform differential transcript usage analysis. We compared results from the sequins dataset to the ground truth, and results of the mouse dataset to a previous short-read study on equivalent samples. Overall, our work shows that transcriptomic analysis of long-read nanopore data using long-read specific preprocessing methods together with short-read differential expression methods and software that are already in wide use can yield meaningful results.

Comprehensive characterization of single-cell full-length isoforms in human and mouse with long-read sequencing.

A modified Chromium 10x droplet-based protocol that subsamples cells for both short-read and long-read (nanopore) sequencing together with a new computational pipeline (FLAMES) is developed to enable isoform discovery, splicing analysis, and mutation detection in single cells. We identify thousands of unannotated isoforms and find conserved functional modules that are enriched for alternative transcript usage in different cell types and species, including ribosome biogenesis and mRNA splicing. Analysis at the transcript level allows data integration with scATAC-seq on individual promoters, improved correlation with protein expression data, and linked mutations known to confer drug resistance to transcriptome heterogeneity.

A guide to creating design matrices for gene expression experiments.

Differential expression analysis of genomic data types, such as RNA-sequencing experiments, use linear models to determine the size and direction of the changes in gene expression. For RNA-sequencing, there are several established software packages for this purpose accompanied with analysis pipelines that are well described. However, there are two crucial steps in the analysis process that can be a stumbling block for many -- the set up an appropriate model via design matrices and the set up of comparisons of interest via contrast matrices. These steps are particularly troublesome because an extensive catalogue for design and contrast matrices does not currently exist. One would usually search for example case studies across different platforms and mix and match the advice from those sources to suit the dataset they have at hand. This article guides the reader through the basics of how to set up design and contrast matrices. We take a practical approach by providing code and graphical representation of each case study, starting with simpler examples (e.g. models with a single explanatory variable) and move onto more complex ones (e.g. interaction models, mixed effects models, higher order time series and cyclical models). Although our work has been written specifically with a limma-style pipeline in mind, most of it is also applicable to other software packages for differential expression analysis, and the ideas covered can be adapted to data analysis of other high-throughput technologies. Where appropriate, we explain the interpretation and differences between models to aid readers in their own model choices. Unnecessary jargon and theory is omitted where possible so that our work is accessible to a wide audience of readers, from beginners to those with experience in genomics data analysis.

Covering all your bases: incorporating intron signal from RNA-seq data.

RNA-seq datasets can contain millions of intron reads per library that are typically removed from downstream analysis. Only reads overlapping annotated exons are considered to be informative since mature mRNA is assumed to be the major component sequenced, especially for poly(A) RNA libraries. In this study, we show that intron reads are informative, and through exploratory data analysis of read coverage that intron signal is representative of both pre-mRNAs and intron retention. We demonstrate how intron reads can be utilized in differential expression analysis using our index method where a unique set of differentially expressed genes can be detected using intron counts. In exploring read coverage, we also developed the superintronic software that quickly and robustly calculates user-defined summary statistics for exonic and intronic regions. Across multiple datasets, superintronic enabled us to identify several genes with distinctly retained introns that had similar coverage levels to that of neighbouring exons. The work and ideas presented in this paper is the first of its kind to consider multiple biological sources for intron reads through exploratory data analysis, minimizing bias in discovery and interpretation of results. Our findings open up possibilities for further methods development for intron reads and RNA-seq data in general.

Germline heterozygous mutations in Nxf1 perturb RNA metabolism and trigger thrombocytopenia and lymphopenia in mice.

In eukaryotic cells, messenger RNA (mRNA) molecules are exported from the nucleus to the cytoplasm, where they are translated. The highly conserved protein nuclear RNA export factor1 (Nxf1) is an important mediator of this process. Although studies in yeast and in human cell lines have shed light on the biochemical mechanisms of Nxf1 function, its contribution to mammalian physiology is less clear. Several groups have identified recurrent NXF1 mutations in chronic lymphocytic leukemia (CLL), placing it alongside several RNA-metabolism factors (including SF3B1, XPO, RPS15) whose dysregulation is thought to contribute to CLL pathogenesis. We report here an allelic series of germline point mutations in murine Nxf1. Mice heterozygous for these loss-of-function Nxf1 mutations exhibit thrombocytopenia and lymphopenia, together with milder hematological defects. This is primarily caused by cell-intrinsic defects in the survival of platelets and peripheral lymphocytes, which are sensitized to intrinsic apoptosis. In contrast, Nxf1 mutations have almost no effect on red blood cell homeostasis. Comparative transcriptome analysis of platelets, lymphocytes, and erythrocytes from Nxf1-mutant mice shows that, in response to impaired Nxf1 function, the cytoplasmic representation of transcripts encoding regulators of RNA metabolism is altered in a unique, lineage-specific way. Thus, blood cell lineages exhibit differential requirements for Nxf1-mediated global mRNA export.

Easy and efficient ensemble gene set testing with EGSEA.

Gene set enrichment analysis is a popular approach for prioritising the biological processes perturbed in genomic datasets. The Bioconductor project hosts over 80 software packages capable of gene set analysis. Most of these packages search for enriched signatures amongst differentially regulated genes to reveal higher level biological themes that may be missed when focusing only on evidence from individual genes. With so many different methods on offer, choosing the best algorithm and visualization approach can be challenging. The EGSEA package solves this problem by combining results from up to 12 prominent gene set testing algorithms to obtain a consensus ranking of biologically relevant results.This workflow demonstrates how EGSEA can extend limma-based differential expression analyses for RNA-seq and microarray data using experiments that profile 3 distinct cell populations important for studying the origins of breast cancer. Following data normalization and set-up of an appropriate linear model for differential expression analysis, EGSEA builds gene signature specific indexes that link a wide range of mouse or human gene set collections obtained from MSigDB, GeneSetDB and KEGG to the gene expression data being investigated. EGSEA is then configured and the ensemble enrichment analysis run, returning an object that can be queried using several S4 methods for ranking gene sets and visualizing results via heatmaps, KEGG pathway views, GO graphs, scatter plots and bar plots. Finally, an HTML report that combines these displays can fast-track the sharing of results with collaborators, and thus expedite downstream biological validation. EGSEA is simple to use and can be easily integrated with existing gene expression analysis pipelines for both human and mouse data.

Identification of quiescent and spatially restricted mammary stem cells that are hormone responsive.

Despite accumulating evidence for a mammary differentiation hierarchy, the basal compartment comprising stem cells remains poorly characterized. Through gene expression profiling of Lgr5+ basal epithelial cells, we identify a new marker, Tetraspanin8 (Tspan8). Fractionation based on Tspan8 and Lgr5 expression uncovered three distinct mammary stem cell (MaSC) subsets in the adult mammary gland. These exist in a largely quiescent state but differ in their reconstituting ability, spatial localization, and their molecular and epigenetic signatures. Interestingly, the deeply quiescent MaSC subset (Lgr5+Tspan8hi) resides within the proximal region throughout life, and has a transcriptome strikingly similar to that of claudin-low tumours. Lgr5+Tspan8hi cells appear to originate from the embryonic mammary primordia before switching to a quiescent state postnatally but can be activated by ovarian hormones. Our findings reveal an unexpected degree of complexity within the adult MaSC compartment and identify a dormant subset poised for activation in response to physiological stimuli.