Comparison of normalization and differential expression analyses using RNA-Seq data from 726 individual Drosophila melanogaster

A comprehensive workflow for optimizing RNA-seq data analysis

BACKGROUND

Current RNA-seq analysis software for RNA-seq data tends to use similar parameters across different species without considering species-specific differences. However, the suitability and accuracy of these tools may vary when analyzing data from different species, such as humans, animals, plants, fungi, and bacteria. For most laboratory researchers lacking a background in information science, determining how to construct an analysis workflow that meets their specific needs from the array of complex analytical tools available poses a significant challenge.

RESULTS

By utilizing RNA-seq data from plants, animals, and fungi, it was observed that different analytical tools demonstrate some variations in performance when applied to different species. A comprehensive experiment was conducted specifically for analyzing plant pathogenic fungal data, focusing on differential gene analysis as the ultimate goal. In this study, 288 pipelines using different tools were applied to analyze five fungal RNA-seq datasets, and the performance of their results was evaluated based on simulation. This led to the establishment of a relatively universal and superior fungal RNA-seq analysis pipeline that can serve as a reference, and certain standards for selecting analysis tools were derived for reference. Additionally, we compared various tools for alternative splicing analysis. The results based on simulated data indicated that rMATS remained the optimal choice, although consideration could be given to supplementing with tools such as SpliceWiz.

CONCLUSION

The experimental results demonstrate that, in comparison to the default software parameter configurations, the analysis combination results after tuning can provide more accurate biological insights. It is beneficial to carefully select suitable analysis software based on the data, rather than indiscriminately choosing tools, in order to achieve high-quality analysis results more efficiently.

The role of DNA methylation reprogramming during sex determination and sex reversal in the Pacific oyster Crassostrea gigas

DNA methylation is instrumental in vertebrate sex reversal. However, the mechanism of DNA methylation regulation regarding sex reversal in invertebrates is unclear. In this study, we used whole genome bisulfite sequencing (WGBS) to map single-base resolution methylation profiles of the Pacific oyster, including female-to-male (FMa-to-FMb) and male-to-female (MFa-to-MFb) sex reversal, as well as sex non-reversed males and females (MMa-to-MMb and FFa-to-FFb). The results showed that global DNA methylation levels increase during female-to-male sex reversals, with a particular increase in the proportion of high methylation levels (mCGs >0.75) and a decrease in the proportion of intermediate methylation levels (0.25 < mCGs <0.75). This increase in DNA methylation was mainly associated with the elevated expression of DNA methylase genes. Genome-wide methylation patterns of females were accurately remodeled to those of males after sex reversal, while the opposite was true for the male-to-female reversal. Those findings directly indicate that alterations in DNA methylation play a significant role in sex reversal in Pacific oysters. Comparative analysis of the DNA methylomes of pre- and post- sex reversal gonadal tissues (FMb-vs-FMa or MFb-vs-MFa) revealed that differentially methylated genes were mainly involved in the biological processes of sex determination or gonadal development. However critical genes such as Dmrt1, Foxl2 and Sox-like, which are involved in the putative sex determination pathway in Pacific oysters, showed almost an absence of methylation modifications, varying greatly from vertebrates. Additionally, comparative analysis of the DNA methylomes of sexual reversal and sex non-reversal (FMa-vs-FFa or MFa-vs-MMa) revealed that heat shock protein genes, such as Hsp68-like and Hsp70B, were important for the occurrence of sex reversal. These findings shed light on the epigenetic mechanisms underlying the maintenance of gonadal plasticity and the reversal of organ architecture in oysters.

Predicting miRNA-disease association via graph attention learning and multiplex adaptive modality fusion

miRNAs are a class of small non-coding RNA molecules that play important roles in gene regulation. They are crucial for maintaining normal cellular functions, and dysregulation or dysfunction of miRNAs which are linked to the onset and advancement of multiple human diseases. Research on miRNAs has unveiled novel avenues in the realm of the diagnosis, treatment, and prevention of human diseases. However, clinical trials pose challenges and drawbacks, such as complexity and time-consuming processes, which create obstacles for many researchers. Graph Attention Network (GAT) has shown excellent performance in handling graph-structured data for tasks such as link prediction. Some studies have successfully applied GAT to miRNA-disease association prediction. However, there are several drawbacks to existing methods. Firstly, most of the previous models rely solely on concatenation operations to merge features of miRNAs and diseases, which results in the deprivation of significant modality-specific information and even the inclusion of redundant information. Secondly, as the number of layers in GAT increases, there is a possibility of excessive smoothing in the feature extraction process, which significantly affects the prediction accuracy. To address these issues and effectively complete miRNA disease prediction tasks, we propose an innovative model called Multiplex Adaptive Modality Fusion Graph Attention Network (MAMFGAT). MAMFGAT utilizes GAT as the main structure for feature aggregation and incorporates a multi-modal adaptive fusion module to extract features from three interconnected networks: the miRNA-disease association network, the miRNA similarity network, and the disease similarity network. It employs adaptive learning and cross-modality contrastive learning to fuse more effective miRNA and disease feature embeddings as well as incorporates multi-modal residual feature fusion to tackle the problem of excessive feature smoothing in GATs. Finally, we employ a Multi-Layer Perceptron (MLP) model that takes the embeddings of miRNA and disease features as input to anticipate the presence of potential miRNA-disease associations. Extensive experimental results provide evidence of the superior performance of MAMFGAT in comparison to other state-of-the-art methods. To validate the significance of various modalities and assess the efficacy of the designed modules, we performed an ablation analysis. Furthermore, MAMFGAT shows outstanding performance in three cancer case studies, indicating that it is a reliable method for studying the association between miRNA and diseases. The implementation of MAMFGAT can be accessed at the following GitHub repository: https://github.com/zixiaojin66/MAMFGAT-master.

Plant Transcriptome Analysis with HISAT-StringTie-Ballgown and TopHat-Cufflinks Pipelines

The development of next-generation sequencing technology has led to a burst of data in a single assay. Management of a large dataset requires high demands on bioinformatic skills and computing resources. Here we present two popular pipelines for RNA-seq data analysis, using open-source software tools HISAT-StringTie-Ballgown and TopHat-Cufflinks. To meet the need of plant scientist, we describe in detail how to perform such comprehensive analysis beginning with raw RNA-seq reads and available reference genome. It allows biologists to align short reads to a reference genome, measure the transcript abundance, and analyze gene differential expression under two or more conditions. We also discuss other RNA-seq tools that are comparable or alternative to this protocol.

Statistical normalization methods in microbiome data with application to microbiome cancer research

Mounting evidence has shown that gut microbiome is associated with various cancers, including gastrointestinal (GI) tract and non-GI tract cancers. But microbiome data have unique characteristics and pose major challenges when using standard statistical methods causing results to be invalid or misleading. Thus, to analyze microbiome data, it not only needs appropriate statistical methods, but also requires microbiome data to be normalized prior to statistical analysis. Here, we first describe the unique characteristics of microbiome data and the challenges in analyzing them (Section 2). Then, we provide an overall review on the available normalization methods of 16S rRNA and shotgun metagenomic data along with examples of their applications in microbiome cancer research (Section 3). In Section 4, we comprehensively investigate how the normalization methods of 16S rRNA and shotgun metagenomic data are evaluated. Finally, we summarize and conclude with remarks on statistical normalization methods (Section 5). Altogether, this review aims to provide a broad and comprehensive view and remarks on the promises and challenges of the statistical normalization methods in microbiome data with microbiome cancer research examples.

Transcriptome software results show significant variation among different commercial pipelines

BACKGROUND

