Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data

Assessing the impact of transcriptomics data analysis pipelines on downstream functional enrichment results

Transcriptomics is widely used to assess the state of biological systems. There are many tools for the different steps, such as normalization, differential expression, and enrichment. While numerous studies have examined the impact of method choices on differential expression results, little attention has been paid to their effects on further downstream functional analysis, which typically provides the basis for interpretation and follow-up experiments. To address this, we introduce FLOP, a comprehensive nextflow-based workflow combining methods to perform end-to-end analyses of transcriptomics data. We illustrate FLOP on datasets ranging from end-stage heart failure patients to cancer cell lines. We discovered effects not noticeable at the gene-level, and observed that not filtering the data had the highest impact on the correlation between pipelines in the gene set space. Moreover, we performed three benchmarks to evaluate the 12 pipelines included in FLOP, and confirmed that filtering is essential in scenarios of expected moderate-to-low biological signal. Overall, our results underscore the impact of carefully evaluating the consequences of the choice of preprocessing methods on downstream enrichment analyses. We envision FLOP as a valuable tool to measure the robustness of functional analyses, ultimately leading to more reliable and conclusive biological findings.

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.

Viromes of Haemaphysalis longicornis reveal different viral abundance and diversity in free and engorged ticks

Haemaphysalis longicornis ticks, commonly found in East Asia, can transmit various pathogenic viruses, including the severe fever with thrombocytopenia syndrome virus (SFTSV) that has caused febrile diseases among humans in Hubei Province. However, understanding of the viromes of H. longicornis was limited, and the prevalence of viruses among H. longicornis ticks in Hubei was not well clarified. This study investigates the viromes of both engorged (fed) and free (unfed) H. longicornis ticks across three mountainous regions in Hubei Province from 2019 to 2020. RNA-sequencing analysis identified viral sequences that were related to 39 reference viruses belonging to unclassified viruses and seven RNA viral families, namely Chuviridae, Nairoviridae, Orthomyxoviridae, Parvoviridae, Phenuiviridae, Rhabdoviridae, and Totiviridae. Viral abundance and diversity in these ticks were analysed, and phylogenetic characteristics of the Henan tick virus (HNTV), Dabieshan tick virus (DBSTV), Okutama tick virus (OKTV), and Jingmen tick virus (JMTV) were elucidated based on their full genomic sequences. Prevalence analysis demonstrated that DBSTV was the most common virus found in individual H. longicornis ticks (12.59%), followed by HNTV (0.35%), whereas JMTV and OKTV were not detected. These results improve our understanding of H. longicornis tick viromes in central China and highlight the role of tick feeding status and geography in shaping the viral community. The findings of new viral strains and their potential impact on public health raise the need to strengthen surveillance efforts for comprehensively assessing their spillover potentials.

Normalization of RNA-Seq data using adaptive trimmed mean with multi-reference

The normalization of RNA sequencing data is a primary step for downstream analysis. The most popular method used for the normalization is the trimmed mean of M values (TMM) and DESeq. The TMM tries to trim away extreme log fold changes of the data to normalize the raw read counts based on the remaining non-deferentially expressed genes. However, the major problem with the TMM is that the values of trimming factor M are heuristic. This paper tries to estimate the adaptive value of M in TMM based on Jaeckel's Estimator, and each sample acts as a reference to find the scale factor of each sample. The presented approach is validated on SEQC, MAQC2, MAQC3, PICKRELL and two simulated datasets with two-group and three-group conditions by varying the percentage of differential expression and the number of replicates. The performance of the present approach is compared with various state-of-the-art methods, and it is better in terms of area under the receiver operating characteristic curve and differential expression.

Physiological Characteristic Changes and Transcriptome Analysis of Maize (Zea mays L.) Roots under Drought Stress

Water deficit is a key limiting factor for limiting yield in maize (Zea mays L.). It is crucial to elucidate the molecular regulatory networks of stress tolerance for genetic enhancement of drought tolerance. The mechanism of drought tolerance of maize was explored by comparing physiological and transcriptomic data under normal conditions and drought treatment at polyethylene glycol- (PEG-) induced drought stress (5%, 10%, 15%, and 20%) in the root during the seedling stage. The content of saccharide, SOD, CAT, and MDA showed an upward trend, proteins showed a downward trend, and the levels of POD first showed an upward trend and then decreased. Compared with the control group, a total of 597, 2748, 6588, and 5410 differentially expressed genes were found at 5%, 10%, 15%, and 20% PEG, respectively, and 354 common DEGs were identified in these comparisons. Some differentially expressed genes were remarkably enriched in the MAPK signaling pathway and plant hormone signal transduction. The 50 transcription factors (TFs) divided into 15 categories were screened from the 354 common DEGs during drought stress. Auxin response factor 10 (ARF10), auxin-responsive protein IAA9 (IAA9), auxin response factor 14 (ARF14), auxin-responsive protein IAA1 (IAA1), auxin-responsive protein IAA27 (IAA27), and 1 ethylene response sensor 2 (ERS2) were upregulated. The two TFs, including bHLH 35 and bHLH 96, involved in the MAPK signal pathway and plant hormones pathway, are significantly upregulated in 5%, 10%, 15%, and 20% PEG stress groups. The present study provides greater insight into the fundamental transcriptome reprogramming of grain crops under drought.

Genomic and phylogenetic profiling of RNA of tick-borne arboviruses in Hainan Island, China

Ticks act as vectors and hosts of numerous arboviruses. Examples of medically important arboviruses include the tick-borne encephalitis virus, Crimean Congo hemorrhagic fever, and severe fever with thrombocytopenia syndrome. Recently, some novel arboviruses have been identified in blood specimens of patients with unexplained fever and a history of tick bites in Inner Mongolia. Consequently, tick-borne viruses are a major focus of infectious disease research. However, the spectrum of tick-borne viruses in subtropical areas of China has yet to be sufficiently characterized. In this study, we collected 855 ticks from canine and bovine hosts in four locations in Hainan Province. The ticks were combined into 18 pools according to genus and location. Viral RNA-sequence libraries were subjected to transcriptome sequencing analysis. Molecular clues from metagenomic analyses were used to classify sequence reads into virus species, genera, or families. The diverse viral reads closely associated with mammals were assigned to 12 viral families and important tick-borne viruses, such as Jingmen, Beiji nairovirus, and Colorado tick fever. Our virome and phylogenetic analyses of the arbovirus strains provide basic data for preventing and controlling human infectious diseases caused by tick-borne viruses in the subtropical areas of China.

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.

VDR promotes testosterone synthesis in mouse Leydig cells via regulation of cholesterol side chain cleavage cytochrome P450 (Cyp11a1) expression

BACKGROUND

The vitamin D receptor (VDR) mediates the pleiotropic biological actions that include osteoporosis, immune responses and androgen synthesis.VDR is widely expressed in testis cells such as Leydig cells, Sertoli cells, and sperm. The levels of steroids are critical for sexual development. In the early stage of steroidogenesis, cholesterol is converted to pregnenolone (precursor of most steroid hormones) by cholesterol side-chain lyase (CYP11A1), which eventually synthesizes the male hormone testosterone.

OBJECTIVE

This study aims to reveal how VDR regulates CYP11A1 expression and affects testosterone synthesis in murine Leydig cells.

METHODS

The levels of VDR, CYP11A1 were determined by quantitative real-time polymerase chain reaction (RT-qPCR) or western blot. Targeted relationship between VDR and Cyp11a1 was evaluated by dual-luciferase reporter assay. The levels of testosterone concentrations in cell culture media serum by enzyme-linked immunosorbent assay (ELISA).

