Comparative exomics of Phalariscultivars under salt stress

RoDEO evaluation on real human data

While simulation studies are usually effective and serve a useful purpose in evaluating the efficacy of algorithms, we found that in the context of DE detection, the RNA-seq simulators themselves are compelled to use some parametric model for transcript abundances. Also, on most public DE data sets, the gold standard is usually the result of some computational method. To the best of our knowledge, there is no public dataset for stable genes. To avoid any bias we used real human data, where the DE genes were experimentally validated by qRT-PCR. These RNA samples are well characterized and have been used in several benchmarking studies and DE method comparisons [1–3]. We evaluated the performance of RoDEO against parametric DE detectors that assume a negative binomial count distribution model (baySeq, edgeR, sSeq) and demonstrate that it outperforms the others on several scenarios (details in the Methods section). Moreover, RoDEO's parameter-free framework is very suitable for the diverse and challenging data sets of Phalaris.

Phalarisreference transcriptome construction

Total RNA was extracted from Phalaris arundinacea cv Venture from five different tissues: leaf and stem of mature (40d) plants under drought stress, stem of mature plant under waterlogging stress, shoot and roots of young (9d old) plants. The samples were pooled, converted into a normalized sequencing library, and sequenced on the Roche 454 FLX platform. In total 494, 477 reads with mean length 389bp were assembled into 18, 682 transcripts with the GS De Novo Assembler software from Roche (which has been designed for the analysis of 454 read data), using the default parameters for cDNA analysis and the adapter-trimmed raw data as input. The sequence data are available at the Sequence Read Archive at NCBI, accession number SRP045256.

All assembled transcriptome contigs were searched against the non-redundant (nr) database at the National Center for Biotechnology Information (NCBI) using BLAST [4]. Gene Ontology annotations for cellular component, biological and molecular processes were predicted for the nr derived BLAST hits using Blast2GO [5]. In addition, InterPro [6] scans were performed to predict protein domain structure of the assembled contigs and Kyoto Encyclopedia of Genes and Genomes (KEGG) [7] pathway assignments were estimated based on enzyme commission (EC) numbers. The distribution of top BLAST hits among annotated plant species are shown in Figure 1. The largest number of sequences match purple false brome (Brachypodium distachyon). The species list reflects the known phylogenetic affiliations of Phalaris, thus further validating its evolutionary relation to these species. Note that large-scale genomic or transcriptomic datasets have not been previously published on Phalaris.

Salt stress experiment

Seeds of three Phalaris arundinacea and one Phalaris aquatica genotypes were obtained from Lantmännen lantbruk (former Svalöf Weibull seed). SWN RF 9901 is a Norwegian accession from the Nordic Seed bank (NGB7147), SW RF 5008 is a SW breeding line containing European and North American material. Palaton is a OECD listed cultivar in Estonia, Finland, and United States. CPI 19315 is a Moroccan accession of P. aquatica, a species closely related to P. arundinacea [8]. Two cultivars (CPI 19315, SWN RF 9901) exhibit a negative response to stress, while the other two (Palaton 2, SW RF 5008) tolerate salt stress having a positive response. Figure 2 illustrates some of their phenotypes during salt stress. CPI and Palaton data were collected 55 days after starting the salt treatment, thus called early stress, while SWN and SW data were collected after 100 days salt stress, called late stress. At sampling time the plants' response to salt of was termed either positive, if the plants looked largely unaffected, or negative, if the plants were severely damaged.

RNA-seq libraries during salt stress and control libraries before salt stress were constructed using the Illumina TruSeq protocol and chemistry (version 2). Libraries were sequenced on Illumina HiSeq2000. The sequence data are available at the Sequence Read Archive at NCBI, accession number SRP045256. The reads were mapped with bwa [9] to the reference transcriptome. In each experiment, 44 - 62 million reads mapped to the reference, with fewer than 10% transcripts having no reads mapping to them.

See Additional File 1 for more details on the salt stress experiment and RNA sequencing.

Differential expression

