Differential gene expression analysis using coexpression and RNA-Seq data

BackgroundDifferential gene expression analysis using RNA-seq data is a popular approach for discovering specific regulation mechanisms under certain environmental settings. Both gene ontology (GO) and KEGG pathway enrichment analysis are major processes for investigating gene groups that participate in common biological responses or possess related functions. However, traditional approaches based on differentially expressed genes only detect a few significant GO terms and pathways, which are frequently insufficient to explain all-inclusive gene regulation mechanisms.MethodsTranscriptomes of survivin (birc5) gene knock-down experimental and wild-type control zebrafish embryos were sequenced and assembled, and a differential expression (DE) gene list was obtained for traditional functional enrichment analysis. In addition to including DE genes with significant fold-change levels, we considered additional associated genes near or overlapped with differentially expressed long noncoding RNAs (DE lncRNAs), which may directly or indirectly activate or inhibit target genes and play important roles in regulation networks. Both the original DE gene list and the additional DE lncRNA-associated genes were combined to perform a comprehensive overrepresentation analysis.ResultsIn this study, a total of 638 DE genes and 616 DE lncRNA-associated genes (lncGenes) were leveraged simultaneously in searching for significant GO terms and KEGG pathways. Compared to the traditional approach of only using a differential expression gene list, the proposed method of employing DE lncRNA-associated genes identified several additional important GO terms and KEGG pathways. In GO enrichment analysis, 60% more GO terms were obtained, and several neuron development functional terms were retrieved as complete annotations. We also observed that additional important pathways such as the FoxO and MAPK signaling pathways were retrieved, which were shown in previous reports to play important roles in apoptosis and neuron development functions regulated by the survivin gene.ConclusionsWe demonstrated that incorporating genes near or overlapped with DE lncRNAs into the DE gene list outperformed the traditional enrichment analysis method for effective biological functional interpretations. These hidden interactions between lncRNAs and target genes could facilitate more comprehensive analyses.

BackgroundRNA sequencing (RNA-seq) is a powerful tool for genome-wide expression profiling of biological samples with the advantage of high-throughput and high resolution. There are many existing algorithms nowadays for quantifying expression levels and detecting differential gene expression, but none of them takes the misaligned reads that are mapped to non-exonic regions into account. We developed a novel algorithm, XBSeq, where a statistical model was established based on the assumption that observed signals are the convolution of true expression signals and sequencing noises. The mapped reads in non-exonic regions are considered as sequencing noises, which follows a Poisson distribution. Given measureable observed and noise signals from RNA-seq data, true expression signals, assuming governed by the negative binomial distribution, can be delineated and thus the accurate detection of differential expressed genes.ResultsWe implemented our novel XBSeq algorithm and evaluated it by using a set of simulated expression datasets under different conditions, using a combination of negative binomial and Poisson distributions with parameters derived from real RNA-seq data. We compared the performance of our method with other commonly used differential expression analysis algorithms. We also evaluated the changes in true and false positive rates with variations in biological replicates, differential fold changes, and expression levels in non-exonic regions. We also tested the algorithm on a set of real RNA-seq data where the common and different detection results from different algorithms were reported.ConclusionsIn this paper, we proposed a novel XBSeq, a differential expression analysis algorithm for RNA-seq data that takes non-exonic mapped reads into consideration. When background noise is at baseline level, the performance of XBSeq and DESeq are mostly equivalent. However, our method surpasses DESeq and other algorithms with the increase of non-exonic mapped reads. Only in very low read count condition XBSeq had a slightly higher false discovery rate, which may be improved by adjusting the background noise effect in this situation. Taken together, by considering non-exonic mapped reads, XBSeq can provide accurate expression measurement and thus detect differential expressed genes even in noisy conditions.

Differential Expression of Genes Associated with Oncogene-Induced Senescence and Senescence Associated Secretory Phenotype in the Absence of Differential Expression of High Molecular Risk Genes and Genes Associated with JAK-STAT Pathway in Sorted Cells of Patients with Polycythemia Vera and Primary Myelofibrosis

Heterosis is defined as the advantage of hybrid performance over its parents in growth and productivity. Previous studies showed that differential gene expression in hybrids and their parents is responsible for the heterosis. Plant hormones, espe-cially GAs are found to be related to heterosis in development of shoot and plant height. However, information on systematic identification and characterization of the differentially expressed genes in shoots of hybrid and their parents are not found, and molecular evidence for the contribution of plant hormone to plant heterosis are also absent. The purpose of this paper was to in-vestigate the differential gene express in stems of hybrid and their parents, and further, to investigate the possible role of genes participating hormone signaling in the formation of heterosis of plant stem. Using a wheat cross combination that exhibited vig-orous heterosis in plant height, the differential gene expression in upper-most internode of wheat hybrids and their parents was analyzed by cDNA-AFLP. Ninety-eight differentially expressing cDNA fragments were obtained, which were sequenced and analyzed by Blast-x. Reverse Northern hybridization was performed in order to test the expression pattern of isolated ESTs in hybrid and their parents. The differentially expressed genes included a cdc2 gene and a PAS1 gene, which were induced by GAs and CTKs respectively and participated in plant cell division. The cdc2 gene was up-regulated in hybrid while the expression of PAS1 in hybrid was near the higher parent. The differentially expressed genes also included kinase participating in signal trans-duction, genes involved in metabolism and genes related to transcription regulations. It can be inferred that numerous genes were differentially expressed in the upper-most internode of hybrids and their parents, including genes participating in cell division and cell expansion, kinase related to signal transduction and so on. Reverse Northern blot analysis indicated that the expression pattern in our cDNA-AFLP analysis is credible. Since it has been much reported that genes participating in cell division and cell expan-sion are involved in plant stem growth and elongation, and genes that participate in gibberellins signals have been proved to be up-regulated in wheat hybrid and contributed to heterosis in wheat plant height, it can be modestly speculated that the differential expression of genes involved in phytohormone signals may be related to heterosis in wheat intenode elongation and plant height.

