Dispersion estimation and its effect on test performance in RNA-seq data analysis: a simulation-based comparison of methods

Perturbations in the neuroactive ligand-receptor interaction and renin angiotensin system pathways are associated with cancer-related cognitive impairment

PURPOSE

This study reports on the results from our data-driven approach that identified perturbations in neuroactive ligand-receptor interaction and renin-angiotensin system (RAS) pathways in oncology patients with and without self-reported cancer-related cognitive impairment (CRCI).

METHODS

In a sample of oncology patients receiving chemotherapy (n = 1343), the Attentional Function Index (AFI) was used to assess CRCI. Patients were grouped into low (AFI score of < 5) versus high (AFI score of > 7.5) levels of cognitive function. Gene expression analyses were done using RNA-seq (n = 185) and microarray (n = 158) technologies. Pathway impact analysis was used to evaluate for perturbations in biological pathways associated with self-reported CRCI.

RESULTS

The combined pathway impact analysis revealed that the neuroactive ligand-receptor interaction and RAS pathways were significantly perturbed between the patients with low versus high AFI scores.

CONCLUSIONS

Findings from this study suggest that in addition to inflammatory pathways, numerous mechanisms may contribute to the underlying mechanisms for the development and/or persistence of self-reported CRCI.

Differentially expressed heterogeneous overdispersion genes testing for count data

The mRNA-seq data analysis is a powerful technology for inferring information from biological systems of interest. Specifically, the sequenced RNA fragments are aligned with genomic reference sequences, and we count the number of sequence fragments corresponding to each gene for each condition. A gene is identified as differentially expressed (DE) if the difference in its count numbers between conditions is statistically significant. Several statistical analysis methods have been developed to detect DE genes based on RNA-seq data. However, the existing methods could suffer decreasing power to identify DE genes arising from overdispersion and limited sample size, where overdispersion refers to the empirical phenomenon that the variance of read counts is larger than the mean of read counts. We propose a new differential expression analysis procedure: heterogeneous overdispersion genes testing (DEHOGT) based on heterogeneous overdispersion modeling and a post-hoc inference procedure. DEHOGT integrates sample information from all conditions and provides a more flexible and adaptive overdispersion modeling for the RNA-seq read count. DEHOGT adopts a gene-wise estimation scheme to enhance the detection power of differentially expressed genes when the number of replicates is limited as long as the number of conditions is large. DEHOGT is tested on the synthetic RNA-seq read count data and outperforms two popular existing methods, DESeq2 and EdgeR, in detecting DE genes. We apply the proposed method to a test dataset using RNAseq data from microglial cells. DEHOGT tends to detect more differently expressed genes potentially related to microglial cells under different stress hormones treatments.

Chemotherapy for pain: reversing inflammatory and neuropathic pain with the anticancer agent mithramycin A

ABSTRACT

The persistence of inflammatory and neuropathic pain is poorly understood. We investigated a novel therapeutic paradigm by targeting gene networks that sustain or reverse persistent pain states. Our prior observations found that Sp1-like transcription factors drive the expression of TRPV1, a pain receptor, that is blocked in vitro by mithramycin A (MTM), an inhibitor of Sp1-like factors. Here, we investigate the ability of MTM to reverse in vivo models of inflammatory and chemotherapy-induced peripheral neuropathy (CIPN) pain and explore MTM's underlying mechanisms. Mithramycin reversed inflammatory heat hyperalgesia induced by complete Freund adjuvant and cisplatin-induced heat and mechanical hypersensitivity. In addition, MTM reversed both short-term and long-term (1 month) oxaliplatin-induced mechanical and cold hypersensitivity, without the rescue of intraepidermal nerve fiber loss. Mithramycin reversed oxaliplatin-induced cold hypersensitivity and oxaliplatin-induced TRPM8 overexpression in dorsal root ganglion (DRG). Evidence across multiple transcriptomic profiling approaches suggest that MTM reverses inflammatory and neuropathic pain through broad transcriptional and alternative splicing regulatory actions. Mithramycin-dependent changes in gene expression following oxaliplatin treatment were largely opposite to and rarely overlapped with changes in gene expression induced by oxaliplatin alone. Notably, RNAseq analysis revealed MTM rescue of oxaliplatin-induced dysregulation of mitochondrial electron transport chain genes that correlated with in vivo reversal of excess reactive oxygen species in DRG neurons. This finding suggests that the mechanism(s) driving persistent pain states such as CIPN are not fixed but are sustained by ongoing modifiable transcription-dependent processes.

Detection of genes with differential expression dispersion unravels the role of autophagy in cancer progression

