Discovering Novel Genes in Non-Model Fly Accessory Glands Using De Novo Nanopore Transcriptomics

BackgroundOxford Nanopore Technologies (ONT) long-read transcriptomes offer many advantages including long reads (>10kbp), end-to-end transcripts, structural variants, isoform-level resolution of genes and expression. However, uptake of ONT transcriptomics is still low, largely due to high error rates (2 to 13%) and reliance on reference databases that are unavailable for many non-model species. Additionally, bioinformatics tools and pipelines for de novo ONT transcriptomics are still in early stages of development. ResultsHere, we use de novo ONT GridION transcriptomics to discover novel genes from the male accessory glands (AG) of a widespread, non-model dung fly, Sepsis punctum. Insect AGs are of particular interest for this as they are hotspots for rapid evolution of novel reproductive genes, and they synthesize seminal fluid proteins that lack homology to any other known proteins. We implement a completely de novo ONT GridION transcriptome pipeline, incorporating quality-filtering and rigorous error-correction procedures, to characterize this novel gene set and to quantify their expression. Specifically, we compare these ONT genes and their expression against de novo lllumina HiSeq transcriptome data. We find 40 high-quality and high-confidence ONT genes that cross-verify against Illumina genes; twenty-six of which are novel and specific to S. punctum. Read count based expression quantification in ONT samples is highly congruent with Illumina’s Transcript per Million (TPM), both in overall pattern and within functional categories. Novel genes account for an average of 81% of total gene expression underscoring their functional importance in S. punctum AGs. Eighty percentage of these genes are secretory in nature, responsible for 74% total gene expression. Notably, median sequence similarities of ONT nucleotide and protein sequences match within-Illumina sequence similarities indicating that our de novo ONT transcriptome pipeline successfully mitigated sequencing errors. ConclusionsThis is the first study to adapt ONT transcriptomics for completely de novo characterization of novel genes in animals. Our study demonstrates that ONT long-reads, constituting a quarter of the number of bases sequenced at less than a third the cost of Illumina reads, can be a resource-friendly and cost-effective solution for end-to-end sequencing of unknown genes even in the absence of a reference database.

Time-Course Transcriptome Analysis Reveals Global Gene Expression Profiling and Dynamic Developmental Signatures across Complete Life Cycle of Bombyx mori

Background: The silkworm (Bombyx mori) is an important lepidopteran model insect worldwide which undergoes a complete metamorphosis developmental process. Although genome sequencing has been long performed, no transcriptome data covering the complete life cycle are available. Methods: Herein, a total of 10 samples were collected consecutively at four different developmental stages, including eggs of 24 h after oviposition (Ed) and eggs of 24 h after artificial egg-hatching (E); larvae from fist to fifth instar (L1–L5); early and late pupa (P4 and P8); and adult moth (M), were subjected to Illumina RNA-Seq and time-course analysis. Results: The summations of the gene expression of the silkworm ten developmental stages show: at Ed stage, eggs develop towards diapause status, the total gene expression level is relatively low; at E stage, after artificial egg-hatching, the expression level improves rapidly; during larval stages from L1–L5, the expression level rises gradually and reaches a peak at L5 stage; during pupae and moth stages, the total gene expression decline stage by stage. The results revealed a dynamical gene expression profile exhibiting significant differential expressions throughout the silkworm life cycle. Moreover, stage-specific key genes were identified at different developmental stages, suggesting their functions mainly characterized in maintaining insect development and immunity homeostasis or driving metamorphosis. GO annotation and KEGG enrichment analysis further revealed the most significantly enriched and fundamentally biological processes during silkworm growth. Conclusion: Collectively, our omics data depicted the first comprehensive landscape of dynamic transcriptome throughout complete developmental processes of B. mori. Our findings also provide valuable references and novel insights into understanding the molecular developmental remodeling events for other Lepidoptera species.

Robust and annotation-free analysis of isoform variation using short-read scRNA-seq data

Although isoform diversity is acknowledged as a fundamental and pervasive aspect of gene expression in higher eukaryotes, it is often omitted from single-cell studies due to quantification challenges inherent to commonly used short-read sequencing technologies. To address this issue, we have developed a suite of computational tools to investigate isoform variation by focusing on splice junction usage patterns, which can often be well characterized in spite of technical difficulties. Our method, which we name scQuint (single-cell quantification of introns), can perform accurate quantification, dimensionality reduction, and differential splicing analysis using short-read, full-length single-cell RNA-seq data. Notably, scQuint does not require transcriptome annotations and is robust to technical artifacts. In applications across diverse mouse tissues from Tabula Muris and the primary motor cortex from the BRAIN Initiative Cell Census Network, we find evidence of strong cell-type-specific isoform variation, complementary to total gene expression, and also identify a large volume of previously unannotated splice junctions. As a community resource, we provide ways to interactively visualize and explore these results, accessible at https://github.com/songlab-cal/scquint-analysis/.

Cleft lip and cleft palate (CL/P) in Esrp1 KO mice is associated with alterations in Wnt signaling and epithelial-mesenchymal crosstalk