Article Figures and data Abstract Editor's evaluation Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Replication of the genome must be coordinated with gene transcription and cellular metabolism, especially following replication stress in the presence of limiting deoxyribonucleotides. The Saccharomyces cerevisiae Rad53 (CHEK2 in mammals) checkpoint kinase plays a major role in cellular responses to DNA replication stress. Cell cycle regulated, genome-wide binding of Rad53 to chromatin was examined. Under replication stress, the kinase bound to sites of active DNA replication initiation and fork progression, but unexpectedly to the promoters of about 20% of genes encoding proteins involved in multiple cellular functions. Rad53 promoter binding correlated with changes in expression of a subset of genes. Rad53 promoter binding to certain genes was influenced by sequence-specific transcription factors and less by checkpoint signaling. However, in checkpoint mutants, untimely activation of late-replicating origins reduces the transcription of nearby genes, with concomitant localization of Rad53 to their gene bodies. We suggest that the Rad53 checkpoint kinase coordinates genome-wide replication and transcription under replication stress conditions. Editor's evaluation The unexpected localization of a cell cycle checkpoint kinase, Rad53, to promoters in response to replication stress suggests that Rad53 may help coordinate transcription in response to disrupted replication. This work will be of interest to those interested in the interplay between genome stability and gene expression. https://doi.org/10.7554/eLife.84320.sa0 Decision letter eLife's review process Introduction Rad53, the homolog of the mammalian CHEK2 kinase present in the budding yeast Saccharomyces cerevisiae, is an essential kinase involved in a multitude of cellular processes. The most well-established function of Rad53 is its role as an effector kinase in response to various sources of DNA damage. To maintain genome stability during S-phase, a DNA replication checkpoint (DRC) is activated in response to replication stress via the sensor kinase Mec1 (the mammalian ATM/ATR), the replication fork protein Mrc1 (Claspin in mammals), and other fork proteins (Lanz et al., 2019; Osborn and Elledge, 2003; Pardo et al., 2017; Paulovich and Hartwell, 1995; Saldivar et al., 2017). A second DNA damage checkpoint mediated by Rad9 (TP53BP1 in mammals) responds to double strand DNA breaks. Both branches converge on the effector kinase Rad53, which triggers a wide range of downstream events, including slowing or halting cell cycle progression, transiently preventing initiation of DNA replication at origins that have not yet replicated, activating DNA repair pathways and elevating the synthesis of deoxyribonucleoside triphosphates (dNTPs) by upregulating genes encoding the ribonucleotide reductase (RNR). The checkpoint signaling also promotes widespread changes in gene expression (Jaehnig et al., 2013; Pardo et al., 2017). Unlike most of the checkpoint genes, both Mec1 and Rad53 kinases are essential for cell viability in unperturbed cells. This can be partly explained by their role in regulating dNTP pools (Desany et al., 1998; Forey et al., 2020; Zhao et al., 2000). However, under the bypass conditions in cells without Sml1, the inhibitor of the RNR enzyme, cells lacking either kinases are viable. Yet, it is important to note that these kinase-null mutants are extremely sick and sensitive to various types of exogenous stress. Furthermore, lacking Rad53 causes a more severe defect than lacking Mec1, implying that Rad53 has activities beyond checkpoint signaling. Consistent with this suggestion, the kinase-deficient mutant rad53K227A lacks checkpoint function but retains growth-associated activity and an additional rad53 mutation was found that is checkpoint proficient but supports cell growth poorly (Gunjan and Verreault, 2003; Hoch et al., 2013; Holzen and Sclafani, 2010; Pellicioli et al., 1999). Rad53 is central to the transcriptional response to DNA damage. It is well established that the Dun1 protein kinase acts downstream of Rad53 to phosphorylate and inactivate the transcriptional repressor Rfx1/Crt1 and thereby upregulate target genes (Huang et al., 1998), including RNR2, RNR3, and RNR4, encoding subunits of RNR. However, the induced expression of RNR1, which encodes the major isoform of the RNR large subunit, is not controlled by the Rfx1 repressor, but by Ixr1 binding to the RNR1 promoter upon genotoxic stress (Tsaponina et al., 2011). This Ixr1-dependent regulation of RNR1 is independent of Dun1 but requires Rad53. Another Rad53-dependent, Dun1-independent regulation of RNR1 involves dissociation of the Nrm1 repressor from MBF following Nrm1 phosphorylation by Rad53 (Travesa et al., 2012). Thus, the Rad53-dependent transcription control is complex, but not yet fully understood. In addition to upregulating the dNTP pools, defects in cells lacking Rad53 can be suppressed by manipulating factors functioning in transcription regulation, cell wall maintenance, proteolysis, and cell cycle control (Desany et al., 1998; Manfrini et al., 2012). Moreover, Rad53 kinase targets and interaction partners found in biochemical and proteomic studies suggest that the kinase is pleiotropic (Gunjan and Verreault, 2003; Jaehnig et al., 2013; Lao et al., 2018; Smolka et al., 2007; Smolka et al., 2006). In this study, while investigating the response of yeast cells to replication stresses caused by depletion of dNTPs, we found that Rad53 localized to more than 20% of gene promoters in the S. cerevisiae genome, suggesting a multifaceted role in coordinating stress responses. Furthermore, we provide evidence that in the absence of the DRC checkpoint pathway, untimely activation of replication from late origins can negatively affect transcription activity of nearby genes. Results Rad53 is recruited to genomic loci other than replication origins in proliferating yeast cells Previous studies focusing on specific genomic features, mainly subsets of DNA replication origins, have detected Rad53 binding to replication origins and replication forks (Can et al., 2019; Dohrmann and Sclafani, 2006). However, given the evidence that Rad53 has functions beyond the control of DNA replication stress, it is possible that Rad53 also functions at additional genomic loci. We therefore investigated the genome-wide distribution of Rad53 by chromatin immunoprecipitation and deep sequencing (ChIP-seq) in proliferating yeast cells (Behrouzi et al., 2016). Since the sml1 null mutation (sml1Δ) allows cells to bypass the requirement for Rad53, or its upstream kinase Mec1, for growth, we also performed ChIP-seq in the sml1Δ mutant and the rad53 null (rad53Δ sml1Δ) as controls for antibody specificity (Figure 1). Figure 1 Download asset Open asset Rad53 is recruited to genomic loci other than replication origins in proliferating yeast cells. (a) Coverage tracks of Rad53 ChIP-seq signals in WT, sml1Δ, rad53Δsml1Δ, and mec1Δ sml1Δ for chromosome III. Asynchronous yeast cultures were processed for ChIP-seq analysis for distribution of Rad53 at genomic loci. The results from two independent experiments are presented. Experiment 1 compares only sml1Δ, rad53Δsml1Δ. (b) Heatmaps and average signals of Rad53 across 2-kb intervals centered on transcription start sites (TSSs) for proliferating WT, sml1Δ, rad53Δsml1Δ, and mec1Δ sml1Δ cells. (c) Heatmaps of ChIP-seq signals across 30-kb intervals centered on all origins annotated in OriDB database (Siow et al., 2012). Visual inspection of the ChIP-seq peaks from normalized coverage tracks suggested that Rad53 bound to sites throughout the genome (Figure 1a). Rad53 binding to the same sites also occurred in the sml1Δ mutant and the mec1Δ sml1Δ mutant, but was absent in the rad53Δ sml1Δ mutant. Computational analysis was performed to generate heatmaps of Rad53 ChIP-seq signal across 2-kb intervals centered on transcription start sites (TSSs) (Figure 1b). The average Rad53 ChIP-seq signal on top of the heatmap show that Rad53 was most concentrated within 500-bp upstream of TSSs in WT, sml1Δ, and mec1Δ sml1Δ mutants. A smaller fraction of genes also have signal within the gene bodies. No Rad53 signal beyond background was detected at TSSs in the rad53Δ sml1Δ cells, demonstrating antibody specificity. Analysis by calculating empirical p values for the ChIP-seq signals 500-bp upstream of TSS found that, out of 6604 genes, 1464, 1983, and 1293 in WT, sml1Δ, and mec1Δ sml1Δ, respectively, exhibited significant Rad53 signal at this region. This represents ~20% or greater of all gene promoters. Heatmap analysis was also performed on Rad53 ChIP-seq signal across 30-kb intervals centered on all origins annotated in the OriDB database (Siow et al., 2012). The average signal across the whole region was higher in WT, sml1Δ, and mec1Δ sml1Δ than in rad53Δ sml1Δ after normalization (Figure 1c). However, the signal was not concentrated at the origins. Since these ChIP-seq datasets are from asynchronous cells, this finding is consistent with the idea that recruitment of Rad53 to replication origins is cell cycle regulated rather than constant binding to origins. Binding of Rad53 to upstream TSS changes with cell cycle stages To gain insight into the dynamics of Rad53 recruitment to genomic loci such as promoters and replication origins, Rad53 ChIP-seq was analyzed in samples from synchronized cell cultures. Cells were arrested in G1 using α-factor and then released into media containing hydroxyurea (HU) to induce replication stress caused by limiting dNTPs. Three stages of synchronous cell cultures were collected and referred to as G1 (for cells arrested in G1), HU45 and HU90 (for cells released from G1 into HU for 45 and 90 min, respectively). These cell samples were then processed for ChIP-seq analysis. Heatmaps of the Rad53 signals at 2-kb intervals centered on all TSSs show a trend of increasing Rad53 binding as cells progress from G1-phase into HU45 or HU90 (Figure 2a), concomitant with increased levels of Rad53 protein in cells treated with HU (Figure 2—figure supplement 1). The increase in Rad53 parallels entry into S-phase, as measured by Orc6 phosphorylation (Figure 2—figure supplement 1a). Similar to the data from asynchronous cell samples, Rad53 signals were most concentrated upstream of TSSs (Figure 2a, left panel). The ChIP-seq signals using antibodies against γ-H2A, which like Rad53 is a target of the sensor kinase Mec1, were not enriched upstream of TSSs (Figure 2a, right panel). In fact, the signals for γ-H2A were lower immediately upstream of TSSs than the surrounding regions, consistent with promoters being histone-free regions. Figure 2 with 4 supplements see all Download asset Open asset Rad53 is recruited to transcription start site (TSS) and the binding changes with the cell cycle stage. Cells were synchronized in G1 phase and released into YPD containing 0.2M hydroxyurea (HU) for 45 and 90min (HU45 and HU90, respectively). (a) Heatmaps and average signals of Rad53 and γ-H2A ChIP-seq signals across 2-kb intervals centered on TSSs for WT cells at stages of G1, HU45, and HU90. (b) Distribution of aggregated peak numbers around TSS using merged Rad53 ChIP-seq peaks from all three stages (G1, HU45, and HU90). Pie charts showing the distribution of Rad53 ChIP-seq peaks in relation to gene features. (c) Rad53 ChIP-seq profiles near RNR1, PCL1, and TOS6 genes in proliferating WT, sml1Δ, rad53Δ sml1Δ, and mec1Δ sml1Δ cells, and WT cells at stages of G1, HU45, and HU90. (d) Lorenz curves for ChIP-seq read counts of Swi6 and Rad53 mapped to gene promoters showing inequality for promoter binding. Further analysis showed that >85% of the Rad53 ChIP-seq peaks overlap with gene promoters (defined as 500-bp upstream and 50-bp downstream of TSS) at all stages investigated. The distribution of aggregated peak numbers around TSSs showed highest count numbers immediately upstream of the TSS (Figure 2b). By visual inspection of normalized ChIP-seq coverage tracks, we noticed that signals at promoters of genes such as RNR1, PCL1, and TOS6 varied depending on the cell cycle stage (Figure 2c, bottom three tracks for WT at stages of G1, HU45, and HU90) while signals at neighboring gene promoters remained largely constant. Thus, the recruitment of Rad53 to a subset of gene promoters is regulated. Because RNR1, PCL1, and TOS6 are known targets of the cell cycle regulator SBF, a sequence-specific transcription factor composed of subunits Swi4 and Swi6, we compared Rad53 ChIP-seq data and a previous ChIP-seq dataset of Swi6 (Park et al., 2013). As indicators for the gene specificity of protein binding, the Gini indices were computed from Lorenz curves of ChIP-seq data for Swi6 and two of our Rad53 replicates (Figure 2d), being 0.763, 0.2918, and 0.2982, respectively. Rad53 has a coverage for many promoters while Swi6 shows substantially high coverage only for a limited number of promoters. Thus, it is likely Rad53 would effect a wider range of genes than the Swi6 regulatory network. Previous studies have found that under certain conditions, regions of the genome are promiscuously present in ChIP-seq studies independent of the antibody used (Park et al., 2013; Teytelman et al., 2013). These regions were termed hyper-ChIPable regions and were found enriched for sequences in and around gene bodies of highly expressed genes. We therefore examined whether these regions were promiscuously present under our conditions. Analysis of ChIP-seq data for Rad53 and γ-H2A with or without the hyper-ChIPable regions observed by Teytelman et al. and Park et al. did not alter the pattern of Rad53 binding to TSSs (Figure 2—figure supplement 2a–d). When we specifically examined the pattern of Rad53 ChIP-seq signals around TSSs of the 296 genes associated with those hyper-ChIPable sequences, we observed not only localization to TSSs, but also enrichment to the gene bodies as previously reported (Park et al., 2013; Teytelman et al., 2013). This pattern is distinct from the pattern of promoter localization (compare Figure 2—figure supplement 2d, e). Furthermore, KEGG analysis of the genes enriched in the studies by Teytelman et al. and Park et al. showed predominantly genes encoding snoRNAs and tRNAs, genes we did not find in the promoter binding for Rad53 (Figure 2—figure supplement 2f). Finally, we did not observe any enrichment at gene promoters when the RAD53 gene was deleted from the strain (Figure 1a, b) or when anti-γ-H2A antibodies were used (Figure 2a). Thus, we suggest that the Rad53 binding observed here is not the same as the promiscuous, non-specific enrichment of hyper-ChIPable regions reported previously. Moreover, Rad53 binding to a subset of promoters is transcription factor dependent (see below). Rad53 is recruited to sites of DNA synthesis independent of checkpoint signaling Since previous studies reported localization of Rad53 to replication origins (Can et al., 2019; Dohrmann and Sclafani, 2006), we also performed heatmap analysis of ChIP-seq signal around replication origins for Rad53, γ-H2A, and Cdc45, a component of active helicase complex and hence the marker for active replication forks (Figure 2—figure supplement 3). The 30-kb window for the heatmap around origins was chosen based on our prior knowledge that the extent of DNA synthesis in WT cells under similar growth conditions at HU90 was about 10 kb (Sheu et al., 2014). In WT cells, the Rad53 signal was present at regions associated with early firing origins but not with late firing origins that are inactive due to inhibition by the HU-induced checkpoint. However, Rad53 signal was also present at late origins in the kinase-deficient rad53K227A and mrc1Δ mutants, both of which allow activation of late origins as a result of the checkpoint defect. Thus, Rad53 was recruited only to replication forks associated with activated origins and the pattern was similar to that of Cdc45. Surprisingly, Rad53 binding to replication forks did not require Mrc1 or its own kinase activity, suggesting checkpoint-independent recruitment of Rad53 to sites of DNA synthesis. Interestingly, γ-H2A was observed at genomic regions surrounding the very late origins in G1-phase in both WT and mutants (Figure 2—figure supplement 3b). It is possible that these γ-H2A signals reflect a low level of ssDNA gaps at these late-replicating regions that were tolerated and carried over from the previous cell cycle, similar to unrepaired post-replication gaps resulting from low level of UV irradiation found in S. pombe G2-phase (Callegari and Kelly, 2006). Identification of genes with differential binding of Rad53 at promoters Rad53 promoter binding was temporally dynamic in a subset of genes, suggesting regulation by cell cycle progression or DNA replication stress. To identify genes with differential or dynamic binding of Rad53 at their promoters (DB genes), we applied residual analysis. The read count difference was investigated for the promoter regions (500-bp window upstream of TSSs) of all genes for Rad53 ChIP-seq. The comparison was done between stages G1 and HU45 of WT samples from two independent experiments (termed CP and TF, Figure 3a). Since each dataset had two biological replicates from stages G1 and HU45, we extracted the top 1000 genes that displayed dynamic binding from the aggregated read coverage and called these genes the top DB genes (Figure 3a). Among the 1000 top DB genes from each set of comparison, 435 genes were identified in both sets (435 top DB overlap). Figure 3 Download asset Open asset Identification of genes with Rad53-binding changes at the promoters. (a) Scatter plots compare the signals in G1 (x-axis) and HU45 (y-axis) at 500-bp intervals upstream of transcription start sites (TSSs) for all genes in WT datasets from two independent experiments, TF and CP. Orange dots indicated the 1000 genes with highest binding changes (top 1000 DB) and satisfying the filter of minimal signal of −0.075 (Maximal = 1). The Venn diagram to the right illustrates the identification of 435 genes that are found in both experiments (435 top DB overlap). (b) Binding changes (y-axis: residuals from analysis of experiment TF) for 435 top DB overlap. (c) Examples of coverage tracks for selected genes show Rad53 signal changes at the indicated gene promoters from G1 to HU45 and HU90. Overall, during the G1- to S-phase transition (HU45), there are more genes with increased Rad53 promoter binding than those with decreased binding: within the list of 435 genes, 337 show increased Rad53 binding at their promoters, while 98 show decreased binding (Figure 3b). These genes include those involved in cell cycle progression (e.g., genes encoding cyclins and regulators of DNA replication) and cell growth (e.g., cell wall maintenance and mating response). Figure 3c shows Rad53 dynamic binding, either up or down, at the promoters of representative genes as cells transitioned from G1 phase to HU45 and HU90 time points. The relationship between Rad53 promoter binding and gene expression To gain insight into the relationship between Rad53 promoter recruitment and gene expression, RNA-seq analysis was performed under the same growth conditions used for the ChIP-seq experiments. RNA-seq replicates of datasets from four yeast strains (WT, rad9Δ, rad53K227A, and mrc1Δ), each with three stages (G1, HU45, and HU90) were analyzed using rank data analysis (Figure 4a). In the hierarchical clustering, cell cycle stage contributes more to similarities in expression than the genotype. In particular, the expression profiles in G1 were very similar among all strains. In HU45, however, two subgroups clustered by genotype were evident: rad9Δ clustered with WT, consistent with Rad9 having no role in the DRC checkpoint branch, while rad53K227A and mrc1Δ clustered together, consistent with Rad53 and Mrc1 functioning together in the DRC response to HU stress. Figure 4 with 1 supplement see all Download asset Open asset Gene expression changes in WT and checkpoint mutants under stress and the tendency of higher Rad53 binding at promoter of genes with significant differential expression. WT, rad9Δ, rad53K227A, and mrc1Δ cells were synchronized in G1-phase and released into YPD containing 0.2M hydroxyurea (HU). Cells at stages of G1, HU45, and HU90 were collected and processed for RNA-seq analysis. (a) Rank data analysis of RNA-seq samples. (b) Bar graph summarizing on the x-axis the number of genes that show statistically significant differential expression (DEGs). The types of pairwise comparison are indicated to the right. Blue bars, downregulated DEGs. Orange bars, upregulated DEGs. (c) Average Rad53 ChIP-seq signal across 2-kb intervals centered on at transcription start site (TSS) for statistically significant DEGs (red) and the insignificant (cyan). Left panel shows the result of DEGs from comparison of G1 and HU45 datasets of WT (DEGs: WT(G1 → HU45)) and right panel shows DEGs from comparison of G1 and HU90 datasets (DEGs: WT(G1 → HU90)). Pairwise comparisons of G1 to HU45 samples showed that, in both WT and rad9Δ cells, ~2300 genes exhibited significant expression changes (differentially expressed genes, DEGs; Figure 4b). The number of DEGs increased to ~3000 when comparing G1 to HU90. Moreover, in both rad53K227A and mrc1Δ mutants, ~2500 DEGs were detected from G1 to HU45, which increased to >3400 in G1 to HU90. The response to cell cycle stage was largely equally distributed between up- and downregulation. Comparison of WT to rad53K227A in the HU45 and HU90 conditions found 517 and 2234 DEGs, respectively, while comparison of WT to mrc1Δ revealed 370 and 1166 DEGs in HU45 and HU90. A WT and rad9Δ comparison at all stages showed only five DEGs, including the deleted RAD9 gene and its marker cassette HIS3. Thus, Rad9 did not contribute much to gene expression changes under HU stress. on ChIP-seq the average Rad53-binding upstream of TSSs was higher in the significant DEGs than in the when comparing G1 to HU45 or to HU90 of the WT datasets (Figure and Figure supplement 1a). These results suggested that Rad53 may a role in the control of gene expression at a subset of expressed genes. Similar were also when comparing samples in the checkpoint mutants (Figure supplement demonstrating that the functions of Rad53 at gene promoters were not checkpoint we performed gene analysis based on in the et al., analysis of significant DEGs in WT(G1 → 10 (Figure Rad53 binding at promoter regions occurred in most (Figure with enriched functions including cell cycle regulation, mating proteolysis, and (Figure the 435 top DB overlap gene (Figure genes were also detected as DEGs in the WT(G1 → The between binding changes and expression changes was significant (Figure left p of Rad53-binding changes against gene expression changes of these revealed a between Rad53-binding and gene expression with a = (Figure left panel). Among this out of genes with decreased Rad53 signal were downregulated in genes with increased Rad53 signals into both up- and and respectively). Further of the by revealed that genes in 1 and exhibited the most significant between Rad53-binding and gene expression changes (Figure = and = respectively). Thus, specific subsets of DEGs in the from G1 → HU exhibited between a in gene expression and Rad53 promoter binding. Figure Download asset Open asset between differential binding of Rad53 at promoter and differential gene expression (a) for significant DEGs in WT(G1 → for Gene level of differential expression in (b) Heatmaps of Rad53 ChIP-seq signal across 2-kb intervals centered on transcription start site (TSS) by the in within each are in to the differential expression level (c) Scatter plots of Rad53-binding changes at the promoter against expression changes for the significant DEGs in the 435 top DB overlap and subgroups in and is indicated on top of each number of genes in the (d) of for between binding changes and expression changes for in not Examples of genes in the are the mutants of gene expression near promiscuously active late origins We the localization of Rad53 in a subset of from the HU45 comparison (Figure and the HU45 comparison (Figure supplement 1a). In these we noticed a in which downregulated genes to have a Rad53 signal not only upstream of the but a signal within gene bodies (Figure and Figure supplement 1b). This pattern was in the mrc1Δ mutant in HU45 and in HU90. The gene localization was also found more transiently in HU45 from rad53K227A cells. This signal pattern was not as in the WT HU45 and HU90 samples. Since Rad53 is also recruited to active origins and with the replication we that these gene signals in the checkpoint mutants were caused by the activation of origins that are inactive in WT, between DNA replication and gene The of the Rad53 localization at gene bodies in this of genes in the rad53K227A mutant is also consistent with the signal pattern at origins (Figure 2—figure supplement Thus, we investigated the relationship between these genes and their replication origins. Figure with 1 supplement see all Download asset Open asset expressed genes are in the mrc1Δ mutant. DEGs from HU45 comparison were (a) for significant DEGs. for Gene level of differential expression in (b) Heatmaps of Rad53 ChIP-seq signal across 2-kb intervals centered on transcription start site (TSS) by the in ChIP-seq signal in WT, rad53K227A, and mrc1Δ mutant cells at stages G1, HU45, and HU90 are (c) of relation for DEGs between each TSS and its is indicated in as well as and and are (d) charts number (top and of and upregulated DEGs as by TSS to within of TSS no origins within of and The from the TSS of each gene to the replication origins, the of gene transcription to the or and the late or see of the types in Materials and was and in the same as in the heatmaps of in Figure (Figure Overall, most of the downregulated genes in 1 of this were very to active origins kb between and and Interestingly, the pattern of for to TSS largely the of the Rad53 signal within the gene bodies (Figure This pattern was not found in the WT To the relationship between replication origins and genes, we the of up- and downregulated DEGs from the HU45 comparison in of the of their relation to the origins (Figure In the we this to the range of between TSS and origins (Figure left this genes kb or more from the active origins were distributed between up- and of gene expression. However, for those genes that are to an active the downregulated genes those genes less than 1 kb from active origins, more than were