We applied our differential expression method RoDEO on the Phalaris salt stress read counts per transcript, with parameters P = 15, I = 100 iterations, R equaling the number of mapped reads per experiment, and DE threshold dist(·, ·) = 1.0 and mode distance ≥ 5 (see Methods for details). Stable transcripts were identified by pairwise comparison of the re-sampled iterations, and taking the intersection of the stable sets from all the iterations (see Additional file 1 for details and results).

The general differential expression landscape for each cultivar is shown in Figure 3. Thousands of genes have the maximum norm distance 1.0 between their character functions in the control and stress sample, and those transcripts are further ranked by their mode distance to refine the DE candidates set to the most interesting transcripts, i.e., those with mode distance ≥ 5.

We call robust positive up-regulated transcripts those that are up-regulated in both positive response cultivars (similarly for down-regulation, and for negative response). These transcripts are more likely to be truly differentially expressed, since they are found DE in both cultivars having the same stress response. Figure 4 shows the overview of the number of transcripts that are up- or down-regulated in each cultivar.

We find 8 robust down-regulated and 19 robust up-regulated transcripts in negative response cultivars. In positive response cultivars, 2 out of 9 robust up-regulated transcripts are shared with the negative response. In addition, there are 2 down-regulated transcripts for positive response. Thus we have in total 10 down-regulated and 26 up-regulated response-specific robust DE transcripts. Annotations for these transcripts are included in Additional file 2.

The annotations for these robust response-specific DE transcripts show several proteins up-regulated under positive response which are widely associated with plant stress proteins, such as sulphur-rich thionin like protein, thiol protease, polyamine oxidase, universal stress protein a like protein, and maize proteinase inhibitor. Interestingly, these proteins do not show up in the plants with negative stress response, indicating that absence of these proteins could be correlated with a poor response to stress. Instead, enzymes from various catabolic processes are up-regulated which are not specifically active in stress response, such as amino acid or carbohydrate metabolism.

In addition, we find robust DE transcripts for early and late stress response. There are 4 and 36 down-regulated transcripts in early and late stress, respectively. Out of 29 for early and 23 for late up-regulated transcripts, 2 are shared between early and late. Thus we have in total 40 down-regulated and 50 up-regulated time-specific robust DE transcripts. Annotations for these transcripts are included in Additional file 2.

The functional distinction between early and late DE up-regulation is less clear than between positive and negative response. In both early and late categories several proteins associated with abiotic stress response appear, such as cbl-interacting protein kinase 1-like protein, metallothionein-like protein type 2, delta-1-pyrroline-5-carboxylate synthase-like (early up-regulated), and e.g. cysteine proteinase inhibitor, sulphur-rich thionin like protein, thiol protease (late up-regulated). There is little overlap between the two categories, indicating that the plant response changes between early and late (and more severe) phase of salt stress. The late phase of salt stress reveals proteins that are more similar to the DE proteins detected in the positive response plants that were better able to cope with salt stress.

Two transcripts are up-regulated in all cultivars' stress vs. control samples. One of them has a close match (89% similarity) to Cytochrome P450 71C4 protein in barley (the other has a weaker match to a hypothetical Aegilopsis tauschii grass protein). Some other members of the Cytochrome P450 gene family have previously been studied w.r.t. salt stress in Arabidopsis [10].

Cultivar comparative exomics

We performed a global comparative exomics analysis of the Phalaris cultivars based on the variations observed in the RNA-seq reads in the stress experiment. Additionally, we looked for SNP variations across cultivars in the candidate transcripts from the DE analysis. Possibly due to the lack of a reference genome, including non-coding regions and promoters, and because we are limited by the sequenced mRNA information, we did not observe any DE loci where both cultivars with positive response would have allele differences from both cultivars with negative response. Globally, we identified a handful of transcripts with such loci, but they did not exhibit differential expression.

We computed pairwise distances between the cultivars, based on 200k loci where at least one cultivar's RNA-seq reads differed from the reference transcriptome, and there were at least 100 reads covering the SNP. The dendograms of the pairwise distances (based on Hamming distances of the alleles per loci, and average linkage) are shown in Figure 5. The phylogenetic trees from using all transcripts, DE transcripts, or stable transcripts tell a similar story. Identifying Phalaris stable genes is discussed in Additional file 1. CPI (P. aquatica), originating from a distant geographical location, is also genetically the most distant from the three North European P. arundinacea genotypes. The tree constructed from the DE transcripts differs slightly from the other trees, differentiating Palaton from SWN and SW. It is encouraging that these inferences can be made based on only a partial reference transcriptome and RNA-seq reads from the cultivars.

