DESeq2: Create Data Objects¶

Objective: Create a DESeqDataSet object¶

Initial note¶

First step is to create a countData and colData object (see ? DESeqDataSet)

countData: for matrix input: a matrix of non-negative integers

colData: for matrix input: a ‘DataFrame’ or ‘data.frame’ with at least a single column. Rows of colData correspond to columns of countData

Load packages¶

[1]:

library(tidyverse)
library(DESeq2)

Registered S3 methods overwritten by 'ggplot2':
  method         from
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang
── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.1     ✔ purrr   0.3.2
✔ tibble  2.1.2     ✔ dplyr   0.8.1
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min


Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: DelayedArray
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following objects are masked from ‘package:Biobase’:

    anyMissing, rowMedians

The following object is masked from ‘package:dplyr’:

    count

Loading required package: BiocParallel

Attaching package: ‘DelayedArray’

The following objects are masked from ‘package:matrixStats’:

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

The following object is masked from ‘package:purrr’:

    simplify

The following objects are masked from ‘package:base’:

    aperm, apply, rowsum

Load the 2019 pilot count objects from the image file¶

[2]:

curdir <- "/home/jovyan/work/scratch/analysis_output"

imgdir <- file.path(curdir, "img")

imgfile <- file.path(imgdir, "pilotcnt2019.RData")

imgfile

attach(imgfile)

tools::md5sum(imgfile)

### List the objects that have been attached
ls(2)

cnt2019 <- cnt2019
mtdf2019 <- mtdf2019

detach(pos = 2)

'/home/jovyan/work/scratch/analysis_output/img/pilotcnt2019.RData'

/home/jovyan/work/scratch/analysis_output/img/pilotcnt2019.RData: '0d1cdfc8bb97baf52ca4f8f33e2c0ae4'

'cnt2019'
'mtdf2019'

Check dimensions of the two objects¶

[3]:

dim(cnt2019)
dim(mtdf2019)

24
8500

24
22

[4]:

mtdf2019 %>% head

A tibble: 6 × 22
Label	sample_year	group	enrich_rep	RNA_sample_num	genotype	condition	libprep_person	enrichment_method	enrichment_short	⋯	i5_primer	i7_primer	library_num	bio_replicate	Nanodrop_260_280	Nanodrop_260_230	Nanodrop_concentration_ng_ul	Bioanalyzer_concentration_ng_ul	RIN_normal_threshold	RIN_lowered_threshold
<chr>	<dbl>	<chr>	<dbl>	<dbl>	<chr>	<chr>	<chr>	<chr>	<chr>	⋯	<chr>	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<dbl>
1_2019_P_M1	2019	P	1	1	WT	pH4	C	mRNA	M	⋯	i501	i701	1	1	2.14	1.52	293	197	N/A	9.8
2_2019_P_M1	2019	P	1	2	WT	pH4	C	mRNA	M	⋯	i502	i701	2	2	2.12	1.79	290	225	N/A	9.9
3_2019_P_M1	2019	P	1	3	WT	pH4	C	mRNA	M	⋯	i503	i701	3	3	2.11	2.49	302	241	N/A	9.9
4_2019_P_M1	2019	P	1	4	WT	pH4	P	mRNA	M	⋯	i504	i701	4	4	2.13	1.15	296	189	N/A	9.7
5_2019_P_M1	2019	P	1	5	WT	pH4	P	mRNA	M	⋯	i505	i701	5	5	2.09	2.42	337	268	10	10.0
6_2019_P_M1	2019	P	1	6	WT	pH4	P	mRNA	M	⋯	i506	i701	6	6	2.08	2.40	319	276	10	10.0

[5]:

cnt2019[,1:5]

A tibble: 24 × 5
Label	CNAG_00001	CNAG_00002	CNAG_00003	CNAG_00004
<chr>	<int>	<int>	<int>	<int>
1_2019_P_M1	0	158	201	904
10_2019_P_M1	0	119	131	513
11_2019_P_M1	0	90	121	573
12_2019_P_M1	0	81	151	533
13_2019_P_M1	0	188	215	474
14_2019_P_M1	0	177	154	440
15_2019_P_M1	0	216	197	425
16_2019_P_M1	0	224	195	548
17_2019_P_M1	0	234	211	517
18_2019_P_M1	0	338	201	464
19_2019_P_M1	0	191	118	458
2_2019_P_M1	0	204	156	902
20_2019_P_M1	0	201	133	405
21_2019_P_M1	0	192	109	339
22_2019_P_M1	0	166	113	406
23_2019_P_M1	0	202	154	442
24_2019_P_M1	0	235	169	523
3_2019_P_M1	0	149	161	941
4_2019_P_M1	0	176	171	795
5_2019_P_M1	0	161	162	849
6_2019_P_M1	0	148	103	688
7_2019_P_M1	0	172	172	768
8_2019_P_M1	0	169	170	744
9_2019_P_M1	0	124	175	659

