DESeq2: Basic Differential Expression (DE) analysis

Objective: Carry out a basic set of DE analyses using DESeq2 and visualize the results

Load packages

[1]:

library(tidyverse)
library(DESeq2)
library(dendextend)
library(RColorBrewer)

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


---------------------
Welcome to dendextend version 1.12.0
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
Or contact: <tal.galili@gmail.com>

        To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------


Attaching package: ‘dendextend’

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

    cutree

Load the 2019 pilot dds object from image file

[2]:


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

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

imgfile

attach(imgfile)

tools::md5sum(imgfile)

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

dds2019 <- dds2019

detach(pos = 2)
'/home/jovyan/work/scratch/analysis_output/img/pilotdds2019.RData'
/home/jovyan/work/scratch/analysis_output/img/pilotdds2019.RData: '6632de5c8a2eed06d8c40c958e6f5d03'
'dds2019'

Check dimensions of the two objects

Inspect object & Slots of an S4 class

Let’s has a look at the object we have created.

[3]:

dds2019

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

see the class of dds object

[4]:

class(dds2019)
'DESeqDataSet'

DESeqDataSet is a S4 object. Recall that a S4 object was taught when introducing bioconductor. Note that S4 objects allow users to wrap up multiple elements into a single variables where each element is called a slot.

[5]:

slotNames(dds2019)
  1. 'design'
  2. 'dispersionFunction'
  3. 'rowRanges'
  4. 'colData'
  5. 'assays'
  6. 'NAMES'
  7. 'elementMetadata'
  8. 'metadata'

The metadata (columnData) is stored in the slot colData

[6]:

dds2019@colData %>% as.data.frame %>% head(3)
A data.frame: 3 × 22
Labelsample_yeargroupenrich_repRNA_sample_numgenotypeconditionlibprep_personenrichment_methodenrichment_short⋯i5_primeri7_primerlibrary_numbio_replicateNanodrop_260_280Nanodrop_260_230Nanodrop_concentration_ng_ulBioanalyzer_concentration_ng_ulRIN_normal_thresholdRIN_lowered_threshold
<chr><dbl><chr><dbl><dbl><fct><fct><chr><chr><chr>⋯<chr><chr><dbl><dbl><dbl><dbl><dbl><dbl><chr><dbl>
1_2019_P_M11_2019_P_M12019P11WTpH4CmRNAM⋯i501i701112.141.52293197N/A9.8
2_2019_P_M12_2019_P_M12019P12WTpH4CmRNAM⋯i502i701222.121.79290225N/A9.9
3_2019_P_M13_2019_P_M12019P13WTpH4CmRNAM⋯i503i701332.112.49302241N/A9.9

The design formula is stored in the slot design. The design holds the R formula which expresses how the counts depend on the variables in colData.

[7]:

dds2019@design

~condition + genotype

The first thing you may want to do is have a look at the raw counts you have imported. The DESeq2::counts function extracts a matrix of counts (with the genes along the rows and samples along the columns). Let us first verify the dimension of this matrix.

[8]:

dim(counts(dds2019))
  1. 8499
  2. 24
[9]:

head(counts(dds2019),3)
A matrix: 3 × 24 of type int
1_2019_P_M12_2019_P_M13_2019_P_M14_2019_P_M15_2019_P_M16_2019_P_M17_2019_P_M18_2019_P_M19_2019_P_M110_2019_P_M1⋯15_2019_P_M116_2019_P_M117_2019_P_M118_2019_P_M119_2019_P_M120_2019_P_M121_2019_P_M122_2019_P_M123_2019_P_M124_2019_P_M1
CNAG_00001 0 0 0 0 0 0 0 0 0 0⋯ 0 0 0 0 0 0 0 0 0 0
CNAG_00002158204149176161148172169124119⋯216224234338191201192166202235
CNAG_00003201156161171162103172170175131⋯197195211201118133109113154169

This slot returns gene specific information (it will be populated later)

[10]:

dds2019@dispersionFunction

function ()
NULL

Estimate Size Factors and Dispersion Parameters

You recall that DESeq requires that we have estimates for sample specific size factors and gene specific dispersion factors. More specifically, recall that DESeq models the count \(K_{ij}\) (gene \(i\), sample \(j\)) as negative binomial with mean \(\mu_{ij}\) and dispersion parameter \(\alpha_i\). Here \(\mu_{ij}=s_j q_{ij}\) where \(\log_2(q_{ij}) = \beta_{0i} + \beta_{1i} z_j\). Here \(s_j\) is the sample \(j\) specific size factor.

Summarize of notation - \(K_{ij}\) denotes the observed number of reads mapped to gene \(i\) for sample \(j\) - \(K_{ij}\) follows a negative binomial distribution with - Mean \(\mu_{ij}\) - Dispersion parameter \(\alpha_i\) - Modelling - \(K_{ij} \sim NB(\mu_{ij}, \alpha_i)\) - \(\mu_{ij} = s_{j}q_{ij}\) - \(s_j\) is sample \(j\) specific normalization constant - \(\log_2(q_{ij}) = \beta_{0i} + \beta_{1i} z_j\)

01 Size Factors

We begin by estimating the size factors \(s_1,\ldots,s_n\):

[11]:

dds2019 <- estimateSizeFactors(dds2019)

Now, compare the dds object to that of before applying the estimateSizeFactors() function. What has changed? What remains unchanged?

[12]:

dds2019

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(23): Label sample_year ... RIN_lowered_threshold
  sizeFactor

Note that there is a sizeFactor added to colData. Let’s look at it more carefully