The mathematical framework

Overview. The primary objective is to compare the expression level of genes in high-throughput sequencing (HTS) data between a pair of biological samples. Due to the un-targeted nature of HTS, it is only natural to assume that the observed expression values of the genes are not independent. Second, the technologies produce data in terms of integer read counts. This is in contrast to earlier technologies which produced real-valued intensity measures. Also the total number of genes being simultaneously assayed earlier were determined by the technology while in the sequencing approaches, these read counts are a result of the application of mapping software pipeline, at whose core reside some sophisticated (parametric) string matching algorithms. Such a setting of read count computation is not naturally conducive to normalization processes. Hence the read count of a gene g is simply divided by the total number of reads (in millions) in the assay and by the length of the gene g (in kilobases) to yield a quantity called the RPKM value of g. Instead of directly working with the expression value of g, we define a character function ${\varphi }_g$ for each gene g. The two most desirable properties of this function are (1) ${\varphi }_g$ depend on the expression values of all the other genes in the assay and (2) ${\varphi }_g$ are scale invariant. We also define another family of distributions, built on these functions (called ψ), to study the overall behavior of the changes captured by the two transcriptomes.

Notation

whose components, ${\varphi }_g$, termed character functions, are defined below. In particular, these functions are scale invariant.

Definition 1 (character functions)

a character function if there exists a ${\left({\delta }_p\right)}_{1\le p\le P}\in {ℝ}_{\ge 0}^P$, satisfying the following conditions:

1 if ${\varphi }_g\left(t\right)={\varphi }_{g^{\prime }}\left(t\right)=p$ then |t _g − t _g′ | ≤ δ _p , and

2 if ${\varphi }_g\left(t\right)<{\varphi }_{g^{\prime }}\left(t\right)$, then t _g < t _g′ ,

for all g, g′ ∈ G, and ${\left(t_g\right)}_{g\in G}\in {ℝ}_{\ge 0}^G$.

Notice that such a character function satisfies the following additional properties.

1 ${\varphi }_g$ is scale-invariant: ${\varphi }_g\left(ct\right)={\varphi }_g\left(t\right)$ for all c > 0 and $t\in {ℝ}_{\ge 0}^G$.

This follows from the first condition of the definition of ${\varphi }_g$.

2 For all triplets, g, g′, g″, if t _g <t _g′ <t _g″ with ${\varphi }_g\left(t\right)={\varphi }_{g^{″}}\left(t\right)$ then ${\varphi }_g\left(t\right)={\varphi }_{g^{\prime }}\left(t\right)={\varphi }_{g^{″}}\left(t\right)$.

This property follows from the second condition of the definition. In fact this property leads to the algorithms for actually computing the scale-invariant maps.

whose components ${\psi }_p\left({\text{Φ}}^b|{\text{Φ}}^a\right)$, termed dispersion functions, are defined below.

Definition 2 (dispersion function)

where q ∈ {1, 2, .., P}.

This intrinsic measure can be used to evaluate the accuracy of the estimation of distributions of random variable ${\text{Φ}}^a$, based on single or multiple replicates. In practice, we have used ${\widetilde{\mathsf{\text{N}}}}_p\left({\text{Φ}}^a\right)$ to evaluate the character functions $\varphi$ (and thus estimate parameter P as well).

Estimating distributions of random variables $\text{Φ}\left(\mathcal{T}\right)$ and $\text{Ψ}\left(\text{Φ}\right)$

Critical to the computations of $\text{Φ}$ and $\text{Ψ}$ functions, is the estimation of the distribution of the random variables X _g of Eqn (1). Our model is influenced by the observation that the total number of genes (L) in the RNA-Seq experiment is based on the biological sample and not the technology and that this also renders the observed abundances of each gene to be not quite independent of each other. This is reflected in the details of the computations discussed here.