The majority of gene expression studies focus on the search for genes whose mean expression is different between two or more populations of samples in the so-called "differential expression analysis" approach. However, a difference in variance in gene expression may also be biologically and physiologically relevant. In the classical statistical model used to analyze RNA-sequencing (RNA-seq) data, the dispersion, which defines the variance, is only considered as a parameter to be estimated prior to identifying a difference in mean expression between conditions of interest. Here, we propose to evaluate four recently published methods, which detect differences in both the mean and dispersion in RNA-seq data. We thoroughly investigated the performance of these methods on simulated datasets and characterized parameter settings to reliably detect genes with a differential expression dispersion. We applied these methods to The Cancer Genome Atlas datasets. Interestingly, among the genes with an increased expression dispersion in tumors and without a change in mean expression, we identified some key cellular functions, most of which were related to catabolism and were overrepresented in most of the analyzed cancers. In particular, our results highlight autophagy, whose role in cancerogenesis is context-dependent, illustrating the potential of the differential dispersion approach to gain new insights into biological processes and to discover new biomarkers.

Perturbations in Neuroinflammatory Pathways Are Associated With a Worst Pain Profile in Oncology Patients Receiving Chemotherapy

Unrelieved pain occurs in 55% of cancer patients. Identification of molecular mechanisms for pain may provide insights into therapeutic targets. Purpose was to evaluate for perturbations in neuroinflammatory pathways between oncology patients with and without severe pain. Worst pain severity was rated using a 0 to 10 numeric rating scale six times over two cycles of chemotherapy. Latent profile analysis was used to identify subgroups of patients with distinct pain profiles. Pathway impact analyses were performed in two independent samples using gene expression data obtained from RNA sequencing (n = 192) and microarray (n = 197) technologies. Fisher's combined probability test was used to identify significantly perturbed pathways between None versus the Severe pain classes. In the RNA sequencing and microarray samples, 62.5% and 56.3% of patients were in the Severe pain class, respectively. Nine perturbed pathways were related to neuroinflammatory mechanisms (i.e., retrograde endocannabinoid signaling, gamma-aminobutyric acid synapse, glutamatergic synapse, Janus kinase-signal transducer and activator of transcription signaling, phagosome, complement and coagulation cascades, cytokine-cytokine receptor interaction, chemokine signaling, calcium signaling). First study to identify perturbations in neuroinflammatory pathways associated with severe pain in oncology outpatients. Findings suggest that complex neuroimmune interactions are involved in the maintenance of chronic pain conditions. Perspective: In this study that compared oncology patients with none versus severe pain, nine perturbed neuroinflammatory pathways were identified. Findings suggest that complex neuroimmune interactions are involved in the maintenance of persistent pain conditions.

Building an embryo: An auxin gene toolkit for zygotic and somatic embryogenesis in Brazilian pine

Many studies in the model species Arabidopsis thaliana characterized genes involved in embryo formation. However, much remains to be learned about the portfolio of genes that are involved in signal transduction and transcriptional regulation during plant embryo development in other species, particularly in an evolutionary context, especially considering that some genes involved in embryo patterning are not exclusive of land plants. This study, used a combination of domain architecture phylostratigraphy and phylogenetic reconstruction to investigate the evolutionary history of embryo patterning and auxin metabolism (EPAM) genes in Viridiplantae. This approach shed light on the co-optation of auxin metabolism and other molecular mechanisms that contributed to the radiation of land plants, and specifically to embryo formation. These results have potential to assist conservation programs, by directing the development of tools for obtaining somatic embryos. In this context, we employed this methodology with critically endangered and non-model species Araucaria angustifolia, the Brazilian pine, which is current focus of conservation efforts using somatic embryogenesis. So far, this approach had little success since somatic embryos fail to completely develop. By profiling the expression of genes that we identified as necessary for the emergence of land-plant embryos, we found striking differences between zygotic and somatic embryos that might explain the developmental arrest and be used to improve A. angustifolia somatic culture.

Cancer-related cognitive impairment is associated with perturbations in inflammatory pathways

