Inference and Sample Size Calculations Based on Statistical Tests in a Negative Binomial Distribution for Differential Gene Expression in RNA-seq Data

The high throughput RNA sequencing (RNA-seq) technology has become the popular method of choice for transcriptomics and the detection of differentially expressed genes. Sample size calculations for RNA-seq experimental design are an important consideration in biological research and clinical trials. Currently, the sample size formulas derived from the Wald and the likelihood ratio statistical tests with a Poisson distribution to model RNA-seq data have been developed. However, since the mean read counts in the real RNA-seq data are not equal to the variance, an extended method to calculate sample sizes based on a negative binomial distribution using an exact test statistic was proposed by Li et al. in 2013. In this study, we alternatively derive five sample size calculation methods based on the negative binomial distribution using the Wald test, the log-transformed Wald test and the log-likelihood ratio test statistics. A comparison of our five methods and an existing method was performed by calculating the sample sizes and the simulated power in different scenarios. We first calculated the sample sizes for testing a single gene using the six methods given a nominal significance level α at 0.05 and 80% power. Then, we calculated the sample sizes for testing multiple genes given a false discovery rate (FDR) at 0.05 and 0.10. The empirical power and true prognostic genes for differential gene expression analysis corresponding to the estimated sample sizes from the six methods are also estimated via the simulation studies. Using the sample size formulas derived from log-transformed and Wald-based tests, we observed smaller sample properties while maintaining the nominal power close to or higher than 80% in all the settings compared to other methods. Moreover, the Wald test based sample size calculation method is easier to compute and faster in an RNA-seq experimental design. Citation: Li X, Cooper NGF, Shyr Y, Wu D, Rouchka EC, et al. (2017) Inference and Sample Size Calculations Based on Statistical Tests in a Negative Binomial Distribution for Differential Gene Expression in RNA-seq Data. J Biom Biostat 8: 332. doi:10.4172/2155-6180.1000332


Introduction
Sample size calculations are a prerequisite in an experimental design for biological research and clinical trials. Recently, highthroughput RNA sequencing (RNA-seq) technology has been widely used for gene expression studies in a variety of applications, such as expression profiling of mRNAs or non-coding RNA [1][2][3], de novo assembly and characterization of transcriptomes [4,5], and the identification of novel alternatively spliced transcript [6,7]. These novel transcripts or differentially expressed genes (DEGs) identified from RNA-seq data may serve as human disease biomarkers or gene signatures for the clinical diagnosis [8][9][10]. With the rapid growth of RNA-seq applications, sample size calculation methods derived from test statistics with an appropriate distribution are important issues to be explored and discussed.
Due to the initial high cost of RNA-seq, sample size, in terms of the number of biological replicates, was not seriously considered as part of the experimental design. As a case in point, one RNA-seq review article [11] documented several RNA-seq studies showing that many had only one or a few biological replicates. While thousands of DEGs were identified within these studies, the lack of biological replicates leads to an absence of knowledge concerning biological variations and may result in a high percentage of false positive genes. Therefore, ignorance of biological variation is the fundamental problem with the analyzed results collected from un-replicated data. A recent paper [12] was the first to point out that conclusions drawn from un-replicated samples can be misleading and unrealistic. Later, another study [13] further addressed design and validation issues due to the lack of biological replicates in RNA-seq data.
One of the key questions in an experimental design is to determine the number of biological replicates needed for differential expression analysis in order to achieve a desired statistical power given a significance level α and an underlying distribution. Since RNA-seq data are read counts, a Poisson distribution is commonly used as the model for identifying DEGs in RNA-seq data [14,15]. Fang and Cui were the first one to derive the sample size calculations based on a Wald statistical test with a Poisson distribution for single gene in RNA-seq.
Like microarray data, an RNA-seq dataset contains thousands of genes to be tested simultaneously and independently for differential expression analysis. A method for the adjustment or correction of p-values is required to control the type I error rate when multiple pairwise comparisons are performed. Instead of setting the critical value α at 0.05 or 0.01 for significance, a much lower critical value α* is required to correct for the inflation of α. The most common method to control the family-wise error rate (FWER) is the Bonferroni correction in which the adjusted p-value is computed via dividing the critical p-value by the total number of comparisons being made. The other widely used method for this multiple correction problem is an FDR correction [20]. Since the Bonferroni correction with a large number of tests is more conservative than the FDR correction, using the Bonferroni correction results in a cost of increasing the probability of producing false negatives and consequently reducing the statistical power. For high dimensional microarray data analysis, an extension of the FDR correction was proposed [21] and is widely used by many researchers. To address similar issues for high-dimensional RNA-seq data, a sample size determination based on the extended FDR correction from microarray data analysis was further proposed [16].
Our study is motivated by exploring sample size calculations using the well-known test statistics (the Wald test and LRT) and a negative binomial distribution to model RNA-seq data. In Section 2, we first define the Wald test, the log-transformed Wald test and the LRT statistics using the negative binomial distribution to model a single gene in RNA-seq data. Then, we derive sample size calculations based on these defined statistical tests. Lastly, we derived sample size calculation methods for testing multiple genes while controlling the FDR [16]. In Section 3, we simulated power for testing single gene and multiple genes corresponding to the sample sizes estimated from our proposed methods and an existing method. The performance of these six methods is compared and evaluated via the required sample sizes and the estimated power. An application of real RNA-seq data to illustrate sample size calculations is presented in Section 4. Finally, we end with a discussion and conclusion in Section 5.