RESULTS

Phylogenetic and motif analysis showed that the Cyp11a1 family had sequence loss, which may have special biological functions during evolution. The results of promoter prediction showed that vitamin D response element (VDRE) existed in the upstream promoter region of murine Cyp11a1. Dual-luciferase assay confirmed that VDR could bind candidate VDREs in upstream region of Cyp11a1, and enhance gene expression. Tissue distribution and localizatio analysis showed that Cyp11a1 was mainly expressed in testis, and dominantly existed in murine Leydig cells. Furthermore, over-expression VDR and CYP11A1 significantly increased testosterone synthesis in mice Leydig cells.

CONCLUSIONS

Active vitamin D3 (VD3) and Vdr interference treatment showed that VD3/VDR had a positive regulatory effect on Cyp11a1 expression and testosterone secretion. VDR promotes testosterone synthesis in male mice by up-regulating Cyp11a1 expression, which played an important role for male reproduction.

Transcriptomics for Clinical and Experimental Biology Research: Hang on a Seq

Sequencing the human genome empowers translational medicine, facilitating transcriptome-wide molecular diagnosis, pathway biology, and drug repositioning. Initially, microarrays are used to study the bulk transcriptome; but now short-read RNA sequencing (RNA-seq) predominates. Positioned as a superior technology, that makes the discovery of novel transcripts routine, most RNA-seq analyses are in fact modeled on the known transcriptome. Limitations of the RNA-seq methodology have emerged, while the design of, and the analysis strategies applied to, arrays have matured. An equitable comparison between these technologies is provided, highlighting advantages that modern arrays hold over RNA-seq. Array protocols more accurately quantify constitutively expressed protein coding genes across tissue replicates, and are more reliable for studying lower expressed genes. Arrays reveal long noncoding RNAs (lncRNA) are neither sparsely nor lower expressed than protein coding genes. Heterogeneous coverage of constitutively expressed genes observed with RNA-seq, undermines the validity and reproducibility of pathway analyses. The factors driving these observations, many of which are relevant to long-read or single-cell sequencing are discussed. As proposed herein, a reappreciation of bulk transcriptomic methods is required, including wider use of the modern high-density array data-to urgently revise existing anatomical RNA reference atlases and assist with more accurate study of lncRNAs.

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.

Current Scenario of Exogenously Induced RNAi for Lepidopteran Agricultural Pest Control: From dsRNA Design to Topical Application

Invasive insects cost the global economy around USD 70 billion per year. Moreover, increasing agricultural insect pests raise concerns about global food security constraining and infestation rising after climate changes. Current agricultural pest management largely relies on plant breeding-with or without transgenes-and chemical pesticides. Both approaches face serious technological obsolescence in the field due to plant resistance breakdown or development of insecticide resistance. The need for new modes of action (MoA) for managing crop health is growing each year, driven by market demands to reduce economic losses and by consumer demand for phytosanitary measures. The disabling of pest genes through sequence-specific expression silencing is a promising tool in the development of environmentally-friendly and safe biopesticides. The specificity conferred by long dsRNA-base solutions helps minimize effects on off-target genes in the insect pest genome and the target gene in non-target organisms (NTOs). In this review, we summarize the status of gene silencing by RNA interference (RNAi) for agricultural control. More specifically, we focus on the engineering, development and application of gene silencing to control Lepidoptera through non-transforming dsRNA technologies. Despite some delivery and stability drawbacks of topical applications, we reviewed works showing convincing proof-of-concept results that point to innovative solutions. Considerations about the regulation of the ongoing research on dsRNA-based pesticides to produce commercialized products for exogenous application are discussed. Academic and industry initiatives have revealed a worthy effort to control Lepidoptera pests with this new mode of action, which provides more sustainable and reliable technologies for field management. New data on the genomics of this taxon may contribute to a future customized target gene portfolio. As a case study, we illustrate how dsRNA and associated methodologies could be applied to control an important lepidopteran coffee pest.

Temporal progress of gene expression analysis with RNA-Seq data: A review on the relationship between computational methods

Analysis of differential gene expression from RNA-seq data has become a standard for several research areas. The steps for the computational analysis include many data types and file formats, and a wide variety of computational tools that can be applied alone or together as pipelines. This paper presents a review of the differential expression analysis pipeline, addressing its steps and the respective objectives, the principal methods available in each step, and their properties, therefore introducing an organized overview to this context. This review aims to address mainly the aspects involved in the differentially expressed gene (DEG) analysis from RNA sequencing data (RNA-seq), considering the computational methods. In addition, a timeline of the computational methods for DEG is shown and discussed, and the relationships existing between the most important computational tools are presented by an interaction network. A discussion on the challenges and gaps in DEG analysis is also highlighted in this review. This paper will serve as a tutorial for new entrants into the field and help established users update their analysis pipelines.

microRNA signatures associated with fetal growth restriction: a systematic review

Placental-origin microRNA (miRNA) profiles can be useful toward early diagnosis and management of fetal growth restriction (FGR) and associated complications. We conducted a systematic review to identify case-control studies that have examined miRNA signatures associated with human FGR. We systematically searched PubMed and ScienceDirect databases for relevant articles and manually searched reference lists of the relevant articles till May 18th, 2021. Of the 2133 studies identified, 21 were included. FGR-associated upregulation of miR-210 and miR-424 and downregulation of a placenta-specific miRNA cluster miRNA located on C19MC (miR-518b, miR-519d) and miR-221-3p was reported by >1 included studies. Analysis of the target genes of these miRNA as well as pathway analysis pointed to the involvement of angiogenesis and growth signaling pathways, such as the phosphatidylinositol 3-kinase- protein kinase B (PI3K-Akt) pathway. Only 3 out of the 21 included studies reported FGR-associated miRNAs in matched placental and maternal blood samples. We conclude that FGR-associated placental miRNAs could be utilized to inform clinical practice towards early diagnosis of FGR, provided enough evidence from studies on matched placental and maternal blood samples become available.Prospective Register of Systematic Reviews (PROSPERO) registration number: CRD42019136762.

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.

Analysis of Gut Microbiome Structure Based on GMPR+Spectrum

The gut microbiome is related to many major human diseases, and it is of great significance to study the structure of the gut microbiome under different conditions. Multivariate statistics or pattern recognition methods were often used to identify different structural patterns in gut microbiome data. However, these methods have some limitations. Minimal hepatic encephalopathy (MHE) datasets were taken as an example. Due to the physical lack or insufficient sampling of the gut microbiome in the sequencing process, the microbiome data contains many zeros. Therefore, the geometric mean of pairwise ratios (GMPR) was used to normalize gut microbiome data, then Spectrum was used to analyze the structure of the gut microbiome, and lastly, the structure of core microflora was compared with Network analysis. GMPR calculates the Intraclass correlation coefficient (ICC), whose reproducibility was significantly better than other normalization methods. In addition, running-time, Normalized Mutual Information (NMI), Davies-Boulding Index (DBI), and Calinski-Harabasz index (CH) of GMPR+Spectrum were far superior to other clustering algorithms such as M3C, iClusterPlus. GMPR+Spectrum can not only perform better but also effectively identify the structural differences of intestinal microbiota in different patients and excavate the unique critical bacteria such as Akkermansia, and Lactobacillus in MHE patients, which may provide a new reference for the study of the gut microbiome in disease.