Cancer-related cognitive impairment (CRCI) is a significant problem for patients receiving chemotherapy. While a growing amount of pre-clinical and clinical evidence suggests that inflammatory mechanisms underlie CRCI, no clinical studies have evaluated for associations between CRCI and changes in gene expression. Therefore, the purpose of this study was to evaluate for differentially expressed genes and perturbed inflammatory pathways across two independent samples of patients with cancer who did and did not report CRCI. The Attentional Function Index (AFI) was the self-report measure used to assess CRCI. AFI scores of <5 and of >7.5 indicate low versus high levels of cognitive function, respectively. Of the 185 patients in Sample 1, 49.2% had an AFI score of <5 and 50.8% had an AFI score of >7.5. Of the 158 patients in Sample 2, 50.6% had an AFI score of <5 and 49.4% had an AFI score of >7.5. Data from 182 patients in Sample 1 were analyzed using RNA-seq. Data from 158 patients in Sample 2 were analyzed using microarray. Twelve KEGG signaling pathways were significantly perturbed between the AFI groups, five of which were signaling pathways related to inflammatory mechanisms (e.g., cytokine-cytokine receptor interaction, tumor necrosis factor signaling). This study is the first to describe perturbations in inflammatory pathways associated with CRCI. Findings highlight the role of cytokines both in terms of cytokine-specific pathways, as well as pathways involved in cytokine production and cytokine activation. These findings have the potential to identify new targets for therapeutics and lead to the development of interventions to improve cognition in patients with cancer.

Alterations in Patterns of Gene Expression and Perturbed Pathways in the Gut-Brain Axis Are Associated With Chemotherapy-Induced Nausea

CONTEXT

Despite current advances in antiemetic treatments, approximately 50% of oncology patients experience chemotherapy-induced nausea (CIN).

OBJECTIVES

The purpose of this study was to evaluate for differentially expressed genes and perturbed pathways associated with the gut-brain axis (GBA) across two independent samples of oncology patients who did and did not experience CIN.

METHODS

Oncology patients (n = 735) completed study questionnaires in the week before their second or third cycle of chemotherapy. CIN occurrence was assessed using the Memorial Symptom Assessment Scale. Gene expression analyses were performed in two independent samples using ribonucleic acid sequencing (Sample 1, n = 357) and microarray (Sample 2, n = 352) methodologies. Fisher's combined probability method was used to determine genes that were differentially expressed and pathways that were perturbed between the two nausea groups across both samples.

RESULTS

CIN was reported by 63.6% of the patients in Sample 1 and 48.9% of the patients in Sample 2. Across the two samples, 703 genes were differentially expressed, and 37 pathways were found to be perturbed between the two CIN groups. We identified nine perturbed pathways that are involved in mechanisms associated with alterations in the GBA (i.e., mucosal inflammation, disruption of gut microbiome).

CONCLUSION

Persistent CIN remains a significant clinical problem. Our study is the first to identify novel GBA-related pathways associated with the occurrence of CIN. Our findings warrant confirmation and suggest directions for future clinical studies to decrease CIN occurrence.

Genomic profiling of Nipah virus using NGS driven RNA-Seq expression data

Nipah virus (NiV) is an ssRNA, enveloped paramyxovirus in the genus Henipaveridae with a case fatality rate >70%. We analyzed the NGS RNA-Seq gene expression data of NiV to detect differentially expressed genes (DEGs) using the statistical R package limma. We used the Cytoscape, Ensembl, and STRING tools to construct the gene-gene interaction tree, phylogenetic gene tree and protein-protein interaction networks towards functional annotation. We identified 2707 DEGs (p-value <0.05) among 54359 NiV genes. The top-up and down-regulated DEGs were EPST1, MX1, IFIT3, RSAD2, OAS1, OASL, CMPK2 and SLFN13, SPAC977.17 using log2FC criteria with optimum threshold 1.0. The top 20 up-regulated gene-gene interaction trees showed no significant association between Nipah and Tularemia virus. Similarly, the top 20 down-regulated genes of neither Ebola nor Tularemia virus showed an association with the Nipah virus. Hence, we document the top-up and down-regulated DEGs for further consideration as biomarkers and candidates for vaccine or drug design against Nipah virus to combat infection.

Transcriptomic response to aquaculture intensification in Nile tilapia

To meet future global demand for fish protein, more fish will need to be farmed using fewer resources, and this will require the selection of nonaggressive individuals that perform well at high densities. Yet, the genetic changes underlying loss of aggression and adaptation to crowding during aquaculture intensification are largely unknown. We examined the transcriptomic response to aggression and crowding in Nile tilapia, one of the oldest and most widespread farmed fish, whose social structure shifts from social hierarchies to shoaling with increasing density. A mirror test was used to quantify aggression and skin darkening (a proxy for stress) of fish reared at low and high densities, and gene expression in the hypothalamus was analysed among the most and least aggressive fish at each density. Fish reared at high density were darker, had larger brains, were less active and less aggressive than those reared at low density and had differentially expressed genes consistent with a reactive stress-coping style and activation of the hypothalamus-pituitary-interrenal (HPI) axis. Differences in gene expression among aggressive fish were accounted for by density and the interaction between density and aggression levels, whereas for nonaggressive fish differences in gene expression were associated with individual variation in skin brightness and social stress. Thus, the response to crowding in Nile tilapia is context dependent and involves different neuroendocrine pathways, depending on social status. Knowledge of genes associated with the response to crowding may pave the way for more efficient fish domestication, based on the selection of nonaggressive individuals with increasing tolerance to chronic stress necessary for aquaculture intensification.

