Development of RnaSeqSampleSize version 2

Shilin Zhao
Department of Biostatistics
Vanderbilt University School of Medicine

2020-06-26

library(knitr)
library(ggplot2)
library(RColorBrewer)
library(reshape2)
library(RnaSeqSampleSize)
library(DESeq2)

#source('D:/source/cqsR/R/summaryTableFunctions.R')
#devtools::source_url("https://raw.github.com/slzhao/cqsR/master/R/summaryTableFunctions.R")
library(cqsR)

Make/Load Data

set.seed(123)
rawData<-matrix(nrow=200,ncol=10,sample(1:1000,2000,replace = TRUE))
row.names(rawData)<-paste0("Gene",1:200)
colnames(rawData)<-paste0("Sample",1:10)

rawDataDesign=data.frame(
                          #Sample=colnames(rawData),
                         #Group1=c(rep("Treat1",3),rep("Treat2",3),rep("Control",4)),
                         Group1=c(rep("Treat",5),rep("Control",5)),
                         Group2=rep(c("Strain1","Strain2"),5))

Estimate distribution of genes by DESeq2

dds <- DESeqDataSetFromMatrix(countData = rawData,
                              colData = rawDataDesign,
                              design = ~ Group1+Group2)
dds <- estimateSizeFactors(dds)

dds <- estimateDispersions(dds) #this not work sometimes

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

dds <- nbinomLRT(dds, reduced = ~ Group2) #interested in Group1, so compare with model with only Group2
mypar<-data.frame(coef(dds),disp=dispersions(dds))
colnames(mypar)[1:3]<-c("b0","b1","b2")

#Only one factor
dds <- DESeqDataSetFromMatrix(countData = rawData,
                              colData = rawDataDesign,
                              design = ~ Group1)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

dds <- nbinomLRT(dds, reduced = ~ 1) 
myparOneFactor<-data.frame(coef(dds),disp=dispersions(dds))
colnames(mypar)[1:3]<-c("b0","b1")

## Warning in colnames(mypar)[1:3] <- c("b0", "b1"): number of items to replace is not a multiple of replacement length

Estimation

source('D:/source/cqsR/R/RnaSeqSampleSize2.R')

## Loading required package: ssanv

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:VariantAnnotation':
## 
##     select

aa<-apply(mypar[1:6,], 1, get.v)