A Robust Immuno-Prognostic Model of Non-Muscle-Invasive Bladder Cancer Indicates Dynamic Interaction in Tumor Immune Microenvironment Contributes to Cancer Progression

Non-muscle-invasive bladder cancer (NMIBC) accounts for more than 70% of urothelial cancer. More than half of NMIBC patients experience recurrence, progression, or metastasis, which essentially reduces life quality and survival time. Identifying the high-risk patients prone to progression remains the primary concern of risk management of NMIBC. In this study, we included 1370 NMIBC transcripts data from nine public datasets, identified nine tumor-infiltrating marker cells highly related to the survival of NMIBC, quantified the cells’ proportion by self-defined differentially expressed signature genes, and established a robust immuno-prognostic model dividing NMIBC patients into low-risk versus high-risk progression groups. Our model implies that the loss of crosstalk between tumor cells and adjacent normal epithelium, along with enriched cell proliferation signals, may facilitate tumor progression. Thus, evaluating tumor progression should consider various components in the tumor immune microenvironment instead of the single marker in a single dimension. Moreover, we also appeal to the necessity of using appropriate meta-analysis methods to integrate the evidence from multiple sources in the feature selection step from large-scale heterogeneous omics data such as our study.

Evaluation of critical data processing steps for reliable prediction of gene co-expression from large collections of RNA-seq data

Motivation

Gene co-expression analysis is an attractive tool for leveraging enormous amounts of public RNA-seq datasets for the prediction of gene functions and regulatory mechanisms. However, the optimal data processing steps for the accurate prediction of gene co-expression from such large datasets remain unclear. Especially the importance of batch effect correction is understudied.

Results

We processed RNA-seq data of 68 human and 76 mouse cell types and tissues using 50 different workflows into 7,200 genome-wide gene co-expression networks. We then conducted a systematic analysis of the factors that result in high-quality co-expression predictions, focusing on normalization, batch effect correction, and measure of correlation. We confirmed the key importance of high sample counts for high-quality predictions. However, choosing a suitable normalization approach and applying batch effect correction can further improve the quality of co-expression estimates, equivalent to a >80% and >40% increase in samples. In larger datasets, batch effect removal was equivalent to a more than doubling of the sample size. Finally, Pearson correlation appears more suitable than Spearman correlation, except for smaller datasets.

Conclusion

A key point for accurate prediction of gene co-expression is the collection of many samples. However, paying attention to data normalization, batch effects, and the measure of correlation can significantly improve the quality of co-expression estimates.

Systems biology under heat stress in Indian cattle

Transcriptome profiling of Vrindavani and Tharparkar cattle (n = 5 each) revealed that more numbers of genes were dysregulated in Vrindavani than in Tharparkar. A contrast in gene expression was observed with 18.9 % of upregulated genes in Vrindavani downregulated in Tharparkar and 17.8% upregulated genes in Tharparkar downregulated in Vrindavani. Functional annotation of genes differentially expressed in Tharparkar and Vrindavani revealed that the systems biology in Tharparkar is moving towards counteracting the effects due to heat stress. Unlike Vrindavani, Tharparkar is not only endowed with higher expression of the scavengers (UBE2G1, UBE2S, and UBE2H) of misfolded proteins but also with protectors (VCP, Serp1, and CALR) of naïve unfolded proteins. Further, higher expression of the antioxidants in Tharparkar enables it to cope up with higher levels of free radicals generated as a result of heat stress. In this study, we found relevant genes dysregulated in Tharparkar in the direction that can counter heat stress.

GeneCloudOmics: A Data Analytic Cloud Platform for High-Throughput Gene Expression Analysis

Gene expression profiling techniques, such as DNA microarray and RNA-Sequencing, have provided significant impact on our understanding of biological systems. They contribute to almost all aspects of biomedical research, including studying developmental biology, host-parasite relationships, disease progression and drug effects. However, the high-throughput data generations present challenges for many wet experimentalists to analyze and take full advantage of such rich and complex data. Here we present GeneCloudOmics, an easy-to-use web server for high-throughput gene expression analysis that extends the functionality of our previous ABioTrans with several new tools, including protein datasets analysis, and a web interface. GeneCloudOmics allows both microarray and RNA-Seq data analysis with a comprehensive range of data analytics tools in one package that no other current standalone software or web-based tool can do. In total, GeneCloudOmics provides the user access to 23 different data analytical and bioinformatics tasks including reads normalization, scatter plots, linear/non-linear correlations, PCA, clustering (hierarchical, k-means, t-SNE, SOM), differential expression analyses, pathway enrichments, evolutionary analyses, pathological analyses, and protein-protein interaction (PPI) identifications. Furthermore, GeneCloudOmics allows the direct import of gene expression data from the NCBI Gene Expression Omnibus database. The user can perform all tasks rapidly through an intuitive graphical user interface that overcomes the hassle of coding, installing tools/packages/libraries and dealing with operating systems compatibility and version issues, complications that make data analysis tasks challenging for biologists. Thus, GeneCloudOmics is a one-stop open-source tool for gene expression data analysis and visualization. It is freely available at http://combio-sifbi.org/GeneCloudOmics.

Constructing a Defined Starter for Multispecies Vinegar Fermentation via Evaluating the Vitality and Dominance of Functional Microbes in Autochthonous Starter

Mature vinegar culture has usually been used as a type of autochthonous starter for rapidly initiate initiating the next batch of acetic acid fermentation (AAF) and maintaining the batch-to-batch uniformity of AAF in the production of traditional cereal vinegar. However, the vitality and dominance of functional microbes in autochthonous starters remain unclear, which hinders further improvement of fermentation yield and production. Here, based on metagenomic (MG), metatranscriptomic (MT), and 16S rRNA gene sequencings, 11 bacterial operational taxonomic units (OTUs) with significant metabolic activity (MT/MG ratio >1) and dominance (relative abundance >1%) were targeted in the autochthonous vinegar starter, all of which were assigned to 4 species (Acetobacter pasteurianus, Lactobacillus acetotolerans, L. helveticus, Acetilactobacillus jinshanensis). Then, we evaluated the successions and interactions of these 11 bacterial OTUs at different AAF stages. Last, a defined starter was constructed with 4 core species isolated from the autochthonous starter (A. pasteurianus, L. acetotolerans, L. helveticus, Ac. jinshanensis). The defined starter culture could rapidly initiate the AAF in a sterile or unsterilized environment and similar dynamics of metabolites (ethanol, titratable acidity, acetic acid, lactic acid, and volatile compounds) and environmental indexes (temperature, pH) of fermentation were observed as compared with that of autochthonous starter (P > 0.05). This work provides a method to construct a defined microbiota from a complex system while preserving its metabolic function. IMPORTANCE Complex microorganisms are beneficial to the flavor formation in natural food fermentation, but they also pose challenges to the mass production of standardized products. It is attractive to construct a defined starter to rapidly initiate fermentation process and significantly improve fermentation yield. This study provides a comprehensive understanding of vital and dominant species in the autochthonous vinegar starter via multi-omics, and designs a defined microbial community for the efficient fermentation of cereal vinegar.

RNA-Seq reveals different responses to drought in Neotropical trees from savannas and seasonally dry forests

BACKGROUND