DREAMSeq: An Improved Method for Analyzing Differentially Expressed Genes in RNA-seq Data

RNA sequencing (RNA-seq) has become a widely used technology for analyzing global gene-expression changes during certain biological processes. It is generally acknowledged that RNA-seq data displays equidispersion and overdispersion characteristics; therefore, most RNA-seq analysis methods were developed based on a negative binomial model capable of capturing both equidispersed and overdispersed data. In this study, we reported that in addition to equidispersion and overdispersion, RNA-seq data also displays underdispersion characteristics that cannot be adequately captured by general RNA-seq analysis methods. Based on a double Poisson model capable of capturing all data characteristics, we developed a new RNA-seq analysis method (DREAMSeq). Comparison of DREAMSeq with five other frequently used RNA-seq analysis methods using simulated datasets showed that its performance was comparable to or exceeded that of other methods in terms of type I error rate, statistical power, receiver operating characteristics (ROC) curve, area under the ROC curve, precision-recall curve, and the ability to detect the number of differentially expressed genes, especially in situations involving underdispersion. These results were validated by quantitative real-time polymerase chain reaction using a real Foxtail dataset. Our findings demonstrated DREAMSeq as a reliable, robust, and powerful new method for RNA-seq data mining. The DREAMSeq R package is available at http://tanglab.hebtu.edu.cn/tanglab/Home/DREAMSeq.

Expression of mitochondrial dysfunction-related genes and pathways in paclitaxel-induced peripheral neuropathy in breast cancer survivors

Background

Paclitaxel is one of the most commonly used drugs to treat breast cancer. Its major dose-limiting toxicity is paclitaxel-induced peripheral neuropathy (PIPN). PIPN persists into survivorship and has a negative impact on patient’s mood, functional status, and quality of life. No interventions are available to treat PIPN. A critical barrier to the development of efficacious interventions is the lack of understanding of the mechanisms that underlie PIPN. Mitochondrial dysfunction has been evaluated in preclinical studies as a hypothesized mechanism for PIPN, but clinical data to support this hypothesis are limited. The purpose of this pilot study was to evaluate for differential gene expression and perturbed pathways between breast cancer survivors with and without PIPN.

Methods

Gene expression in peripheral blood was assayed using RNA-seq. Differentially expressed genes (DEG) and pathways associated with mitochondrial dysfunction were identified between survivors who received paclitaxel and did (n = 25) and did not (n = 25) develop PIPN.

Results

Breast cancer survivors with PIPN were significantly older; more likely to be unemployed; reported lower alcohol use; had a higher body mass index and poorer functional status; and had a higher number of lower extremity sites with loss of light touch, cold, and pain sensations and higher vibration thresholds. No between-group differences were found in the cumulative dose of paclitaxel received or in the percentage of patients who had a dose reduction or delay due to PIPN. Five DEGs and nine perturbed pathways were associated with mitochondrial dysfunction related to oxidative stress, iron homeostasis, mitochondrial fission, apoptosis, and autophagy.

Conclusions

This study is the first to provide molecular evidence that a number of mitochondrial dysfunction mechanisms identified in preclinical models of various types of neuropathic pain including chemotherapy-induced peripheral neuropathy are found in breast cancer survivors with persistent PIPN and suggest genes for validation and as potential therapeutic targets.

Fully Bayesian analysis of RNA-seq counts for the detection of gene expression heterosis

Heterosis, or hybrid vigor, is the enhancement of the phenotype of hybrid progeny relative to their inbred parents. Heterosis is extensively used in agriculture, and the underlying mechanisms are unclear. To investigate the molecular basis of phenotypic heterosis, researchers search tens of thousands of genes for heterosis with respect to expression in the transcriptome. Difficulty arises in the assessment of heterosis due to composite null hypotheses and non-uniform distributions for p-values under these null hypotheses. Thus, we develop a general hierarchical model for count data and a fully Bayesian analysis in which an efficient parallelized Markov chain Monte Carlo algorithm ameliorates the computational burden. We use our method to detect gene expression heterosis in a two-hybrid plant-breeding scenario, both in a real RNA-seq maize dataset and in simulation studies. In the simulation studies, we show our method has well-calibrated posterior probabilities and credible intervals when the model assumed in analysis matches the model used to simulate the data. Although model misspecification can adversely affect calibration, the methodology is still able to accurately rank genes. Finally, we show that hyperparameter posteriors are extremely narrow and an empirical Bayes (eBayes) approach based on posterior means from the fully Bayesian analysis provides virtually equivalent posterior probabilities, credible intervals, and gene rankings relative to the fully Bayesian solution. This evidence of equivalence provides support for the use of eBayes procedures in RNA-seq data analysis if accurate hyperparameter estimates can be obtained.