ObjectiveThe aim of this study was to investigate differentially expressed genes and their functions in the hippocampus and striatum after heroin administration in cynomolgus macaques of different ages.MethodsCynomolgus monkeys were divided by age as follows: 1 year (A1, n = 2); 3 to 4 years (A2, n = 2); 6 to 8 years (A3, n = 2); and older than 11 years (A4, n = 2). After heroin was injected intramuscularly into the monkeys (0.6 mg/kg), we performed large-scale transcriptome profiling in the hippocampus (H) and striatum (S) using RNA sequencing technology. Some genes were validated with real-time quantitative PCR.ResultsIn the hippocampus, the gene expression of A1H was similar to that of A4H, while the gene expression of A2H was similar to that of A3H. Genes associated with the mitogen-activated protein kinase signaling pathway (STMN1, FGF14, and MAPT) and γ-aminobutyric acid-ergic synapses (GABBR2 and GAD1) were differentially expressed among control and heroin-treated animals. Differential gene expression between A1S and A4S was the least significant, while differential gene expression between A3S and A2S was the most significant. Genes associated with the neurotrophin signaling pathway (NTRK1 and NGFR), autophagy (ATG5), and dopaminergic synapses (AKT1) in the striatum were differentially expressed among control and heroin-treated animals.ConclusionThese results suggest that even a single heroin exposure can cause differential gene expression in the hippocampus and striatum of nonhuman primates at different ages.

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