Water is one of the main limiting factors for plant growth and crop productivity. Plants constantly monitor water availability and can rapidly adjust their metabolism by altering gene expression. This leads to phenotypic plasticity, which aids rapid adaptation to climate changes. Here, we address phenotypic plasticity under drought stress by analyzing differentially expressed genes (DEG) in four phylogenetically related neotropical Bignoniaceae tree species: two from savanna, Handroanthus ochraceus and Tabebuia aurea, and two from seasonally dry tropical forests (SDTF), Handroanthus impetiginosus and Handroanthus serratifolius. To the best of our knowledge, this is the first report of an RNA-Seq study comparing tree species from seasonally dry tropical forest and savanna ecosystems.

RESULTS

Using a completely randomized block design with 4 species × 2 treatments (drought and wet) × 3 blocks (24 plants) and an RNA-seq approach, we detected a higher number of DEGs between treatments for the SDTF species H. serratifolius (3153 up-regulated and 2821 down-regulated under drought) and H. impetiginosus (332 and 207), than for the savanna species. H. ochraceus showed the lowest number of DEGs, with only five up and nine down-regulated genes, while T. aurea exhibited 242 up- and 96 down-regulated genes. The number of shared DEGs among species was not related to habitat of origin or phylogenetic relationship, since both T. aurea and H impetiginosus shared a similar number of DEGs with H. serratifolius. All four species shared a low number of enriched gene ontology (GO) terms and, in general, exhibited different mechanisms of response to water deficit. We also found 175 down-regulated and 255 up-regulated transcription factors from several families, indicating the importance of these master regulators in drought response.

CONCLUSION

Our findings show that phylogenetically related species may respond differently at gene expression level to drought stress. Savanna species seem to be less responsive to drought at the transcriptional level, likely due to morphological and anatomical adaptations to seasonal drought. The species with the largest geographic range and widest edaphic-climatic niche, H. serratifolius, was the most responsive, exhibiting the highest number of DEG and up- and down-regulated transcription factors (TF).

Viromes and surveys of RNA viruses in camel-derived ticks revealing transmission patterns of novel tick-borne viral pathogens in Kenya

Tick-borne viruses (TBVs) capable of transmitting between ticks and hosts have been increasingly recognized as a global public health concern. In this study, Hyalomma ticks and serum samples from camels were collected using recorded sampling correlations in eastern Kenya. Viromes of pooled ticks were profiled by metagenomic sequencing, revealing a diverse community of viruses related to at least 11 families. Five highly abundant viruses, including three novel viruses (Iftin tick virus, Mbalambala tick virus [MATV], and Bangali torovirus [BanToV]) and new strains of previously identified viruses (Bole tick virus 4 [BLTV4] and Liman tick virus [LMTV]), were characterized in terms of genome sequences, organizations, and phylogeny, and their molecular prevalence was investigated in individual ticks. Moreover, viremia and antibody responses to these viruses have been investigated in camels. MATV, BLTV4, LMTV, and BanToV were identified as viral pathogens that can potentially cause zoonotic diseases. The transmission patterns of these viruses were summarized, suggesting three different types according to the sampling relationships between viral RNA-positive ticks and camels positive for viral RNA and/or antibodies. They also revealed the frequent transmission of BanToV and limited but effective transmission of other viruses between ticks and camels. Furthermore, follow-up surveys on TBVs from tick, animal, and human samples with definite sampling relationships are suggested. The findings revealed substantial threats from the emerging TBVs and may guide the prevention and control of TBV-related zoonotic diseases in Kenya and in other African countries.

Genome-Wide Identification and Transcriptional Expression Profiles of Transcription Factor WRKY in Common Walnut (Juglans regia L.)

The transcription factor WRKY is widely distributed in the plant kingdom, playing a significant role in plant growth, development and response to stresses. Walnut is an economically important temperate tree species valued for both its edible nuts and high-quality wood, and its response to various stresses is an important factor that determines the quality of its fruit. However, in walnut trees themselves, information about the WRKY gene family remains scarce. In this paper, we perform a comprehensive study of the WRKY gene family in walnut. In total, we identified 103 WRKY genes in the common walnut that are clustered into 4 groups and distributed on 14 chromosomes. The conserved domains all contained a WRKY domain, and motif 2 was observed in most WRKYs, suggesting a high degree of conservation and similar functions within each subfamily. However, gene structure was significantly differentiated between different subfamilies. Synteny analysis indicates that there were 56 gene pairs in J. regia and A. thaliana, 76 in J. regia and J. mandshurica, 75 in J. regia and J. microcarpa, 76 in J. regia and P. trichocarpa, and 33 in J. regia and Q. robur, indicating that the WRKY gene family may come from a common ancestor. GO and KEGG enrichment analysis showed that the WRKY gene family was involved in resistance traits and the plant-pathogen interaction pathway. In anthracnose-resistant F26 fruits (AR) and anthracnose-susceptible F423 fruits (AS), transcriptome and qPCR analysis results showed that JrWRKY83, JrWRKY73 and JrWRKY74 were expressed significantly more highly in resistant cultivars, indicating that these three genes may be important contributors to stress resistance in walnut trees. Furthermore, we investigate how these three genes potentially target miRNAs and interact with proteins. JrWRKY73 was target by the miR156 family, including 12 miRNAs; this miRNA family targets WRKY genes to enhance plant defense. JrWRKY73 also interacted with the resistance gene AtMPK6, showing that it may play a crucial role in walnut defense.

Assessment of reference genes at six different developmental stages of Schistosoma mansoni for quantitative RT-PCR

Reverse-transcription quantitative real-time polymerase chain reaction (RT-qPCR) is the most used, fast, and reproducible method to confirm large-scale gene expression data. The use of stable reference genes for the normalization of RT-qPCR assays is recognized worldwide. No systematic study for selecting appropriate reference genes for usage in RT-qPCR experiments comparing gene expression levels at different Schistosoma mansoni life-cycle stages has been performed. Most studies rely on genes commonly used in other organisms, such as actin, tubulin, and GAPDH. Therefore, the present study focused on identifying reference genes suitable for RT-qPCR assays across six S. mansoni developmental stages. The expression levels of 25 novel candidates that we selected based on the analysis of public RNA-Seq datasets, along with eight commonly used reference genes, were systematically tested by RT-qPCR across six developmental stages of S. mansoni (eggs, miracidia, cercariae, schistosomula, adult males and adult females). The stability of genes was evaluated with geNorm, NormFinder and RefFinder algorithms. The least stable candidate reference genes tested were actin, tubulin and GAPDH. The two most stable reference genes suitable for RT-qPCR normalization were Smp_101310 (Histone H4 transcription factor) and Smp_196510 (Ubiquitin recognition factor in ER-associated degradation protein 1). Performance of these two genes as normalizers was successfully evaluated with females maintained unpaired or paired to males in culture for 8 days, or with worm pairs exposed for 16 days to double-stranded RNAs to silence a protein-coding gene. This study provides reliable reference genes for RT-qPCR analysis using samples from six different S. mansoni life-cycle stages.

Distinct peripheral blood monocyte and neutrophil transcriptional programs following intracerebral hemorrhage and different etiologies of ischemic stroke