Recent stimulant use and leukocyte gene expression in methamphetamine users with treated HIV infection

Stimulant use may accelerate HIV disease progression through biological and behavioral pathways. However, scant research with treated HIV-positive persons has examined stimulant-associated alterations in pathophysiologic processes relevant to HIV pathogenesis. In a sample of 55 HIV-positive, methamphetamine-using sexual minority men with a viral load less than 200 copies/mL, we conducted RNA sequencing to examine patterns of leukocyte gene expression in participants who had a urine sample that was reactive for stimulants (n = 27) as compared to those who tested non-reactive (n = 28). Results indicated differential expression of 32 genes and perturbation of 168 pathways in recent stimulant users. We observed statistically significant differential expression of single genes previously associated with HIV latency, cell cycle regulation, and immune activation in recent stimulant users (false discovery rate p < 0.10). Pathway analyses indicated enrichment for genes associated with inflammation, innate immune activation, neuroendocrine hormone regulation, and neurotransmitter synthesis. Recent stimulant users displayed concurrent elevations in plasma levels of tumor necrosis factor - alpha (TNF-α) but not interleukin 6 (IL-6). Further research is needed to examine the bio-behavioral mechanisms whereby stimulant use may contribute to HIV persistence and disease progression.

Minority stress and leukocyte gene expression in sexual minority men living with treated HIV infection

Sexual minority (i.e., non-heterosexual) individuals experience poorer mental and physical health, accounted for in part by the additional burden of sexual minority stress occurring from being situated in a culture favoring heteronormativity. Informed by previous research, the purpose of this study was to identify the relationship between sexual minority stress and leukocyte gene expression related to inflammation, cancer, immune function, and cardiovascular function. Sexual minority men living with HIV who were on anti-retroviral medication, had viral load < 200 copies/mL, and had biologically confirmed, recent methamphetamine use completed minority stress measures and submitted blood samples for RNA sequencing on leukocytes. Differential gene expression and pathway analyses were conducted comparing those with clinically elevated minority stress (n = 18) and those who did not meet the clinical cutoff (n = 20), covarying reactive urine toxicology results for very recent stimulant use. In total, 90 differentially expressed genes and 138 gene set pathways evidencing 2-directional perturbation were observed at false discovery rate (FDR) < 0.10. Of these, 41 of the differentially expressed genes and 35 of the 2-directionally perturbed pathways were identified as functionally related to hypothesized mechanisms of inflammation, cancer, immune function, and cardiovascular function. The neuroactive-ligand receptor pathway (implicated in cancer development) was identified using signaling pathway impact analysis. Our results suggest several potential biological pathways for future work investigating the relationship between sexual minority stress and health.

Biomarker Identification from RNA-Seq Data using a Robust Statistical Approach

Biomarker identification by differentially expressed genes (DEGs) using RNA-sequencing technology is an important task to characterize the transcriptomics data. This is possible with the advancement of next-generation sequencing technology (NGS). There are a number of statistical techniques to identify DEGs from high-dimensional RNA-seq count data with different groups or conditions such as edgeR, SAMSeq, voom-limma, etc. However, these methods produce high false positives and low accuracy in presence of outliers. We describe a robust t-statistic method to overcome these drawbacks using both simulated and real RNA-seq datasets. The model performance with 61.2%, 35.2%, 21.6%, 6.9%, 74.5%, 78.4%, 93.1%, 35.2% sensitivity, specificity, MER, FDR, AUC, ACC, PPV, and NPV, respectively at 20% outliers is reported. We identified 409 DE genes with p-values<0.05 using robust t-test in HIV viremic vs avirmeic state real dataset. There are 28 up-regulated genes and 381 down-regulated genes estimated by log2 fold change (FC) approach at threshold value 1.5. The up-regulated genes form three clusters and it is found that 11 genes are highly associated in HIV- 1/AIDS. Protein-protein interaction (PPI) of up-regulated genes using STRING database found 21 genes with strong association among themselves. Thus, the identification of potential biomarkers from RNA-seq dataset using a robust t-statistical model is demonstrated.

RNA-Seq Analysis of IL-1B and IL-36 Responses in Epidermal Keratinocytes Identifies a Shared MyD88-Dependent Gene Signature