ABSTRACTCleft lip is one of the most highly prevalent birth defects in human patients. However, there remain a limited number of mouse models of cleft lip and thus much work is needed to further characterize genes and mechanisms that lead to this disorder. It is well established that crosstalk between epithelial and mesenchymal cells underlies formation of the face and palate, yet the basic molecular events mediating this crosstalk are still poorly understood. We previously demonstrated that mice with ablation of the epithelial-specific splicing factor Esrp1 have fully penetrant bilateral CL/P. In this study we further investigated the mechanisms by which ablation of Esrp1 leads to cleft lip as well as cleft palate. These studies included a detailed analysis of the changes in splicing and total gene expression in embryonic ectoderm during formation of the face as well as gene expression changes in adjacent mesenchyme. We identified altered expression in components of pathways previously implicated in cleft lip and/or palate, including numerous components of the Wnt signaling pathway. These findings illustrate that maintenance of an Esrp1 regulated epithelial splicing program is essential for face development through regulation of key signaling pathways.

Exploring the Diversity and Novelty of Toxin Genes in Naja sumatrana, the Equatorial Spitting Cobra from Malaysia through De Novo Venom-Gland Transcriptomics

The equatorial spitting cobra, Naja sumatrana, is a distinct species of medically important venomous snakes, listed as WHO Category 1 in Southeast Asia. The diversity of its venom genes has not been comprehensively examined, although a few toxin sequences annotated to Naja sputatrix were reported previously through cloning studies. To investigate this species venom genes’ diversity, de novo venom-gland transcriptomics of N. sumatrana from West Malaysia was conducted using next-generation sequencing technology. Genes encoding toxins represented only 60 of the 55,396 transcripts, but were highly expressed, contributing to 79.22% of total gene expression (by total FPKM) in the venom-glands. The toxin transcripts belong to 21 families, and 29 transcripts were further identified as full-length. Three-finger toxins (3FTx) composed of long, short, and non-conventional groups, constituted the majority of toxin transcripts (91.11% of total toxin FPKM), followed by phospholipase A2 (PLA2, 7.42%)—which are putatively pro-inflammatory and cytotoxic. The remaining transcripts in the 19 families were expressed at extremely low levels. Presumably, these toxins were associated with ancillary functions. Our findings unveil the diverse toxin genes unique to N. sumatrana, and provide insights into the pathophysiology of N. sumatrana envenoming.

Hierarchical Analysis of Multi-mapping RNA-Seq Reads Improves the Accuracy of Allele-specific Expression

Allele-specific expression (ASE) refers to the differential abundance of the allelic copies of a transcript. Direct RNA sequencing (RNA-Seq) can provide quantitative estimates of ASE for genes with transcribed polymorphisms. However, estimating ASE is challenging due to ambiguities in read alignment. Current approaches do not account for the hierarchy of multiple read alignments to genes, isoforms, and alleles. We have developed EMASE (Expectation-Maximization for Allele Specific Expression), an integrated approach to estimate total gene expression, ASE, and isoform usage based on hierarchical allocation of multi-mapping reads. In simulations, EMASE outperforms standard ASE estimation methods. We apply EMASE to RNA-Seq data from F1 hybrid mice where we observe widespread ASE associated with cis-acting polymorphisms and a small number of parent-of-origin effects at known imprinted genes. The EMASE software is freely available under GNU license at https://github.com/churchill-lab/emase and it can be adapted to other sequencing applications.

Inferring intrinsic and extrinsic noise from a dual fluorescent reporter

Dual fluorescent reporter constructs, which measure gene expression from two identical promoters within the same cell, allow total gene expression noise to be decomposed into an extrinsic component, roughly associated with cell-to-cell fluctuations in cellular component concentrations, and intrinsic noise, roughly associated with inherent stochasticity of the biochemical reactions involved in gene expression [1]. A recent paper by Fu and Pachter presented frequentist statistical estimators for intrinsic and extrinsic noise using data from dual reporters [2]. For comparison, I here present results of a Bayesian analysis of this problem. I show that the orthodox estimators suffer from pathologies such as predicting negative values for a manifestly non-negative quantity, i.e. variance, and show that the Bayesian estimators do not suffer from such pathologies. In addition, I show that the Bayesian analysis automatically identifies that optimal estimates of intrinsic and extrinsic noise depend on a subtle combination of two statistics of the data, allowing for accuracies that are up to twice the accuracy of the orthodox estimators in some parameter regimes.I hope up this little worked out example contrasting orthodox statistical analysis based on ad hoc estimators with estimators resulting from a Bayesian analysis, will be educational for others in the field. I distribute a Mathematica Notebook with this paper that allows users to easily reproduce all results and figures of the paper.

Evaluating intra- and inter-individual variation in the human placental transcriptome