Understanding cell-specific transcriptome responses following intracerebral hemorrhage (ICH) and ischemic stroke (IS) will improve knowledge of the immune response to brain injury. Transcriptomic profiles of 141 samples from 48 subjects with ICH, different IS etiologies, and vascular risk factor controls were characterized using RNA-seq in isolated neutrophils, monocytes and whole blood. In both IS and ICH, monocyte genes were down-regulated, whereas neutrophil gene expression changes were generally up-regulated. The monocyte down-regulated response to ICH included innate, adaptive immune, dendritic, NK cell and atherosclerosis signaling. Neutrophil responses to ICH included tRNA charging, mitochondrial dysfunction, and ER stress pathways. Common monocyte and neutrophil responses to ICH included interferon signaling, neuroinflammation, death receptor signaling, and NFAT pathways. Suppressed monocyte responses to IS included interferon and dendritic cell maturation signaling, phagosome formation, and IL-15 signaling. Activated neutrophil responses to IS included oxidative phosphorylation, mTOR, BMP, growth factor signaling, and calpain proteases-mediated blood-brain barrier (BBB) dysfunction. Common monocyte and neutrophil responses to IS included JAK1, JAK3, STAT3, and thrombopoietin signaling. Cell-type and cause-specific approaches will assist the search for future IS and ICH biomarkers and treatments.

Constrained Standardization of Count Data from Massive Parallel Sequencing

In high-throughput omics disciplines like transcriptomics, researchers face a need to assess the quality of an experiment prior to an in-depth statistical analysis. To efficiently analyze such voluminous collections of data, researchers need triage methods that are both quick and easy to use. Such a normalization method for relative quantitation, CONSTANd, was recently introduced for isobarically-labeled mass spectra in proteomics. It transforms the data matrix of abundances through an iterative, convergent process enforcing three constraints: (I) identical column sums; (II) each row sum is fixed (across matrices) and (III) identical to all other row sums. In this study, we investigate whether CONSTANd is suitable for count data from massively parallel sequencing, by qualitatively comparing its results to those of DESeq2. Further, we propose an adjustment of the method so that it may be applied to identically balanced but differently sized experiments for joint analysis. We find that CONSTANd can process large data sets at well over 1 million count records per second whilst mitigating unwanted systematic bias and thus quickly uncovering the underlying biological structure when combined with a PCA plot or hierarchical clustering. Moreover, it allows joint analysis of data sets obtained from different batches, with different protocols and from different labs but without exploiting information from the experimental setup other than the delineation of samples into identically processed sets (IPSs). CONSTANd's simplicity and applicability to proteomics as well as transcriptomics data make it an interesting candidate for integration in multi-omics workflows.

Robustness of differential gene expression analysis of RNA-seq

RNA-sequencing (RNA-seq) is a relatively new technology that lacks standardisation. RNA-seq can be used for Differential Gene Expression (DGE) analysis, however, no consensus exists as to which methodology ensures robust and reproducible results. Indeed, it is broadly acknowledged that DGE methods provide disparate results. Despite obstacles, RNA-seq assays are in advanced development for clinical use but further optimisation will be needed. Herein, five DGE models (DESeq2, voom + limma, edgeR, EBSeq, NOISeq) for gene-level detection were investigated for robustness to sequencing alterations using a controlled analysis of fixed count matrices. Two breast cancer datasets were analysed with full and reduced sample sizes. DGE model robustness was compared between filtering regimes and for different expression levels (high, low) using unbiased metrics. Test sensitivity estimated as relative False Discovery Rate (FDR), concordance between model outputs and comparisons of a ’population’ of slopes of relative FDRs across different library sizes, generated using linear regressions, were examined. Patterns of relative DGE model robustness proved dataset-agnostic and reliable for drawing conclusions when sample sizes were sufficiently large. Overall, the non-parametric method NOISeq was the most robust followed by edgeR, voom, EBSeq and DESeq2. Our rigorous appraisal provides information for method selection for molecular diagnostics. Metrics may prove useful towards improving the standardisation of RNA-seq for precision medicine.

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.

Identification of transcriptional subtypes in lung adenocarcinoma and squamous cell carcinoma through integrative analysis of microarray and RNA sequencing data

Classification of tumors into subtypes can inform personalized approaches to treatment including the choice of targeted therapies. The two most common lung cancer histological subtypes, lung adenocarcinoma and lung squamous cell carcinoma, have been previously divided into transcriptional subtypes using microarray data, and corresponding signatures were subsequently used to classify RNA-seq data. Cross-platform unsupervised classification facilitates the identification of robust transcriptional subtypes by combining vast amounts of publicly available microarray and RNA-seq data. However, cross-platform classification is challenging because of intrinsic differences in data generated using the two gene expression profiling technologies. In this report, we show that robust gene expression subtypes can be identified in integrated data representing over 3500 normal and tumor lung samples profiled using two widely used platforms, Affymetrix HG-U133 Plus 2.0 Array and Illumina HiSeq RNA sequencing. We tested and analyzed consensus clustering for 384 combinations of data processing methods. The agreement between subtypes identified in single-platform and cross-platform normalized data was then evaluated using a variety of statistics. Results show that unsupervised learning can be achieved with combined microarray and RNA-seq data using selected preprocessing, cross-platform normalization, and unsupervised feature selection methods. Our analysis confirmed three lung adenocarcinoma transcriptional subtypes, but only two consistent subtypes in squamous cell carcinoma, as opposed to four subtypes previously identified. Further analysis showed that tumor subtypes were associated with distinct patterns of genomic alterations in genes coding for therapeutic targets. Importantly, by integrating quantitative proteomics data, we were able to identify tumor subtype biomarkers that effectively classify samples on the basis of both gene and protein expression. This study provides the basis for further integrative data analysis across gene and protein expression profiling platforms.

Recommendations for the introduction of metagenomic next-generation sequencing in clinical virology, part II: bioinformatic analysis and reporting

Metagenomic next-generation sequencing (mNGS) is an untargeted technique for determination of microbial DNA/RNA sequences in a variety of sample types from patients with infectious syndromes. mNGS is still in its early stages of broader translation into clinical applications. To further support the development, implementation, optimization and standardization of mNGS procedures for virus diagnostics, the European Society for Clinical Virology (ESCV) Network on Next-Generation Sequencing (ENNGS) has been established. The aim of ENNGS is to bring together professionals involved in mNGS for viral diagnostics to share methodologies and experiences, and to develop application guidelines. Following the ENNGS publication Recommendations for the introduction of mNGS in clinical virology, part I: wet lab procedure in this journal, the current manuscript aims to provide practical recommendations for the bioinformatic analysis of mNGS data and reporting of results to clinicians.

Gene and Transposable Element Expression Evolution Following Recent and Past Polyploidy Events in Spartina (Poaceae)

Gene expression dynamics is a key component of polyploid evolution, varying in nature, intensity, and temporal scales, most particularly in allopolyploids, where two or more sub-genomes from differentiated parental species and different repeat contents are merged. Here, we investigated transcriptome evolution at different evolutionary time scales among tetraploid, hexaploid, and neododecaploid Spartina species (Poaceae, Chloridoideae) that successively diverged in the last 6-10 my, at the origin of differential phenotypic and ecological traits. Of particular interest are the recent (19th century) hybridizations between the two hexaploids Spartina alterniflora (2n = 6x = 62) and S. maritima (2n = 6x = 60) that resulted in two sterile F1 hybrids: Spartina × townsendii (2n = 6x = 62) in England and Spartina × neyrautii (2n = 6x = 62) in France. Whole genome duplication of S. × townsendii gave rise to the invasive neo-allododecaploid species Spartina anglica (2n = 12x = 124). New transcriptome assemblies and annotations for tetraploids and the enrichment of previously published reference transcriptomes for hexaploids and the allododecaploid allowed identifying 42,423 clusters of orthologs and distinguishing 21 transcribed transposable element (TE) lineages across the seven investigated Spartina species. In 4x and 6x mesopolyploids, gene and TE expression changes were consistent with phylogenetic relationships and divergence, revealing weak expression differences in the tetraploid sister species Spartina bakeri and Spartina versicolor (<2 my divergence time) compared to marked transcriptome divergence between the hexaploids S. alterniflora and S. maritima that diverged 2-4 mya. Differentially expressed genes were involved in glycolysis, post-transcriptional protein modifications, epidermis development, biosynthesis of carotenoids. Most detected TE lineages (except SINE elements) were found more expressed in hexaploids than in tetraploids, in line with their abundance in the corresponding genomes. Comparatively, an astonishing (52%) expression repatterning and deviation from parental additivity were observed following recent reticulate evolution (involving the F1 hybrids and the neo-allododecaploid S. anglica), with various patterns of biased homoeologous gene expression, including genes involved in epigenetic regulation. Downregulation of TEs was observed in both hybrids and accentuated in the neo-allopolyploid. Our results reinforce the view that allopolyploidy represents springboards to new regulatory patterns, offering to worldwide invasive species, such as S. anglica, the opportunity to colonize stressful and fluctuating environments on saltmarshes.