IL-36 cytokines have recently emerged as mediators of inflammation in autoimmune conditions including psoriasis vulgaris (PsV) and generalized pustular psoriasis (GPP). This study used RNA-seq to profile the transcriptome of primary epidermal keratinocytes (KCs) treated with IL-1B, IL-36A, IL-36B, or IL-36G. We identified some early IL-1B-specific responses (8 h posttreatment), but nearly all late IL-1B responses were replicated by IL-36 cytokines (24 h posttreatment). Type I and II interferon genes exhibited time-dependent response patterns, with early induction (8 h) followed by no response or repression (24 h). Altogether, we identified 225 differentially expressed genes (DEGs) with shared responses to all 4 cytokines at both time points (8 and 24 h). These involved upregulation of ligands (IL1A, IL1B, and IL36G) and activating proteases (CTSS) but also upregulation of inhibitors such as IL1RN and IL36RN. Shared IL-1B/IL-36 DEGs overlapped significantly with genes altered in PsV and GPP skin lesions, as well as genes near GWAS loci linked to autoimmune and autoinflammatory diseases (e.g., PsV, psoriatic arthritis, inflammatory bowel disease, and primary biliary cholangitis). Inactivation of MyD88 adapter protein using CRISPR/Cas9 completely abolished expression responses of such DEGs to IL-1B and IL-36G stimulation. These results provide a global view of IL-1B and IL-36 expression responses in epidermal KCs with fine-scale characterization of time-dependent and cytokine-specific response patterns. Our findings support an important role for IL-1B and IL-36 in autoimmune or autoinflammatory conditions and show that MyD88 adaptor protein mediates shared IL-1B/IL-36 responses.

RNA-seq identifies a diminished differentiation gene signature in primary monolayer keratinocytes grown from lesional and uninvolved psoriatic skin

Keratinocyte (KC) hyper-proliferation and epidermal thickening are characteristic features of psoriasis lesions, but the specific contributions of KCs to plaque formation are not fully understood. This study used RNA-seq to investigate the transcriptome of primary monolayer KC cultures grown from lesional (PP) and non-lesional (PN) biopsies of psoriasis patients and control subjects (NN). Whole skin biopsies from the same subjects were evaluated concurrently. RNA-seq analysis of whole skin identified a larger number of psoriasis-increased differentially expressed genes (DEGs), but analysis of KC cultures identified more PP- and PN-decreased DEGs. These latter DEG sets overlapped more strongly with genes near loci identified by psoriasis genome-wide association studies and were enriched for genes associated with epidermal differentiation. Consistent with this, the frequency of AP-1 motifs was elevated in regions upstream of PN-KC-decreased DEGs. A subset of these genes belonged to the same co-expression module, mapped to the epidermal differentiation complex, and exhibited differentiation-dependent expression. These findings demonstrate a decreased differentiation gene signature in PP/PN-KCs that had not been identified by pre-genomic studies of patient-derived monolayers. This may reflect intrinsic defects limiting psoriatic KC differentiation capacity, which may contribute to compromised barrier function in normal-appearing uninvolved psoriatic skin.

Imiquimod has strain-dependent effects in mice and does not uniquely model human psoriasis

BACKGROUND

Imiquimod (IMQ) produces a cutaneous phenotype in mice frequently studied as an acute model of human psoriasis. Whether this phenotype depends on strain or sex has never been systematically investigated on a large scale. Such effects, however, could lead to conflicts among studies, while further impacting study outcomes and efforts to translate research findings.

METHODS

RNA-seq was used to evaluate the psoriasiform phenotype elicited by 6 days of Aldara (5% IMQ) treatment in both sexes of seven mouse strains (C57BL/6 J (B6), BALB/cJ, CD1, DBA/1 J, FVB/NJ, 129X1/SvJ, and MOLF/EiJ).

RESULTS

In most strains, IMQ altered gene expression in a manner consistent with human psoriasis, partly due to innate immune activation and decreased homeostatic gene expression. The response of MOLF males was aberrant, however, with decreased expression of differentiation-associated genes (elevated in other strains). Key aspects of the IMQ response differed between the two most commonly studied strains (BALB/c and B6). Compared with BALB/c, the B6 phenotype showed increased expression of genes associated with DNA replication, IL-17A stimulation, and activated CD8+ T cells, but decreased expression of genes associated with interferon signaling and CD4+ T cells. Although IMQ-induced expression shifts mirrored psoriasis, responses in BALB/c, 129/SvJ, DBA, and MOLF mice were more consistent with other human skin conditions (e.g., wounds or infections). IMQ responses in B6 mice were most consistent with human psoriasis and best replicated expression patterns specific to psoriasis lesions.

CONCLUSIONS

These findings demonstrate strain-dependent aspects of IMQ dermatitis in mice. We have shown that IMQ does not uniquely model psoriasis but in fact triggers a core set of pathways active in diverse skin diseases. Nonetheless, our findings suggest that B6 mice provide a better background than other strains for modeling psoriasis disease mechanisms.