## [1] "x= 1  f(x)= 0.999995041944263"
## [1] "x= 4097  f(x)= 0.478734666504519"
## [1] "x= 8193  f(x)= 0.161319421569493"
## [1] "x= 12289  f(x)= 0.0490743034992552"
## [1] "x= 16385  f(x)= 0.0141844831980985"
## [1] "x= 20481  f(x)= 0.00397507534421078"
## [1] "x= 24577  f(x)= 0.0010910157830575"
## [1] "x= 28673  f(x)= 0.000294914853570686"
## [1] "x= 32769  f(x)= 7.87267519890094e-05"
## [1] "x= 36865  f(x)= 2.07296064889025e-05"
## [1] "x= 40961  f(x)= 5.31656230374118e-06"
## [1] "x= 45057  f(x)= 1.25111597659618e-06"
## [1] "x= 49153  f(x)= 1.85337130419376e-07"
## [1] "x= 53249  f(x)= -9.26430553695727e-08"
## [1] "x= 51201  f(x)= 1.24907899437545e-09"
## [1] "x= 52225  f(x)= -5.35314724364523e-08"
## [1] "x= 51713  f(x)= -2.84419878808225e-08"
## [1] "x= 51457  f(x)= -1.4220955604954e-08"
## [1] "x= 51329  f(x)= -6.64868771212495e-09"
## [1] "x= 51265  f(x)= -2.74135038697818e-09"
## [1] "x= 51233  f(x)= -7.56631514782099e-10"
## [1] "x= 51217  f(x)= 2.4358599574487e-10"
## [1] "x= 51225  f(x)= -2.57180466212458e-10"
## [1] "x= 51221  f(x)= -6.96187771779984e-12"
## [1] "x= 51219  f(x)= 1.18270871385199e-10"
## [1] "x= 51220  f(x)= 5.56442033036466e-11"
## [1] "x= 1  f(x)= 0.999999809994343"
## [1] "x= 4097  f(x)= 0.90318763978619"
## [1] "x= 8193  f(x)= 0.595536811294492"
## [1] "x= 12289  f(x)= 0.311102252208315"
## [1] "x= 16385  f(x)= 0.140319382958395"
## [1] "x= 20481  f(x)= 0.0574105299711482"
## [1] "x= 24577  f(x)= 0.0219307481056441"
## [1] "x= 28673  f(x)= 0.00796328265377202"
## [1] "x= 32769  f(x)= 0.00278118431921175"
## [1] "x= 36865  f(x)= 0.00094185221933813"
## [1] "x= 40961  f(x)= 0.000311037187923665"
## [1] "x= 45057  f(x)= 0.00010054131333095"
## [1] "x= 49153  f(x)= 3.1855297094869e-05"
## [1] "x= 53249  f(x)= 9.85596900663679e-06"
## [1] "x= 57345  f(x)= 2.91979055538144e-06"
## [1] "x= 61441  f(x)= 7.62126878493179e-07"
## [1] "x= 65537  f(x)= 9.87124793894778e-08"
## [1] "x= 69633  f(x)= -1.03198793777398e-07"
## [1] "x= 67585  f(x)= -3.15352746226606e-08"
## [1] "x= 66561  f(x)= 2.39496059040216e-08"
## [1] "x= 67073  f(x)= -5.85932034098446e-09"
## [1] "x= 66817  f(x)= 8.48952406322165e-09"
## [1] "x= 66945  f(x)= 1.18129686107694e-09"
## [1] "x= 67009  f(x)= -2.37184586635987e-09"
## [1] "x= 66977  f(x)= -6.03559288263058e-10"
## [1] "x= 66961  f(x)= 2.86787981323565e-10"
## [1] "x= 66969  f(x)= -1.58904650201857e-10"
## [1] "x= 66965  f(x)= 6.38117660273177e-11"
## [1] "x= 66967  f(x)= -4.75788981630118e-11"
## [1] "x= 66966  f(x)= 8.10831756308481e-12"
## [1] "x= 1  f(x)= 0.999999527321979"
## [1] "x= 4097  f(x)= 0.707508862054515"
## [1] "x= 8193  f(x)= 0.371472671565136"
## [1] "x= 12289  f(x)= 0.173556598543287"
## [1] "x= 16385  f(x)= 0.0761305010583297"
## [1] "x= 20481  f(x)= 0.0320969036313218"
## [1] "x= 24577  f(x)= 0.0131684701485702"
## [1] "x= 28673  f(x)= 0.00529622355055984"
## [1] "x= 32769  f(x)= 0.00209797302156031"
## [1] "x= 36865  f(x)= 0.00082112787379463"
## [1] "x= 40961  f(x)= 0.000318224513229391"
## [1] "x= 45057  f(x)= 0.000122276506236353"
## [1] "x= 49153  f(x)= 4.65986910253992e-05"
## [1] "x= 53249  f(x)= 1.75844201002487e-05"
## [1] "x= 57345  f(x)= 6.52928452358511e-06"
## [1] "x= 61441  f(x)= 2.33929345692734e-06"
## [1] "x= 65537  f(x)= 7.58526018037758e-07"
## [1] "x= 69633  f(x)= 1.645353888636e-07"
## [1] "x= 73729  f(x)= -5.78745795091379e-08"
## [1] "x= 71681  f(x)= 2.65083628745251e-08"
## [1] "x= 72705  f(x)= -2.08517790183082e-08"
## [1] "x= 72193  f(x)= 1.37288564666834e-09"
## [1] "x= 72449  f(x)= -1.00813259368617e-08"
## [1] "x= 72321  f(x)= -4.44233191015322e-09"
## [1] "x= 72257  f(x)= -1.55709021764054e-09"
## [1] "x= 72225  f(x)= -9.77370191100733e-11"
## [1] "x= 72209  f(x)= 6.36160224130349e-10"
## [1] "x= 72217  f(x)= 2.68858757492249e-10"
## [1] "x= 72221  f(x)= 8.54727424692532e-11"
## [1] "x= 72223  f(x)= -6.15415944305634e-12"
## [1] "x= 72222  f(x)= 3.96537849127198e-11"
## [1] "x= 1  f(x)= 0.999999805458011"
## [1] "x= 4097  f(x)= 0.955322737258444"
## [1] "x= 8193  f(x)= 0.842276133098205"
## [1] "x= 12289  f(x)= 0.702850014875448"
## [1] "x= 16385  f(x)= 0.564392500363751"
## [1] "x= 20481  f(x)= 0.440659693565679"
## [1] "x= 24577  f(x)= 0.336790133342676"
## [1] "x= 28673  f(x)= 0.253134406869171"
## [1] "x= 32769  f(x)= 0.187716162946779"
## [1] "x= 36865  f(x)= 0.137675361392116"
## [1] "x= 40961  f(x)= 0.10004708874203"
## [1] "x= 45057  f(x)= 0.0721367923436979"
## [1] "x= 49153  f(x)= 0.0516648019922219"
## [1] "x= 53249  f(x)= 0.0367878564744426"
## [1] "x= 57345  f(x)= 0.026061559106629"
## [1] "x= 61441  f(x)= 0.0183798436474313"
## [1] "x= 65537  f(x)= 0.0129105587867046"
## [1] "x= 69633  f(x)= 0.00903633444799859"
## [1] "x= 73729  f(x)= 0.00630432276429235"
## [1] "x= 77825  f(x)= 0.00438546915792309"
## [1] "x= 81921  f(x)= 0.00304256173928727"
## [1] "x= 86017  f(x)= 0.00210575264624476"
## [1] "x= 90113  f(x)= 0.00145413857389112"
## [1] "x= 94209  f(x)= 0.00100209590057958"
## [1] "x= 98305  f(x)= 0.000689258795567823"
## [1] "x= 102401  f(x)= 0.00047323848032356"
## [1] "x= 106497  f(x)= 0.000324375888396881"
## [1] "x= 110593  f(x)= 0.000221985446777775"
## [1] "x= 114689  f(x)= 0.000151681980589344"
## [1] "x= 118785  f(x)= 0.000103488126499289"
## [1] "x= 122881  f(x)= 7.05003769103456e-05"
## [1] "x= 126977  f(x)= 4.7952612822133e-05"
## [1] "x= 131073  f(x)= 3.25610237515428e-05"
## [1] "x= 135169  f(x)= 2.20673310158224e-05"
## [1] "x= 139265  f(x)= 1.49212061101821e-05"
## [1] "x= 143361  f(x)= 1.00600469664412e-05"
## [1] "x= 147457  f(x)= 6.75663214275451e-06"
## [1] "x= 151553  f(x)= 4.51396289191212e-06"
## [1] "x= 155649  f(x)= 2.9928236348124e-06"
## [1] "x= 159745  f(x)= 1.96197395501867e-06"
## [1] "x= 163841  f(x)= 1.26396065436581e-06"
## [1] "x= 167937  f(x)= 7.91688906439115e-07"
## [1] "x= 172033  f(x)= 4.72390500873266e-07"
## [1] "x= 176129  f(x)= 2.56668982843113e-07"
## [1] "x= 180225  f(x)= 1.11023687196938e-07"
## [1] "x= 184321  f(x)= 1.27541085613263e-08"
## [1] "x= 188417  f(x)= -5.35093098013147e-08"
## [1] "x= 186369  f(x)= -2.36337666537211e-08"
## [1] "x= 185345  f(x)= -6.3357166716905e-09"
## [1] "x= 184833  f(x)= 2.97409671012155e-09"
## [1] "x= 185089  f(x)= -1.73815128682151e-09"
## [1] "x= 184961  f(x)= 6.03460930427361e-10"
## [1] "x= 185025  f(x)= -5.70950898124111e-10"
## [1] "x= 184993  f(x)= 1.53508186019584e-11"
## [1] "x= 185009  f(x)= -2.78025742541966e-10"
## [1] "x= 185001  f(x)= -1.31393930950049e-10"
## [1] "x= 184997  f(x)= -5.80356788337048e-11"
## [1] "x= 184995  f(x)= -2.13459614576005e-11"
## [1] "x= 184994  f(x)= -2.99845434820436e-12"
## [1] "x= 1  f(x)= 0.999999804308077"
## [1] "x= 4097  f(x)= 0.863000993595377"
## [1] "x= 8193  f(x)= 0.589821365838866"
## [1] "x= 12289  f(x)= 0.35450497209389"
## [1] "x= 16385  f(x)= 0.197238886166511"
## [1] "x= 20481  f(x)= 0.104279078572677"
## [1] "x= 24577  f(x)= 0.0531755922369667"
## [1] "x= 28673  f(x)= 0.0263970169873115"
## [1] "x= 32769  f(x)= 0.0128345440209529"
## [1] "x= 36865  f(x)= 0.00613808307229837"
## [1] "x= 40961  f(x)= 0.00289628091752136"
## [1] "x= 45057  f(x)= 0.00135141464964625"
## [1] "x= 49153  f(x)= 0.000624620333055607"
## [1] "x= 53249  f(x)= 0.000286332983020645"
## [1] "x= 57345  f(x)= 0.000130293993516251"
## [1] "x= 61441  f(x)= 5.88742970168938e-05"
## [1] "x= 65537  f(x)= 2.64036306524227e-05"
## [1] "x= 69633  f(x)= 1.17272633188103e-05"
## [1] "x= 73729  f(x)= 5.12791470689971e-06"
## [1] "x= 77825  f(x)= 2.17405916090785e-06"
## [1] "x= 81921  f(x)= 8.57337615784109e-07"
## [1] "x= 86017  f(x)= 2.72555336233719e-07"
## [1] "x= 90113  f(x)= 1.37081647631819e-08"
## [1] "x= 94209  f(x)= -1.00520208321644e-07"
## [1] "x= 92161  f(x)= -5.49477752384696e-08"
## [1] "x= 91137  f(x)= -2.41217029071113e-08"
## [1] "x= 90625  f(x)= -6.17367112072193e-09"
## [1] "x= 90369  f(x)= 3.51305305668804e-09"
## [1] "x= 90497  f(x)= -1.39224944091e-09"
## [1] "x= 90433  f(x)= 1.04471888154417e-09"
## [1] "x= 90465  f(x)= -1.77661087259406e-10"
## [1] "x= 90449  f(x)= 4.32551842080902e-10"
## [1] "x= 90457  f(x)= 1.27201502309199e-10"
## [1] "x= 90461  f(x)= -2.52907127346774e-11"
## [1] "x= 90459  f(x)= 5.09401586612344e-11"
## [1] "x= 90460  f(x)= 1.28209146893304e-11"
## [1] "x= 1  f(x)= 0.999991010906521"
## [1] "x= 4097  f(x)= 0.74144546324518"
## [1] "x= 8193  f(x)= 0.488695152458006"
## [1] "x= 12289  f(x)= 0.309502055430197"
## [1] "x= 16385  f(x)= 0.191775868167439"
## [1] "x= 20481  f(x)= 0.117175539046921"
## [1] "x= 24577  f(x)= 0.0708954134124996"
## [1] "x= 28673  f(x)= 0.0425832552115229"
## [1] "x= 32769  f(x)= 0.025434184188511"
## [1] "x= 36865  f(x)= 0.0151234579325277"
## [1] "x= 40961  f(x)= 0.0089597506516546"
## [1] "x= 45057  f(x)= 0.00529197650522541"
## [1] "x= 49153  f(x)= 0.00311759343046159"
## [1] "x= 53249  f(x)= 0.0018325568742401"
## [1] "x= 57345  f(x)= 0.00107511301834022"
## [1] "x= 61441  f(x)= 0.000629656471837324"
## [1] "x= 65537  f(x)= 0.000368192710333557"
## [1] "x= 69633  f(x)= 0.000214986807308055"
## [1] "x= 73729  f(x)= 0.000125350283303191"
## [1] "x= 77825  f(x)= 7.29766106047129e-05"
## [1] "x= 81921  f(x)= 4.24118501249757e-05"
## [1] "x= 86017  f(x)= 2.45937626140581e-05"
## [1] "x= 90113  f(x)= 1.42166148286515e-05"
## [1] "x= 94209  f(x)= 8.1783754853507e-06"
## [1] "x= 98305  f(x)= 4.66769270535405e-06"
## [1] "x= 102401  f(x)= 2.62806393168263e-06"
## [1] "x= 106497  f(x)= 1.44389203767284e-06"
## [1] "x= 110593  f(x)= 7.5681527086659e-07"
## [1] "x= 114689  f(x)= 3.5839360327739e-07"
## [1] "x= 118785  f(x)= 1.27481921662799e-07"
## [1] "x= 122881  f(x)= -6.27933822715722e-09"
## [1] "x= 120833  f(x)= 5.15245027200475e-08"
## [1] "x= 121857  f(x)= 2.06518037589225e-08"
## [1] "x= 122369  f(x)= 6.72655029171424e-09"
## [1] "x= 122625  f(x)= 1.12570338304904e-10"
## [1] "x= 122753  f(x)= -3.11067159345074e-09"
## [1] "x= 122689  f(x)= -1.50593082298864e-09"
## [1] "x= 122657  f(x)= -6.98407630959545e-10"
## [1] "x= 122641  f(x)= -2.9335141437195e-10"
## [1] "x= 122633  f(x)= -9.0498845370886e-11"
## [1] "x= 122629  f(x)= 1.10086552039019e-11"
## [1] "x= 122631  f(x)= -3.97518660952276e-11"
## [1] "x= 122630  f(x)= -1.43732984242317e-11"