A comprehensive analysis of tumor microenvironment-related genes in colon cancer

BACKGROUND

The development and progression of colon cancer are significantly affected by the tumor microenvironment, which has attracted much attention. The goal of our study was primarily to find out all possible tumor microenvironment-related genes in colon cancer.

METHOD

This study quantified the immune and stromal landscape using the ESTIMATION algorithm using the gene expression matrix obtained from the UCSC Xena database. Dysregulated genes were harvested using the limma R package, and relevant pathways and biofunctions were identified using enrichment analysis. A least absolute shrinkage and selection operator (LASSO) regression was used to select the pivotal genes from the DEGs. Then, survival analysis was performed to determine the hub genes and a prognostic model was constructed by these hub genes with (or) TNM stage. Besides, associations between hub gene expressions and immune cell infiltration were assessed.

RESULTS

A total of 725 DEGs were identified. Most of the results of the enrichment analysis were immune-related items. 13 genes were selected as the hub genes and a moderate-to-strong positive correlation between most hub genes and several immune cells were observed. Besides, the prognostic value of the hub genes were comparable to TNM staging.

CONCLUSIONS

Our study provides a better understanding of how interactions between the 13 immune-prognostic hub genes and immune cells in the tumor microenvironment affect biological processes in colon cancer. These genes exhibit an equivalent ability to TNM staging in prognosis prediction. They are particularly expected to become novel prognostic biomarkers and targets of immunotherapies for colon cancer.

MODEL-BASED FEATURE SELECTION AND CLUSTERING OF RNA-SEQ DATA FOR UNSUPERVISED SUBTYPE DISCOVERY

Clustering is a form of unsupervised learning that aims to uncover latent groups within data based on similarity across a set of features. A common application of this in biomedical research is in delineating novel cancer subtypes from patient gene expression data, given a set of informative genes. However, it is typically unknown a priori what genes may be informative in discriminating between clusters, and what the optimal number of clusters are. Few methods exist for performing unsupervised clustering of RNA-seq samples, and none currently adjust for between-sample global normalization factors, select cluster-discriminatory genes, or account for potential confounding variables during clustering. To address these issues, we propose the Feature Selection and Clustering of RNA-seq (FSCseq): a model-based clustering algorithm that utilizes a finite mixture of regression (FMR) model and the quadratic penalty method with a Smoothly-Clipped Absolute Deviation (SCAD) penalty. The maximization is done by a penalized Classification EM algorithm, allowing us to include normalization factors and confounders in our modeling framework. Given the fitted model, our framework allows for subtype prediction in new patients via posterior probabilities of cluster membership, even in the presence of batch effects. Based on simulations and real data analysis, we show the advantages of our method relative to competing approaches.

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.

Formula Diet Alters the Ileal Metagenome and Transcriptome at Weaning and during the Postweaning Period in a Porcine Model

Exclusive breastfeeding impacts the intestinal microbiome and is associated with a better immune function than is seen with milk formula (MF) feeding in infants and yet with mechanisms poorly defined. The porcine model was used to evaluate the impact of MF on ileum microbial communities and gene expression relative to human milk (HM)-fed piglets. Fifty-two Dutch Landrace male piglets were fed an isocaloric diet of either HM (n = 26) or MF (n = 26) from day 2 through day 21 of age and weaned to a solid diet until day 51. Eleven piglets from each group were euthanized at day 21, while the remaining piglets (HM, n = 15; MF, n = 15) were euthanized at day 51 to collect ileal epithelium (EP) scrapings and ileal (IL) tissues. The epithelial mucosa was subjected to shotgun metagenome sequencing, and EP and IL tissues were used for transcriptome analysis. On day 21, transcriptome data revealed that the levels of pathways involved in inflammation and apoptosis were significantly higher in MF piglets than in HM piglets, whereas the levels of tight junctions and pathogen detection systems were lower in MF piglets than in HM piglets. The MF impacts on the small intestine were maintained over the postweaning period (day 51) as indicated by higher levels of Dialister invisus bacteria and higher levels of expression of genes associated with inflammation and apoptosis pathways relative to HM group. The current study demonstrated that MF might impact local intestinal inflammation, apoptosis, and tight junctions and might suppress pathogen recognition in the small intestine compared with HM.IMPORTANCE Exclusive human milk (HM) breastfeeding for the first 6 months of age in infants is recommended to improve health outcomes during early life and beyond. When women are unable to provide sufficient HM, milk formula (MF) is often recommended as a complementary or alternative source of nutrition. Previous studies in piglets demonstrated that MF alters the gut microbiome and induces inflammatory cytokine production. The links between MF feeding, gut microbiome, and inflammation status are unclear due to challenges associated with the collection of intestinal samples from human infants. The current report provides the first insight into MF-microbiome-inflammation connections in the small intestine compared with HM feeding using a porcine model. The present results showed that, compared with HM, MF might impact immune function through the induction of ileal inflammation, apoptosis, and tight junction disruptions and likely compromised immune defense against pathogen detection in the small intestine relative to piglets that were fed HM.

Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols

In recent years, RNA-sequencing (RNA-seq) has emerged as a powerful technology for transcriptome profiling. For a given gene, the number of mapped reads is not only dependent on its expression level and gene length, but also the sequencing depth. To normalize these dependencies, RPKM (reads per kilobase of transcript per million reads mapped) and TPM (transcripts per million) are used to measure gene or transcript expression levels. A common misconception is that RPKM and TPM values are already normalized, and thus should be comparable across samples or RNA-seq projects. However, RPKM and TPM represent the relative abundance of a transcript among a population of sequenced transcripts, and therefore depend on the composition of the RNA population in a sample. Quite often, it is reasonable to assume that total RNA concentration and distributions are very close across compared samples. Nevertheless, the sequenced RNA repertoires may differ significantly under different experimental conditions and/or across sequencing protocols; thus, the proportion of gene expression is not directly comparable in such cases. In this review, we illustrate typical scenarios in which RPKM and TPM are misused, unintentionally, and hope to raise scientists' awareness of this issue when comparing them across samples or different sequencing protocols.

PenDA, a rank-based method for personalized differential analysis: Application to lung cancer