Let $c=\left(c_{g_1},c_{g_2},\dots ,c_{g_L}\right)$ be the observed abundances in the experiment. Without loss of generality, let c _g be the number of reads (read-counts) of gene g, and Σ _g∈G c _g = M. To estimate the distribution we use a Bernoulli process with parameters c _g /M, for g ∈ G, and M. This is a natural extension of the method for the single replicate. Let ${\left(c_g^1\right)}_{g\in G},{\left(c_g^2\right)}_{g\in G},\dots ,{\left(c_g^K\right)}_{g\in G}$ be the observed abundances in the K experiments under identical conditions. As before, let $c_g^k$ be the number of reads (read-counts) of gene g in the kth replicate, for 1 ≤ k ≤ K. Let $M^k={\sum }_g{c}_g^k$. We use K Bernoulli processes each with parameters $c_g^k/M^k$, for g ∈ G, and M ^k , as in the single replicate case. Let $c_g={\sum }_{k=1}^K{c}_g^k$ and let $M={\sum }_g{c}_g$. We estimate the probability, p _g , of gene g to be c _g /M.

Thus, modeling the experiment, with K ≥ 1 replicates, as a Bernoulli process, $\mathcal{T}={\left(X_g\right)}_{g\in G}$ of Eqn (1) follows a multinomial distribution with parameters M and (p _g ) _g∈G . A closed form of the distribution of ${\varphi }_g$ or ${\psi }_g$ is not straightforward to estimate and moreover, the argument to ${\varphi }_g$ is an L-tuple $t\in {ℝ}_{\ge 0}^G$. Hence we simulate the Bernoulli process in our implementation. Then based on Defn (1), for a fixed P, ${\varphi }_g$ is evaluated for each g using the L-tuple, t _i , of trial i. Similarly, based on Defn (2), ${\psi }_p$ is evaluated for each p using the L-tuples of the two experiments. The details are discussed below.

Simulating the character functions ${\varphi }_g$

Simulating the Bernoulli process for one experiment with n genes is described below, for some number of iterations I (e.g., I = 10³), and some number of reads R (e.g., R = 10⁶), using a fixed P (e.g., P = 20). Consider iteration i:

1 Repeat for each read r ∈ R: Assign the read to gene g ∈ G with probability $p_g=\mathsf{\text{log}}{c}_g/{\sum }_{i=1}^n\mathsf{\text{log}}{c}_i$.

2 Order genes according to decreasing assigned read count. Assign ${\varphi }_g^i=1$ for genes g with no reads.

3 Apply linear (regression) segmentation with P − 1 segments to the cumulative sum of the ordered read counts; assign ${\varphi }_g^i=P$ for genes g in the first segment, P − 1 for the second segment, and so forth. Least-squares linear segmentation can be performed in time O(n²P), see [12]. In practice when n is large one may want to subsample the genes, e.g., only use every 10th gene for the segmentation.

The values ${\varphi }_g^i$ are stored in each iteration i and together make up the distribution ${\varphi }_g$. The values of the dispersion function ψ _p are computed from ${\varphi }_g$. In practice, to compute $Pr\left({\Psi }_p\right({\text{Φ}}^b\left|{\text{Φ}}^a\right)=q)$, one counts the number of pairs of iterations (i, i′) ∈ (1...I, 1...I) where ${\varphi }_g^{a_i}=p$ and ${\varphi }_g^{b_{i^{\prime }}}=q$, where ${\varphi }_g^{a_i}$ denotes the bin of gene g in iteration i in dataset a. The probability Pr is the frequency of pairs $\left({\varphi }_g^{a_i}=p,{\varphi }_g^{b_{i^{\prime }}}=q\right)$ among all the pairs $\left({\varphi }_g^{a_i}=p,{\varphi }_g^{b_{i^{\prime }}}=p^{\prime }\right)|p^{\prime }\in 1\dots P$.