BackgroundDifferential expression and feature selection analyses are essential steps for the development of accurate diagnostic/prognostic classifiers of complicated human diseases using transcriptomics data. These steps are particularly challenging due to the curse of dimensionality and the presence of technical and biological noise. A promising strategy for overcoming these challenges is the incorporation of pre-existing transcriptomics data in the identification of differentially expressed (DE) genes. This approach has the potential to improve the quality of selected genes, increase classification performance, and enhance biological interpretability. While a number of methods have been developed that use pre-existing data for differential expression analysis, existing methods do not leverage the identities of experimental conditions to create a robust metric for identifying DE genes.ResultsIn this study, we propose a novel differential expression and feature selection method—GEOlimma—which combines pre-existing microarray data from the Gene Expression Omnibus (GEO) with the widely-applied Limma method for differential expression analysis. We first quantify differential gene expression across 2481 pairwise comparisons from 602 curated GEO Datasets, and we convert differential expression frequencies to DE prior probabilities. Genes with high DE prior probabilities show enrichment in cell growth and death, signal transduction, and cancer-related biological pathways, while genes with low prior probabilities were enriched in sensory system pathways. We then applied GEOlimma to four differential expression comparisons within two human disease datasets and performed differential expression, feature selection, and supervised classification analyses. Our results suggest that use of GEOlimma provides greater experimental power to detect DE genes compared to Limma, due to its increased effective sample size. Furthermore, in a supervised classification analysis using GEOlimma as a feature selection method, we observed similar or better classification performance than Limma given small, noisy subsets of an asthma dataset.ConclusionsOur results demonstrate that GEOlimma is a more effective method for differential gene expression and feature selection analyses compared to the standard Limma method. Due to its focus on gene-level differential expression, GEOlimma also has the potential to be applied to other high-throughput biological datasets.

