Stability of methods for differential expression analysis of RNA-seq data

Flowering-associated gene expression and metabolic characteristics in adzuki bean (Vigna angularis L.) with different short-day induction periods

Background

The adzuki bean is a typical short-day plant and an important grain crop that is widely used due to its high nutritional and medicinal value. The adzuki bean flowering time is affected by multiple environmental factors, particularly the photoperiod. Adjusting the day length can induce flower synchronization in adzuki bean and accelerate the breeding process. In this study, we used RNA sequencing analysis to determine the effects of different day lengths on gene expression and metabolic characteristics related to adzuki bean flowering time.

Methods

'Tangshan hong xiao dou' was used as the experimental material in this study and field experiments were conducted in 2022 using a randomized block design with three treatments: short-day induction periods of 5 d (SD-5d), 10 d (SD-10d), and 15 d (SD-15d).

Results

A total of 5,939 differentially expressed genes (DEGs) were identified, of which 38.09% were up-regulated and 23.81% were down-regulated. Gene ontology enrichment analysis was performed on the target genes to identify common functions related to photosystems I and II. Kyoto Encyclopedia of Genes and Genomes enrichment analysis identified two pathways involved in the antenna protein and circadian rhythm. Furthermore, florescence was promoted by down-regulating genes in the circadian rhythm pathway through the blue light metabolic pathway; whereas, antenna proteins promoted flowering by enhancing the reception of light signals and accelerating electron transport. In these two metabolic pathways, the number of DEGs was the greatest between the SD-5d VS SD-15d groups. Real-time reverse transcription‒quantitative polymerase chain reaction analysis results of eight DEGs were consistent with the sequencing results. Thus, the sequencing results were accurate and reliable and eight genes were identified as candidates for the regulation of short-day induction at the adzuki bean seedling stage.

Conclusions

Short-day induction was able to down-regulate the expression of genes related to flowering according to the circadian rhythm and up-regulate the expression of certain genes in the antenna protein pathway. The results provide a theoretical reference for the molecular mechanism of short-day induction and multi-level information for future functional studies to verify the key genes regulating adzuki bean flowering.

Critical Differential Expression Assessment for Individual Bulk RNA-Seq Projects

Finding the right balance of quality and quantity can be important, and it is essential that project quality does not drop below the level where important main conclusions are missed or misstated. We use knock-out and over-expression studies as a simplification to test recovery of a known causal gene in RNA-Seq cell line experiments. When single-end RNA-Seq reads are aligned with STAR and quantified with htseq-count, we found potential value in testing the use of the Generalized Linear Model (GLM) implementation of edgeR with robust dispersion estimation more frequently for either single-variate or multi-variate 2-group comparisons (with the possibility of defining criteria less stringent than |fold-change| > 1.5 and FDR < 0.05). When considering a limited number of patient sample comparisons with larger sample size, there might be some decreased variability between methods (except for DESeq1). However, at the same time, the ranking of the gene identified using immunohistochemistry (for ER/PR/HER2 in breast cancer samples from The Cancer Genome Atlas) showed as possible shift in performance compared to the cell line comparisons, potentially highlighting utility for standard statistical tests and/or limma-based analysis with larger sample sizes. If this continues to be true in additional studies and comparisons, then that could be consistent with the possibility that it may be important to allocate time for potential methods troubleshooting for genomics projects. Analysis of public data presented in this study does not consider all experimental designs, and presentation of downstream analysis is limited. So, any estimate from this simplification would be an underestimation of the true need for some methods testing for every project. Additionally, this set of independent cell line experiments has a limitation in being able to determine the frequency of missing a highly important gene if the problem is rare (such as 10% or lower). For example, if there was an assumption that only one method can be tested for "initial" analysis, then it is not completely clear to the extent that using edgeR-robust might perform better than DESeq2 in the cell line experiments. Importantly, we do not wish to cause undue concern, and we believe that it should often be possible to define a gene expression differential expression workflow that is suitable for some purposes for many samples. Nevertheless, at the same time, we provide a variety of measures that we believe emphasize the need to critically assess every individual project and maximize confidence in published results.

Association of Chromosome 17 Aneuploidy, TP53 Deletion, Expression and Its rs1042522 Variant with Multiple Myeloma Risk and Response to Thalidomide/Bortezomib Treatment