Test of parallel Estimation

library("future.apply")

## Warning: package 'future.apply' was built under R version 3.6.1

## Loading required package: future

## Warning: package 'future' was built under R version 3.6.1

## 
## Attaching package: 'future'

## The following object is masked from 'package:survival':
## 
##     cluster

## The following object is masked from 'package:SummarizedExperiment':
## 
##     values

## The following object is masked from 'package:GenomicRanges':
## 
##     values

## The following object is masked from 'package:IRanges':
## 
##     values

## The following object is masked from 'package:S4Vectors':
## 
##     values

plan(multiprocess) ## Run in parallel on local computer

aaParallel <- future_apply(mypar[1:6,],1, get.v)

## [1] "x= 1  f(x)= 0.999995041944263"
## [1] "x= 4097  f(x)= 0.478734666504519"
## [1] "x= 8193  f(x)= 0.161319421569493"
## [1] "x= 12289  f(x)= 0.0490743034992552"
## [1] "x= 16385  f(x)= 0.0141844831980985"
## [1] "x= 20481  f(x)= 0.00397507534421078"
## [1] "x= 24577  f(x)= 0.0010910157830575"
## [1] "x= 28673  f(x)= 0.000294914853570686"
## [1] "x= 32769  f(x)= 7.87267519890094e-05"
## [1] "x= 36865  f(x)= 2.07296064889025e-05"
## [1] "x= 40961  f(x)= 5.31656230374118e-06"
## [1] "x= 45057  f(x)= 1.25111597659618e-06"
## [1] "x= 49153  f(x)= 1.85337130419376e-07"
## [1] "x= 53249  f(x)= -9.26430553695727e-08"
## [1] "x= 51201  f(x)= 1.24907899437545e-09"
## [1] "x= 52225  f(x)= -5.35314724364523e-08"
## [1] "x= 51713  f(x)= -2.84419878808225e-08"
## [1] "x= 51457  f(x)= -1.4220955604954e-08"
## [1] "x= 51329  f(x)= -6.64868771212495e-09"
## [1] "x= 51265  f(x)= -2.74135038697818e-09"
## [1] "x= 51233  f(x)= -7.56631514782099e-10"
## [1] "x= 51217  f(x)= 2.4358599574487e-10"
## [1] "x= 51225  f(x)= -2.57180466212458e-10"
## [1] "x= 51221  f(x)= -6.96187771779984e-12"
## [1] "x= 51219  f(x)= 1.18270871385199e-10"
## [1] "x= 51220  f(x)= 5.56442033036466e-11"
## [1] "x= 1  f(x)= 0.999999809994343"
## [1] "x= 4097  f(x)= 0.90318763978619"
## [1] "x= 8193  f(x)= 0.595536811294492"
## [1] "x= 12289  f(x)= 0.311102252208315"
## [1] "x= 16385  f(x)= 0.140319382958395"
## [1] "x= 20481  f(x)= 0.0574105299711482"
## [1] "x= 24577  f(x)= 0.0219307481056441"
## [1] "x= 28673  f(x)= 0.00796328265377202"
## [1] "x= 32769  f(x)= 0.00278118431921175"
## [1] "x= 36865  f(x)= 0.00094185221933813"
## [1] "x= 40961  f(x)= 0.000311037187923665"
## [1] "x= 45057  f(x)= 0.00010054131333095"
## [1] "x= 49153  f(x)= 3.1855297094869e-05"
## [1] "x= 53249  f(x)= 9.85596900663679e-06"
## [1] "x= 57345  f(x)= 2.91979055538144e-06"
## [1] "x= 61441  f(x)= 7.62126878493179e-07"
## [1] "x= 65537  f(x)= 9.87124793894778e-08"
## [1] "x= 69633  f(x)= -1.03198793777398e-07"
## [1] "x= 67585  f(x)= -3.15352746226606e-08"
## [1] "x= 66561  f(x)= 2.39496059040216e-08"
## [1] "x= 67073  f(x)= -5.85932034098446e-09"
## [1] "x= 66817  f(x)= 8.48952406322165e-09"
## [1] "x= 66945  f(x)= 1.18129686107694e-09"
## [1] "x= 67009  f(x)= -2.37184586635987e-09"
## [1] "x= 66977  f(x)= -6.03559288263058e-10"
## [1] "x= 66961  f(x)= 2.86787981323565e-10"
## [1] "x= 66969  f(x)= -1.58904650201857e-10"
## [1] "x= 66965  f(x)= 6.38117660273177e-11"
## [1] "x= 66967  f(x)= -4.75788981630118e-11"
## [1] "x= 66966  f(x)= 8.10831756308481e-12"
## [1] "x= 1  f(x)= 0.999999527321979"
## [1] "x= 4097  f(x)= 0.707508862054515"
## [1] "x= 8193  f(x)= 0.371472671565136"
## [1] "x= 12289  f(x)= 0.173556598543287"
## [1] "x= 16385  f(x)= 0.0761305010583297"
## [1] "x= 20481  f(x)= 0.0320969036313218"
## [1] "x= 24577  f(x)= 0.0131684701485702"
## [1] "x= 28673  f(x)= 0.00529622355055984"
## [1] "x= 32769  f(x)= 0.00209797302156031"
## [1] "x= 36865  f(x)= 0.00082112787379463"
## [1] "x= 40961  f(x)= 0.000318224513229391"
## [1] "x= 45057  f(x)= 0.000122276506236353"
## [1] "x= 49153  f(x)= 4.65986910253992e-05"
## [1] "x= 53249  f(x)= 1.75844201002487e-05"
## [1] "x= 57345  f(x)= 6.52928452358511e-06"
## [1] "x= 61441  f(x)= 2.33929345692734e-06"
## [1] "x= 65537  f(x)= 7.58526018037758e-07"
## [1] "x= 69633  f(x)= 1.645353888636e-07"
## [1] "x= 73729  f(x)= -5.78745795091379e-08"
## [1] "x= 71681  f(x)= 2.65083628745251e-08"
## [1] "x= 72705  f(x)= -2.08517790183082e-08"
## [1] "x= 72193  f(x)= 1.37288564666834e-09"
## [1] "x= 72449  f(x)= -1.00813259368617e-08"
## [1] "x= 72321  f(x)= -4.44233191015322e-09"
## [1] "x= 72257  f(x)= -1.55709021764054e-09"
## [1] "x= 72225  f(x)= -9.77370191100733e-11"
## [1] "x= 72209  f(x)= 6.36160224130349e-10"
## [1] "x= 72217  f(x)= 2.68858757492249e-10"
## [1] "x= 72221  f(x)= 8.54727424692532e-11"
## [1] "x= 72223  f(x)= -6.15415944305634e-12"
## [1] "x= 72222  f(x)= 3.96537849127198e-11"
## [1] "x= 1  f(x)= 0.999999805458011"
## [1] "x= 4097  f(x)= 0.955322737258444"
## [1] "x= 8193  f(x)= 0.842276133098205"
## [1] "x= 12289  f(x)= 0.702850014875448"
## [1] "x= 16385  f(x)= 0.564392500363751"
## [1] "x= 20481  f(x)= 0.440659693565679"
## [1] "x= 24577  f(x)= 0.336790133342676"
## [1] "x= 28673  f(x)= 0.253134406869171"
## [1] "x= 32769  f(x)= 0.187716162946779"
## [1] "x= 36865  f(x)= 0.137675361392116"
## [1] "x= 40961  f(x)= 0.10004708874203"
## [1] "x= 45057  f(x)= 0.0721367923436979"
## [1] "x= 49153  f(x)= 0.0516648019922219"
## [1] "x= 53249  f(x)= 0.0367878564744426"
## [1] "x= 57345  f(x)= 0.026061559106629"
## [1] "x= 61441  f(x)= 0.0183798436474313"
## [1] "x= 65537  f(x)= 0.0129105587867046"
## [1] "x= 69633  f(x)= 0.00903633444799859"
## [1] "x= 73729  f(x)= 0.00630432276429235"
## [1] "x= 77825  f(x)= 0.00438546915792309"
## [1] "x= 81921  f(x)= 0.00304256173928727"
## [1] "x= 86017  f(x)= 0.00210575264624476"
## [1] "x= 90113  f(x)= 0.00145413857389112"
## [1] "x= 94209  f(x)= 0.00100209590057958"
## [1] "x= 98305  f(x)= 0.000689258795567823"
## [1] "x= 102401  f(x)= 0.00047323848032356"
## [1] "x= 106497  f(x)= 0.000324375888396881"
## [1] "x= 110593  f(x)= 0.000221985446777775"
## [1] "x= 114689  f(x)= 0.000151681980589344"
## [1] "x= 118785  f(x)= 0.000103488126499289"
## [1] "x= 122881  f(x)= 7.05003769103456e-05"
## [1] "x= 126977  f(x)= 4.7952612822133e-05"
## [1] "x= 131073  f(x)= 3.25610237515428e-05"
## [1] "x= 135169  f(x)= 2.20673310158224e-05"
## [1] "x= 139265  f(x)= 1.49212061101821e-05"
## [1] "x= 143361  f(x)= 1.00600469664412e-05"
## [1] "x= 147457  f(x)= 6.75663214275451e-06"
## [1] "x= 151553  f(x)= 4.51396289191212e-06"
## [1] "x= 155649  f(x)= 2.9928236348124e-06"
## [1] "x= 159745  f(x)= 1.96197395501867e-06"
## [1] "x= 163841  f(x)= 1.26396065436581e-06"
## [1] "x= 167937  f(x)= 7.91688906439115e-07"
## [1] "x= 172033  f(x)= 4.72390500873266e-07"
## [1] "x= 176129  f(x)= 2.56668982843113e-07"
## [1] "x= 180225  f(x)= 1.11023687196938e-07"
## [1] "x= 184321  f(x)= 1.27541085613263e-08"
## [1] "x= 188417  f(x)= -5.35093098013147e-08"
## [1] "x= 186369  f(x)= -2.36337666537211e-08"
## [1] "x= 185345  f(x)= -6.3357166716905e-09"
## [1] "x= 184833  f(x)= 2.97409671012155e-09"
## [1] "x= 185089  f(x)= -1.73815128682151e-09"
## [1] "x= 184961  f(x)= 6.03460930427361e-10"
## [1] "x= 185025  f(x)= -5.70950898124111e-10"
## [1] "x= 184993  f(x)= 1.53508186019584e-11"
## [1] "x= 185009  f(x)= -2.78025742541966e-10"
## [1] "x= 185001  f(x)= -1.31393930950049e-10"
## [1] "x= 184997  f(x)= -5.80356788337048e-11"
## [1] "x= 184995  f(x)= -2.13459614576005e-11"
## [1] "x= 184994  f(x)= -2.99845434820436e-12"
## [1] "x= 1  f(x)= 0.999999804308077"
## [1] "x= 4097  f(x)= 0.863000993595377"
## [1] "x= 8193  f(x)= 0.589821365838866"
## [1] "x= 12289  f(x)= 0.35450497209389"
## [1] "x= 16385  f(x)= 0.197238886166511"
## [1] "x= 20481  f(x)= 0.104279078572677"
## [1] "x= 24577  f(x)= 0.0531755922369667"
## [1] "x= 28673  f(x)= 0.0263970169873115"
## [1] "x= 32769  f(x)= 0.0128345440209529"
## [1] "x= 36865  f(x)= 0.00613808307229837"
## [1] "x= 40961  f(x)= 0.00289628091752136"
## [1] "x= 45057  f(x)= 0.00135141464964625"
## [1] "x= 49153  f(x)= 0.000624620333055607"
## [1] "x= 53249  f(x)= 0.000286332983020645"
## [1] "x= 57345  f(x)= 0.000130293993516251"
## [1] "x= 61441  f(x)= 5.88742970168938e-05"
## [1] "x= 65537  f(x)= 2.64036306524227e-05"
## [1] "x= 69633  f(x)= 1.17272633188103e-05"
## [1] "x= 73729  f(x)= 5.12791470689971e-06"
## [1] "x= 77825  f(x)= 2.17405916090785e-06"
## [1] "x= 81921  f(x)= 8.57337615784109e-07"
## [1] "x= 86017  f(x)= 2.72555336233719e-07"
## [1] "x= 90113  f(x)= 1.37081647631819e-08"
## [1] "x= 94209  f(x)= -1.00520208321644e-07"
## [1] "x= 92161  f(x)= -5.49477752384696e-08"
## [1] "x= 91137  f(x)= -2.41217029071113e-08"
## [1] "x= 90625  f(x)= -6.17367112072193e-09"
## [1] "x= 90369  f(x)= 3.51305305668804e-09"
## [1] "x= 90497  f(x)= -1.39224944091e-09"
## [1] "x= 90433  f(x)= 1.04471888154417e-09"
## [1] "x= 90465  f(x)= -1.77661087259406e-10"
## [1] "x= 90449  f(x)= 4.32551842080902e-10"
## [1] "x= 90457  f(x)= 1.27201502309199e-10"
## [1] "x= 90461  f(x)= -2.52907127346774e-11"
## [1] "x= 90459  f(x)= 5.09401586612344e-11"
## [1] "x= 90460  f(x)= 1.28209146893304e-11"
## [1] "x= 1  f(x)= 0.999991010906521"
## [1] "x= 4097  f(x)= 0.74144546324518"
## [1] "x= 8193  f(x)= 0.488695152458006"
## [1] "x= 12289  f(x)= 0.309502055430197"
## [1] "x= 16385  f(x)= 0.191775868167439"
## [1] "x= 20481  f(x)= 0.117175539046921"
## [1] "x= 24577  f(x)= 0.0708954134124996"
## [1] "x= 28673  f(x)= 0.0425832552115229"
## [1] "x= 32769  f(x)= 0.025434184188511"
## [1] "x= 36865  f(x)= 0.0151234579325277"
## [1] "x= 40961  f(x)= 0.0089597506516546"
## [1] "x= 45057  f(x)= 0.00529197650522541"
## [1] "x= 49153  f(x)= 0.00311759343046159"
## [1] "x= 53249  f(x)= 0.0018325568742401"
## [1] "x= 57345  f(x)= 0.00107511301834022"
## [1] "x= 61441  f(x)= 0.000629656471837324"
## [1] "x= 65537  f(x)= 0.000368192710333557"
## [1] "x= 69633  f(x)= 0.000214986807308055"
## [1] "x= 73729  f(x)= 0.000125350283303191"
## [1] "x= 77825  f(x)= 7.29766106047129e-05"
## [1] "x= 81921  f(x)= 4.24118501249757e-05"
## [1] "x= 86017  f(x)= 2.45937626140581e-05"
## [1] "x= 90113  f(x)= 1.42166148286515e-05"
## [1] "x= 94209  f(x)= 8.1783754853507e-06"
## [1] "x= 98305  f(x)= 4.66769270535405e-06"
## [1] "x= 102401  f(x)= 2.62806393168263e-06"
## [1] "x= 106497  f(x)= 1.44389203767284e-06"
## [1] "x= 110593  f(x)= 7.5681527086659e-07"
## [1] "x= 114689  f(x)= 3.5839360327739e-07"
## [1] "x= 118785  f(x)= 1.27481921662799e-07"
## [1] "x= 122881  f(x)= -6.27933822715722e-09"
## [1] "x= 120833  f(x)= 5.15245027200475e-08"
## [1] "x= 121857  f(x)= 2.06518037589225e-08"
## [1] "x= 122369  f(x)= 6.72655029171424e-09"
## [1] "x= 122625  f(x)= 1.12570338304904e-10"
## [1] "x= 122753  f(x)= -3.11067159345074e-09"
## [1] "x= 122689  f(x)= -1.50593082298864e-09"
## [1] "x= 122657  f(x)= -6.98407630959545e-10"
## [1] "x= 122641  f(x)= -2.9335141437195e-10"
## [1] "x= 122633  f(x)= -9.0498845370886e-11"
## [1] "x= 122629  f(x)= 1.10086552039019e-11"
## [1] "x= 122631  f(x)= -3.97518660952276e-11"
## [1] "x= 122630  f(x)= -1.43732984242317e-11"