MotivationThe sophisticated lfcShrink() function implemented in the DESeq2 package for differential gene expression analysis aims to reduce noise from low read count and/or highly variable genes in bulk RNA-sequencing data, thus circumventing the need for arbitrary filtering thresholds. However, difficulties can arise when analysing clinical data with multiple biologically-relevant groupings. In particular, changing the reference group can dramatically alter the ranking of differentially expressed genes, instead of merely ‘mirroring’ the up- and down-regulated genes in reciprocal comparisons.ResultsHere, we present mirrorCheck, an R package to automate methodical lfcShrink() usage and data visualization for quality control and data-driven decision-making during analysis.Availability and implementationThe source code, including documentation, is available on github at https://github.com/kescull/mirrorCheck.

BackgroundSocial dominance is important for the reproductive success of males in many species. In the black-faced blenny (Tripterygion delaisi) during the reproductive season, some males change color and invest in nest making and defending a territory, whereas others do not change color and ‘sneak’ reproductions when females lay their eggs. Using RNAseq, we profiled differential gene expression between the brains of territorial males, sneaker males, and females to study the molecular signatures of male dimorphism.ResultsWe found that more genes were differentially expressed between the two male phenotypes than between males and females, suggesting that during the reproductive period phenotypic plasticity is a more important factor in differential gene expression than sexual dimorphism. The territorial male overexpresses genes related to synaptic plasticity and the sneaker male overexpresses genes involved in differentiation and development.ConclusionsPreviously suggested candidate genes for social dominance in the context of alternative mating strategies seem to be predominantly species-specific. We present a list of novel genes which are differentially expressed in Tripterygion delaisi. This is the first genome-wide study for a molecular non-model species in the context of alternative mating strategies and provides essential information for further studies investigating the molecular basis of social dominance.