Create columnData object¶

[6]:

# columnData --- metadata
mtdf2019 %>%
    DataFrame ->
        columnData2019

### Add the labels as rownames
rownames(columnData2019) <- columnData2019[["Label"]]

columnData2019[, c("Label", "genotype", "condition")] %>% head

DataFrame with 6 rows and 3 columns
                  Label    genotype   condition
            <character> <character> <character>
1_2019_P_M1 1_2019_P_M1          WT         pH4
2_2019_P_M1 2_2019_P_M1          WT         pH4
3_2019_P_M1 3_2019_P_M1          WT         pH4
4_2019_P_M1 4_2019_P_M1          WT         pH4
5_2019_P_M1 5_2019_P_M1          WT         pH4
6_2019_P_M1 6_2019_P_M1          WT         pH4

[7]:

### Note that libraries are across rows and libraries across columns
### DESeq2 requires that the matrix is transposed as that the gene names become row names
cnt2019[1:4,1:5]

A tibble: 4 × 5
Label	CNAG_00001	CNAG_00002	CNAG_00003	CNAG_00004
<chr>	<int>	<int>	<int>	<int>
1_2019_P_M1	0	158	201	904
10_2019_P_M1	0	119	131	513
11_2019_P_M1	0	90	121	573
12_2019_P_M1	0	81	151	533

[8]:

### Transpose the count matrix (so that libraries are across the columns and genes across rows)
### Note that as.matrix() converts the tibble to a matrix object
cnt2019 %>%
    gather(key = gene, value = value, 2:ncol(cnt2019)) %>%
            spread_(key = names(cnt2019)[1],value = 'value') %>%
                column_to_rownames("gene") %>%
                    as.matrix ->
                        countData2019

countData2019[1:5, 1:6]

A matrix: 5 × 6 of type int
	1_2019_P_M1	10_2019_P_M1	11_2019_P_M1	12_2019_P_M1	13_2019_P_M1	14_2019_P_M1
CNAG_00001	0	0	0	0	0	0
CNAG_00002	158	119	90	81	188	177
CNAG_00003	201	131	121	151	215	154
CNAG_00004	904	513	573	533	474	440
CNAG_00005	22	24	18	20	25	13

Reorder the columns of the count matrix according to the order of Label in columnData¶

[9]:

### The two sets coincide
setequal(columnData2019[["Label"]], colnames(countData2019))

TRUE

[10]:

### but they are identical as they follow different orders
identical(columnData2019[["Label"]], colnames(countData2019))

FALSE

Reorder the columns of countData¶

[11]:

countData2019 <- countData2019[,columnData2019[["Label"]]]

[12]:

countData2019[1:4,1:5]

A matrix: 4 × 5 of type int
	1_2019_P_M1	2_2019_P_M1	3_2019_P_M1	4_2019_P_M1	5_2019_P_M1
CNAG_00001	0	0	0	0	0
CNAG_00002	158	204	149	176	161
CNAG_00003	201	156	161	171	162
CNAG_00004	904	902	941	795	849

make sure that labels match¶

[13]:

### The two variables coincide
setequal(columnData2019[["Label"]], colnames(countData2019))
### but they are identical as they follow different orders
identical(columnData2019[["Label"]], colnames(countData2019))

TRUE

Make DESeq object on the basis of the counts¶

The design option allows you to specify an additive or a multiplicitive model

Additive model

[14]:

dds_add <- DESeqDataSetFromMatrix(
    countData2019,                      # Count matrix
    columnData2019,                     # metadata
    ~ condition + genotype)             # design formula

Warning message in DESeqDataSet(se, design = design, ignoreRank):
“some variables in design formula are characters, converting to factors”

Inspect object

[15]:

dds_add

class: DESeqDataSet
dim: 8499 24
metadata(1): version
assays(1): counts
rownames(8499): CNAG_00001 CNAG_00002 ... large_MTrRNA small_MTrRNA
rowData names(0):
colnames(24): 1_2019_P_M1 2_2019_P_M1 ... 23_2019_P_M1 24_2019_P_M1
colData names(22): Label sample_year ... RIN_normal_threshold
  RIN_lowered_threshold

[16]:

slotNames(dds_add)

'design'
'dispersionFunction'
'rowRanges'
'colData'
'assays'
'NAMES'
'elementMetadata'
'metadata'

Check design

[17]:

dds_add@design

~condition + genotype

Check column data

[18]:

dds_add@colData

DataFrame with 24 rows and 22 columns
                    Label sample_year       group enrich_rep RNA_sample_num
              <character>   <numeric> <character>  <numeric>      <numeric>