Reproduce Real Examples in Paper

#This matches with Table 2
exemplary.data=betaToData(Betas=c(1.61,0.77),phi=0.25,sampleSize=10,error=10^(-7))
dataToPower(exemplary.data,alpha=0.05)

## [1] 0.4294029

exemplary.data=betaToData(Betas=c(1.61,0.77),phi=0.25,sampleSize=20,error=10^(-7))
dataToPower(exemplary.data,alpha=0.05)

## [1] 0.7122415

#3.1 Simulation studies, Table 3
exemplary.data=betaToData(Betas=c(1.61,log(1.5)),phi=0.25,sampleSize=86,error=10^(-7))
dataToPower(exemplary.data)

## [1] 0.8116428

est_power_fdr(exemplary.data,fdr=0.01,m=10000,m1=10000*0.01)

## [1] 0.1366154

#3.2.1 Transcript regulation data set
exemplary.data=betaToData(Betas=c(log(1.67),log(4)),phi=0.6513,sampleSize=20,error=10^(-7))
dataToPower(exemplary.data)

## [1] 0.7860664

est_power_fdr(exemplary.data,fdr=0.1,m=13267,m1=13267*0.1)

## [1] 0.5686321

est_power_fdr(exemplary.data,fdr=0.1,m=13267,m1=13267*0.01)

## [1] 0.2858171

#Table 5
exemplary.data=betaToData(Betas=c(log(2.5),log(1.1)),phi=1,sampleSize=1576,error=10^(-7))
dataToPower(exemplary.data)

## [1] 0.3628025

est_power_fdr(exemplary.data,fdr=0.01,m=13267,m1=13267*0.1)

## [1] 0.04344502

est_power_fdr(exemplary.data,fdr=0.1,m=13267,m1=13267*0.1)

## [1] 0.1656871

exemplary.data=betaToData(Betas=c(log(2.5),log(2)),phi=0.6513,sampleSize=70,error=10^(-7))
dataToPower(exemplary.data)

## [1] 0.8195874

est_power_fdr(exemplary.data,fdr=0.1,m=13267)

## [1] 0.3282692