``` > dds # (before estimateSizeFactors) class: DESeqDataSet 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

dds # (after estimateSizeFactors) 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(23): Label sample_year … RIN_lowered_threshold sizeFactor

You can also get the size factors directly

[13]:

sizeFactors(dds2019)
1_2019_P_M1
1.5916393221245
2_2019_P_M1
1.65871103315141
3_2019_P_M1
1.51289455728943
4_2019_P_M1
1.41966861598247
5_2019_P_M1
1.43738473025559
6_2019_P_M1
1.23710002282889
7_2019_P_M1
1.40148291158442
8_2019_P_M1
1.44217254922035
9_2019_P_M1
1.17444094436861
10_2019_P_M1
1.07636468283386
11_2019_P_M1
1.04216270302148
12_2019_P_M1
0.971238282275506
13_2019_P_M1
0.770986843971703
14_2019_P_M1
0.752382049912615
15_2019_P_M1
0.717679610028228
16_2019_P_M1
0.897840290244175
17_2019_P_M1
0.89364164432414
18_2019_P_M1
0.931618100592125
19_2019_P_M1
0.749435710766038
20_2019_P_M1
0.775571527912074
21_2019_P_M1
0.583832529967166
22_2019_P_M1
0.634399759403834
23_2019_P_M1
0.755186085659437
24_2019_P_M1
0.808836405192349

It is preferable to limit the number of decimal places. Next show the size factors rounded to 3 decimal places

[14]:

round(sizeFactors(dds2019),3)
1_2019_P_M1
1.592
2_2019_P_M1
1.659
3_2019_P_M1
1.513
4_2019_P_M1
1.42
5_2019_P_M1
1.437
6_2019_P_M1
1.237
7_2019_P_M1
1.401
8_2019_P_M1
1.442
9_2019_P_M1
1.174
10_2019_P_M1
1.076
11_2019_P_M1
1.042
12_2019_P_M1
0.971
13_2019_P_M1
0.771
14_2019_P_M1
0.752
15_2019_P_M1
0.718
16_2019_P_M1
0.898
17_2019_P_M1
0.894
18_2019_P_M1
0.932
19_2019_P_M1
0.749
20_2019_P_M1
0.776
21_2019_P_M1
0.584
22_2019_P_M1
0.634
23_2019_P_M1
0.755
24_2019_P_M1
0.809

Summarize size factors

[15]:

summary(sizeFactors(dds2019))

  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.5838  0.7670  0.9514  1.0515  1.4060  1.6587

Do you see a trend?

[16]:

sizeFactors(dds2019) %>%
    as.data.frame %>%
        rownames_to_column %>%
            mutate(libnum = as.integer(str_remove(rowname, "_2019_P_M1"))) -> mydf

colnames(mydf)[2] <- "sizefac"

mydf
A data.frame: 24 × 3
rownamesizefaclibnum
<chr><dbl><int>
1_2019_P_M1 1.5916393 1
2_2019_P_M1 1.6587110 2
3_2019_P_M1 1.5128946 3
4_2019_P_M1 1.4196686 4
5_2019_P_M1 1.4373847 5
6_2019_P_M1 1.2371000 6
7_2019_P_M1 1.4014829 7
8_2019_P_M1 1.4421725 8
9_2019_P_M1 1.1744409 9
10_2019_P_M11.076364710
11_2019_P_M11.042162711
12_2019_P_M10.971238312
13_2019_P_M10.770986813
14_2019_P_M10.752382014
15_2019_P_M10.717679615
16_2019_P_M10.897840316
17_2019_P_M10.893641617
18_2019_P_M10.931618118
19_2019_P_M10.749435719
20_2019_P_M10.775571520
21_2019_P_M10.583832521
22_2019_P_M10.634399822
23_2019_P_M10.755186123
24_2019_P_M10.808836424
 
[17]:

ggplot(mydf, aes(x = libnum, y = sizefac)) + geom_point()
../_images/DESeq2_02_deseq_DE_38_0.png

Now that the size factors have been estimated, we can get “normalized” counts (DESeq2 normalizes against size factor)

[18]:

# original counts for libraries 1 and 24
counts(dds2019)[1:5,c(1,24)]

# normalized count
counts(dds2019, normalize = TRUE)[1:5, c(1,24)]

# Size factor

sizeFactors(dds2019)[c(1,24)]
A matrix: 5 × 2 of type int
1_2019_P_M124_2019_P_M1
CNAG_00001 0 0
CNAG_00002158235
CNAG_00003201169
CNAG_00004904523
CNAG_00005 22 40
A matrix: 5 × 2 of type dbl
1_2019_P_M124_2019_P_M1
CNAG_00001 0.00000 0.00000
CNAG_00002 99.26872290.54083
CNAG_00003126.28489208.94213
CNAG_00004567.96787646.60789
CNAG_00005 13.82223 49.45376
1_2019_P_M1
1.5916393221245
24_2019_P_M1
0.808836405192349
[19]:

# normalized manually using size factors for library 1
counts(dds2019)[1:5, 1] / sizeFactors(dds2019)[1]
CNAG_00001
0
CNAG_00002
99.2687211252757
CNAG_00003
126.284892064433
CNAG_00004
567.9678727674
CNAG_00005
13.822226992127
[20]:

# normalized manually using size factors for library 24
counts(dds2019)[1:5, 24] / sizeFactors(dds2019)[24]
CNAG_00001
0
CNAG_00002
290.54082938331
CNAG_00003
208.942128365019
CNAG_00004
646.607888372217
CNAG_00005
49.4537581929038

How do you get the raw counts for gene “GeneID: CNAG_05845”?

[21]:

counts(dds2019, normalize = TRUE)["CNAG_05845",]
1_2019_P_M1
1753.53797886484
2_2019_P_M1
1483.68217900155
3_2019_P_M1
1395.99946990619
4_2019_P_M1
1526.41255544016
5_2019_P_M1
1540.99313383281
6_2019_P_M1
1662.7596492127
7_2019_P_M1
1483.42871883441
8_2019_P_M1
1396.50418466832
9_2019_P_M1
1611.83072599508
10_2019_P_M1
1532.00864567435
11_2019_P_M1
1625.46596139805
12_2019_P_M1
1693.71412764532
13_2019_P_M1
1629.08097566201
14_2019_P_M1
1600.25083019968
15_2019_P_M1
1574.51874654145
16_2019_P_M1
1551.50087953968
17_2019_P_M1
1732.23798357198
18_2019_P_M1
1523.15632242235
19_2019_P_M1
1822.70471019287
20_2019_P_M1
1799.96293025118
21_2019_P_M1
1872.1121963887
22_2019_P_M1
1692.93885137862
23_2019_P_M1
1832.66088488836
24_2019_P_M1
1791.46239053794

02 Dispersion Parameters

Next, we get the dispersion factors \(\alpha_1,\ldots,\alpha_{m}\)

[22]:

dds2019 <- estimateDispersions(dds2019)

gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
Now inspect the dds object again and note that the rowRanges slot has extra information (“metadata column names(0):” before versus “column names(9): baseMean baseVar … dispOutlier dispMAP”) - before: - metadata column names(0): - after:
- column names(9): baseMean baseVar ...
[23]:

dds2019

class: DESeqDataSet
dim: 8499 24
metadata(1): version
assays(2): counts mu
rownames(8499): CNAG_00001 CNAG_00002 ... large_MTrRNA small_MTrRNA
rowData names(10): baseMean baseVar ... dispOutlier dispMAP
colnames(24): 1_2019_P_M1 2_2019_P_M1 ... 23_2019_P_M1 24_2019_P_M1
colData names(23): Label sample_year ... RIN_lowered_threshold
  sizeFactor

Can you notice the difference?

> dds (before dispersion)
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(23): Label sample_year ... RIN_lowered_threshold
  sizeFactor

> dds (after dispersion)
class: DESeqDataSet
dim: 8499 24
metadata(1): version
assays(2): counts mu
rownames(8499): CNAG_00001 CNAG_00002 ... large_MTrRNA small_MTrRNA
rowData names(10): baseMean baseVar ... dispOutlier dispMAP
colnames(24): 1_2019_P_M1 2_2019_P_M1 ... 23_2019_P_M1 24_2019_P_M1
colData names(23): Label sample_year ... RIN_lowered_threshold
  sizeFactor

Note that the dispersionfunction slot is now populated

[24]:

dds2019@dispersionFunction

structure(function (q)
coefs[1] + coefs[2]/q, coefficients = c(asymptDisp = 0.0141803487671824,
extraPois = 0.87630661036897), fitType = "parametric", varLogDispEsts = 0.902466256228483, dispPriorVar = 0.802549300169356)

We can extract the gene specific dispersion factors using dispersions(). Note that there will be one number per gene. We look at the first four genes (rounded to 4 decimal places)

[25]:

alphas <- dispersions(dds2019)

Verify that the number of dispersion factors equals the number of genes

[26]:

# number of disperion factors
length(alphas)
8499
[27]:

round(alphas[1:4], 4)
  1. <NA>
  2. 0.0132
  3. 0.0231
  4. 0.0044

Extract the metadata using mcols() for the first four genes

Terms

Description

baseMean

mean of normalized counts for all samples

baseVar

variance of normalized counts for all samples

allZero

all counts for a gene are zero

dispGeneEst

gene-wise estimates of dispersion

dispFit

fitted values of dispersion

dispersion

final estimate of dispersion

dispIter

number of iterations

dispOut

dispersion flagged as outlier

dispMAP

maximum a posteriori estimate

 
[28]:

mcols(dds2019)[1:4,] %>% as.data.frame
A data.frame: 4 × 10
baseMeanbaseVarallZerodispGeneEstdispGeneIterdispFitdispersiondispIterdispOutlierdispMAP
<dbl><dbl><lgl><dbl><dbl><dbl><dbl><dbl><lgl><dbl>
CNAG_00001 0.0000 0.000 TRUE NANA NA NANA NA NA
CNAG_00002192.45648133.183FALSE0.012090962 90.018733620.013202310 9FALSE0.013202310
CNAG_00003164.99083089.489FALSE0.023820456100.019491590.023065791 7FALSE0.023065791
CNAG_00004571.05471819.185FALSE0.003107148 50.015714890.00443058612FALSE0.004430586

Summarize the dispersion factors using a box plot (may want to log transform)

[29]:

boxplot(log(dispersions(dds2019)))
../_images/DESeq2_02_deseq_DE_60_0.png

Differential Expression Analysis

We can now conduct a differential expression analysis using the DESeq() function. Keep in mind that to get to this step, we first estimated the size factors and then the dispersion parameters.

[30]:

### Carry out DE analysis
ddsDE <- DESeq(dds2019)

using pre-existing size factors
estimating dispersions
found already estimated dispersions, replacing these
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

[31]:

### Look at object
ddsDE

class: DESeqDataSet
dim: 8499 24
metadata(1): version
assays(4): counts mu H cooks
rownames(8499): CNAG_00001 CNAG_00002 ... large_MTrRNA small_MTrRNA
rowData names(26): baseMean baseVar ... deviance maxCooks
colnames(24): 1_2019_P_M1 2_2019_P_M1 ... 23_2019_P_M1 24_2019_P_M1
colData names(23): Label sample_year ... RIN_lowered_threshold
  sizeFactor

[32]:

### Look at some of the results
results(ddsDE)

log2 fold change (MLE): genotype WT vs sre1d
Wald test p-value: genotype WT vs sre1d
DataFrame with 8499 rows and 6 columns
                        baseMean     log2FoldChange              lfcSE
                       <numeric>          <numeric>          <numeric>
CNAG_00001                     0                 NA                 NA
CNAG_00002      192.456384076694 0.0459290105698189 0.0812849125012254
CNAG_00003      164.990819436254 0.0432979308487266  0.101045100070341
CNAG_00004      571.054651467718 0.0347715410374543 0.0464312203030827
CNAG_00005      28.7106195205257 -0.377056805602812  0.192347436642309
...                          ...                ...                ...
ENSRNA049551942                0                 NA                 NA
ENSRNA049551964                0                 NA                 NA
ENSRNA049551993                0                 NA                 NA
large_MTrRNA    5333.97356461621 -0.396337460272526   0.13434088145007
small_MTrRNA    1705.75472334291 -0.436802837852531  0.295538981353034
                             stat              pvalue                padj
                        <numeric>           <numeric>           <numeric>
CNAG_00001                     NA                  NA                  NA
CNAG_00002      0.565037337883909   0.572048368119426   0.675469542016723
CNAG_00003      0.428501043777338   0.668286373790347   0.756450803050288
CNAG_00004      0.748882773497678   0.453927864041203   0.566263547566551
CNAG_00005      -1.96029025488907  0.0499618744055703   0.094333851813714
...                           ...                 ...                 ...
ENSRNA049551942                NA                  NA                  NA
ENSRNA049551964                NA                  NA                  NA
ENSRNA049551993                NA                  NA                  NA
large_MTrRNA    -2.95023715785154 0.00317530091247218 0.00896807960414441
small_MTrRNA     -1.4779872213566   0.139411198621898   0.221714577729148

Note that currently, the model we have is an additive model, which does not include the interaction term of Media and Strain

Look at some of the results (tidy version)

[33]:

results(ddsDE, tidy = TRUE)
A data.frame: 8499 × 7
rowbaseMeanlog2FoldChangelfcSEstatpvaluepadj
<chr><dbl><dbl><dbl><dbl><dbl><dbl>
CNAG_00001 0.00000 NA NA NA NA NA
CNAG_00002 192.45638 0.04592901060.08128491 0.5650373385.720484e-016.754695e-01
CNAG_00003 164.99082 0.04329793080.10104510 0.4285010446.682864e-017.564508e-01
CNAG_00004 571.05465 0.03477154100.04643122 0.7488827734.539279e-015.662635e-01
CNAG_00005 28.71062-0.37705680560.19234744 -1.9602902554.996187e-029.433385e-02
CNAG_000063147.39964-0.13136969270.04955489 -2.6509934578.025540e-031.995982e-02
CNAG_000071484.87545-0.25112143440.08051906 -3.1187824551.816000e-035.492473e-03
CNAG_00008 938.01932-0.00957755680.04760749 -0.2011775048.405598e-018.873725e-01
CNAG_00009 305.19881-0.18664588700.09952871 -1.8752969226.075189e-021.116044e-01
CNAG_00010 989.04227 0.29910115060.06322330 4.7308687972.235610e-061.474309e-05
CNAG_000111001.75238-0.01506299410.08611946 -0.1749081338.611518e-019.032619e-01
CNAG_00012 773.99385 0.19474879020.05765064 3.3780853377.299243e-042.477973e-03
CNAG_00013 485.75319-0.00028230690.07949652 -0.0035511859.971666e-019.974113e-01
CNAG_00014 180.53988 0.04893072500.06795183 0.7200795534.714760e-015.825377e-01
CNAG_00015 318.87538-0.52256276610.04948141-10.5607897344.528359e-266.591188e-24
CNAG_00016 901.25567-0.17439720990.04458282 -3.9117585499.162649e-053.953666e-04
CNAG_00017 184.65031 0.25011413670.08576066 2.9164204443.540730e-039.829867e-03
CNAG_00018 876.12980 0.05336718550.03737492 1.4278876781.533242e-012.399204e-01
CNAG_000191208.90542-0.03761997930.05412226 -0.6950924854.869973e-015.980888e-01
CNAG_00020 968.39030 0.16334678170.06202149 2.6337124988.445697e-032.079157e-02
CNAG_00021 539.92676 0.34559598600.05392145 6.4092492931.462379e-102.337226e-09
CNAG_00022 746.84601-0.03380024950.06149652 -0.5496286335.825741e-016.837382e-01
CNAG_00023 278.87826-1.46009376560.26611238 -5.4867562794.093813e-083.977196e-07
CNAG_000247703.04647-0.10749591950.04323925 -2.4860728151.291615e-023.001984e-02
CNAG_000252021.99612-0.09663889950.09609375 -1.0056730593.145729e-014.265652e-01
CNAG_000262237.77134-0.22179505750.06090899 -3.6414173392.711412e-041.034194e-03
CNAG_00027 425.60765 0.12268413780.06051990 2.0271700384.264503e-028.258484e-02
CNAG_00028 381.48552-0.10529428930.05837561 -1.8037376927.127244e-021.274834e-01
CNAG_00029 52.66523-0.06987920340.10076540 -0.6934841354.880058e-015.990298e-01
CNAG_00030 76.57328-0.32881123370.12296567 -2.6740084517.495056e-031.881006e-02
⋮⋮⋮⋮⋮⋮⋮
ENSRNA049550980 0.000 NA NA NA NA NA
ENSRNA049551037 0.000 NA NA NA NA NA
ENSRNA049551063 0.000 NA NA NA NA NA
ENSRNA049551074 0.000 NA NA NA NA NA
ENSRNA049551126 0.000 NA NA NA NA NA
ENSRNA049551197 0.000 NA NA NA NA NA
ENSRNA049551244 0.000 NA NA NA NA NA
ENSRNA049551269 0.000 NA NA NA NA NA
ENSRNA049551298 0.000 NA NA NA NA NA
ENSRNA049551337 0.000 NA NA NA NA NA
ENSRNA049551365 0.000 NA NA NA NA NA
ENSRNA049551391 0.000 NA NA NA NA NA
ENSRNA049551484 0.000 NA NA NA NA NA
ENSRNA049551550 0.000 NA NA NA NA NA
ENSRNA049551574 0.000 NA NA NA NA NA
ENSRNA049551598 0.000 NA NA NA NA NA
ENSRNA049551636 0.000 NA NA NA NA NA
ENSRNA049551673 0.000 NA NA NA NA NA
ENSRNA049551717 0.000 NA NA NA NA NA
ENSRNA049551745 0.000 NA NA NA NA NA
ENSRNA049551774 0.000 NA NA NA NA NA
ENSRNA049551798 0.000 NA NA NA NA NA
ENSRNA049551814 0.000 NA NA NA NA NA
ENSRNA049551862 0.000 NA NA NA NA NA
ENSRNA049551899 0.000 NA NA NA NA NA
ENSRNA049551942 0.000 NA NA NA NA NA
ENSRNA049551964 0.000 NA NA NA NA NA
ENSRNA049551993 0.000 NA NA NA NA NA
large_MTrRNA 5333.974-0.39633750.1343409-2.9502370.0031753010.00896808
small_MTrRNA 1705.755-0.43680280.2955390-1.4779870.1394111990.22171458

We can get the results for the differential expression analysis using results(). Here, we can compare two group of samples specified by the contrast. (If not, the default contrast would be the last term in your additive model design(dds)).

[34]:

# DE with respect to condition
myres_condition4v8 <- results(ddsDE, contrast = c("condition", "pH4", "pH8"))
myres_condition4v8

log2 fold change (MLE): condition pH4 vs pH8
Wald test p-value: condition pH4 vs pH8
DataFrame with 8499 rows and 6 columns
                        baseMean      log2FoldChange              lfcSE
                       <numeric>           <numeric>          <numeric>
CNAG_00001                     0                  NA                 NA
CNAG_00002      192.456384076694   -1.34351851614176 0.0813854454817818
CNAG_00003      164.990819436254  -0.831905577724638  0.101018400723969
CNAG_00004      571.054651467718 -0.0871138906737344 0.0465166816365353
CNAG_00005      28.7106195205257   -1.29275886849841  0.192729463777133
...                          ...                 ...                ...
ENSRNA049551942                0                  NA                 NA
ENSRNA049551964                0                  NA                 NA
ENSRNA049551993                0                  NA                 NA
large_MTrRNA    5333.97356461621  -0.948564278537808  0.134340883270906
small_MTrRNA    1705.75472334291  -0.691472418621256  0.295538979277226
                             stat               pvalue                 padj
                        <numeric>            <numeric>            <numeric>
CNAG_00001                     NA                   NA                   NA
CNAG_00002      -16.5080931631996 3.20862403991611e-61 2.87085560366149e-60
CNAG_00003      -8.23518855735805 1.79275046151211e-16 4.79576928512805e-16
CNAG_00004      -1.87274516601186   0.0611035831980367    0.074895534834165
CNAG_00005      -6.70763485334717 1.97804137561036e-11 4.32020773113613e-11
...                           ...                  ...                  ...
ENSRNA049551942                NA                   NA                   NA
ENSRNA049551964                NA                   NA                   NA
ENSRNA049551993                NA                   NA                   NA
large_MTrRNA    -7.06087570248421 1.65456549165796e-12 3.75246614426935e-12
small_MTrRNA    -2.33969955608675   0.0192992582107572   0.0250291573072207

[35]:

# DE with respect to condition (flip order)
myres_condition8v4 <- results(ddsDE, contrast = c("condition", "pH8", "pH4"))
myres_condition8v4

log2 fold change (MLE): condition pH8 vs pH4
Wald test p-value: condition pH8 vs pH4
DataFrame with 8499 rows and 6 columns
                        baseMean     log2FoldChange              lfcSE
                       <numeric>          <numeric>          <numeric>
CNAG_00001                     0                 NA                 NA
CNAG_00002      192.456384076694   1.34351851614176 0.0813854454817818
CNAG_00003      164.990819436254  0.831905577724638  0.101018400723969
CNAG_00004      571.054651467718 0.0871138906737344 0.0465166816365353
CNAG_00005      28.7106195205257   1.29275886849841  0.192729463777133
...                          ...                ...                ...
ENSRNA049551942                0                 NA                 NA
ENSRNA049551964                0                 NA                 NA
ENSRNA049551993                0                 NA                 NA
large_MTrRNA    5333.97356461621  0.948564278537808  0.134340883270906
small_MTrRNA    1705.75472334291  0.691472418621256  0.295538979277226
                            stat               pvalue                 padj
                       <numeric>            <numeric>            <numeric>
CNAG_00001                    NA                   NA                   NA
CNAG_00002      16.5080931631996 3.20862403991611e-61 2.87085560366149e-60
CNAG_00003      8.23518855735805 1.79275046151211e-16 4.79576928512805e-16
CNAG_00004      1.87274516601186   0.0611035831980367    0.074895534834165
CNAG_00005      6.70763485334717 1.97804137561036e-11 4.32020773113613e-11
...                          ...                  ...                  ...
ENSRNA049551942               NA                   NA                   NA
ENSRNA049551964               NA                   NA                   NA
ENSRNA049551993               NA                   NA                   NA
large_MTrRNA    7.06087570248421 1.65456549165796e-12 3.75246614426935e-12
small_MTrRNA    2.33969955608675   0.0192992582107572   0.0250291573072207

[36]:

### DE with respect to genotype
myres_strainvWT <- results(ddsDE, contrast = c("genotype", "sre1d", "WT"))
myres_strainvWT

log2 fold change (MLE): genotype sre1d vs WT
Wald test p-value: genotype sre1d vs WT
DataFrame with 8499 rows and 6 columns
                        baseMean      log2FoldChange              lfcSE
                       <numeric>           <numeric>          <numeric>
CNAG_00001                     0                  NA                 NA
CNAG_00002      192.456384076694 -0.0459290105698189 0.0812849125012254
CNAG_00003      164.990819436254 -0.0432979308487266  0.101045100070341
CNAG_00004      571.054651467718 -0.0347715410374543 0.0464312203030827
CNAG_00005      28.7106195205257   0.377056805602812  0.192347436642309
...                          ...                 ...                ...
ENSRNA049551942                0                  NA                 NA
ENSRNA049551964                0                  NA                 NA
ENSRNA049551993                0                  NA                 NA
large_MTrRNA    5333.97356461621   0.396337460272526   0.13434088145007
small_MTrRNA    1705.75472334291   0.436802837852531  0.295538981353034
                              stat              pvalue                padj
                         <numeric>           <numeric>           <numeric>
CNAG_00001                      NA                  NA                  NA
CNAG_00002      -0.565037337883909   0.572048368119426   0.675469542016723
CNAG_00003      -0.428501043777338   0.668286373790347   0.756450803050288
CNAG_00004      -0.748882773497678   0.453927864041203   0.566263547566551
CNAG_00005        1.96029025488907  0.0499618744055703   0.094333851813714
...                            ...                 ...                 ...
ENSRNA049551942                 NA                  NA                  NA
ENSRNA049551964                 NA                  NA                  NA
ENSRNA049551993                 NA                  NA                  NA
large_MTrRNA      2.95023715785154 0.00317530091247218 0.00896807960414441
small_MTrRNA       1.4779872213566   0.139411198621898   0.221714577729148

Let’s look at the results for the first four genes

[37]:

### Tidy the results
myres_condition8v4 <- results(ddsDE, contrast = c("condition", "pH8", "pH4"), tidy = TRUE)
myres_condition8v4
A data.frame: 8499 × 7
rowbaseMeanlog2FoldChangelfcSEstatpvaluepadj
<chr><dbl><dbl><dbl><dbl><dbl><dbl>
CNAG_00001 0.00000 NA NA NA NA NA
CNAG_00002 192.45638 1.3435185160.08138545 16.5080932 3.208624e-61 2.870856e-60
CNAG_00003 164.99082 0.8319055780.10101840 8.2351886 1.792750e-16 4.795769e-16
CNAG_00004 571.05465 0.0871138910.04651668 1.8727452 6.110358e-02 7.489553e-02
CNAG_00005 28.71062 1.2927588680.19272946 6.7076349 1.978041e-11 4.320208e-11
CNAG_000063147.39964-0.1393277750.04956000 -2.8112951 4.934250e-03 6.756102e-03
CNAG_000071484.87545-0.7479203280.08053624 -9.2867544 1.590647e-20 4.927921e-20
CNAG_00008 938.01932-0.0791101960.04765929 -1.6599115 9.693226e-02 1.159517e-01
CNAG_00009 305.19881-0.5588380660.09968395 -5.6060989 2.069377e-08 3.933651e-08
CNAG_00010 989.04227 0.5378936680.06321490 8.5089694 1.754861e-17 4.868574e-17
CNAG_000111001.75238-2.1488943670.08634115-24.88841429.932146e-1372.998405e-135
CNAG_00012 773.99385 0.1614128170.05766095 2.7993437 5.120661e-03 6.994889e-03
CNAG_00013 485.75319-0.7229664970.07965037 -9.0767500 1.118652e-19 3.374959e-19
CNAG_00014 180.53988-0.2958865600.06870554 -4.3065895 1.657908e-05 2.684468e-05
CNAG_00015 318.87538-0.3285995880.05019105 -6.5469752 5.871409e-11 1.252823e-10
CNAG_00016 901.25567 0.3127175960.04460378 7.0110118 2.366008e-12 5.340717e-12
CNAG_00017 184.65031 0.3166832430.08573971 3.6935425 2.211516e-04 3.341873e-04
CNAG_00018 876.12980 0.1552110010.03742605 4.1471376 3.366577e-05 5.348074e-05
CNAG_000191208.90542 0.7471236790.05412045 13.8048310 2.383308e-43 1.414883e-42
CNAG_00020 968.39030-0.4701641250.06206999 -7.5747418 3.598425e-14 8.737196e-14
CNAG_00021 539.92676 0.0079220360.05395599 0.1468240 8.832709e-01 8.975865e-01
CNAG_00022 746.84601 0.0202990740.06152583 0.3299277 7.414546e-01 7.689054e-01
CNAG_00023 278.87826 0.2320639630.26610452 0.8720782 3.831657e-01 4.210277e-01
CNAG_000247703.04647 0.9970433840.04323929 23.05873811.202071e-1172.760023e-116
CNAG_000252021.99612 0.8810554030.09609371 9.1687104 4.787115e-20 1.463607e-19
CNAG_000262237.77134 0.1674335120.06091140 2.7488043 5.981308e-03 8.132384e-03
CNAG_00027 425.60765-0.9296469980.06112216-15.2096555 3.051445e-52 2.273522e-51
CNAG_00028 381.48552-0.8174229400.05910049-13.8310684 1.655433e-43 9.892547e-43
CNAG_00029 52.66523-3.3160308810.16631275-19.9385244 1.885448e-88 2.729713e-87
CNAG_00030 76.57328-1.1036474980.12592110 -8.7645955 1.874484e-18 5.414216e-18
⋮⋮⋮⋮⋮⋮⋮
ENSRNA049550980 0.000 NA NA NA NA NA
ENSRNA049551037 0.000 NA NA NA NA NA
ENSRNA049551063 0.000 NA NA NA NA NA
ENSRNA049551074 0.000 NA NA NA NA NA
ENSRNA049551126 0.000 NA NA NA NA NA
ENSRNA049551197 0.000 NA NA NA NA NA
ENSRNA049551244 0.000 NA NA NA NA NA
ENSRNA049551269 0.000 NA NA NA NA NA
ENSRNA049551298 0.000 NA NA NA NA NA
ENSRNA049551337 0.000 NA NA NA NA NA
ENSRNA049551365 0.000 NA NA NA NA NA
ENSRNA049551391 0.000 NA NA NA NA NA
ENSRNA049551484 0.000 NA NA NA NA NA
ENSRNA049551550 0.000 NA NA NA NA NA
ENSRNA049551574 0.000 NA NA NA NA NA
ENSRNA049551598 0.000 NA NA NA NA NA
ENSRNA049551636 0.000 NA NA NA NA NA
ENSRNA049551673 0.000 NA NA NA NA NA
ENSRNA049551717 0.000 NA NA NA NA NA
ENSRNA049551745 0.000 NA NA NA NA NA
ENSRNA049551774 0.000 NA NA NA NA NA
ENSRNA049551798 0.000 NA NA NA NA NA
ENSRNA049551814 0.000 NA NA NA NA NA
ENSRNA049551862 0.000 NA NA NA NA NA
ENSRNA049551899 0.000 NA NA NA NA NA
ENSRNA049551942 0.000 NA NA NA NA NA
ENSRNA049551964 0.000 NA NA NA NA NA
ENSRNA049551993 0.000 NA NA NA NA NA
large_MTrRNA 5333.9740.94856430.13434097.0608761.654565e-123.752466e-12
small_MTrRNA 1705.7550.69147240.29553902.3397001.929926e-022.502916e-02
 
[38]:

### Tidy the results for DE with respect to condition
### Results are sorted in ascending order by adjusted p-value
### Here ph4 is the reference level
### log2FC > 0 suggests that higher pH (pH8) is associated with increased expression
### log2FC < 0 suggests that higher pH (pH8) is associated with lower expression
myres_condition8v4 <- results(ddsDE, contrast = c("condition", "pH8", "pH4"), tidy = TRUE)

myres_condition8v4 %>%
    arrange(desc(-padj)) %>%
        head(10)
A data.frame: 10 × 7
rowbaseMeanlog2FoldChangelfcSEstatpvaluepadj
<chr><dbl><dbl><dbl><dbl><dbl><dbl>
CNAG_00275 1072.9219-4.0481060.10667356-37.9485400
CNAG_00309 539.1312-3.3315240.08444635-39.4513700
CNAG_00409 1093.4775-5.0306550.07452918-67.4991300
CNAG_0053115371.8856 3.5636910.05663035 62.9290000
CNAG_00601 506.4752 4.6615770.09625726 48.4283200
CNAG_0089712255.8281 4.8569240.07970883 60.9333200
CNAG_0127249334.5209 1.9610090.04521777 43.3681200
CNAG_01275 2534.1015 2.2730880.04888985 46.4940600
CNAG_01344 3528.6160-1.5821220.04028591-39.2723400
CNAG_01713 9652.4535 2.0164860.05310869 37.9690400

Visualize DE effect

Looking at the results for these two genes:

  • The estimated log2FC for CNAG_00275 is negative. We will verify visually that ph8, compared to pH4, is associated with lower expression

  • The estimated log2FC for CNAG_00531 is positive. We will verify visually that ph8, compared to pH4, is associated with higher expression

[39]:

results(ddsDE, tidy = TRUE) %>%
    filter(row %in% c("CNAG_00275","CNAG_00531"))
A data.frame: 2 × 7
rowbaseMeanlog2FoldChangelfcSEstatpvaluepadj
<chr><dbl><dbl><dbl><dbl><dbl><dbl>
CNAG_00275 1072.922-0.32684900.10536149-3.1021670.00192109430.005769654
CNAG_0053115371.886 0.19024240.05662427 3.3597330.00078017940.002625616
 
[40]:

### This dot plot verify visually that exposure to ph8, compared to pH4,  is associated with lower expression
plotCounts(dds2019, "CNAG_00275", intgroup = "condition")
../_images/DESeq2_02_deseq_DE_78_0.png 
[41]:

### This dot plot verify visually that exposure to ph8, compared to pH4,  is associated with higher expression
plotCounts(dds2019, "CNAG_00531", intgroup = "condition")
../_images/DESeq2_02_deseq_DE_79_0.png

Volcano plot

[42]:

### Volcano plot for con effect
ggplot(results(ddsDE, contrast = c("condition", "pH4", "pH8"), tidy = TRUE),
       aes(x = log2FoldChange, y = -log10(padj))) + geom_point()

Warning message:
“Removed 348 rows containing missing values (geom_point).”
../_images/DESeq2_02_deseq_DE_81_1.png 
[43]:

### Genotype Effect
ggplot(results(ddsDE, contrast = c("genotype", "sre1d", "WT"), tidy = TRUE),
       aes(x = log2FoldChange, y = -log10(padj))) + geom_point()

Warning message:
“Removed 348 rows containing missing values (geom_point).”
../_images/DESeq2_02_deseq_DE_82_1.png

Clustering

Regularized log transformation

The regularized log transform can be obtained using the rlog() function. Note that an important argument for this function is blind (TRUE by default). The default “blinds” the normalization to the design. This is very important so as to not bias the analyses (e.g. class discovery)

[44]:

rld <- rlog(dds2019, blind = TRUE)

Dendrogram of samples: showing strain & media of each sample

Hierarchical clustering using rlog transformation

[45]:

options(repr.plot.width = 9, repr.plot.height = 5)
dists <- dist(t(assay(rld)))
plot(hclust(dists))
../_images/DESeq2_02_deseq_DE_88_0.png

Store the dendrogram of samples using hierarchical clustering

[46]:

assay(rld) %>%
    t() %>%
    dist %>%
    hclust(method = "complete") %>%
    as.dendrogram ->
    mydend

Dendrogram of samples: showing strain of each sample

[47]:


dendplot <- function(mydend, columndata, labvar, colvar, pchvar) {
    cols <- factor(columndata[[colvar]][order.dendrogram(mydend)])
    collab <- brewer.pal(max(3,nlevels(cols)),"Set1")[cols]
    pchs <- factor(columndata[[pchvar]][order.dendrogram(mydend)])
    pchlab <- seq_len(nlevels(pchs))[pchs]
    lablab <- columndata[[labvar]][order.dendrogram(mydend)]

    mydend %>%
        set("labels_cex",1) %>%
        set("labels_col",collab) %>%
        set("leaves_pch",pchlab) %>%
        set("labels", lablab)
}



[48]:

options(repr.plot.width = 9, repr.plot.height = 5)
dendplot(mydend, dds2019@colData,
         "genotype",    # variable that show in label
         "genotype",    # variable that define color
         "condition") %>% # variable that define shape of points
    plot
../_images/DESeq2_02_deseq_DE_93_0.png

Dendrogram of samples: showing media of each sample

Customize presentation

[49]:

### Merge gene expression with meta data
myDEplotData <- function(mydds, geneid, mergelab) {
    counts(mydds, normalize = TRUE) %>%
        as_tibble(rownames="gene") %>%
        filter(gene == geneid) %>%
        gather(Label, geneexp, -gene) %>%
        select(-gene) -> genedat

    colData(mydds) %>%
        as.data.frame %>%
        as_tibble %>%
        full_join(genedat, by = mergelab) -> genedat

    return(genedat)
}

myDEplotData(dds2019, "CNAG_00003", "Label")[,c("Label", "genotype", "condition" , "geneexp")]
A tibble: 24 × 4
Labelgenotypeconditiongeneexp
<chr><fct><fct><dbl>
1_2019_P_M1 WT pH4126.28489
2_2019_P_M1 WT pH4 94.04893
3_2019_P_M1 WT pH4106.41852
4_2019_P_M1 WT pH4120.45064
5_2019_P_M1 WT pH4112.70469
6_2019_P_M1 WT pH4 83.25923
7_2019_P_M1 sre1dpH4122.72715
8_2019_P_M1 sre1dpH4117.87771
9_2019_P_M1 sre1dpH4149.00707
10_2019_P_M1sre1dpH4121.70596
11_2019_P_M1sre1dpH4116.10471
12_2019_P_M1sre1dpH4155.47163
13_2019_P_M1WT pH8278.86338
14_2019_P_M1WT pH8204.68325
15_2019_P_M1WT pH8274.49575
16_2019_P_M1WT pH8217.18785
17_2019_P_M1WT pH8236.11254
18_2019_P_M1WT pH8215.75364
19_2019_P_M1sre1dpH8157.45180
20_2019_P_M1sre1dpH8171.48644
21_2019_P_M1sre1dpH8186.69737
22_2019_P_M1sre1dpH8178.12113
23_2019_P_M1sre1dpH8203.92325
24_2019_P_M1sre1dpH8208.94213
 
[50]:

### Basic function

myDEplot <- function(mydds, geneid, mergelab) {
    mydat <- myDEplotData(mydds, geneid, mergelab)
    ggplot(mydat, aes(x = condition, y = geneexp))+ geom_point()
}



[51]:

### Allow for grouping by any factor in dataframe

myDEplot <- function(mydds, geneid, grpvar, mergelab) {
    mydat <- myDEplotData(mydds, geneid, mergelab)
    ggplot(mydat, aes_string(x=grpvar, y = "geneexp"))+ geom_point()
}

myDEplot(dds2019, "CNAG_00003", "genotype", "Label")
myDEplot(dds2019, "CNAG_00003", "condition", "Label")
../_images/DESeq2_02_deseq_DE_98_0.png
../_images/DESeq2_02_deseq_DE_98_1.png 
[52]:

### Add color

myDEplot <- function(mydds, geneid, grpvar, mergelab) {
    mydat <- myDEplotData(mydds, geneid, mergelab)
    ggplot(mydat, aes_string(x=grpvar, y = "geneexp", col = grpvar))+ geom_point()
}
myDEplot(dds2019, "CNAG_00003", "condition", "Label")
../_images/DESeq2_02_deseq_DE_99_0.png 
[53]:

### Alow for coloring with respect to another factor
myDEplot <- function(mydds, geneid, grpvar, colvar, mergelab) {
    mydat <- myDEplotData(mydds, geneid, mergelab)
    ggplot(mydat, aes_string(x=grpvar, y = "geneexp", col = colvar))+ geom_point()
}
myDEplot(dds2019, "CNAG_00003", "condition", "genotype", "Label")
../_images/DESeq2_02_deseq_DE_100_0.png 
[54]:

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] RColorBrewer_1.1-2          dendextend_1.12.0
 [3] DESeq2_1.24.0               SummarizedExperiment_1.14.0
 [5] DelayedArray_0.10.0         BiocParallel_1.18.0
 [7] matrixStats_0.54.0          Biobase_2.44.0
 [9] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0
[11] IRanges_2.18.1              S4Vectors_0.22.0
[13] BiocGenerics_0.30.0         forcats_0.4.0
[15] stringr_1.4.0               dplyr_0.8.1
[17] purrr_0.3.2                 readr_1.3.1
[19] tidyr_0.8.3                 tibble_2.1.2
[21] ggplot2_3.1.1               tidyverse_1.2.1

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