Derivation of test statistics
In an RNA-seq experiment, the data contains thousands of genes (g=1,…,G) with different number of reads for each sample mapped to the reference genome. Since the total number of reads among samples is different, the distribution of the gene in the sample with the same condition is not identical. A normalization factor called the size factor is used to model RNA-seq data with a negative binomial distribution. For simplicity, the following statistical tests and consequent sample size calculation methods are based on a single gene tests for DEG analysis.
For a single gene in RNA-seq data, suppose that, for each condition i (i=0,1), the observations X ij (j=1,…,n i ), are independent and identically follow a negative binomial distribution as X ij ~ NB (s ij γ i , ϕ) [17,22]. Under this setting, γ i is the true gene expression level in condition i, s ij is a size factor to normalize the raw read for the total number of reads mapped in the sample j, and ϕ is a dispersion parameter with the assumption ϕ i = ϕ.
Thus, the summation of reads per gene per condition For detection of a differentially expressed gene from RNA-seq data, the ratio γ 1 /γ 0 typically represents the fold change (ρ=γ 1 /γ 0 ). If the fold change equals one, we can say that this gene is not differentially expressed. Therefore, we are interested in making an inference about the ratio using the Wald statistics and the likelihood ratio methods for sample size calculations.
For testing the hypothesis about the fold change in regards to the DEGs, it is equivalent to test the hypothesis, Since γ i (i=0,1) and ф are unknown parameters, we use the following two sample estimates for these parameters under the negative binomial distribution in RNA-seq data [23].
Setting the first derivative to zero, we obtain the unrestricted MLEs with respect to x s γ = . Subsequently, we can obtain the MLE of φ by solving the following equation: Several mathematical optimization methods, such as the Newton-Raphson method, can be used to estimate φ for equation (4). Since there is no closing form to estimate the dispersion ϕ, we derived the sample size formula based on a constant value of ϕ estimated from the data.

Constrained maximum likelihood estimate (CMLE) for γ i and ф
The parameters are estimated under the null hypothesis H 0 : γ 0 =γ 1 .
Using partial derivatives of the equation (5) Setting equation (6) to zero, we obtain: Thus, setting the derivative of ф to the zero and 0 0 and 1 t t ρ = =  , the CMLE  a for ф is estimated by solving the following equation: and Since there is no closing form to estimate the dispersion ϕ from MLE and CMLE, we derived the sample size formula based on a fixed and constant value. For simplicity, we set the dispersion in MLE and CMLE to be equal with a combination of fixed dispersion (0.1, 0.5 and 1).

Wald statistical test and log-transformed Wald statistical test
Wald statistical test: The Wald statistical test is an asymptotic test based on the normal approximation, which utilizes the large-sample properties of the MLE. Following procedures from the studies [24][25][26] for comparing two independent Poisson rates with unequal sample frames, we derived the Wald's inference procedures using the properties of ˆi γ and φ estimated from MLE and ˆi γ and φ from CMLE with a negative binomial distribution in two conditions ( ) , where X 0 and X 1 in two conditions are assumed to be independent. For simplicity, we set φ = =  a a with a constant for the sample size and power analysis.
The null hypothesis H 0 in equation (1) is equivalent to H 0 : γ 1 -γ 0 = 0, and consequently we make inferences based on the quantity . In this case, the variance of T is where 1 0 w s s = , is the ratio of total size factors between the two conditions.
Similarly, for using the parameters estimated from CMLE. By substituting ( )