1_2019_P_M1   1_2019_P_M1        2019           P          1              1
2_2019_P_M1   2_2019_P_M1        2019           P          1              2
3_2019_P_M1   3_2019_P_M1        2019           P          1              3
4_2019_P_M1   4_2019_P_M1        2019           P          1              4
5_2019_P_M1   5_2019_P_M1        2019           P          1              5
...                   ...         ...         ...        ...            ...
20_2019_P_M1 20_2019_P_M1        2019           P          1             20
21_2019_P_M1 21_2019_P_M1        2019           P          1             21
22_2019_P_M1 22_2019_P_M1        2019           P          1             22
23_2019_P_M1 23_2019_P_M1        2019           P          1             23
24_2019_P_M1 24_2019_P_M1        2019           P          1             24
             genotype condition libprep_person enrichment_method
             <factor>  <factor>    <character>       <character>
1_2019_P_M1        WT       pH4              C              mRNA
2_2019_P_M1        WT       pH4              C              mRNA
3_2019_P_M1        WT       pH4              C              mRNA
4_2019_P_M1        WT       pH4              P              mRNA
5_2019_P_M1        WT       pH4              P              mRNA
...               ...       ...            ...               ...
20_2019_P_M1    sre1d       pH8              C              mRNA
21_2019_P_M1    sre1d       pH8              C              mRNA
22_2019_P_M1    sre1d       pH8              P              mRNA
23_2019_P_M1    sre1d       pH8              P              mRNA
24_2019_P_M1    sre1d       pH8              P              mRNA
             enrichment_short    i7_index    i5_index   i5_primer   i7_primer
                  <character> <character> <character> <character> <character>
1_2019_P_M1                 M    ATTACTCG    AGGCTATA        i501        i701
2_2019_P_M1                 M    ATTACTCG    GCCTCTAT        i502        i701
3_2019_P_M1                 M    ATTACTCG    AGGATAGG        i503        i701
4_2019_P_M1                 M    ATTACTCG    TCAGAGCC        i504        i701
5_2019_P_M1                 M    ATTACTCG    CTTCGCCT        i505        i701
...                       ...         ...         ...         ...         ...
20_2019_P_M1                M    CGCTCATT    TCAGAGCC        i504        i703
21_2019_P_M1                M    CGCTCATT    CTTCGCCT        i505        i703
22_2019_P_M1                M    CGCTCATT    TAAGATTA        i506        i703
23_2019_P_M1                M    CGCTCATT    ACGTCCTG        i507        i703
24_2019_P_M1                M    CGCTCATT    GTCAGTAC        i508        i703
             library_num bio_replicate Nanodrop_260_280 Nanodrop_260_230
               <numeric>     <numeric>        <numeric>        <numeric>
1_2019_P_M1            1             1             2.14             1.52
2_2019_P_M1            2             2             2.12             1.79
3_2019_P_M1            3             3             2.11             2.49
4_2019_P_M1            4             4             2.13             1.15
5_2019_P_M1            5             5             2.09             2.42
...                  ...           ...              ...              ...
20_2019_P_M1          20             2             2.14             2.39
21_2019_P_M1          21             3             2.01             2.22
22_2019_P_M1          22             4             2.22             2.33
23_2019_P_M1          23             5             2.13             2.56
24_2019_P_M1          24             6             2.11             2.16
             Nanodrop_concentration_ng_ul Bioanalyzer_concentration_ng_ul
                                <numeric>                       <numeric>
1_2019_P_M1                           293                             197
2_2019_P_M1                           290                             225
3_2019_P_M1                           302                             241
4_2019_P_M1                           296                             189
5_2019_P_M1                           337                             268
...                                   ...                             ...
20_2019_P_M1                          172                             131
21_2019_P_M1                          397                             295
22_2019_P_M1                          567                             499
23_2019_P_M1                          254                             206
24_2019_P_M1                          215                             336
             RIN_normal_threshold RIN_lowered_threshold
                      <character>             <numeric>
1_2019_P_M1                   N/A                   9.8
2_2019_P_M1                   N/A                   9.9
3_2019_P_M1                   N/A                   9.9
4_2019_P_M1                   N/A                   9.7
5_2019_P_M1                    10                    10
...                           ...                   ...
20_2019_P_M1                  N/A                    10
21_2019_P_M1                  N/A                    10
22_2019_P_M1                   10                    10
23_2019_P_M1                  N/A                    10
24_2019_P_M1                   10                    10

Get count matrix

[19]:

counts(dds_add)[1:10,1:10]

A matrix: 10 × 10 of type int
	1_2019_P_M1	2_2019_P_M1	3_2019_P_M1	4_2019_P_M1	5_2019_P_M1	6_2019_P_M1	7_2019_P_M1	8_2019_P_M1	9_2019_P_M1	10_2019_P_M1