Background: Gene expression variation is a phenotypic trait of particular interest as it represents the initial link between genotype and other phenotypes. Analyzing how such variation apportions among and within groups allows for the evaluation of how genetic and environmental factors influence such traits. It also provides opportunities to identify genes and pathways that may have been influenced by non-neutral processes. Here we use a population genetics framework and next generation sequencing to evaluate how gene expression variation is apportioned among four human groups in a natural biological tissue, the placenta. Results: We estimate that on average, 33.2%, 58.9% and 7.8% of the placental transcriptome is explained by variation within individuals, among individuals and among human groups, respectively. Additionally, when technical and biological traits are included in models of gene expression they account for roughly 2% of total gene expression variation. Notably, the variation that is significantly different among groups is enriched in biological pathways associated with immune response, cell signaling and metabolism. Many biological traits demonstrated correlated changes in expression in numerous pathways of potential interest to clinicians and evolutionary biologists. Finally, we estimate that the majority of the human placental transcriptome (65% of expressed genes) exhibits expression profiles consistent with neutrality; the remainder are consistent with stabilizing selection (26%), directional selection (4.9%), or diversifying selection (4.8%). Conclusion: We apportion placental gene expression variation into individual, population and biological trait factors and identify how each influence the transcriptome. Additionally, we advance methods to associate expression profiles with different forms of selection.

Alternative Splicing Is a Frequent Event and Impacts Clinical Outcome in Myeloma: A Large RNA-Seq Data Analysis of Newly-Diagnosed Myeloma Patients

Alternative splicing is an important post-translational change that alters gene function. Misregulation of alternative splicing has been implicated in number of disease processes including cancer. Here we have analyzed alternative splicing in myeloma using high throughput RNA-seq. Our analytic pipeline for RNA-seq data used in this investigation not only provides information on expression levels for genes, but also provides information on the expression of known splice variants of genes (isoforms), and can identify novel exon level events across individuals (i.e. exon skipping events). We conducted a study of 328 newly-diagnosed patients with multiple myeloma treated homogeneously with novel agent combination containting lenalidomide, bortezomib and dexamethsone with or without high-dose melphalan followed by lenalidomide maintenance in the IFM/DFCI study. RNA isolated from purified CD138+ MM cells collected at the time of diagnosis and from 18 normal donor plasma cells were processed by RNA-seq (100 million paired end reads on Illumina HiSeq) and analyzed using a custom computational and statistical pipeline. Following read alignment to hg19, we utilized RSEM to quantify both gene-level and isoform-level expression of known ENSEMBL transcripts. We then implemented a novel testing approach based on compositional regression to discover genes that show significant isoform switching between the 328 MM samples and 18 Normal Plasma Cell (NPC) samples from healthy donors. Using various programs and their modifications, we also identified novel alternative splicing events, such as exon skipping and mutually exclusive exon usage, among others. Patient data for MM characteristics, cytogenetic and FISH as well as clinical survival outcomes were also analyzed and correlated with genomic data. We observed over 600 genes showing significant changes in relative isoform abundances (isoform switching) between MM and normal samples. A number of previously characterized genes including MYCL1 (adj. p = 0.0014) and CCND3 (adj. p = 0.0013), and MAP kinase-related genes (MAP3K8, MAPKAPK2, MAPKAPK3, MAP4K4) exhibited significant isoform switching compared to normal, in addition to some not well characterized genes. Genes showing the greatest magnitude of isoform switching include MEFV (adj. p = 2.7 x 10-5), showing a two fold change in the relative major isoform abundance compared to normal, and has been previously shown to have a role in lymphoid neoplasms. We applied hierarchical clustering to the isoforms showing significant changes in isoform-switching and identified 4 distinct clusters, which are currently being investigated for correlation with clinical subtypes of MM. Exon level analyses of alternative splicing events, such as exon skipping, are currently underway. Clinical data including MM characteristics, cytogenetics, FISH and survival outcomes was available for a subset of 265 patients. We found that 109 genes showed significant isoform switching between t(4;14) and non-t(4;14) patients, such as CD44 (adj. p =1.8 x 10-6) and WHSC1 (adj. p =5.1 x 10-28). Comparing del17p (28 in total) and non del17p patients, we found no significant splicing changes after multiple testing adjustment. Of these genes, only a subset (40%) were shown to be differentially expressed in terms of total gene expression, suggesting the importance of examining alternative splicing events in addition to total gene expression. With respect to treatment response, we compared the expression of gene isoforms between patients achieving complete response (CR) versus others and identified 38 isoforms associated with response to treatment (adj. p value < 0.05), with SEPT9, SLC2A5, and UBX6 having the strongest associations (adj. p-value < 3 x 10-4). Using a univariate cox regression model, 4 spliced isoforms relating to 3 genes were identified as having significant correlation with event-free survival (EFS) (FDR-adjusted cox p value < 0.05). We are in the process of now integrating the gene expression data with altered splicing data to develop an integrated survival model. In summary, this study highlights the significant frequency, biological and clinical importance of alternative splicing in MM and points to the need for evaluation of not only the expression level of genes but also post-translational modifications. The genes identified here are important targets for therapy as well as possible immune modulation. Disclosures Moreau: Celgene Corporation: Honoraria, Membership on an entity's Board of Directors or advisory committees.