In PEER-paper they write:

> When RNA-seq estimates are used for transcript abundance, we recommend using DESeq to estimate library sizes and variance-stabilize expression data sets.

And I got confused by DESeq2 and DESeq. In DESeq (varianceStabili... depends on the estimateDispersions):

```r
cdsBlind = estimateDispersions( cds, method="blind" )
vsd = varianceStabilizingTransformation( cdsBlind )
```

In DESeq2:

    vsd <- varianceStabilizingTransformation(dds)


I am getting the same results from `varianceStabilizingTransformation()` if I use DESeq functions as `DESeq()`, `estimateSizeFactors()`, `estimateDispersions()` before or not.

```r
dds = DESeqDataSetFromMatrix(countData=histone_m, colData=Design, design=~condition)
pik = DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

vst=varianceStabilizingTransformation(pik)
head(assay(vst))
      WEN1     WEN2     WNN1     WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555

vst=varianceStabilizingTransformation(dds)
head(assay(vst))
      WEN1     WEN2     WNN1     WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555

cds=estimateSizeFactors(dds)
cds = estimateDispersions( cds)
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

vst=varianceStabilizingTransformation(cds)
head(assay(vst))
      WEN1     WEN2     WNN1     WNN2
1 8.568387 8.206238 7.843942 8.086705
2 8.962147 8.688945 8.320699 8.750248
3 8.135942 7.859804 7.843942 7.878520
4 7.321410 6.863988 7.165757 6.908507
5 6.474284 6.379227 6.420262 6.351544
6 6.561254 6.468225 6.555806 6.861555
```

So, it seems it does not matter if I use DESeq function before or not. Does it make sense? Also, once I estimate SizeFactors and varianceStabilizingTransformation, how should I plug it into PEER?
From my understanding, there is negligible difference between the Variance Stabilizing Transformation functions of DESeq and DESeq2. I would suggest you use DESeq2 because it is the most up-to-date version.

```
library(peer)

# assign vst data to variable called expr_set
expr_set = assay(vst)    

# the PEER tutorial recommends k == 25% of "number of samples"
k = ncol(expr_set) * .25

# create empty peer model
model = PEER()

# set number of "hidden factors" searched for
PEER_setNk(model,k)

# PEER ask NxG matrix, where N=samples and G=genes
PEER_setPhenoMean(model,as.matrix(t(expr_set)))

PEER_setAdd_mean(model, TRUE)

# These can be changed to the tolerance you desire
PEER_setTolerance(model, .001)
PEER_setVarTolerance(model, 0.0005)

# usually doesn't go over 200 iterations        
PEER_setNmax_iterations(model, 1000)

# run peer
PEER_update(model)
```