An example We compute the different distributions for a dataset with 10,000 genes and their (observed) abundances c. Figure 6 (a) shows ${\varphi }_g$ for all the genes g (as defined in Eqn 3) and Figure 6 (b) shows the ${\Psi }_p\left({\text{Φ}}^b\right|{\text{Φ}}^a)$ and ${\Psi }_p\left({\text{Φ}}^a\right|{\text{Φ}}^b)$, for each p (according to Eqn 5), while (c) summarizes the Ψ distributions as ${\widetilde{\mathsf{\text{O}}}}_p\left(\text{Φ}\right)$, ${\widetilde{\mathsf{\text{U}}}}_p\left(\text{Φ}\right)$, and ${\widetilde{\mathsf{\text{N}}}}_p\left(\text{Φ}\right)$ (according to Eqns 6-8).

Comparison across a pair ${\mathcal{T}}^a$ and ${\mathcal{T}}^b$

Instead of directly studying ${\mathcal{T}}^a$ and ${\mathcal{T}}^b$, we observe instead the behavior of the scale-invariant maps $\text{Φ}\left({\mathcal{T}}^a\right)$ and $\text{Φ}\left({\mathcal{T}}^b\right)$. We also distinguish two classes of genes: one whose members are differentially expressed (DE) the others whose members are stable across the two samples a and b with sufficient confidence (based on a threshold on some measure). The definition of stable genes in this context, and detecting them via the combinatorial problem of longest increasing subsequence is discussed in Additional file 1. Unfortunately, these two which appear complementary in their characteristics do not exactly partition the genes G into two because of the need to meet the confidence thresholds.

Differentially expressed (DE) genes

Character functions ${\varphi }_g\left({\mathcal{T}}^a\right)$ and ${\varphi }_g\left({\mathcal{T}}^b\right)$ are used to estimate the DE genes. We define a real-valued distance, dis(·, ·) ≥ 0, between the two functions as discussed below. For some fixed δ > 0, if $dis\left({\varphi }_g\left({\mathcal{T}}^a\right),{\varphi }_g\left({\mathcal{T}}^b\right)\right)>\delta$, then g is estimated to be differentially expressed. It is important to note that the mathematical framework discussed here is agnostic to the underlying models in the diverse read count simulators prevalent in literature.

To define the amount of differential expression between datasets a and b, we rank each gene according to the differences in their distributions ${\varphi }_g\left({\mathcal{T}}^a\right)$ and ${\varphi }_g\left({\mathcal{T}}^b\right)$. The ranking process is described below.

The genes are ordered by decreasing maximum norm, then by decreasing mode distance, and the genes at the top of the list are the most likely DE candidates.

Comparison with existing methods

We compare RoDEO against existing methods baySeq [13], edgeR [14], and sSeq [15] for identifying differentially expressed genes. In an effort to avoid biases often associated with generating simulated data using the exact distribution assumed by the methods (negative binomial in this case), we only used real data. Specifically, we used read counts and quantitative expression data from the MAQC human tissue samples, for which DE genes have been experimentally validated by qRT-PCR [1, 2]. The aim here is to detect differentially expressed genes between human reference RNA and brain RNA samples.

The methods were run on all 33k genes for which there were reads sequenced with Ion Torrent Proton technology (MAQC PRO dataset), while the DE results were evaluated only on those 17k genes with available qRT-PCR quantification. We used log-fold change cutoff > 1.5 to define DE genes and < 1.0 to define non-DE genes (other genes remained unassigned). Each method was run with default parameters (if any), RoDEO was run with P = 20 bins, I = 100 iterations, and R = 10⁷ reads. Each compared method provides a list of genes ordered by extent of their DE, and this list is used in the evaluation.

Figure 7 shows the DE methods' performance on the this MAQC PRO data consisting of two treatments, each with three replicates. Results on a single-replicate case using the first replicate only are shown in Figure 8. The four panels illustrate (i) For the number of top DE candidates on the x-axis, how many are false positives, (ii) False positive rate vs. true positive rate, including AUC (Area Under the Curve) measurement, (iii) For the top up to x = 500 true DE genes, how many are included in the method's DE list of size x (ideally x = y), (iv) AUC when varying the threshold for calling DE genes from the qRT-PCR results (threshold 1.5 is used in panels i-iii).

RoDEO outperforms all the other DE methods on this dataset, based on the AUC measurement at varying levels of the PCR DE threshold (panel iv), both on the multiple replicate and single replicate cases.