The design formula in the resulting DESeqDataSet is just the sum of all variables in designmat from the ThinData object (except the intercept term). You should change this design formula if you want to study other models.

ThinDataToDESeqDataSet(obj)

Arguments

obj

A ThinData S3 object. This is generally output by either thin_diff, thin_2group, thin_lib, thin_gene, or thin_all.

Value

A DESeqDataSet S4 object. This will allow you to insert the simulated data directly into DESeq2.

Author

David Gerard

Examples

# \donttest{
## Generate simulated data and modify using thin_diff().
## In practice, you would use real data, not simulated.
set.seed(1)
n <- 10
p <- 1000
Z <- matrix(abs(rnorm(n, sd = 4)))
alpha <- matrix(abs(rnorm(p, sd = 1)))
mat <- round(2^(alpha %*% t(Z) + abs(matrix(rnorm(n * p, sd = 5),
                                            nrow = p,
                                            ncol = n))))
design_perm <- cbind(rep(c(0, 1), length.out = n), runif(n))
coef_perm   <- matrix(rnorm(p * ncol(design_perm), sd = 6), nrow = p)
design_obs  <- matrix(rnorm(n), ncol = 1)
target_cor <- matrix(c(0.9, 0))
thout <- thin_diff(mat            = mat,
                   design_perm    = design_perm,
                   coef_perm      = coef_perm,
                   target_cor     = target_cor,
                   design_obs     = design_obs,
                   permute_method = "hungarian")

## Convert ThinData object to DESeqDataSet object.
seobj <- ThinDataToDESeqDataSet(thout)
#> renaming the first element in assays to 'counts'
#>   the design formula contains one or more numeric variables with integer values,
#>   specifying a model with increasing fold change for higher values.
#>   did you mean for this to be a factor? if so, first convert
#>   this variable to a factor using the factor() function
class(seobj)
#> [1] "DESeqDataSet"
#> attr(,"package")
#> [1] "DESeq2"

## The "O1" variable in the colData corresponds to design_obs.
## The "P1" and "P2" variables in colData correspond to design_perm.
seobj
#> class: DESeqDataSet 
#> dim: 1000 10 
#> metadata(1): version
#> assays(1): counts
#> rownames(1000): gene1 gene2 ... gene999 gene1000
#> rowData names(2): true_P1 true_P2
#> colnames(10): sample1 sample2 ... sample9 sample10
#> colData names(3): O1 P1 P2
# }