The hopes of precision medicine rely on our capacity to measure various high-throughput genomic information of a patient and to integrate them for personalized diagnosis and adapted treatment. Reaching these ambitious objectives will require the development of efficient tools for the detection of molecular defects at the individual level. Here, we propose a novel method, PenDA, to perform Personalized Differential Analysis at the scale of a single sample. PenDA is based on the local ordering of gene expressions within individual cases and infers the deregulation status of genes in a sample of interest compared to a reference dataset. Based on realistic simulations of RNA-seq data of tumors, we showed that PenDA outcompetes existing approaches with very high specificity and sensitivity and is robust to normalization effects. Applying the method to lung cancer cohorts, we observed that deregulated genes in tumors exhibit a cancer-type-specific commitment towards up- or down-regulation. Based on the individual information of deregulation given by PenDA, we were able to define two new molecular histologies for lung adenocarcinoma cancers strongly correlated to survival. In particular, we identified 37 biomarkers whose up-regulation lead to bad prognosis and that we validated on two independent cohorts. PenDA provides a robust, generic tool to extract personalized deregulation patterns that can then be used for the discovery of therapeutic targets and for personalized diagnosis. An open-access, user-friendly R package is available at https://github.com/bcm-uga/penda.

The hopes of precision medicine rely on our capacity to measure individual molecular information for personalized diagnosis and treatment. These challenging perspectives will be only possible with the development of efficient methodological tools to identify patient-specific molecular defects from the many precise molecular information that one can access at the single-individual, single tissue or even single-cell levels. Such methods will provide a better understanding of disease-specific biological mechanisms and will promote the development of personalized therapeutic strategies. Here we describe a novel method, named PenDA, to perform differential analysis of gene expression at the individual level. Based on a realistic benchmark of simulated tumors, we demonstrated that PenDA reaches very high efficiency in detecting sample-specific deregulated genes. We then applied the method to two large cohorts associated with lung cancer. A detailed statistical analysis of the results allowed to isolate genes with specific deregulation patterns, like genes that are up-regulated in all tumors or genes that are expressed but never deregulated in any tumors. Given their specificities, these genes are likely to be of interest in therapeutic research. In particular, we were able to identified 37 new biomarkers associated to bad prognosis that we validated on two independent cohorts.

A Zipf-plot based normalization method for high-throughput RNA-seq data

Normalization is crucial in RNA-seq data analyses. Due to the existence of excessive zeros and a large number of small measures, it is challenging to find reliable linear rescaling normalization parameters. We propose a Zipf plot based normalization method (ZN) assuming that all gene profiles have similar upper tail behaviors in their expression distributions. The new normalization method uses global information of all genes in the same profile without gene-level expression alteration. It doesn’t require the majority of genes to be not differentially expressed (DE), and can be applied to data where the majority of genes are weakly or not expressed. Two normalization schemes are implemented with ZN: a linear rescaling scheme and a non-linear transformation scheme. The linear rescaling scheme can be applied alone or together with the non-linear normalization scheme. The performance of ZN is benchmarked against five popular linear normalization methods for RNA-seq data. Results show that the linear rescaling normalization scheme by itself works well and is robust. The non-linear normalization scheme can further improve the normalization outcomes and is optional if the Zipf plots show parallel patterns.

Homoeologous gene expression and co-expression network analyses and evolutionary inference in allopolyploids

Polyploidy is a widespread phenomenon throughout eukaryotes. Due to the coexistence of duplicated genomes, polyploids offer unique challenges for estimating gene expression levels, which is essential for understanding the massive and various forms of transcriptomic responses accompanying polyploidy. Although previous studies have explored the bioinformatics of polyploid transcriptomic profiling, the causes and consequences of inaccurate quantification of transcripts from duplicated gene copies have not been addressed. Using transcriptomic data from the cotton genus (Gossypium) as an example, we present an analytical workflow to evaluate a variety of bioinformatic method choices at different stages of RNA-seq analysis, from homoeolog expression quantification to downstream analysis used to infer key phenomena of polyploid expression evolution. In general, EAGLE-RC and GSNAP-PolyCat outperform other quantification pipelines tested, and their derived expression dataset best represents the expected homoeolog expression and co-expression divergence. The performance of co-expression network analysis was less affected by homoeolog quantification than by network construction methods, where weighted networks outperformed binary networks. By examining the extent and consequences of homoeolog read ambiguity, we illuminate the potential artifacts that may affect our understanding of duplicate gene expression, including an overestimation of homoeolog co-regulation and the incorrect inference of subgenome asymmetry in network topology. Taken together, our work points to a set of reasonable practices that we hope are broadly applicable to the evolutionary exploration of polyploids.

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.

Genome-Wide Analysis of Cyclophilin Proteins in 21 Oomycetes

Cyclophilins (CYPs), a highly-conserved family of proteins, belong to a subgroup of immunophilins. Ubiquitous in eukaryotes and prokaryotes, CYPs have peptidyl-prolyl cis–trans isomerase (PPIase) activity and have been implicated as virulence factors in plant pathogenesis by oomycetes. We identified 16 CYP orthogroups from 21 diverse oomycetes. Each species was found to encode 15 to 35 CYP genes. Three of these orthogroups contained proteins with signal peptides at the N-terminal end, suggesting a role in secretion. Multidomain analysis revealed five conserved motifs of the CYP domain of oomycetes shared with other eukaryotic PPIases. Expression analysis of CYP proteins in different asexual life stages of the hemibiotrophic Phytophthora infestans and the biotrophic Plasmopara halstedii demonstrated distinct expression profiles between life stages. In addition to providing detailed comparative information on the CYPs in multiple oomycetes, this study identified candidate CYP effectors that could be the foundation for future studies of virulence.

A protocol to evaluate RNA sequencing normalization methods

BACKGROUND

RNA sequencing technologies have allowed researchers to gain a better understanding of how the transcriptome affects disease. However, sequencing technologies often unintentionally introduce experimental error into RNA sequencing data. To counteract this, normalization methods are standardly applied with the intent of reducing the non-biologically derived variability inherent in transcriptomic measurements. However, the comparative efficacy of the various normalization techniques has not been tested in a standardized manner. Here we propose tests that evaluate numerous normalization techniques and applied them to a large-scale standard data set. These tests comprise a protocol that allows researchers to measure the amount of non-biological variability which is present in any data set after normalization has been performed, a crucial step to assessing the biological validity of data following normalization.

RESULTS

In this study we present two tests to assess the validity of normalization methods applied to a large-scale data set collected for systematic evaluation purposes. We tested various RNASeq normalization procedures and concluded that transcripts per million (TPM) was the best performing normalization method based on its preservation of biological signal as compared to the other methods tested.

CONCLUSION

Normalization is of vital importance to accurately interpret the results of genomic and transcriptomic experiments. More work, however, needs to be performed to optimize normalization methods for RNASeq data. The present effort helps pave the way for more systematic evaluations of normalization methods across different platforms. With our proposed schema researchers can evaluate their own or future normalization methods to further improve the field of RNASeq normalization.

An expanded landscape of human long noncoding RNA

Long noncoding RNAs (lncRNAs) are emerging as key regulators of multiple essential biological processes involved in physiology and pathology. By analyzing the largest compendium of 14,166 samples from normal and tumor tissues, we significantly expand the landscape of human long noncoding RNA with a high-quality atlas: RefLnc (Reference catalog of LncRNA). Powered by comprehensive annotation across multiple sources, RefLnc helps to pinpoint 275 novel intergenic lncRNAs correlated with sex, age or race as well as 369 novel ones associated with patient survival, clinical stage, tumor metastasis or recurrence. Integrated in a user-friendly online portal, the expanded catalog of human lncRNAs provides a valuable resource for investigating lncRNA function in both human biology and cancer development.