Multiple myeloma (MM) is a multifactorial genetic disorder caused by interactive effects of environmental and genetic factors. The proper locus of the TP53 gene (17p13.1) and its protein is essential in genomic stability. The most common variant of the TP53 gene-p.P72R (rs1042522)-shows functional variation. The aim of our study was a complex analysis of the TP53 p.P72R variant and TP53 gene expression in relation to chromosomal changes of the TP53 gene locus, as well as MM risk and outcome. Genomic DNA from 129 newly diagnosed MM patients was analyzed by methods of automated DNA sequencing (for TP53 variant analysis) and cIg-FISH (for chromosomal aberrations analysis). RNA was used in real-time PCR to determine the TP53 expression. In MM patients, the TP53 variant was not in Hardy-Weinberg equilibrium. The RR genotype was associated with lower MM risk (OR = 0.44, p = 0.004). A higher number of plasma cells was found in patients with RR genotype in comparison to those with PP + PR genotypes (36.74% vs. 28.30%, p = 0.02). A higher expression of the TP53 gene was observed in PP + PR genotypes vs. RR homozygote (p < 0.001), in smokers vs. non-smokers (p = 0.02). A positive Pearson's correlation was found between the TP53 expression level and the number of plasma cells (r = 0.26, p = 0.04). The presence of chromosome 17 aberrations with or without TP53 locus did not affect the MM risk and outcome. Similar results were observed in the case of TP53 gene expression and the p.P72R variant.

A field-wide assessment of differential expression profiling by high-throughput sequencing reveals widespread bias

We assess inferential quality in the field of differential expression profiling by high-throughput sequencing (HT-seq) based on analysis of datasets submitted from 2008 to 2020 to the NCBI GEO data repository. We take advantage of the parallel differential expression testing over thousands of genes, whereby each experiment leads to a large set of p-values, the distribution of which can indicate the validity of assumptions behind the test. From a well-behaved p-value set π0, the fraction of genes that are not differentially expressed can be estimated. We found that only 25% of experiments resulted in theoretically expected p-value histogram shapes, although there is a marked improvement over time. Uniform p-value histogram shapes, indicative of <100 actual effects, were extremely few. Furthermore, although many HT-seq workflows assume that most genes are not differentially expressed, 37% of experiments have π0-s of less than 0.5, as if most genes changed their expression level. Most HT-seq experiments have very small sample sizes and are expected to be underpowered. Nevertheless, the estimated π0-s do not have the expected association with N, suggesting widespread problems of experiments with controlling false discovery rate (FDR). Both the fractions of different p-value histogram types and the π0 values are strongly associated with the differential expression analysis program used by the original authors. While we could double the proportion of theoretically expected p-value distributions by removing low-count features from the analysis, this treatment did not remove the association with the analysis program. Taken together, our results indicate widespread bias in the differential expression profiling field and the unreliability of statistical methods used to analyze HT-seq data.

Prediction of RNA-interacting residues in a protein using CNN and evolutionary profile