Evaluation of logistic regression models and effect of covariates for case-control study in RNA-Seq analysis

Background

Next generation sequencing provides a count of RNA molecules in the form of short reads, yielding discrete, often highly non-normally distributed gene expression measurements. Although Negative Binomial (NB) regression has been generally accepted in the analysis of RNA sequencing (RNA-Seq) data, its appropriateness has not been exhaustively evaluated. We explore logistic regression as an alternative method for RNA-Seq studies designed to compare cases and controls, where disease status is modeled as a function of RNA-Seq reads using simulated and Huntington disease data. We evaluate the effect of adjusting for covariates that have an unknown relationship with gene expression. Finally, we incorporate the data adaptive method in order to compare false positive rates.

Results

When the sample size is small or the expression levels of a gene are highly dispersed, the NB regression shows inflated Type-I error rates but the Classical logistic and Bayes logistic (BL) regressions are conservative. Firth’s logistic (FL) regression performs well or is slightly conservative. Large sample size and low dispersion generally make Type-I error rates of all methods close to nominal alpha levels of 0.05 and 0.01. However, Type-I error rates are controlled after applying the data adaptive method. The NB, BL, and FL regressions gain increased power with large sample size, large log2 fold-change, and low dispersion. The FL regression has comparable power to NB regression.

Conclusions

We conclude that implementing the data adaptive method appropriately controls Type-I error rates in RNA-Seq analysis. Firth’s logistic regression provides a concise statistical inference process and reduces spurious associations from inaccurately estimated dispersion parameters in the negative binomial framework.

Electronic supplementary material

The online version of this article (doi:10.1186/s12859-017-1498-y) contains supplementary material, which is available to authorized users.

NBLDA: negative binomial linear discriminant analysis for RNA-Seq data

BACKGROUND

RNA-sequencing (RNA-Seq) has become a powerful technology to characterize gene expression profiles because it is more accurate and comprehensive than microarrays. Although statistical methods that have been developed for microarray data can be applied to RNA-Seq data, they are not ideal due to the discrete nature of RNA-Seq data. The Poisson distribution and negative binomial distribution are commonly used to model count data. Recently, Witten (Annals Appl Stat 5:2493-2518, 2011) proposed a Poisson linear discriminant analysis for RNA-Seq data. The Poisson assumption may not be as appropriate as the negative binomial distribution when biological replicates are available and in the presence of overdispersion (i.e., when the variance is larger than or equal to the mean). However, it is more complicated to model negative binomial variables because they involve a dispersion parameter that needs to be estimated.

RESULTS

In this paper, we propose a negative binomial linear discriminant analysis for RNA-Seq data. By Bayes' rule, we construct the classifier by fitting a negative binomial model, and propose some plug-in rules to estimate the unknown parameters in the classifier. The relationship between the negative binomial classifier and the Poisson classifier is explored, with a numerical investigation of the impact of dispersion on the discriminant score. Simulation results show the superiority of our proposed method. We also analyze two real RNA-Seq data sets to demonstrate the advantages of our method in real-world applications.

CONCLUSIONS

We have developed a new classifier using the negative binomial model for RNA-seq data classification. Our simulation results show that our proposed classifier has a better performance than existing works. The proposed classifier can serve as an effective tool for classifying RNA-seq data. Based on the comparison results, we have provided some guidelines for scientists to decide which method should be used in the discriminant analysis of RNA-Seq data. R code is available at http://www.comp.hkbu.edu.hk/~xwan/NBLDA.R or https://github.com/yangchadam/NBLDA.

Multiple functional linear model for association analysis of RNA-seq with imaging

Emerging integrative analysis of genomic and anatomical imaging data which has not been well developed, provides invaluable information for the holistic discovery of the genomic structure of disease and has the potential to open a new avenue for discovering novel disease susceptibility genes which cannot be identified if they are analyzed separately. A key issue to the success of imaging and genomic data analysis is how to reduce their dimensions. Most previous methods for imaging information extraction and RNA-seq data reduction do not explore imaging spatial information and often ignore gene expression variation at the genomic positional level. To overcome these limitations, we extend functional principle component analysis from one dimension to two dimensions (2DFPCA) for representing imaging data and develop a multiple functional linear model (MFLM) in which functional principal scores of images are taken as multiple quantitative traits and RNA-seq profile across a gene is taken as a function predictor for assessing the association of gene expression with images. The developed method has been applied to image and RNA-seq data of ovarian cancer and kidney renal clear cell carcinoma (KIRC) studies. We identified 24 and 84 genes whose expressions were associated with imaging variations in ovarian cancer and KIRC studies, respectively. Our results showed that many significantly associated genes with images were not differentially expressed, but revealed their morphological and metabolic functions. The results also demonstrated that the peaks of the estimated regression coefficient function in the MFLM often allowed the discovery of splicing sites and multiple isoforms of gene expressions.