BackgroundMicroarray studies using in vitro cultures of synchronized, blood-stage Plasmodium falciparum malaria parasites have revealed a ‘just-in-time’ cascade of gene expression with some indication that these transcriptional patterns remain stable even in the presence of external stressors. However, direct analysis of transcription in P. falciparum blood-stage parasites obtained from the blood of infected patients suggests that parasite gene expression may be modulated by factors present in the in vivo environment of the host. The aim of this study was to examine changes in gene expression of the rodent malaria parasite, Plasmodium yoelii 17X, while varying the in vivo setting of replication.MethodsUsing P. yoelii 17X parasites replicating in vivo, differential gene expression in parasites isolated from individual mice, from independent infections, during ascending, peak and descending parasitaemia and in the presence and absence of host antibody responses was examined using P. yoelii DNA microarrays. A genome-wide analysis to identify coordinated changes in groups of genes associated with specific biological pathways was a primary focus, although an analysis of the expression patterns of two multi-gene families in P. yoelii, the yir and pyst-a families, was also completed.ResultsAcross experimental conditions, transcription was surprisingly stable with little evidence for distinct transcriptional states or for consistent changes in specific pathways. Differential gene expression was greatest when comparing differences due to parasite load and/or host cell availability. However, the number of differentially expressed genes was generally low. Of genes that were differentially expressed, many involved biologically diverse pathways. There was little to no differential expression of members of the yir and pyst-a multigene families that encode polymorphic proteins associated with the membrane of infected erythrocytes. However, a relatively large number of these genes were expressed during blood-stage infection regardless of experimental condition.ConclusionsTaken together, these results indicate that 1) P. yoelii gene expression remains stable in the presence of a changing host environment, and 2) concurrent expression of a large number of the polymorphic yir and pyst-a genes, rather than differential expression in response to specific host factors, may in itself limit the effectiveness of host immune responses.