We have been documenting the biological responses to low levels of radiation (natural background) and very low level radiation (below background), and thus these studies are testing mild external stimuli to which we would expect relatively mild biological responses. We recently published a transcriptome software comparison study based on RNA-Seqs from a below background radiation treatment of two model organisms, E. coli and C. elegans (Thawng and Smith, BMC Genomics 23:452, 2022). We reported DNAstar-D (Deseq2 in the DNAstar software pipeline) to be the more conservative, realistic tool for differential gene expression compared to other transcriptome software packages (CLC, Partek and DNAstar-E (using edgeR). Here we report two follow-up studies (one with a new model organism, Aedes aegypti and another software package (Azenta) on transcriptome responses from varying dose rates using three different sources of natural radiation.

RESULTS

When E. coli was exposed to varying levels of K40, we again found that the DNAstar-D pipeline yielded a more conservative number of DEGs and a lower fold-difference than the CLC pipeline and DNAstar-E run in parallel. After a 30 read minimum cutoff criterion was applied to the data, the number of significant DEGs ranged from 0 to 81 with DNAstar-D, while the number of significant DEGs ranged from 4 to 117 and 14 to 139 using DNAstar-E and the CLC pipelines, respectively. In terms of the extent of expression, the highest foldchange DEG was observed in DNAstar-E with 19.7-fold followed by 12.5-fold in CLC and 4.3-fold in DNAstar-D. In a recently completed study with Ae. Aegypti and using another software package (Azenta), we analyzed the RNA-Seq response to similar sources of low-level radiation and again found the DNAstar-D pipeline to give the more conservative number and fold-expression of DEGs compared to other softwares. The number of significant DEGs ranged 31-221 in Azenta and 31 to 237 in CLC, 19-252 in DNAstar-E and 0-67 in DNAStar-D. The highest fold-change of DEGs were found in CLC (1,350.9-fold), with DNAstar-E (5.9 -fold) and Azenta (5.5-fold) intermediate, and the lowest levels of expression (4-fold) found in DNAstar-D.

CONCLUSIONS

This study once again highlights the importance of choosing appropriate software for transcriptome analysis. Using three different biological models (bacteria, nematode and mosquito) in four different studies testing very low levels of radiation (Van Voorhies et al., Front Public Health 8:581796, 2020; Thawng and Smith, BMC Genomics 23:452, 2022; current study), the CLC software package resulted in what appears to be an exaggerated gene expression response in terms of numbers of DEGs and extent of expression. Setting a 30-read cutoff diminishes this exaggerated response in most of the software tested. We have further affirmed that DNAstar-Deseq2 gives a more conservative transcriptome expression pattern which appears more suitable for studies expecting subtle gene expression patterns.

Ploidy evolution in a wild yeast is linked to an interaction between cell type and metabolism

Ploidy is an evolutionarily labile trait, and its variation across the tree of life has profound impacts on evolutionary trajectories and life histories. The immediate consequences and molecular causes of ploidy variation on organismal fitness are frequently less clear, although extreme mating type skews in some fungi hint at links between cell type and adaptive traits. Here, we report an unusual recurrent ploidy reduction in replicate populations of the budding yeast Saccharomyces eubayanus experimentally evolved for improvement of a key metabolic trait, the ability to use maltose as a carbon source. We find that haploids have a substantial, but conditional, fitness advantage in the absence of other genetic variation. Using engineered genotypes that decouple the effects of ploidy and cell type, we show that increased fitness is primarily due to the distinct transcriptional program deployed by haploid-like cell types, with a significant but smaller contribution from absolute ploidy. The link between cell-type specification and the carbon metabolism adaptation can be traced to the noncanonical regulation of a maltose transporter by a haploid-specific gene. This study provides novel mechanistic insight into the molecular basis of an environment-cell type fitness interaction and illustrates how selection on traits unexpectedly linked to ploidy states or cell types can drive karyotypic evolution in fungi.

Long noncoding RNA study: Genome-wide approaches

Long noncoding RNAs (lncRNAs) have been confirmed to play a crucial role in various biological processes across several species. Though many efforts have been devoted to the expansion of the lncRNAs landscape, much about lncRNAs is still unknown due to their great complexity. The development of high-throughput technologies and the constantly improved bioinformatic methods have resulted in a rapid expansion of lncRNA research and relevant databases. In this review, we introduced genome-wide research of lncRNAs in three parts: (i) novel lncRNA identification by high-throughput sequencing and computational pipelines; (ii) functional characterization of lncRNAs by expression atlas profiling, genome-scale screening, and the research of cancer-related lncRNAs; (iii) mechanism research by large-scale experimental technologies and computational analysis. Besides, primary experimental methods and bioinformatic pipelines related to these three parts are summarized. This review aimed to provide a comprehensive and systemic overview of lncRNA genome-wide research strategies and indicate a genome-wide lncRNA research system.

Change in energy-consuming strategy, nucleolar metabolism and physical defense in Macrobrachium rosenbergii after acute and chronic polystyrene nanoparticles exposure

The COVID-19 pandemic has further intensified plastic pollution due to the escalated use of single-use gloves and masks, consequently leading to the widespread presence of microplastics (MPs) and nanoplastics (NPs) in major rivers and lakes worldwide. Macrobrachium rosenbergii has become an important experimental subject due to its ecological role and environmental sensitivity. In this study, we sought to comprehend the ramifications of NPs on the widely-distributed freshwater prawn, M rosenbergii, by conducting a detailed analysis of its responses to NPs after both 96 h and 30 days of exposure. The transcriptome analysis revealed 918 differentially expressed unigenes (DEGs) after 30 days of NPs exposure (356 upregulated, 562 downregulated) and 2376 DEGs after 96 h of NPs exposure (1541 upregulated, 835 downregulated). The results of DEGs expression indicated that acute NPs exposure enhanced carbohydrate transport and metabolism, fostering chitin and extracellular matrix processes. In contrast, chronic NPs exposure induced nucleolar stress in M. rosenbergii, impeding ribosome development and mRNA maturation while showing no significant changes in glucose metabolism. Our findings underscore the M. rosenbergii distinct coping mechanisms during acute and chronic NPs exposure, elucidating its vital adaptive strategies. These results contribute to our understanding of the ecological implications of NPs pollution and its impact on aquatic animals.

Variability in donor leukocyte counts confound the use of common RNA sequencing data normalization strategies in transcriptomic biomarker studies performed with whole blood

Gene expression data generated from whole blood via next generation sequencing is frequently used in studies aimed at identifying mRNA-based biomarker panels with utility for diagnosis or monitoring of human disease. These investigations often employ data normalization techniques more typically used for analysis of data originating from solid tissues, which largely operate under the general assumption that specimens have similar transcriptome composition. However, this assumption may be violated when working with data generated from whole blood, which is more cellularly dynamic, leading to potential confounds. In this study, we used next generation sequencing in combination with flow cytometry to assess the influence of donor leukocyte counts on the transcriptional composition of whole blood specimens sampled from a cohort of 138 human subjects, and then subsequently examined the effect of four frequently used data normalization approaches on our ability to detect inter-specimen biological variance, using the flow cytometry data to benchmark each specimens true cellular and molecular identity. Whole blood samples originating from donors with differing leukocyte counts exhibited dramatic differences in both genome-wide distributions of transcript abundance and gene-level expression patterns. Consequently, three of the normalization strategies we tested, including median ratio (MRN), trimmed mean of m-values (TMM), and quantile normalization, noticeably masked the true biological structure of the data and impaired our ability to detect true interspecimen differences in mRNA levels. The only strategy that improved our ability to detect true biological variance was simple scaling of read counts by sequencing depth, which unlike the aforementioned approaches, makes no assumptions regarding transcriptome composition.

Viruses in Laboratory Drosophila and Their Impact on Host Gene Expression

Drosophila melanogaster has one of the best characterized antiviral immune responses among invertebrates. However, relatively few easily transmitted natural virus isolates are available, and so many Drosophila experiments have been performed using artificial infection routes and artificial host-virus combinations. These may not reflect natural infections, especially for subtle phenotypes such as gene expression. Here, to explore the laboratory virus community and to better understand how natural virus infections induce changes in gene expression, we have analysed seven publicly available D. melanogaster transcriptomic sequencing datasets that were originally sequenced for projects unrelated to virus infection. We have found ten known viruses-including five that have not been experimentally isolated-but no previously unknown viruses. Our analysis of host gene expression revealed that numerous genes were differentially expressed in flies that were naturally infected with a virus. For example, flies infected with nora virus showed patterns of gene expression consistent with intestinal vacuolization and possible host repair via the upd3 JAK/STAT pathway. We also found marked sex differences in virus-induced differential gene expression. Our results show that natural virus infection in laboratory Drosophila does indeed induce detectable changes in gene expression, suggesting that this may form an important background condition for experimental studies in the laboratory.

De novo transcriptome analysis of Dysoxylum binectariferum to unravel the biosynthesis of pharmaceutically relevant specialized metabolites

The tropical tree, D. binectariferum, is a prominent source of chromone alkaloid rohitukine, which is used in the semi-syntheses of anticancer molecules such as flavopiridol and P-276-00. The biosynthetic pathway of rohitukine or its derivatives is currently unknown in plants. Here, we explored chromone alkaloid biosynthesis in D. binectariferum through targeted transcriptome sequencing. Illumina sequencing of leaves and roots of a year-old D. binectariferum seedling generated, 42.43 and 38.74 million paired-end short reads, respectively. Quality filtering and de novo assembly of the transcriptome generated 274,970 contigs and 126,788 unigenes with an N50 contig length of 1560 bp. The assembly generated 117,619 translated unigene protein sequences and 51,598 non-redundant sequences. Nearly 80% of these non-redundant sequences were annotated to publicly available protein and nucleotide databases, suggesting the completeness and effectiveness of the transcriptome assembly. Using the assembly, we identified a chalcone synthase (CHS) and three type III polyketide synthases (PKS-III; non-CHS type) that are likely to be involved in the biosynthesis of chromone ring/noreugenin moiety of rohitukine. We also identified key enzymes like lysine decarboxylase in the piperidine pathway that make the piperidine moiety of rohitukine. Besides these, the upstream enzymes in flavonoid biosynthesis like phenylalanine ammonia-lyase (PAL), trans-cinnamate 4-hydroxylase (C4H),4-coumarate-CoA ligase (4CL), and chalcone isomerase (CHI) have also been identified. Also, terpene synthases that are likely to be involved in the biosynthesis of various terpenoid scaffolds have been identified. Together, the D. binectariferum transcriptome resource forms a basis for further exploration of biosynthetic pathways of these valuable compounds through functional validation of the candidate genes and metabolic engineering in heterologous hosts. Additionally, the transcriptome dataset generated will serve as an important resource for research on functional genomics and enzyme discovery in D. binectariferum and comparative analysis with other Meliaceae family members.

Transcriptomic Insight into the Pollen Tube Growth of Olea europaea L. subsp. europaea Reveals Reprogramming and Pollen-Specific Genes Including New Transcription Factors

The pollen tube is a key innovation of land plants that is essential for successful fertilisation. Its development and growth have been profusely studied in model organisms, but in spite of the economic impact of olive trees, little is known regarding the genome-wide events underlying pollen hydration and growth in this species. To fill this gap, triplicate mRNA samples at 0, 1, 3, and 6 h of in vitro germination of olive cultivar Picual pollen were analysed by RNA-seq. A bioinformatics R workflow called RSeqFlow was developed contemplating the best practices described in the literature, covering from expression data filtering to differential expression and clustering, to finally propose hub genes. The resulting olive pollen transcriptome consisted of 22,418 reliable transcripts, where 5364 were differentially expressed, out of which 173 have no orthologue in plants and up to 3 of them might be pollen-specific transcription factors. Functional enrichment revealed a deep transcriptional reprogramming in mature olive pollen that is also dependent on protein stability and turnover to allow pollen tube emergence, with many hub genes related to heat shock proteins and F-box-containing proteins. Reprogramming extends to the first 3 h of growth, including processes consistent with studies performed in other plant species, such as global down-regulation of biosynthetic processes, vesicle/organelle trafficking and cytoskeleton remodelling. In the last stages, growth should be maintained from persistent transcripts. Mature pollen is equipped with transcripts to successfully cope with adverse environments, even though the in vitro growth seems to induce several stress responses. Finally, pollen-specific transcription factors were proposed as probable drivers of pollen germination in olive trees, which also shows an overall increased number of pollen-specific gene isoforms relative to other plants.

Integrating massive RNA-seq data to elucidate transcriptome dynamics in Drosophila melanogaster

The volume of ribonucleic acid (RNA)-seq data has increased exponentially, providing numerous new insights into various biological processes. However, due to significant practical challenges, such as data heterogeneity, it is still difficult to ensure the quality of these data when integrated. Although some quality control methods have been developed, sample consistency is rarely considered and these methods are susceptible to artificial factors. Here, we developed MassiveQC, an unsupervised machine learning-based approach, to automatically download and filter large-scale high-throughput data. In addition to the read quality used in other tools, MassiveQC also uses the alignment and expression quality as model features. Meanwhile, it is user-friendly since the cutoff is generated from self-reporting and is applicable to multimodal data. To explore its value, we applied MassiveQC to Drosophila RNA-seq data and generated a comprehensive transcriptome atlas across 28 tissues from embryogenesis to adulthood. We systematically characterized fly gene expression dynamics and found that genes with high expression dynamics were likely to be evolutionarily young and expressed at late developmental stages, exhibiting high nonsynonymous substitution rates and low phenotypic severity, and they were involved in simple regulatory programs. We also discovered that human and Drosophila had strong positive correlations in gene expression in orthologous organs, revealing the great potential of the Drosophila system for studying human development and disease.

Pathogen detection in RNA-seq data with Pathonoia

BACKGROUND

Bacterial and viral infections may cause or exacerbate various human diseases and to detect microbes in tissue, one method of choice is RNA sequencing. The detection of specific microbes using RNA sequencing offers good sensitivity and specificity, but untargeted approaches suffer from high false positive rates and a lack of sensitivity for lowly abundant organisms.

RESULTS

We introduce Pathonoia, an algorithm that detects viruses and bacteria in RNA sequencing data with high precision and recall. Pathonoia first applies an established k-mer based method for species identification and then aggregates this evidence over all reads in a sample. In addition, we provide an easy-to-use analysis framework that highlights potential microbe-host interactions by correlating the microbial to the host gene expression. Pathonoia outperforms state-of-the-art methods in microbial detection specificity, both on in silico and real datasets.

CONCLUSION

Two case studies in human liver and brain show how Pathonoia can support novel hypotheses on microbial infection exacerbating disease. The Python package for Pathonoia sample analysis and a guided analysis Jupyter notebook for bulk RNAseq datasets are available on GitHub.

RNA-seq preprocessing and sample size considerations for gene network inference

Background

Gene network inference (GNI) methods have the potential to reveal functional relationships between different genes and their products. Most GNI algorithms have been developed for microarray gene expression datasets and their application to RNA-seq data is relatively recent. As the characteristics of RNA-seq data are different from microarray data, it is an unanswered question what preprocessing methods for RNA-seq data should be applied prior to GNI to attain optimal performance, or what the required sample size for RNA-seq data is to obtain reliable GNI estimates.

Results

We ran 9144 analysis of 7 different RNA-seq datasets to evaluate 300 different preprocessing combinations that include data transformations, normalizations and association estimators. We found that there was no single best performing preprocessing combination but that there were several good ones. The performance varied widely over various datasets, which emphasized the importance of choosing an appropriate preprocessing configuration before GNI. Two preprocessing combinations appeared promising in general: First, Log-2 TPM (transcript per million) with Variance-stabilizing transformation (VST) and Pearson Correlation Coefficient (PCC) association estimator. Second, raw RNA-seq count data with PCC. Along with these two, we also identified 18 other good preprocessing combinations. Any of these algorithms might perform best in different datasets. Therefore, the GNI performances of these approaches should be measured on any new dataset to select the best performing one for it. In terms of the required biological sample size of RNA-seq data, we found that between 30 to 85 samples were required to generate reliable GNI estimates.

Conclusions

This study provides practical recommendations on default choices for data preprocessing prior to GNI analysis of RNA-seq data to obtain optimal performance results.

A transcriptome software comparison for the analyses of treatments expected to give subtle gene expression responses

BACKGROUND

In this comparative study we evaluate the performance of four software tools: DNAstar-D (DESeq2), DNAstar-E (edgeR), CLC Genomics and Partek Flow for identification of differentially expressed genes (DEGs) using a transcriptome of E. coli. The RNA-seq data are from the effect of below-background radiation 5.5 nGy total dose (0.2nGy/hr) on E. coli grown shielded from natural radiation 655 m below ground in a pre-World War II steel vault. The gene expression response to three supplemented sources of radiation designed to mimic natural background, 1952 - 5720 nGy in total dose (71-208 nGy/hr), are compared to this "radiation-deprived" treatment. In addition, RNA-seq data of Caenorhabditis elegans nematode from similar radiation treatments was analyzed by three of the software packages.

RESULTS

In E. coli, the four software programs identified one of the supplementary sources of radiation (KCl) to evoke about 5 times more transcribed genes than the minus-radiation treatment (69-114 differentially expressed genes, DEGs), and so the rest of the analyses used this KCl vs "Minus" comparison. After imposing a 30-read minimum cutoff, one of the DNAStar options shared two of the three steps (mapping, normalization, and statistic) with Partek Flow (they both used median of ratios to normalize and the DESeq2 statistical package), and these two programs identified the highest number of DEGs in common with each other (53). In contrast, when the programs used different approaches in each of the three steps, between 31 and 40 DEGs were found in common. Regarding the extent of expression differences, three of the four programs gave high fold-change results (15-178 fold), but one (DNAstar's DESeq2) resulted in more conservative fold-changes (1.5-3.5). In a parallel study comparing three qPCR commercial validation software programs, these programs also gave variable results as to which genes were significantly regulated. Similarly, the C. elegans analysis showed exaggerated fold-changes in CLC and DNAstar's edgeR while DNAstar-D was more conservative.

CONCLUSIONS

Regarding the extent of expression (fold-change), and considering the subtlety of the very low level radiation treatments, in E. coli three of the four programs gave what we consider exaggerated fold-change results (15 - 178 fold), but one (DNAstar's DESeq2) gave more realistic fold-changes (1.5-3.5). When RT-qPCR validation comparisons to transcriptome results were carried out, they supported the more conservative DNAstar-D's expression results. When another model organism's (nematode) response to these radiation differences was similarly analyzed, DNAstar-D also resulted in the most conservative expression patterns. Therefore, we would propose DESeq2 ("DNAstar-D") as an appropriate software tool for differential gene expression studies for treatments expected to give subtle transcriptome responses.

OBIF: an omics-based interaction framework to reveal molecular drivers of synergy

Bioactive molecule library screening may empirically identify effective combination therapies, but molecular mechanisms underlying favorable drug–drug interactions often remain unclear, precluding further rational design. In the absence of an accepted systems theory to interrogate synergistic responses, we introduce Omics-Based Interaction Framework (OBIF) to reveal molecular drivers of synergy through integration of statistical and biological interactions in synergistic biological responses. OBIF performs full factorial analysis of feature expression data from single versus dual exposures to identify molecular clusters that reveal synergy-mediating pathways, functions and regulators. As a practical demonstration, OBIF analyzed transcriptomic and proteomic data of a dyad of immunostimulatory molecules that induces synergistic protection against influenza A and revealed unanticipated NF-κB/AP-1 cooperation that is required for antiviral protection. To demonstrate generalizability, OBIF analyzed data from a diverse array of Omics platforms and experimental conditions, successfully identifying the molecular clusters driving their synergistic responses. Hence, unlike existing synergy quantification and prediction methods, OBIF is a phenotype-driven systems model that supports multiplatform interrogation of synergy mechanisms.

Systematic Identification of Microproteins during the Development of Drosophila melanogaster

Short open reading frame-encoded peptides (SEPs) are microproteins with less than 100 amino acids that play an essential role in the growth and development of organisms. There are plenty of short open reading frames in Drosophila melanogaster that potentially code polypeptides. We chose 11 time points during the life cycle of Drosophila to investigate microproteins, particularly those related to development. Finally, we identified a total of 410 microproteins, of which 27 were noncoding RNA-encoded proteins. Of the 410 microproteins, 74 were expressed in all stages from embryo to adults, whereas 300 microproteins were only found in one or two time points. Approximately, one-third of the microproteins were not reported previously and 44 were obtained from de novo sequencing, validated by synthetic peptides. These microproteins are related to the main bioprocesses of growth and development, such as multicellular organism reproduction, postmating behavior, and oviposition. Over half of the microproteins have predicted functional domains and are conserved across species, suggesting that these microproteins have critical functions in fly development. This work enriches the D. melanogaster proteome and provides a significant data resource for growth and development research.

cdev: a ground-truth based measure to evaluate RNA-seq normalization performance

Normalization of RNA-seq data has been an active area of research since the problem was first recognized a decade ago. Despite the active development of new normalizers, their performance measures have been given little attention. To evaluate normalizers, researchers have been relying on ad hoc measures, most of which are either qualitative, potentially biased, or easily confounded by parametric choices of downstream analysis. We propose a metric called condition-number based deviation, or cdev, to quantify normalization success. cdev measures how much an expression matrix differs from another. If a ground truth normalization is given, cdev can then be used to evaluate the performance of normalizers. To establish experimental ground truth, we compiled an extensive set of public RNA-seq assays with external spike-ins. This data collection, together with cdev, provides a valuable toolset for benchmarking new and existing normalization methods.

CAMPAREE: a robust and configurable RNA expression simulator

Background

The accurate interpretation of RNA-Seq data presents a moving target as scientists continue to introduce new experimental techniques and analysis algorithms. Simulated datasets are an invaluable tool to accurately assess the performance of RNA-Seq analysis methods. However, existing RNA-Seq simulators focus on modeling the technical biases and artifacts of sequencing, rather than on simulating the original RNA samples. A first step in simulating RNA-Seq is to simulate RNA.

Results

To fill this need, we developed the Configurable And Modular Program Allowing RNA Expression Emulation (CAMPAREE), a simulator using empirical data to simulate diploid RNA samples at the level of individual molecules. We demonstrated CAMPAREE’s use for generating idealized coverage plots from real data, and for adding the ability to generate allele-specific data to existing RNA-Seq simulators that do not natively support this feature.

Conclusions

Separating input sample modeling from library preparation/sequencing offers added flexibility for both users and developers to mix-and-match different sample and sequencing simulators to suit their specific needs. Furthermore, the ability to maintain sample and sequencing simulators independently provides greater agility to incorporate new biological findings about transcriptomics and new developments in sequencing technologies. Additionally, by simulating at the level of individual molecules, CAMPAREE has the potential to model molecules transcribed from the same genes as a heterogeneous population of transcripts with different states of degradation and processing (splicing, editing, etc.). CAMPAREE was developed in Python, is open source, and freely available at https://github.com/itmat/CAMPAREE.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12864-021-07934-2.

ANPELA: analysis and performance assessment of the label-free quantification workflow for metaproteomic studies

Label-free quantification (LFQ) with a specific and sequentially integrated workflow of acquisition technique, quantification tool and processing method has emerged as the popular technique employed in metaproteomic research to provide a comprehensive landscape of the adaptive response of microbes to external stimuli and their interactions with other organisms or host cells. The performance of a specific LFQ workflow is highly dependent on the studied data. Hence, it is essential to discover the most appropriate one for a specific data set. However, it is challenging to perform such discovery due to the large number of possible workflows and the multifaceted nature of the evaluation criteria. Herein, a web server ANPELA (https://idrblab.org/anpela/) was developed and validated as the first tool enabling performance assessment of whole LFQ workflow (collective assessment by five well-established criteria with distinct underlying theories), and it enabled the identification of the optimal LFQ workflow(s) by a comprehensive performance ranking. ANPELA not only automatically detects the diverse formats of data generated by all quantification tools but also provides the most complete set of processing methods among the available web servers and stand-alone tools. Systematic validation using metaproteomic benchmarks revealed ANPELA's capabilities in 1 discovering well-performing workflow(s), (2) enabling assessment from multiple perspectives and (3) validating LFQ accuracy using spiked proteins. ANPELA has a unique ability to evaluate the performance of whole LFQ workflow and enables the discovery of the optimal LFQs by the comprehensive performance ranking of all 560 workflows. Therefore, it has great potential for applications in metaproteomic and other studies requiring LFQ techniques, as many features are shared among proteomic studies.

Effect of high variation in transcript expression on identifying differentially expressed genes in RNA-seq analysis

Great efforts have been made on the algorithms that deal with RNA-seq data to enhance the accuracy and efficiency of differential expression (DE) analysis. However, no consensus has been reached on the proper threshold values of fold change and adjusted p-value for filtering differentially expressed genes (DEGs). It is generally believed that the more stringent the filtering threshold, the more reliable the result of a DE analysis. Nevertheless, by analyzing the impact of both adjusted p-value and fold change thresholds on DE analyses, with RNA-seq data obtained for three different cancer types from the Cancer Genome Atlas (TCGA) database, we found that, for a given sample size, the reproducibility of DE results became poorer when more stringent thresholds were applied. No matter which threshold level was applied, the overlap rates of DEGs were generally lower for small sample sizes than for large sample sizes. The raw read count analysis demonstrated that the transcript expression of the same gene in different samples, whether in tumor groups or in normal groups, showed high variations, which resulted in a drastic fluctuation in fold change values and adjustedp-values when different sets of samples were used. Overall, more stringent thresholds did not yield more reliable DEGs due to high variations in transcript expression; the reliability of DEGs obtained with small sample sizes was more susceptible to these variations. Therefore, less stringent thresholds are recommended for screening DEGs. Moreover, large sample sizes should be considered in RNA-seq experimental designs to reduce the interfering effect of variations in transcript expression on DEG identification.

Human Placental Trophoblasts Infected by Listeria monocytogenes Undergo a Pro-Inflammatory Switch Associated With Poor Pregnancy Outcomes

The placenta controls the growth of the fetus and ensures its immune protection. Key to these functions, the syncytiotrophoblast (SYN) is a syncytium formed by fusion of underlying mononuclear trophoblasts. The SYN covers the placental surface and is bathed in maternal blood to mediate nutritional and waste exchanges between the mother and fetus. The bacterial pathogen Listeria monocytogenes breaches the trophoblast barrier and infects the placental/fetal unit resulting in poor pregnancy outcomes. In this work, we analyzed the L. monocytogenes intracellular lifecycle in primary human trophoblasts. In accordance with previous studies, we found that the SYN is 20-fold more resistant to infection compared to mononuclear trophoblasts, forming a protective barrier to infection at the maternal interface. We show for the first time that this is due to a significant reduction in L. monocytogenes uptake by the SYN rather than inhibition of the bacterial intracellular division or motility. We here report the first transcriptomic analysis of L. monocytogenes-infected trophoblasts (RNA sequencing). Pathway analysis showed that infection upregulated TLR2, NOD-like, and cytosolic DNA sensing pathways, as well as downstream pro-inflammatory circuitry (NF-κB, AP-1, IRF4, IRF7) leading to the production of mediators known to elicit the recruitment and activation of maternal leukocytes (IL8, IL6, TNFα, MIP-1). Signature genes associated with poor pregnancy outcomes were also upregulated upon infection. Measuring the release of 54 inflammatory mediators confirmed the transcriptomic data and revealed sustained production of tolerogenic factors (IL-27, IL-10, IL-1RA, TSLP) despite infection. Both the SYN and mononuclear trophoblasts produced cytokines, but surprisingly, some cytokines were predominantly produced by the SYN (IL-8, IL-6) or by non-fused trophoblasts (TNFα). Collectively, our data support that trophoblasts act as placental gatekeepers that limit and detect L. monocytogenes infection resulting in a pro-inflammatory response, which may contribute to the poor pregnancy outcomes if the pathogen persists.

Metatranscriptomics for the Human Microbiome and Microbial Community Functional Profiling

Shotgun metatranscriptomics (MTX) is an increasingly practical way to survey microbial community gene function and regulation at scale. This review begins by summarizing the motivations for community transcriptomics and the history of the field. We then explore the principles, best practices, and challenges of contemporary MTX workflows: beginning with laboratory methods for isolation and sequencing of community RNA, followed by informatics methods for quantifying RNA features, and finally statistical methods for detecting differential expression in a community context. In thesecond half of the review, we survey important biological findings from the MTX literature, drawing examples from the human microbiome, other (nonhuman) host-associated microbiomes, and the environment. Across these examples, MTX methods prove invaluable for probing microbe-microbe and host-microbe interactions, the dynamics of energy harvest and chemical cycling, and responses to environmental stresses. We conclude with a review of open challenges in the MTX field, including making assays and analyses more robust, accessible, and adaptable to new technologies; deciphering roles for millions of uncharacterized microbial transcripts; and solving applied problems such as biomarker discovery and development of microbial therapeutics.

Statistical approaches for differential expression analysis in metatranscriptomics

Motivation

Metatranscriptomics (MTX) has become an increasingly practical way to profile the functional activity of microbial communities in situ. However, MTX remains underutilized due to experimental and computational limitations. The latter are complicated by non-independent changes in both RNA transcript levels and their underlying genomic DNA copies (as microbes simultaneously change their overall abundance in the population and regulate individual transcripts), genetic plasticity (as whole loci are frequently gained and lost in microbial lineages) and measurement compositionality and zero-inflation. Here, we present a systematic evaluation of and recommendations for differential expression (DE) analysis in MTX.

Results

We designed and assessed six statistical models for DE discovery in MTX that incorporate different combinations of DNA and RNA normalization and assumptions about the underlying changes of gene copies or species abundance within communities. We evaluated these models on multiple simulated and real multi-omic datasets. Models adjusting transcripts relative to their encoding gene copies as a covariate were significantly more accurate in identifying DE from MTX in both simulated and real datasets. Moreover, we show that when paired DNA measurements (metagenomic data) are not available, models normalizing MTX measurements within-species while also adjusting for total-species RNA balance sensitivity, specificity and interpretability of DE detection, as does filtering likely technical zeros. The efficiency and accuracy of these models pave the way for more effective MTX-based DE discovery in microbial communities.

Availability and implementation

The analysis code and synthetic datasets used in this evaluation are available online at http://huttenhower.sph.harvard.edu/mtx2021.

Supplementary information

Supplementary data are available at Bioinformatics online.

SeqNet: An R Package for Generating Gene-Gene Networks and Simulating RNA-Seq Data

Gene expression data provide an abundant resource for inferring connections in gene regulatory networks. While methodologies developed for this task have shown success, a challenge remains in comparing the performance among methods. Gold-standard datasets are scarce and limited in use. And while tools for simulating expression data are available, they are not designed to resemble the data obtained from RNA-seq experiments. SeqNet is an R package that provides tools for generating a rich variety of gene network structures and simulating RNA-seq data from them. This produces in silico RNA-seq data for benchmarking and assessing gene network inference methods. The package is available on CRAN and on GitHub at https://github.com/tgrimes/SeqNet.

TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-seq Data from the NCI Patient-Derived Models Repository

Background

In order to correctly decode phenotypic information from RNA-sequencing (RNA-seq) data, careful selection of the RNA-seq quantification measure is critical for inter-sample comparisons and for downstream analyses, such as differential gene expression between two or more conditions. Several methods have been proposed and continue to be used. However, a consensus has not been reached regarding the best gene expression quantification method for RNA-seq data analysis.

Methods

In the present study, we used replicate samples from each of 20 patient-derived xenograft (PDX) models spanning 15 tumor types, for a total of 61 human tumor xenograft samples available through the NCI patient-derived model repository (PDMR). We compared the reproducibility across replicate samples based on TPM (transcripts per million), FPKM (fragments per kilobase of transcript per million fragments mapped), and normalized counts using coefficient of variation, intraclass correlation coefficient, and cluster analysis.

Results

Our results revealed that hierarchical clustering on normalized count data tended to group replicate samples from the same PDX model together more accurately than TPM and FPKM data. Furthermore, normalized count data were observed to have the lowest median coefficient of variation (CV), and highest intraclass correlation (ICC) values across all replicate samples from the same model and for the same gene across all PDX models compared to TPM and FPKM data.

Conclusion

We provided compelling evidence for a preferred quantification measure to conduct downstream analyses of PDX RNA-seq data. To our knowledge, this is the first comparative study of RNA-seq data quantification measures conducted on PDX models, which are known to be inherently more variable than cell line models. Our findings are consistent with what others have shown for human tumors and cell lines and add further support to the thesis that normalized counts are the best choice for the analysis of RNA-seq data across samples.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12967-021-02936-w.

A systematic comparison of normalization methods for eQTL analysis

Expression quantitative trait loci (eQTL) analysis has been widely used in interpreting disease-associated loci through correlating genetic variant loci with the expression of specific genes. RNA-sequencing (RNA-Seq), which can quantify gene expression at the genome-wide level, is often used in eQTL identification. Since different normalization methods of gene expression have substantial impacts on RNA-seq downstream analysis, it is of great necessity to systematically compare the effects of these methods on eQTL identification. Here, by using RNA-seq and genotype data of four different cancers in The Cancer Genome Atlas (TCGA) database, we comprehensively evaluated the effect of eight commonly used normalization methods on eQTL identification. Our results showed that the application of different methods could cause 20-30% differences in the final results of eQTL identification. Among these methods, COUNT, Median of Ratio (MED) and Trimmed Mean of M-values (TMM) generated similar results for identifying eQTLs, while Fragments Per Kilobase Million (FPKM) or RANK produced more differential results compared with other methods. Based on the accuracy and receiver operating characteristic (ROC) curve, the TMM method was found to be the optimal method for normalizing gene expression data in eQTLs analysis. In addition, we also evaluated the performance of different pairwise combinations of these methods. As a result, compared with single normalization methods, the combination of methods can not only identify more cis-eQTLs, but also improve the performance of the ROC curve. Overall, this study provides a comprehensive comparison of normalization methods for identifying eQTLs from RNA-seq data, and proposes some practical recommendations for diverse scenarios.

Probing into the Effects of Grapevine Leafroll-Associated Viruses on the Physiology, Fruit Quality and Gene Expression of Grapes

Grapevine leafroll is one of the most widespread and highly destructive grapevine diseases that is responsible for great economic losses to the grape and wine industries throughout the world. Six distinct viruses have been implicated in this disease complex. They belong to three genera, all in the family Closteroviridae. For the sake of convenience, these viruses are named as grapevine leafroll-associated viruses (GLRaV-1, -2, -3, -4, -7, and -13). However, their etiological role in the disease has yet to be established. Furthermore, how infections with each GLRaV induce the characteristic disease symptoms remains unresolved. Here, we first provide a brief overview on each of these GLRaVs with a focus on genome structure, expression strategies and gene functions, where available. We then provide a review on the effects of GLRaV infection on the physiology, fruit quality, fruit chemical composition, and gene expression of grapevine based on the limited information so far reported in the literature. We outline key methodologies that have been used to study how GLRaV infections alter gene expression in the grapevine host at the transcriptomic level. Finally, we present a working model as an initial attempt to explain how infections with GLRaVs lead to the characteristic symptoms of grapevine leafroll disease: leaf discoloration and downward rolling. It is our hope that this review will serve as a starting point for grapevine virology and the related research community to tackle this vastly important and yet virtually uncharted territory in virus-host interactions involving woody and perennial fruit crops.

The overexpression of DNA repair genes in invasive ductal and lobular breast carcinomas: Insights on individual variations and precision medicine

In the era of precision medicine, analyzing the transcriptomic profile of patients is essential to tailor the appropriate therapy. In this study, we explored transcriptional differences between two invasive breast cancer subtypes; infiltrating ductal carcinoma (IDC) and lobular carcinoma (LC) using RNA-Seq data deposited in the TCGA-BRCA project. We revealed 3854 differentially expressed genes between normal ductal tissues and IDC. In addition, IDC to LC comparison resulted in 663 differentially expressed genes. We then focused on DNA repair genes because of their known effects on patients' response to therapy and resistance. We here report that 36 DNA repair genes are overexpressed in a significant number of both IDC and LC patients' samples. Despite the upregulation in a significant number of samples, we observed a noticeable variation in the expression levels of the repair genes across patients of the same cancer subtype. The same trend is valid for the expression of miRNAs, where remarkable variations between patients' samples of the same cancer subtype are also observed. These individual variations could lie behind the differential response of patients to treatment. The future of cancer diagnostics and therapy will inevitably depend on high-throughput genomic and transcriptomic data analysis. However, we propose that performing analysis on individual patients rather than a big set of patients' samples will be necessary to ensure that the best treatment is determined, and therapy resistance is reduced.

Advances in Current Diabetes Proteomics: From the Perspectives of Label- free Quantification and Biomarker Selection

BACKGROUND

Due to its prevalence and negative impacts on both the economy and society, the diabetes mellitus (DM) has emerged as a worldwide concern. In light of this, the label-free quantification (LFQ) proteomics and diabetic marker selection methods have been applied to elucidate the underlying mechanisms associated with insulin resistance, explore novel protein biomarkers, and discover innovative therapeutic protein targets.

OBJECTIVE

The purpose of this manuscript is to review and analyze the recent computational advances and development of label-free quantification and diabetic marker selection in diabetes proteomics.

METHODS

Web of Science database, PubMed database and Google Scholar were utilized for searching label-free quantification, computational advances, feature selection and diabetes proteomics.

RESULTS

In this study, we systematically review the computational advances of label-free quantification and diabetic marker selection methods which were applied to get the understanding of DM pathological mechanisms. Firstly, different popular quantification measurements and proteomic quantification software tools which have been applied to the diabetes studies are comprehensively discussed. Secondly, a number of popular manipulation methods including transformation, pretreatment (centering, scaling, and normalization), missing value imputation methods and a variety of popular feature selection techniques applied to diabetes proteomic data are overviewed with objective evaluation on their advantages and disadvantages. Finally, the guidelines for the efficient use of the computationbased LFQ technology and feature selection methods in diabetes proteomics are proposed.

CONCLUSION

In summary, this review provides guidelines for researchers who will engage in proteomics biomarker discovery and by properly applying these proteomic computational advances, more reliable therapeutic targets will be found in the field of diabetes mellitus.

The discovery, distribution, and diversity of DNA viruses associated with Drosophila melanogaster in Europe

Drosophila melanogaster is an important model for antiviral immunity in arthropods, but very few DNA viruses have been described from the family Drosophilidae. This deficiency limits our opportunity to use natural host-pathogen combinations in experimental studies, and may bias our understanding of the Drosophila virome. Here, we report fourteen DNA viruses detected in a metagenomic analysis of 6668 pool-sequenced Drosophila, sampled from forty-seven European locations between 2014 and 2016. These include three new nudiviruses, a new and divergent entomopoxvirus, a virus related to Leptopilina boulardi filamentous virus, and a virus related to Musca domestica salivary gland hypertrophy virus. We also find an endogenous genomic copy of galbut virus, a double-stranded RNA partitivirus, segregating at very low frequency. Remarkably, we find that Drosophila Vesanto virus, a small DNA virus previously described as a bidnavirus, may be composed of up to twelve segments and thus represent a new lineage of segmented DNA viruses. Two of the DNA viruses, Drosophila Kallithea nudivirus and Drosophila Vesanto virus are relatively common, found in 2 per cent or more of wild flies. The others are rare, with many likely to be represented by a single infected fly. We find that virus prevalence in Europe reflects the prevalence seen in publicly available datasets, with Drosophila Kallithea nudivirus and Drosophila Vesanto virus the only ones commonly detectable in public data from wild-caught flies and large population cages, and the other viruses being rare or absent. These analyses suggest that DNA viruses are at lower prevalence than RNA viruses in D.melanogaster, and may be less likely to persist in laboratory cultures. Our findings go some way to redressing an earlier bias toward RNA virus studies in Drosophila, and lay the foundation needed to harness the power of Drosophila as a model system for the study of DNA viruses.

Identification of Genes Contributing to a Long Circadian Period in Drosophila Melanogaster

The endogenous circadian period of animals and humans is typically very close to 24 h. Individuals with much longer circadian periods have been observed, however, and in the case of humans, these deviations have health implications. Previously, we observed a line of Drosophila with a very long average period of 31.3 h for locomotor activity behavior. Preliminary mapping indicated that the long period did not map to known canonical clock genes but instead mapped to multiple chromosomes. Using RNA-Seq, we surveyed the whole transcriptome of fly heads from this line across time and compared it with a wild-type control. A three-way generalized linear model revealed that approximately two-thirds of the genes were expressed differentially among the two genotypes, while only one quarter of the genes varied across time. Using these results, we applied algorithms to search for genes that oscillated over 24 h, identifying genes not previously known to cycle. We identified 166 differentially expressed genes that overlapped with a previous Genome-wide Association Study (GWAS) of circadian behavior, strongly implicating them in the long-period phenotype. We tested mutations in 45 of these genes for their effect on the circadian period. Mutations in Alk, alph, CG10089, CG42540, CG6034, Kairos (CG6123), CG8768, klg, Lar, sick, and tinc had significant effects on the circadian period, with seven of these mutations increasing the circadian period of locomotor activity behavior. Genetic rescue of mutant Kairos restored the circadian period to wild-type levels, suggesting it has a critical role in determining period length in constant darkness.

Systematic comparison and assessment of RNA-seq procedures for gene expression quantitative analysis

RNA-seq is currently considered the most powerful, robust and adaptable technique for measuring gene expression and transcription activation at genome-wide level. As the analysis of RNA-seq data is complex, it has prompted a large amount of research on algorithms and methods. This has resulted in a substantial increase in the number of options available at each step of the analysis. Consequently, there is no clear consensus about the most appropriate algorithms and pipelines that should be used to analyse RNA-seq data. In the present study, 192 pipelines using alternative methods were applied to 18 samples from two human cell lines and the performance of the results was evaluated. Raw gene expression signal was quantified by non-parametric statistics to measure precision and accuracy. Differential gene expression performance was estimated by testing 17 differential expression methods. The procedures were validated by qRT-PCR in the same samples. This study weighs up the advantages and disadvantages of the tested algorithms and pipelines providing a comprehensive guide to the different methods and procedures applied to the analysis of RNA-seq data, both for the quantification of the raw expression signal and for the differential gene expression.

Transcriptome profiling at osmotic and ionic phases of salt stress response in bread wheat uncovers trait-specific candidate genes

BACKGROUND

Bread wheat is one of the most important crops for the human diet, but the increasing soil salinization is causing yield reductions worldwide. Improving salt stress tolerance in wheat requires the elucidation of the mechanistic basis of plant response to this abiotic stress factor. Although several studies have been performed to analyze wheat adaptation to salt stress, there are still some gaps to fully understand the molecular mechanisms from initial signal perception to the onset of responsive tolerance pathways. The main objective of this study is to exploit the dynamic salt stress transcriptome in underlying QTL regions to uncover candidate genes controlling salt stress tolerance in bread wheat. The massive analysis of 3'-ends sequencing protocol was used to analyze leave samples at osmotic and ionic phases. Afterward, stress-responsive genes overlapping QTL for salt stress-related traits in two mapping populations were identified.

RESULTS

Among the over-represented salt-responsive gene categories, the early up-regulation of calcium-binding and cell wall synthesis genes found in the tolerant genotype are presumably strategies to cope with the salt-related osmotic stress. On the other hand, the down-regulation of photosynthesis-related and calcium-binding genes, and the increased oxidative stress response in the susceptible genotype are linked with the greater photosynthesis inhibition at the osmotic phase. The specific up-regulation of some ABC transporters and Na+/Ca2+ exchangers in the tolerant genotype at the ionic stage indicates their involvement in mechanisms of sodium exclusion and homeostasis. Moreover, genes related to protein synthesis and breakdown were identified at both stress phases. Based on the linkage disequilibrium blocks, salt-responsive genes within QTL intervals were identified as potential components operating in pathways leading to salt stress tolerance. Furthermore, this study conferred evidence of novel regions with transcription in bread wheat.

CONCLUSION

The dynamic transcriptome analysis allowed the comparison of osmotic and ionic phases of the salt stress response and gave insights into key molecular mechanisms involved in the salt stress adaptation of contrasting bread wheat genotypes. The leveraging of the highly contiguous chromosome-level reference genome sequence assembly facilitated the QTL dissection by targeting novel candidate genes for salt tolerance.

Integrative, normalization-insusceptible statistical analysis of RNA-Seq data, with improved differential expression and unbiased downstream functional analysis

The study of differential gene expression patterns through RNA-Seq comprises a routine task in the daily lives of molecular bioscientists, who produce vast amounts of data requiring proper management and analysis. Despite widespread use, there are still no widely accepted golden standards for the normalization and statistical analysis of RNA-Seq data, and critical biases, such as gene lengths and problems in the detection of certain types of molecules, remain largely unaddressed. Stimulated by these unmet needs and the lack of in-depth research into the potential of combinatorial methods to enhance the analysis of differential gene expression, we had previously introduced the PANDORA P-value combination algorithm while presenting evidence for PANDORA's superior performance in optimizing the tradeoff between precision and sensitivity. In this article, we present the next generation of the algorithm along with a more in-depth investigation of its capabilities to effectively analyze RNA-Seq data. In particular, we show that PANDORA-reported lists of differentially expressed genes are unaffected by biases introduced by different normalization methods, while, at the same time, they comprise a reliable input option for downstream pathway analysis. Additionally, PANDORA outperforms other methods in detecting differential expression patterns in certain transcript types, including long non-coding RNAs.

Comparison of Normalization Methods for Analysis of TempO-Seq Targeted RNA Sequencing Data

Analysis of bulk RNA sequencing (RNA-Seq) data is a valuable tool to understand transcription at the genome scale. Targeted sequencing of RNA has emerged as a practical means of assessing the majority of the transcriptomic space with less reliance on large resources for consumables and bioinformatics. TempO-Seq is a templated, multiplexed RNA-Seq platform that interrogates a panel of sentinel genes representative of genome-wide transcription. Nuances of the technology require proper preprocessing of the data. Various methods have been proposed and compared for normalizing bulk RNA-Seq data, but there has been little to no investigation of how the methods perform on TempO-Seq data. We simulated count data into two groups (treated vs. untreated) at seven-fold change (FC) levels (including no change) using control samples from human HepaRG cells run on TempO-Seq and normalized the data using seven normalization methods. Upper Quartile (UQ) performed the best with regard to maintaining FC levels as detected by a limma contrast between treated vs. untreated groups. For all FC levels, specificity of the UQ normalization was greater than 0.84 and sensitivity greater than 0.90 except for the no change and +1.5 levels. Furthermore, K-means clustering of the simulated genes normalized by UQ agreed the most with the FC assignments [adjusted Rand index (ARI) = 0.67]. Despite having an assumption of the majority of genes being unchanged, the DESeq2 scaling factors normalization method performed reasonably well as did simple normalization procedures counts per million (CPM) and total counts (TCs). These results suggest that for two class comparisons of TempO-Seq data, UQ, CPM, TC, or DESeq2 normalization should provide reasonably reliable results at absolute FC levels ≥2.0. These findings will help guide researchers to normalize TempO-Seq gene expression data for more reliable results.

SAGD: a comprehensive sex-associated gene database from transcriptomes

Many animal species present sex differences. Sex-associated genes (SAGs), which have female-biased or male-biased expression, have major influences on the remarkable sex differences in important traits such as growth, reproduction, disease resistance and behaviors. However, the SAGs resulting in the vast majority of phenotypic sex differences are still unknown. To provide a useful resource for the functional study of SAGs, we manually curated public RNA-seq datasets with paired female and male biological replicates from the same condition and systematically re-analyzed the datasets using standardized methods. We identified 27,793 female-biased SAGs and 64,043 male-biased SAGs from 2,828 samples of 21 species, including human, chimpanzee, macaque, mouse, rat, cow, horse, chicken, zebrafish, seven fly species and five worm species. All these data were cataloged into SAGD, a user-friendly database of SAGs (http://bioinfo.life.hust.edu.cn/SAGD) where users can browse SAGs by gene, species, drug and dataset. In SAGD, the expression, annotation, targeting drugs, homologs, ontology and related RNA-seq datasets of SAGs are provided to help researchers to explore their functions and potential applications in agriculture and human health.

Benchmarking RNA-seq differential expression analysis methods using spike-in and simulation data

Benchmarking RNA-seq differential expression analysis methods using spike-in and simulated RNA-seq data has often yielded inconsistent results. The spike-in data, which were generated from the same bulk RNA sample, only represent technical variability, making the test results less reliable. We compared the performance of 12 differential expression analysis methods for RNA-seq data, including recent variants in widely used software packages, using both RNA spike-in and simulation data for negative binomial (NB) model. Performance of edgeR, DESeq2, and ROTS was particularly different between the two benchmark tests. Then, each method was tested under most extensive simulation conditions especially demonstrating the large impacts of proportion, dispersion, and balance of differentially expressed (DE) genes. DESeq2, a robust version of edgeR (edgeR.rb), voom with TMM normalization (voom.tmm) and sample weights (voom.sw) showed an overall good performance regardless of presence of outliers and proportion of DE genes. The performance of RNA-seq DE gene analysis methods substantially depended on the benchmark used. Based on the simulation results, suitable methods were suggested under various test conditions.

The Multifunctional Long-Distance Movement Protein of Pea Enation Mosaic Virus 2 Protects Viral and Host Transcripts from Nonsense-Mediated Decay

The nonsense-mediated decay (NMD) pathway presents a challenge for RNA viruses with termination codons that precede extended 3' untranslated regions (UTRs). The umbravirus Pea enation mosaic virus 2 (PEMV2) is a nonsegmented, positive-sense RNA virus with an unusually long 3' UTR that is susceptible to NMD. To establish a systemic infection, the PEMV2 long-distance movement protein p26 was previously shown to both stabilize viral RNAs and bind them for transport through the plant's vascular system. The current study demonstrated that p26 protects both viral and nonviral messenger RNAs from NMD. Although p26 localizes to both the cytoplasm and nucleolus, p26 exerts its anti-NMD effects exclusively in the cytoplasm independently of long-distance movement. Using a transcriptome-wide approach in the model plant Nicotiana benthamiana, p26 protected a subset of cellular NMD target transcripts, particularly those containing long, structured, GC-rich 3' UTRs. Furthermore, transcriptome sequencing (RNA-seq) revealed that the NMD pathway is highly dysfunctional during PEMV2 infection, with 1,820 (48%) of NMD targets increasing in abundance. Widespread changes in the host transcriptome are common during plant RNA virus infections, and these results suggest that, in at least some instances, virus-mediated NMD inhibition may be a major contributing factor.IMPORTANCE Nonsense-mediated decay (NMD) represents an RNA regulatory pathway that degrades both natural and faulty messenger RNAs with long 3' untranslated regions. NMD targets diverse families of RNA viruses, requiring that viruses counteract the NMD pathway for successful amplification in host cells. A protein required for long-distance movement of Pea enation mosaic virus 2 (PEMV2) is shown to also protect both viral and host mRNAs from NMD. RNA-seq analyses of the Nicotiana benthamiana transcriptome revealed that PEMV2 infection significantly impairs the host NMD pathway. RNA viruses routinely induce large-scale changes in host gene expression, and, like PEMV2, may use NMD inhibition to alter the host transcriptome in an effort to increase virus amplification.

Evaluation of Seven Different RNA-Seq Alignment Tools Based on Experimental Data from the Model Plant Arabidopsis thaliana

Quantification of gene expression is crucial to connect genome sequences with phenotypic and physiological data. RNA-Sequencing (RNA-Seq) has taken a prominent role in the study of transcriptomic reactions of plants to various environmental and genetic perturbations. However, comparative tests of different tools for RNA-Seq read mapping and quantification have been mainly performed on data from animals or humans, which necessarily neglect, for example, the large genetic variability among natural accessions within plant species. Here, we compared seven computational tools for their ability to map and quantify Illumina single-end reads from the Arabidopsis thaliana accessions Columbia-0 (Col-0) and N14. Between 92.4% and 99.5% of all reads were mapped to the reference genome or transcriptome and the raw count distributions obtained from the different mappers were highly correlated. Using the software DESeq2 to determine differential gene expression (DGE) between plants exposed to 20 °C or 4 °C from these read counts showed a large pairwise overlap between the mappers. Interestingly, when the commercial CLC software was used with its own DGE module instead of DESeq2, strongly diverging results were obtained. All tested mappers provided highly similar results for mapping Illumina reads of two polymorphic Arabidopsis accessions to the reference genome or transcriptome and for the determination of DGE when the same software was used for processing.

DFseq: Distribution-Free Method to Detect Differential Gene Expression for RNA-Sequencing Data

Many current RNA-sequencing data analysis methods compare expressions one gene at a time, taking little consideration of the correlations among genes. In this study, we propose a method to convert such an one-dimensional comparison approach into a two-dimensional evaluation of the ratio of standard deviations (SD) of two constructed random variables. This method allows the identification of differentially expressed genes while controlling a preset significance level conditional on the read count mean-variance relationship. Meanwhile, correlations among genes are naturally accommodated due to the clustering of genes with similar distribution in the proposed σ-σ plot. The proposed distribution-free method is designated as DFseq, because it does not depend on a parametric distribution to fit read count. As a result, compared with parametric methods, DFseq can effectively handle genes with a bimodal-like distribution and/or genes with excessive 0 read counts, as well as genes with outlying observations. Besides, DFseq is an ideal platform for comparing performance of different differential gene expression detection methods.

Choice of library size normalization and statistical methods for differential gene expression analysis in balanced two-group comparisons for RNA-seq studies

Background

High-throughput RNA sequencing (RNA-seq) has evolved as an important analytical tool in molecular biology. Although the utility and importance of this technique have grown, uncertainties regarding the proper analysis of RNA-seq data remain. Of primary concern, there is no consensus regarding which normalization and statistical methods are the most appropriate for analyzing this data. The lack of standardized analytical methods leads to uncertainties in data interpretation and study reproducibility, especially with studies reporting high false discovery rates. In this study, we compared a recently developed normalization method, UQ-pgQ2, with three of the most frequently used alternatives including RLE (relative log estimate), TMM (Trimmed-mean M values) and UQ (upper quartile normalization) in the analysis of RNA-seq data. We evaluated the performance of these methods for gene-level differential expression analysis by considering the factors, including: 1) normalization combined with the choice of a Wald test from DESeq2 and an exact test/QL (Quasi-likelihood) F-Test from edgeR; 2) sample sizes in two balanced two-group comparisons; and 3) sequencing read depths.

Results

Using the MAQC RNA-seq datasets with small sample replicates, we found that UQ-pgQ2 normalization combined with an exact test can achieve better performance in term of power and specificity in differential gene expression analysis. However, using an intra-group analysis of false positives from real and simulated data, we found that a Wald test performs better than an exact test when the number of sample replicates is large and that a QL F-test performs the best given sample sizes of 5, 10 and 15 for any normalization. The RLE, TMM and UQ methods performed similarly given a desired sample size.

Conclusion

We found the UQ-pgQ2 method combined with an exact test/QL F-test is the best choice in order to control false positives when the sample size is small. When the sample size is large, UQ-pgQ2 with a QL F-test is a better choice for the type I error control in an intra-group analysis. We observed read depths have a minimal impact for differential gene expression analysis based on the simulated data.

Single-Cell Transcriptomes Distinguish Stem Cell State Changes and Lineage Specification Programs in Early Mammary Gland Development

The mammary gland consists of cells with gene expression patterns reflecting their cellular origins, function, and spatiotemporal context. However, knowledge of developmental kinetics and mechanisms of lineage specification is lacking. We address this significant knowledge gap by generating a single-cell transcriptome atlas encompassing embryonic, postnatal, and adult mouse mammary development. From these data, we map the chronology of transcriptionally and epigenetically distinct cell states and distinguish fetal mammary stem cells (fMaSCs) from their precursors and progeny. fMaSCs show balanced co-expression of factors associated with discrete adult lineages and a metabolic gene signature that subsides during maturation but reemerges in some human breast cancers and metastases. These data provide a useful resource for illuminating mammary cell heterogeneity, the kinetics of differentiation, and developmental correlates of tumorigenesis.

Single-cell RNA sequencing of developing mouse mammary epithelia reveals the timing of lineage specification. Giraddi et al. find that fetal mammary stem cells co-express factors that define distinct lineages in their progeny and bear functionally relevant metabolic program signatures that change with differentiation and are resurrected in human breast cancers and metastases.

A graph-based algorithm for RNA-seq data normalization

The use of RNA-sequencing has garnered much attention in recent years for characterizing and understanding various biological systems. However, it remains a major challenge to gain insights from a large number of RNA-seq experiments collectively, due to the normalization problem. Normalization has been challenging due to an inherent circularity, requiring that RNA-seq data be normalized before any pattern of differential (or non-differential) expression can be ascertained; meanwhile, the prior knowledge of non-differential transcripts is crucial to the normalization process. Some methods have successfully overcome this problem by the assumption that most transcripts are not differentially expressed. However, when RNA-seq profiles become more abundant and heterogeneous, this assumption fails to hold, leading to erroneous normalization. We present a normalization procedure that does not rely on this assumption, nor prior knowledge about the reference transcripts. This algorithm is based on a graph constructed from intrinsic correlations among RNA-seq transcripts and seeks to identify a set of densely connected vertices as references. Application of this algorithm on our synthesized validation data showed that it could recover the reference transcripts with high precision, thus resulting in high-quality normalization. On a realistic data set from the ENCODE project, this algorithm gave good results and could finish in a reasonable time. These preliminary results imply that we may be able to break the long persisting circularity problem in RNA-seq normalization.

The role and robustness of the Gini coefficient as an unbiased tool for the selection of Gini genes for normalising expression profiling data

We recently introduced the Gini coefficient (GC) for assessing the expression variation of a particular gene in a dataset, as a means of selecting improved reference genes over the cohort ('housekeeping genes') typically used for normalisation in expression profiling studies. Those genes (transcripts) that we determined to be useable as reference genes differed greatly from previous suggestions based on hypothesis-driven approaches. A limitation of this initial study is that a single (albeit large) dataset was employed for both tissues and cell lines. We here extend this analysis to encompass seven other large datasets. Although their absolute values differ a little, the Gini values and median expression levels of the various genes are well correlated with each other between the various cell line datasets, implying that our original choice of the more ubiquitously expressed low-Gini-coefficient genes was indeed sound. In tissues, the Gini values and median expression levels of genes showed a greater variation, with the GC of genes changing with the number and types of tissues in the data sets. In all data sets, regardless of whether this was derived from tissues or cell lines, we also show that the GC is a robust measure of gene expression stability. Using the GC as a measure of expression stability we illustrate its utility to find tissue- and cell line-optimised housekeeping genes without any prior bias, that again include only a small number of previously reported housekeeping genes. We also independently confirmed this experimentally using RT-qPCR with 40 candidate GC genes in a panel of 10 cell lines. These were termed the Gini Genes. In many cases, the variation in the expression levels of classical reference genes is really quite huge (e.g. 44 fold for GAPDH in one data set), suggesting that the cure (of using them as normalising genes) may in some cases be worse than the disease (of not doing so). We recommend the present data-driven approach for the selection of reference genes by using the easy-to-calculate and robust GC.

Isolating Male Meiocytes from Maize and Wheat for "-Omics" Analyses

Genome-wide gene expression studies have become a routine approach due to the advances in sequencing technologies, their ease of use, and increasing affordability. Simultaneous investigation of small RNA expression adds further valuable information but is not adopted as widely yet. Both RNA-seq and small RNA-seq benefit from the use of specific cell types. Here, we describe a protocol for the isolation of male meiotic cells from maize or wheat plants, along with the application of downstream RNA sequencing, extendable to other -omics approaches.

Understanding sequencing data as compositions: an outlook and review

Motivation

Although seldom acknowledged explicitly, count data generated by sequencing platforms exist as compositions for which the abundance of each component (e.g. gene or transcript) is only coherently interpretable relative to other components within that sample. This property arises from the assay technology itself, whereby the number of counts recorded for each sample is constrained by an arbitrary total sum (i.e. library size). Consequently, sequencing data, as compositional data, exist in a non-Euclidean space that, without normalization or transformation, renders invalid many conventional analyses, including distance measures, correlation coefficients and multivariate statistical models.

Results

The purpose of this review is to summarize the principles of compositional data analysis (CoDA), provide evidence for why sequencing data are compositional, discuss compositionally valid methods available for analyzing sequencing data, and highlight future directions with regard to this field of study.

Supplementary information

Supplementary data are available at Bioinformatics online.

Machine Learning Classifiers for Endometriosis Using Transcriptomics and Methylomics Data

Endometriosis is a complex and common gynecological disorder yet a poorly understood disease affecting about 176 million women worldwide and causing significant impact on their quality of life and economic burden. Neither a definitive clinical symptom nor a minimally invasive diagnostic method is available, thus leading to an average of 4 to 11 years of diagnostic latency. Discovery of relevant biological patterns from microarray expression or next generation sequencing (NGS) data has been advanced over the last several decades by applying various machine learning tools. We performed machine learning analysis using 38 RNA-seq and 80 enrichment-based DNA methylation (MBD-seq) datasets. We experimented how well various supervised machine learning methods such as decision tree, partial least squares discriminant analysis (PLSDA), support vector machine, and random forest perform in classifying endometriosis from the control samples trained on both transcriptomics and methylomics data. The assessment was done from two different perspectives for improving classification performances: a) implication of three different normalization techniques and b) implication of differential analysis using the generalized linear model (GLM). Several candidate biomarker genes were identified by multiple machine learning experiments including NOTCH3, SNAPC2, B4GALNT1, SMAP2, DDB2, GTF3C5, and PTOV1 from the transcriptomics data analysis and TRPM6, RASSF2, TNIP2, RP3-522J7.6, FGD3, and MFSD14B from the methylomics data analysis. We concluded that an appropriate machine learning diagnostic pipeline for endometriosis should use TMM normalization for transcriptomics data, and quantile or voom normalization for methylomics data, GLM for feature space reduction and classification performance maximization.