Log-transformation of Wald's statistical test:
To test the null hypothesis, it is also equivalent to test The logarithmic transformation is usually adopted for skewness correction and variance stabilization as suggested by these studies [16,24]. The statistical inference on the quantity is ( ) ln X s correspondingly also have an asymptotically normal by the Delta Method and Slutsky's theorem. Thus, the variance of U is 2 a . U/S U can be used for testing H 0 , when We note that the equations defined in (12) and (13) do not exist when X 0 =0 orX 1 =0. In this case, X 0 or X 1 was adjusted to 0.5 [24,25].

Generalized likelihood ratio test (GLRT)
The GLRT statistic is defined as the ratio of the maximum value of the likelihood function under the restriction of the null hypothesis to the maximum likelihood function under the unrestricted parameter space. For a vector of parameters θ∈Θ, the GLRT forH 0 : θ∈Θ 0 versusH 1 : θ∈Θ 1 is expressed as Where ( ) in equation (14) is obtained using the MLE of θ, where a approximately follows a 2 1 χ distribution, the p-value is approximately: The p-value in equation (16) is further adjusted by the FDR correction when multiple genes are used for the data analysis. Combining these together, the parameter estimates based on the MLEs from two assumptions and the following test statistics are summarized in Table 1.

Sample size calculation for a single gene
In order to calculate the sample size, a power function needs to be constructed. The power of a test is the probability that the null hypothesis is rejected when the alternative hypothesis is true. We derived the sample size under the specified power 1-β and the significance level α with an assumption of a balanced design experiment between conditions (i.e., n 0 =n 1 =n) and one-sided statistical test under Derivation of sample size based on the Wald test statistics: The details of derivation for each formula per test statistic are described in the Appendix: Derivation of sample size calculation in Supplementary Information. Briefly, given the parameters (u 0 and ϕ), the sample size formula based on the Wald test statistics (Z w1 ,Z lw1 ,Z w2 and Z lw2 ) are defined as where 0 u  is estimated from CMLE under H 0 in equation (9) and u 0 is the mean read counts under H 1 or the assumed true mean read counts in the control condition. Under the alternative hypothesis (17)(18)(19)(20) are also true.