BackgroundHuman monocytotropic ehrlichiosis is an emerging life-threatening zoonosis caused by obligately intracellular bacterium, Ehrlichia chaffeensis. E. chaffeensis is transmitted by the lone star tick, Amblyomma americanum, and replicates in mononuclear phagocytes in mammalian hosts. Differences in the E. chaffeensis transcriptome in mammalian and arthropod hosts are unknown. Thus, we determined host-specific E. chaffeensis gene expression in human monocyte (THP-1) and in Amblyomma and Ixodes tick cell lines (AAE2 and ISE6) using a whole genome microarray.Methodology/Principal FindingsThe majority (∼80%) of E. chaffeensis genes were expressed during infection in human and tick cells. There were few differences observed in E. chaffeensis gene expression between the vector Amblyomma and non-vector Ixodes tick cells, but extensive host-specific and differential gene expression profiles were detected between human and tick cells, including higher transcriptional activity in tick cells and identification of gene subsets that were differentially expressed in the two hosts. Differentially and host-specifically expressed ehrlichial genes encoded major immunoreactive tandem repeat proteins (TRP), the outer membrane protein (OMP-1) family, and hypothetical proteins that were 30–80 amino acids in length. Consistent with previous observations, high expression of p28 and OMP-1B genes was detected in human and tick cells, respectively. Notably, E. chaffeensis genes encoding TRP32 and TRP47 were highly upregulated in the human monocytes and expressed as proteins; however, although TRP transcripts were expressed in tick cells, the proteins were not detected in whole cell lysates demonstrating that TRP expression was post transcriptionally regulated.Conclusions/Significance Ehrlichia gene expression is highly active in tick cells, and differential gene expression among a wide variety of host-pathogen associated genes occurs. Furthermore, we demonstrate that genes associated with host-pathogen interactions are differentially expressed and regulated by post transcriptional mechanisms.

Background. A common research goal in transcriptome projects is to find genes that are differentially expressed in different phenotype classes. Biologists might wish to validate such gene candidates experimentally, or use them for downstream systems biology analysis. Producing a coherent differential gene expression analysis from RNA-seq count data requires an understanding of how numerous sources of variation such as the replicate size, the hypothesized biological effect size, and the specific method for making differential expression calls interact. We believe an explicit demonstration of such interactions in real RNA-seq data sets is of practical interest to biologists.Results. Using two large public RNA-seq data sets—one representing strong, and another mild, biological effect size—we simulated different replicate size scenarios, and tested the performance of several commonly-used methods for calling differentially expressed genes in each of them. We found that, when biological effect size was mild, RNA-seq experiments should focus on experimental validation of differentially expressed gene candidates. Importantly, at least triplicates must be used, and the differentially expressed genes should be called using methods with high positive predictive value (PPV), such as NOISeq or GFOLD. In contrast, when biological effect size was strong, differentially expressed genes mined from unreplicated experiments using NOISeq, ASC and GFOLD had between 30 to 50% mean PPV, an increase of more than 30-fold compared to the cases of mild biological effect size. Among methods with good PPV performance, having triplicates or more substantially improved mean PPV to over 90% for GFOLD, 60% for DESeq2, 50% for NOISeq, and 30% for edgeR. At a replicate size of six, we found DESeq2 and edgeR to be reasonable methods for calling differentially expressed genes at systems level analysis, as their PPV and sensitivity trade-off were superior to the other methods’.Conclusion. When biological effect size is weak, systems level investigation is not possible using RNAseq data, and no meaningful result can be obtained in unreplicated experiments. Nonetheless, NOISeq or GFOLD may yield limited numbers of gene candidates with good validation potential, when triplicates or more are available. When biological effect size is strong, NOISeq and GFOLD are effective tools for detecting differentially expressed genes in unreplicated RNA-seq experiments for qPCR validation. When triplicates or more are available, GFOLD is a sharp tool for identifying high confidence differentially expressed genes for targeted qPCR validation; for downstream systems level analysis, combined results from DESeq2 and edgeR are useful.