The Impact of Normalization Methods on RNA-Seq Data Analysis

High-throughput sequencing technologies, such as the Illumina Hi-seq, are powerful new tools for investigating a wide range of biological and medical problems. Massive and complex data sets produced by the sequencers create a need for development of statistical and computational methods that can tackle the analysis and management of data. The data normalization is one of the most crucial steps of data processing and this process must be carefully considered as it has a profound effect on the results of the analysis. In this work, we focus on a comprehensive comparison of five normalization methods related to sequencing depth, widely used for transcriptome sequencing (RNA-seq) data, and their impact on the results of gene expression analysis. Based on this study, we suggest a universal workflow that can be applied for the selection of the optimal normalization procedure for any particular data set. The described workflow includes calculation of the bias and variance values for the control genes, sensitivity and specificity of the methods, and classification errors as well as generation of the diagnostic plots. Combining the above information facilitates the selection of the most appropriate normalization method for the studied data sets and determines which methods can be used interchangeably.

An Iterative Leave-One-Out Approach to Outlier Detection in RNA-Seq Data

The discrete data structure and large sequencing depth of RNA sequencing (RNA-seq) experiments can often generate outlier read counts in one or more RNA samples within a homogeneous group. Thus, how to identify and manage outlier observations in RNA-seq data is an emerging topic of interest. One of the main objectives in these research efforts is to develop statistical methodology that effectively balances the impact of outlier observations and achieves maximal power for statistical testing. To reach that goal, strengthening the accuracy of outlier detection is an important precursor. Current outlier detection algorithms for RNA-seq data are executed within a testing framework and may be sensitive to sparse data and heavy-tailed distributions. Therefore, we propose a univariate algorithm that utilizes a probabilistic approach to measure the deviation between an observation and the distribution generating the remaining data and implement it within in an iterative leave-one-out design strategy. Analyses of real and simulated RNA-seq data show that the proposed methodology has higher outlier detection rates for both non-normalized and normalized negative binomial distributed data.

The level of residual dispersion variation and the power of differential expression tests for RNA-Seq data

RNA-Sequencing (RNA-Seq) has been widely adopted for quantifying gene expression changes in comparative transcriptome analysis. For detecting differentially expressed genes, a variety of statistical methods based on the negative binomial (NB) distribution have been proposed. These methods differ in the ways they handle the NB nuisance parameters (i.e., the dispersion parameters associated with each gene) to save power, such as by using a dispersion model to exploit an apparent relationship between the dispersion parameter and the NB mean. Presumably, dispersion models with fewer parameters will result in greater power if the models are correct, but will produce misleading conclusions if not. This paper investigates this power and robustness trade-off by assessing rates of identifying true differential expression using the various methods under realistic assumptions about NB dispersion parameters. Our results indicate that the relative performances of the different methods are closely related to the level of dispersion variation unexplained by the dispersion model. We propose a simple statistic to quantify the level of residual dispersion variation from a fitted dispersion model and show that the magnitude of this statistic gives hints about whether and how much we can gain statistical power by a dispersion-modeling approach.

Integrative RNA-seq and microarray data analysis reveals GC content and gene length biases in the psoriasis transcriptome

Gene expression profiling of psoriasis has driven research advances and may soon provide the basis for clinical applications. For expression profiling studies, RNA-seq is now a competitive technology, but RNA-seq results may differ from those obtained by microarray. We therefore compared findings obtained by RNA-seq with those from eight microarray studies of psoriasis. RNA-seq and microarray datasets identified similar numbers of differentially expressed genes (DEGs), with certain genes uniquely identified by each technology. Correspondence between platforms and the balance of increased to decreased DEGs was influenced by mRNA abundance, GC content, and gene length. Weakly expressed genes, genes with low GC content, and long genes were all biased toward decreased expression in psoriasis lesions. The strength of these trends differed among array datasets, most likely due to variations in RNA quality. Gene length bias was by far the strongest trend and was evident in all datasets regardless of the expression profiling technology. The effect was due to differences between lesional and uninvolved skin with respect to the genome-wide correlation between gene length and gene expression, which was consistently more negative in psoriasis lesions. These findings demonstrate the complementary nature of RNA-seq and microarray technology and show that integrative analysis of both data types can provide a richer view of the transcriptome than strict reliance on a single method alone. Our results also highlight factors affecting correspondence between technologies, and we have established that gene length is a major determinant of differential expression in psoriasis lesions.