Derivation of sample size based on the likelihood ratio test (LRT) statistic:
For the LRT with a negative binomial distribution, it is difficult to derive a closed-form expression of the power function. We used the method [19] to calculate the power of the LRT given a p-value from the equation (16) under LRT. This method originally borrowed a concept from this study [27] to calculate the power. Given a p-value based on the observed joint probability P (X 0 =x 0 , X 1 =x 1 ), the power under the assumption can be expressed as Pr n u w f nw u f nu where X 0 and X 1 are independent, f(u, ф) is the probability mass function of the negative binomial distribution with mean u and dispersion ф, α is the level of significance, and I(⋅) is the indicator function of p-value. Thus, given a nominal power 1-β, the power of the test can be represented as the function of the sample size n in the form of: Therefore, the required sample size n to attain the nominal power 1-β at a significance level α can then be computed by solving equation ( Generalized Likelihood ratio test CMLE/MLE Z lrt γ 0 and γ 1 are true gene expression between two conditions.

Sample size calculation with controlling FDR for testing multiple genes
In an RNA-seq experiment, thousands of genes need to be tested simultaneously for DEGs between conditions. In this case, the sample size calculation for a single gene derived above needs to be further adjusted due to the multiple testing problems. In this section we derive sample size calculations by incorporating FDR controlling based on the statistical tests described in the previous sections. The details of controlling FDR have been given in the study [19]. Briefly, FDR (f) is defined as where m 0 is the number of true null hypotheses, t 1 =E(M 1 ) is the expected number of true rejections, M 0 is the number of false discoveries, M is the total number of genes declared significant, M 1 =M-M 0 and f is the control FDR at a specified level.
By solving equation (23) with respect to α, the marginal type I error level α * for the expected number of true rejections t 1 at a given FDR (f) is Replacing α with α * in equation (24) in equations (17)(18)(19)(20), the corresponding sample size calculation formulas corrected by FDR at level f are, respectively, ( ) Similarly, replacing α with α * for the LRT statistic in equation (22), we obtain the function with respect to n as ( ) Thus, by solving (29) via a numerical approach, the sample size for controlling FDR at level f can be obtained.

Simulation Studies and Comparison of Results
The proposed sample size formulas are derived from the likelihood function based on large sample theory with an approximate normal distribution. The simulation studies include two parts. In the first part, we calculated sample size based on testing single gene from the different formulas. In the second part, we calculated sample size based on testing the multiple genes using FDR adjusted significance α * level. The parameter settings in our simulation studies are based on empirical data sets. A comparison of the simulated power for the sample sizes using different methods was performed.

Sample size calculations and power estimation based on testing a single gene
The purpose of this study is to compare the performance of sample size calculations with the estimated power from our formula with the method based on the exact test in the public study [19]. We set the following inputs based on a single gene. Let the type I error rate α=0.05, the power1-β=0.8, the ratio of total size factors between two condition w=1 and 1.2, the mean counts of gene g in control condition u 0 =1, 5 or 10, the dispersion ф=0.1, 0.5 or 1, and the fold changes ρ=1.5, 2, 3 or 4. Since the read depth across samples in RNA-seq data is usually close to each other, we choose w=1.2 instead of w=2 [19]. For each combination of these designed settings, at first, we used our derived formulas in equations (17-20 and 22) to calculate the required sample size, respectively. Then, we calculated the sample size based on the exact test computed using the R codes with the same input settings. Moreover, for each designed setting, we generated 5000 simulations from independent negative binomial distributions based on the calculated sample size n given the dispersion ф with different mean counts. For the control condition (i=0), we used R to generate random samples given the mean u 0 and ф. For the treatment condition (i=1), we generated random samples given mean u 1 =wρu 0 and ф. The test statistics in equations (10-13 and 14-16) were applied to each simulation sample and the empirical power was obtained as the proportion of simulation samples for which H 0 is rejected with the nominal type I error α=0.05. The results are shown in which reports the estimated sample

Sample size calculation based on testing multiple genes via FDR-controlling method
In this study, we evaluated the performance of the sample size methods based on testing the multiple genes via FDR-controlling method rather than the type I error α, which is widely used in the RNA-seq analysis. We set m=10,000, m 1 =100, m 0 =m -m 1 and want to detect the expected number of true DEG t 1 =80 and the actual power corresponding to the nominal power of 1-β=80%. We also set u 0g =1,5 or 10, ρ g =1.5,2,3 or 4 and ф g =0.1,0.5 or 1.0. With these settings, the new α * =4.25 × 10 -4 was obtained from the equation (24) at a desired FDR (f=0.05). Then, we calculated the sample size by substituting α * and power into the equations (25-29) and the published method using the exact test. For each designed setting, we also conducted 5000 simulations from an independent negative binomial distribution. The number of true DEGs was counted using p-values ≤ α * which is much smaller than the nominal type I error rate 0.05. The empirical power was obtained as the proportion of simulation samples for which H 0 is rejected with the nominal type I error α * =4.25 × 10 -4 . The expected number of true DEGs under α * can be estimated via the multiplication of the estimated power with the total number of true DEGs. The results in Table 2 report the estimated sample size with the empirical power given in parentheses under the cases w=1and 1.2.
Given the sample size obtained from n lw2 , an empirical power for other methods is also computed from 5,000 simulations (Table 3). In order to compare the performance of these methods the paired Wilcoxon signed-rank test is used to test the simulated power in Table 3 for the statistical significance. The results are in Table 4.

Identify the pattern of the sample sizes changing with different values of u 0 , ϕ and ρ for different methods at α and α * levels
First, we identified patterns of sample size changes with different values of u 0 , ϕ and ρ for different methods. Although Table 2 shows similar pattern of sample sizes obtained from the six methods, we found that the required sample sizes n lw1 and n lw2 derived from the log-transformed Wald tests are the smallest compared with the other methods. Figure 1 illustrates n lw2 varying with the values of u 0 when other parameters, such as ρ ∈ (1.5, 2, 3, 4), ф ∈ (.1, .5, 1) and w ∈ (1, 1.2) are fixed. Given the nominal power (1-β=0.8) and α=0.05 or FDR=0.05, as expected, n lw2 decreases as u 0 increases under a fixed ρ and ф. This indicates that for a lowly expressed gene, a larger sample size is required to achieve a detection power of DEGs between two conditions. For a fixed u 0 and ρ, n lw2 increases as ϕ increases (Figure 1). This is also expected because a larger ϕ indicates higher variation of the genes across conditions. Furthermore, for a fixed ϕ and u 0 , n lw2 decreases as ρ increases. This result indicates a smaller n is required for a larger difference of mean read counts between two conditions or vice versa. For the same setting of parameters, we found the n lw2 in Figure 1A and 1C with an equal size factor (w=1) across conditions is slightly larger than the unequal size factor in Figure 1B (w=1). Under the same settings, as expected, the sample sizes in Table 2 with FDR=0.05 are larger than those with α=0.05, indicating that a larger n is required for detecting DEGs while testing thousands of genes simultaneously ( Figure  1C and 1D) compared with testing a single gene (Figure1A and 1B). Table S1 genes in  Table 2 Next we compared the sample sizes (n) estimated from our five derived methods (n w1 , n w2 , n lw1 , n lw2 , n lrt ) and the public method (n exact ) [19]. Figure 2 illustrates that n decreases as increases for all methods given w=1, power=0.8, FDR=0.05, ф ∈ (0.1,0.5,1) and ρ ∈ (1.5,2). Given ρ=1.5, the sample sizes from all methods are getting close to each other when ϕ is 0.1 for small biological variation (Figure 2A-2C). As ϕ increases, the difference between the sample sizes for all methods becomes much larger. Similar patterns were observed given ρ=2 ( Figure 2D-2F). With regards to the n, we noticed that the sample sizes calculated from n lw1 and n lw2 methods (red in Figure 2) are the smallest while maintaining the nominal power close to or above 80%. Among the other four methods, no one performed better than others in all scenarios.

Comparison of the sample calculations from different methods based on testing a single genes in multiple
In addition, the empirical power in parentheses of Table S1 and Table 2 was calculated from simulations with the size of 5,000 corresponding to the estimated n for all methods. The results show almost all of the methods are close to or higher than the desired power. Although the sample sizes calculated using n lw1 and n lw2 are the smallest, we cannot arbitrarily conclude that these two methods are the best because the corresponding empirical powers are varied corresponding to the sample sizes for each method ( Table 2).
For a better comparison with the same settings and a fixed and small n estimated from the log-transformed Wald method (n lw2 ), we observed that n lw1 and n lw2 consistently achieve a better power close to the nominal power 80% or higher in all scenarios compared to other methods (Table  3). We also observed that n lrt and n exact perform similarly and both of them achieve a higher power than n w1 when the fold change is great than 2. However, n w1 achieves a better power than n lrt and n exact when the fold change is at ρ ≤ 2 (Table 3). Table 4 from a paired Wilcoxon ranked test indicates that the empirical power from n lw2 is statistically significant from that achieved using other methods. Table 4 also shows that n lw1 performed significantly better than others four methods. Among these four methods, no one performed better than the others in all scenarios.

Application
Sample size calculation based on RNA-seq data in human breast cancer: To identify DEGs between two conditions, we explored a real human breast cancer dataset to calculate the sample size. Forty Estrogen receptor positive (ER+) and HER2 negative breast cancer primary tumors and 29 uninvolved breast tissue sample that were adjacent to ER+ primary tumors in .fastq format were downloaded from NCBI GEO (series ID GSE58135). The raw sequencing files were mapped to the human hg19 reference genome using tophat2 (v2.0.13) with bowtie version (2.2.3.0). The mapped counts for each gene per sample were then extracted using HTSeq-scripts-count (version 2.7). There are a total of 57,773 genes extracted. After filtering genes with the mean read counts less than one in two groups, 35,112 genes were left. These samples were loaded into edgeR to estimate common dispersion and size factors. With the aid of edgeR, the normalization factors called size factors are estimated using the "RLE" scaling factor n w1 , n w2 , n lw1 and n lw2 are our methods and n exact is the public method. n w1 , n w2 , n lw1 and n lw2 are our methods and n exact is the public method. method [17] and the estimated ratio of total size factors of the samples in each condition (w) is 1.1. The common dispersion ϕ is estimated as 0.48 using the CMLE method [22].
We assumed the top 400 of 35,122 genes (1.1%) are prognostic and have the largest fold changes. The minimum of average read counts among these genes in the control group served as pilot data was estimated as u 0 =1 [16]. In addition, the sample sizes were estimated using u 0 =5,10 and 20. Suppose we want to set the nominal power to 80%, which indicates we want to identify 320 or more of the prognostic genes. Under the control FDR at f=0.10 and 80% power, we can set m=35,112, m 1 =400, m 0 =m-m 1 and t 1 =320. The parameters ρ g (g=1,…, 35, 112) were assumed to be unknown. Let the mean counts in the control condition u 0g =1, 5, 10, or 20 and fold change ρ g =1.5 or 2, 3 and 4 with common dispersion ф g =0.48. With these settings, the new α * =9.992 × 10 -4 was obtained from equation (24) at a desired FDR(f=0.10). Then, we calculated the sample size by substituting α * and the power into equations (25-29) for each method (Table 5).  : Sample size (n) for testing multiple genes is required using the six methods n w1 , n w2 , n lw1 , n lw2 , n lrt and n exact ). n varies with mean counts in control condition (u 0 ) given the combination of fold change (ρ) and dispersion (ϕ) at w=1, 80% power and 0.05 FDR. n lw1 and n lw2 methods in red and blue respectively require smaller n. Figures 2A and 2D Table 5 reports the samples sizes under different scenarios including various minimum mean read counts in control condition (1, 5, 10 and 20) and desired fold changes (1.5, 2, 3 and 4) while controlling the FDR at 0.10. We found that the original RNA-seq experiment [28] with a minimum sample size 31 in each condition can detect more than 80% of the prognostic genes at the FDR (f=0.10) and u 0 =1 if the desired fold change is 3 or more. Moreover, with a minimum sample size 42 in each condition, we found that it can detect more than 80% of the prognostic genes at the FDR (f=0.10) and u 0 =5 if the desired fold change is 2 or more.