Introduction: Colon cancer is the fourth most common cancer in the United States and the third leading cause of cancer-related death. High throughput genomic sequencing has led to a number of significant advances in tumor biology and in identifying novel signaling molecules such as long non-coding RNA (lncRNA). The aim of this study was to identify differentially expressed lncRNA in colon cancer from a large RNA-sequencing (RNA-seq) data set. Methods: The raw RNA-seq files of 398 patients with colon cancer were downloaded from The Cancer Genome Atlas (TCGA). Sequencing files were aligned using STAR (Spliced Transcripts Alignment to a Reference) to the most recent genome annotation. Using a subset of patients with paired colon cancer and normal colon epithelium RNA-seq data (n=40), an exploratory binomial regression model was used to calculate differential RNA expression. The most differentially expressed lncRNA were identified from the exploratory analysis and verified by comparing the larger colon cancer RNA-seq data set (n=358) with the normal colon epithelium RNA-seq data set (n=40). Results: 33,514 genes were identified from the comparative analysis, and using differential expression cut off values of > +1.5 or < -1.5 log fold change and a false discovery rate of <0.05; 543 were upregulated and 1822 were downregulated. Within this dysregulated group, 60 lncRNAs were identified. Using the larger data set (n=358 vs n=40), 41/60 lncRNA remained differentially expressed, of which 15 were downregulated and 26 were upregulated. Twenty-four of these lncRNAs have not previously been described in colon cancer, 12 of which are upregulated and 12 of which are downregulated (Table 1). Conclusions: This analysis of RNA-seq data from TCGA has identified dysregulated lncRNAs which have not previously been described in human colon cancer. These lncRNAs may have significant roles in colon cancer tumor signaling. Table 1.Differential expression of lncRNA not described in colon cancerlncRNALog Fold ChangeP-valueDescription in other cancersKRT16P14.1041.36E-21Increased in lung squamous cell carcinomaBOK-AS13.7784.98E-09Increased in oral squamous cell carcinomaCELP3.1151.30E-09Not availableCLDN10-AS12.3944.61E-12Not availablePOU6F2-AS12.2678.25E-09Not availableSLC7A11-AS12.1931.60E-13Decreased in gastric adenocarcinomaCSAG42.1083.59E-10Not availableLUCAT12.0224.21E-11Increased in Ovarian cancer, renal clear cell carcinoma, head and neck squamous cell carcinoma, Non-small cell cancer, glioma, osteosarcoma, esophageal squamous cell carcinomaVPS9D1-AS11.9795.27E-21Decreased in gastric adenocarcinoma, Increased in non-small cell lung cancerLEF1-AS11.8464.48E-13Increased in glioblastomaCASC81.6485.64E-06Single nucleotide polymorphisms in colorectal adenocarcinomaLINC004911.5903.11E-04Not availableLINC00648-1.5021.02E-15Not availableLINC00163-1.6751.03E-22Not availableSATB2-AS1-1.8074.18E-09Increased in osteosarcomaLINC00702-2.1424.81E-18Not availableTRHDE-AS1-2.2623.44E-24Not availableLINC00908-2.3082.76E-103Not availableLINC00483-2.4012.87E-39Increased in gastric adenocarcinomaLINC00461-2.4401.42E-41Increased in glioma, multiple myelomaFENDRR-2.6872.12E-58Decreased in breast cancer, prostate cancer, osteosarcoma, gastric adenocarcinoma,ARHGEF26-AS1-2.7148.55E-32Not availableADAMTS9-AS2-2.9041.78E-74Increased in lung cancer salivary adenoid cystic carcinoma, Decreased in glioma, gastric adenocarcinomaMT1JP-4.4601.08E-114Decreased in retinoblastoma, Decreased in gastric adenocarcinoma Citation Format: Stephen J. O'Brien, Theodore Kalbfleisch, Sudhir Srivastava, Shesh Rai, Susan Galandiuk. Differential expression of long non-coding RNA in colon adenocarcinoma RNA-sequence data set [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2019; 2019 Mar 29-Apr 3; Atlanta, GA. Philadelphia (PA): AACR; Cancer Res 2019;79(13 Suppl):Abstract nr 1817.

BackgroundDifferential gene expression is important to understand the biological differences between healthy and diseased states. Two common sources of differential gene expression data are microarray studies and the biomedical literature.MethodsWith the aid of text mining and gene expression analysis we have examined the comparative properties of these two sources of differential gene expression data.ResultsThe literature shows a preference for reporting genes associated to higher fold changes in microarray data, rather than genes that are simply significantly differentially expressed. Thus, the resemblance between the literature and microarray data increases when the fold-change threshold for microarray data is increased. Moreover, the literature has a reporting preference for differentially expressed genes that (1) are overexpressed rather than underexpressed; (2) are overexpressed in multiple diseases; and (3) are popular in the biomedical literature at large. Additionally, the degree to which diseases are similar depends on whether microarray data or the literature is used to compare them. Finally, vaguely-qualified reports of differential expression magnitudes in the literature have only small correlation with microarray fold-change data.ConclusionsReporting biases of differential gene expression in the literature can be affecting our appreciation of disease biology and of the degree of similarity that actually exists between different diseases.