This paper describes a method Pprint2, which is an improved version of Pprint developed for predicting RNA-interacting residues in a protein. Training and independent/validation datasets used in this study comprises of 545 and 161 non-redundant RNA-binding proteins, respectively. All models were trained on training dataset and evaluated on the validation dataset. The preliminary analysis reveals that positively charged amino acids such as H, R and K, are more prominent in the RNA-interacting residues. Initially, machine learning based models have been developed using binary profile and obtain maximum area under curve (AUC) 0.68 on validation dataset. The performance of this model improved significantly from AUC 0.68 to 0.76, when evolutionary profile is used instead of binary profile. The performance of our evolutionary profile-based model improved further from AUC 0.76 to 0.82, when convolutional neural network has been used for developing model. Our final model based on convolutional neural network using evolutionary information achieved AUC 0.82 with Matthews correlation coefficient of 0.49 on the validation dataset. Our best model outperforms existing methods when evaluated on the independent/validation dataset. A user-friendly standalone software and web-based server named 'Pprint2' has been developed for predicting RNA-interacting residues (https://webs.iiitd.edu.in/raghava/pprint2 and https://github.com/raghavagps/pprint2).

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.

Gene expression profiling reveals candidate genes for defining spider silk gland types

Molecular studies of the secretory glands involved in spider silk production have revealed candidate genes for silk synthesis and a complicated history of spider silk gene evolution. However, differential gene expression profiles of the multiple silk gland types within an individual orb-web weaving spider are lacking. Each of these gland types produces a functionally distinct silk type. Comparison of gene expression among spider silk gland types would provide insight into the genes that define silk glands generally from non-silk gland tissues, and the genes that define silk glands from each other. Here, we perform 3' tag digital gene expression profiling of the seven silk gland types of the silver garden orb weaver Argiope argentata. Five of these gland types produce silks that are non-adhesive fibers, one silk includes both fibers and glue-like adhesives, and one silk is exclusively glue-like. We identify 1275 highly expressed, significantly upregulated, and tissue specific silk gland specific transcripts (SSTs). These SSTs include seven types of spider silk protein encoding genes known as spidroin genes. We find that the fiber-producing major ampullate and minor ampullate silk glands have more similar expression profiles than any other pair of glands. We also find that a subset of the SSTs is enriched for transmembrane transport and oxidoreductases, and that these transcripts highlight differences and similarities among the major ampullate, minor ampullate, and aggregate silk glands. Furthermore, we show that the wet glue-producing aggregate glands have the most unique SSTs, but still share some SSTs with fiber producing glands. Aciniform glands were the only gland type to share a majority of SSTs with other silk gland types, supporting previous hypotheses that duplication of aciniform glands and subsequent divergence of the duplicates gave rise to the multiple silk gland types within an individual spider.

Application of miRNA-seq in neuropsychiatry: A methodological perspective

MiRNAs are emerging as key molecules to study neuropsychiatric diseases. However, despite the large number of methodologies and software for miRNA-seq analyses, there is little supporting literature for researchers in this area. This review focuses on evaluating how miRNA-seq has been used to study neuropsychiatric diseases to date, analyzing both the main findings discovered and the bioinformatics workflows and tools used from a methodological perspective. The objective of this review is two-fold: first, to evaluate current miRNA-seq procedures used in neuropsychiatry; and second, to offer comprehensive information that can serve as a guide to new researchers in bioinformatics. After conducting a systematic search (from 2016 to June 30, 2020) of articles using miRNA-seq in neuropsychiatry, we have seen that it has already been used for different types of studies in three main categories: diagnosis, prognosis, and mechanism. We carefully analyzed the bioinformatics workflows of each study, observing a high degree of variability with respect to the tools and methods used and several methodological complexities that are identified and discussed in this review.

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.

Effect of ladostigil treatment of aging rats on gene expression in four brain areas associated with regulation of memory

Episodic and spatial memory decline in aging and are controlled by the hippocampus, perirhinal, frontal and parietal cortices and the connections between them. Ladostigil, a drug with antioxidant and anti-inflammatory activity, was shown to prevent the loss of episodic and spatial memory in aging rats. To better understand the molecular effects of aging and ladostigil on these brain regions we characterized the changes in gene expression using RNA-sequencing technology in rats aged 6 and 22 months. We found that the changes induced by aging and chronic ladostigil treatment were brain region specific. In the hippocampus, frontal and perirhinal cortex, ladostigil decreased the overexpression of genes regulating calcium homeostasis, ion channels and those adversely affecting synaptic function. In the parietal cortex, ladostigil increased the expression of several genes that provide neurotrophic support, while reducing that of pro-apoptotic genes and those encoding pro-inflammatory cytokines and their receptors. Ladostigil also decreased the expression of axonal growth inhibitors and those impairing mitochondrial function. Together, these actions could explain the protection by ladostigil against age-related memory decline.

consensusDE: an R package for assessing consensus of multiple RNA-seq algorithms with RUV correction

Extensive evaluation of RNA-seq methods have demonstrated that no single algorithm consistently outperforms all others. Removal of unwanted variation (RUV) has also been proposed as a method for stabilizing differential expression (DE) results. Despite this, it remains a challenge to run multiple RNA-seq algorithms to identify significant differences common to multiple algorithms, whilst also integrating and assessing the impact of RUV into all algorithms. consensusDE was developed to automate the process of identifying significant DE by combining the results from multiple algorithms with minimal user input and with the option to automatically integrate RUV. consensusDE only requires a table describing the sample groups, a directory containing BAM files or preprocessed count tables and an optional transcript database for annotation. It supports merging of technical replicates, paired analyses and outputs a compendium of plots to guide the user in subsequent analyses. Herein, we assess the ability of RUV to improve DE stability when combined with multiple algorithms and between algorithms, through application to real and simulated data. We find that, although RUV increased fold change stability between algorithms, it demonstrated improved FDR in a setting of low replication for the intersect, the effect was algorithm specific and diminished with increased replication, reinforcing increased replication for recovery of true DE genes. We finish by offering some rules and considerations for the application of RUV in a consensus-based setting. consensusDE is freely available, implemented in R and available as a Bioconductor package, under the GPL-3 license, along with a comprehensive vignette describing functionality: http://bioconductor.org/packages/consensusDE/.