Deciphering links between bacterial interactions and spatial organization in multispecies biofilms

Environmental microbes frequently live in multispecies biofilms where mutualistic relationships and co-evolution may occur, defining spatial organization for member species and overall community functions. In this context, intrinsic properties emerging from microbial interactions, such as efficient organization optimizing growth and activities in multispecies biofilms, may become the object of fitness selection. However, little is known on the nature of underlying interspecies interactions during establishment of a predictable spatial organization within multispecies biofilms. We present a comparative metatranscriptomic analysis of bacterial strains residing in triple-species and four-species biofilms, aiming at deciphering molecular mechanisms underpinning bacterial interactions responsible of the remarkably enhanced biomass production and associated typical spatial organization they display. Metatranscriptomic profiles concurred with changes in micro-site occupation in response to the addition/removal of a single species, being driven by both cooperation, competition, and facilitation processes. We conclude that the enhanced biomass production of the four-species biofilm is an intrinsic community property emerging from finely tuned space optimization achieved through concerted antagonistic and mutualistic interactions, where each species occupies a defined micro-site favoring its own growth. Our results further illustrate how molecular mechanisms can be better interpreted when supported by visual imaging of actual microscopic spatial organization, and we propose phenotypic adaptation selected by social interactions as molecular mechanisms stabilizing microbial communities.

Evaluation of commonly used analysis strategies for epigenome- and transcriptome-wide association studies through replication of large-scale population studies

BACKGROUND

A large number of analysis strategies are available for DNA methylation (DNAm) array and RNA-seq datasets, but it is unclear which strategies are best to use. We compare commonly used strategies and report how they influence results in large cohort studies.

RESULTS

We tested the associations of DNAm and RNA expression with age, BMI, and smoking in four different cohorts (n = ~ 2900). By comparing strategies against the base model on the number and percentage of replicated CpGs for DNAm analyses or genes for RNA-seq analyses in a leave-one-out cohort replication approach, we find the choice of the normalization method and statistical test does not strongly influence the results for DNAm array data. However, adjusting for cell counts or hidden confounders substantially decreases the number of replicated CpGs for age and increases the number of replicated CpGs for BMI and smoking. For RNA-seq data, the choice of the normalization method, gene expression inclusion threshold, and statistical test does not strongly influence the results. Including five principal components or excluding correction of technical covariates or cell counts decreases the number of replicated genes.

CONCLUSIONS

Results were not influenced by the normalization method or statistical test. However, the correction method for cell counts, technical covariates, principal components, and/or hidden confounders does influence the results.

Recurrent functional misinterpretation of RNA-seq data caused by sample-specific gene length bias

Data normalization is a critical step in RNA sequencing (RNA-seq) analysis, aiming to remove systematic effects from the data to ensure that technical biases have minimal impact on the results. Analyzing numerous RNA-seq datasets, we detected a prevalent sample-specific length effect that leads to a strong association between gene length and fold-change estimates between samples. This stochastic sample-specific effect is not corrected by common normalization methods, including reads per kilobase of transcript length per million reads (RPKM), Trimmed Mean of M values (TMM), relative log expression (RLE), and quantile and upper-quartile normalization. Importantly, we demonstrate that this bias causes recurrent false positive calls by gene-set enrichment analysis (GSEA) methods, thereby leading to frequent functional misinterpretation of the data. Gene sets characterized by markedly short genes (e.g., ribosomal protein genes) or long genes (e.g., extracellular matrix genes) are particularly prone to such false calls. This sample-specific length bias is effectively removed by the conditional quantile normalization (cqn) and EDASeq methods, which allow the integration of gene length as a sample-specific covariate. Consequently, using these normalization methods led to substantial reduction in GSEA false results while retaining true ones. In addition, we found that application of gene-set tests that take into account gene–gene correlations attenuates false positive rates caused by the length bias, but statistical power is reduced as well. Our results advocate the inspection and correction of sample-specific length biases as default steps in RNA-seq analysis pipelines and reiterate the need to account for intergene correlations when performing gene-set enrichment tests to lessen false interpretation of transcriptomic data.

Analysis of numerous RNA-seq datasets reveals a recurrent sample-specific length bias that causes frequent false positive calls by gene-set enrichment analyses, leading to functional misinterpretation of the data. Its removal requires methods that allow the integration of gene length as sample-specific covariate.

GENAVi: a shiny web application for gene expression normalization, analysis and visualization

BACKGROUND

The development of next generation sequencing (NGS) methods led to a rapid rise in the generation of large genomic datasets, but the development of user-friendly tools to analyze and visualize these datasets has not developed at the same pace. This presents a two-fold challenge to biologists; the expertise to select an appropriate data analysis pipeline, and the need for bioinformatics or programming skills to apply this pipeline. The development of graphical user interface (GUI) applications hosted on web-based servers such as Shiny can make complex workflows accessible across operating systems and internet browsers to those without programming knowledge.

RESULTS

We have developed GENAVi (Gene Expression Normalization Analysis and Visualization) to provide a user-friendly interface for normalization and differential expression analysis (DEA) of human or mouse feature count level RNA-Seq data. GENAVi is a GUI based tool that combines Bioconductor packages in a format for scientists without bioinformatics expertise. We provide a panel of 20 cell lines commonly used for the study of breast and ovarian cancer within GENAVi as a foundation for users to bring their own data to the application. Users can visualize expression across samples, cluster samples based on gene expression or correlation, calculate and plot the results of principal components analysis, perform DEA and gene set enrichment and produce plots for each of these analyses. To allow scalability for large datasets we have provided local install via three methods. We improve on available tools by offering a range of normalization methods and a simple to use interface that provides clear and complete session reporting and for reproducible analysis.

CONCLUSION

The development of tools using a GUI makes them practical and accessible to scientists without bioinformatics expertise, or access to a data analyst with relevant skills. While several GUI based tools are currently available for RNA-Seq analysis we improve on these existing tools. This user-friendly application provides a convenient platform for the normalization, analysis and visualization of gene expression data for scientists without bioinformatics expertise.

Improved cellulase production in recombinant Saccharomyces cerevisiae by disrupting the cell wall protein-encoding gene CWP2

Budding yeast Saccharomyces cerevisiae has been widely used for heterologous protein production. However, low protein production titer and secretion levels continue to challenge its practical applications. The yeast cell wall plays important roles in yeast cell growth and environmental responses. Nevertheless, the effects of yeast cell wall proteins on heterologous protein production and secretion remain unclear. CWP2 encodes a mannoprotein that is the major component of the yeast cell wall. So far, studies on its function have been very limited. Here we show that CWP2 disruption improved extracellular cellobiohydrolase activity by 85.9%. A calcofluor white hypersensitivity assay revealed increased sensitivity of the mutant compared to the parental strain, indicating impaired cell wall integrity. However, no changes were observed in normal cell growth or growth stressed by tunicamycin and dithiothreitol, suggesting that the unfolded protein response pathway was not affected by the gene disruption. Comparative transcriptome analysis revealed changes in multiple genes involved in cell wall structure, biosynthesis, and cell wall integrity induced by CWP2 disruption, suggesting a pivotal role of Cwp2p in yeast cell wall organization. Notably, CWP2 disruption also led to elevated transcription of a large number of genes involved in ribosome biogenesis, which indicated that CWP2 is not only in yeast cell wall biosynthesis, but also in protein translation. This work reveals novel insights into the functions of CWP2 and also presents a new strategy to increase heterologous protein production in yeast strains by manipulating cell wall-related proteins.