CNAG_00001	0	0	0	0	0	0	0	0	0	0
CNAG_00002	158	204	149	176	161	148	172	169	124	119
CNAG_00003	201	156	161	171	162	103	172	170	175	131
CNAG_00004	904	902	941	795	849	688	768	744	659	513
CNAG_00005	22	33	12	15	26	13	29	33	21	24
CNAG_00006	5964	4854	4362	4489	4368	4171	4859	4267	4239	3712
CNAG_00007	3119	3496	2628	2437	2498	2594	2505	2383	2086	2021
CNAG_00008	1481	1744	1602	1391	1433	1183	1313	1389	981	934
CNAG_00009	494	750	541	436	502	522	490	470	404	433
CNAG_00010	1527	1613	1564	1319	1286	1020	949	1227	757	682

Change design: multiplicative model

[20]:

dds_mult <- DESeqDataSetFromMatrix(
    countData2019,                       # Count matrix
    columnData2019,                      # metadata
    ~ condition + genotype + condition:genotype) # design formula

Warning message in DESeqDataSet(se, design = design, ignoreRank):
“some variables in design formula are characters, converting to factors”

In the following demonstration, we will use the additive model. The multiplicitive model will be illustrated in the appendix below.

[21]:

dds2019 <- dds_add

[22]:

curdir <- "/home/jovyan/work/scratch/analysis_output"
imgdir <- file.path(curdir, "img")

imgfile <- file.path(imgdir, "pilotdds2019.RData")

imgfile

'/home/jovyan/work/scratch/analysis_output/img/pilotdds2019.RData'

[23]:

save(dds2019, file = imgfile)
tools::md5sum(imgfile)

/home/jovyan/work/scratch/analysis_output/img/pilotdds2019.RData: '6632de5c8a2eed06d8c40c958e6f5d03'

[24]:

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS:   /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] DESeq2_1.24.0               SummarizedExperiment_1.14.0
 [3] DelayedArray_0.10.0         BiocParallel_1.18.0
 [5] matrixStats_0.54.0          Biobase_2.44.0
 [7] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0
 [9] IRanges_2.18.1              S4Vectors_0.22.0
[11] BiocGenerics_0.30.0         forcats_0.4.0
[13] stringr_1.4.0               dplyr_0.8.1
[15] purrr_0.3.2                 readr_1.3.1
[17] tidyr_0.8.3                 tibble_2.1.2
[19] ggplot2_3.1.1               tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] nlme_3.1-139           bitops_1.0-6           bit64_0.9-7
 [4] lubridate_1.7.4        RColorBrewer_1.1-2     httr_1.4.0
 [7] repr_1.0.1             tools_3.6.0            backports_1.1.4
[10] R6_2.4.0               rpart_4.1-15           DBI_1.0.0
[13] Hmisc_4.2-0            lazyeval_0.2.2         colorspace_1.4-1
[16] nnet_7.3-12            withr_2.1.2            tidyselect_0.2.5
[19] gridExtra_2.3          bit_1.1-14             compiler_3.6.0
[22] cli_1.1.0              rvest_0.3.4            htmlTable_1.13.1
[25] xml2_1.2.0             scales_1.0.0           checkmate_1.9.3
[28] genefilter_1.66.0      pbdZMQ_0.3-3           digest_0.6.19
[31] foreign_0.8-71         XVector_0.24.0         base64enc_0.1-3
[34] pkgconfig_2.0.2        htmltools_0.3.6        htmlwidgets_1.3
[37] rlang_0.3.4            readxl_1.3.1           RSQLite_2.1.1
[40] rstudioapi_0.10        generics_0.0.2         jsonlite_1.6
[43] acepack_1.4.1          RCurl_1.95-4.12        magrittr_1.5
[46] GenomeInfoDbData_1.2.1 Formula_1.2-3          Matrix_1.2-17
[49] Rcpp_1.0.1             IRkernel_1.0.1         munsell_0.5.0
[52] stringi_1.4.3          zlibbioc_1.30.0        plyr_1.8.4
[55] blob_1.1.1             grid_3.6.0             crayon_1.3.4
[58] lattice_0.20-38        IRdisplay_0.7.0        haven_2.1.0
[61] splines_3.6.0          annotate_1.62.0        hms_0.4.2
[64] locfit_1.5-9.1         zeallot_0.1.0          knitr_1.23
[67] pillar_1.4.1           uuid_0.1-2             geneplotter_1.62.0
[70] XML_3.98-1.19          glue_1.3.1             evaluate_0.14
[73] latticeExtra_0.6-28    data.table_1.12.2      modelr_0.1.4
[76] vctrs_0.1.0            cellranger_1.1.0       gtable_0.3.0
[79] assertthat_0.2.1       xfun_0.7               xtable_1.8-4
[82] broom_0.5.2            survival_2.44-1.1      memoise_1.1.0
[85] AnnotationDbi_1.46.0   cluster_2.0.8