Discussion and Conclusion
In this study, five methods (n w1 , n w2 , n lw1 , n lw2 , n lrt ) were derived to calculate sample sizes using the Wald test and LRT statistics based on a negative binomial distribution for modeling an RNA-seq experiment. The parameters are estimated using the MLE and CMLE methods. Since the dispersion estimated from MLE has no closing form, it is difficult to derive the sample size formula. Therefore, all of the methods are based on a fixed and constant dispersion. A log-transformed approach corresponding to the modified Z w1 and Z w2 was used to derive two other test statistics (Z lw1 and Z lw2 ). For all these statistical tests as well as the exact test [19] that are used to derive the sample size calculation formulas, gene expression levels are assumed to be independent in each sample. Although this assumption might not hold in reality, it is widely used in RNA-seq as well as in microRNA data analysis. In this study, we assume equal sample size in the two conditions to derive the sample size formula. The derived formula for sample size calculations can be easily extended to the unequal sample sizes by setting n 1 =kn 0 . In our simulation study, we set the ratio of total size factors in two conditions as 1 and 1.2 instead of w=2 in the study [19]. In reality, the read depths of RNA-seq samples generated from the same run are very close to each other across conditions. Therefore, we think w=1.2 or close to 1 is more common than w=2. Furthermore, in our simulation and application studies, the minimum sample sizes required to achieve a nominal power of 80% with a predefined FDR (f=0.05 or 0.10) are usually larger than those in an RNA-seq experiment due to the real costs. In such a situation, we can increase the read depth per sample to indirectly increase the mean of read counts u 0 in the control condition. Thus, the required sample sizes can be decreased correspondingly.
Among the methods we evaluated, the simulation results show that n lw2 from the log transformed Wald test with the parameters estimated from CMLE is the best method because a smaller sample size is required for designing an RNA-seq experiment while achieving a power close to or higher than 80% at a pre-defined FDR=0.05. The second best method is n lw1 with the parameters estimated from unrestricted MLEs. We also found that n lrt and n exact methods perform better than n w1 method based on the estimated power given the genes with a fold change >2. However, n w1 achieve a better power than n lrt and n exact given a fold change ≤2. In summary, n w1 , n w2 , n lrt and n exact methods varied with different scenarios. Finally, since the log-transformed sample size calculation methods are more robust, simpler and require less time, we hope our tables can help and benefit for researchers and scientists in the design of RNA-seq experiments.  The dispersion ф=0.48 and the ratio of the size factor w=1.1 are estimated using edgeR package particularly for RNA-seq data.