Lorena Pantano Harvard TH Chan School of Public Health, Boston, US

``````library(DEGreport)
data(humanGender)``````

## General QC figures from DE analysis

We are going to do a differential expression analysis with edgeR/DESeq2. We have an object that is coming from the edgeR package. It contains a gene count matrix for 85 TSI HapMap individuals, and the gender information. With that, we are going to apply the `glmFit` function or DESeq2 to get genes differentially expressed between males and females.

``````library(DESeq2)
idx <- c(1:10, 75:85)
dds <- DESeqDataSetFromMatrix(assays(humanGender)[[1]][1:1000, idx],
colData(humanGender)[idx,], design=~group)
dds <- DESeq(dds)
res <- results(dds)``````

We need to extract the experiment design data.frame where the condition is Male or Female.

``````counts <- counts(dds, normalized = TRUE)
design <- as.data.frame(colData(dds))``````

### Size factor QC

A main assumption in library size factor calculation of edgeR and DESeq2 (and others) is that the majority of genes remain unchanged. Plotting the distribution of gene ratios between each gene and the average gene can show how true this is. Not super useful for many samples because the plot becomes crowed.

``degCheckFactors(counts[, 1:6])``

### Mean-Variance QC plots

p-value distribution gives an idea on how well you model is capturing the input data and as well whether it could be some problem for some set of genes. In general, you expect to have a flat distribution with peaks at 0 and 1. In this case, we add the mean count information to check if any set of genes are enriched in any specific p-value range.

Variation (dispersion) and average expression relationship shouldnâ€™t be a factor among the differentially expressed genes. When plotting average mean and standard deviation, significant genes should be randomly distributed.

In this case, it would be good to look at the ones that are totally outside the expected correlation.

You can put this tree plots together using `degQC`.

``degQC(counts, design[["group"]], pvalue = res[["pvalue"]])``

### Covariates effect on count data

Another important analysis to do if you have covariates is to calculate the correlation between PCs from PCA analysis to different variables you may think are affecting the gene expression. This is a toy example of how the function works with raw data, where clearly library size correlates with some of the PCs.

``````resCov <- degCovariates(log2(counts(dds)+0.5),
colData(dds))``````

### Covariates correlation with metrics

Also, the correlation among covariates and metrics from the analysis can be tested. This is useful when the study has multiple variables, like in clinical trials. The following code will return a correlation table, and plot the correlation heatmap for all the covariates and metrics in a table.

``cor <- degCorCov(colData(dds))``

``names(cor)``
``## [1] "cor"    "corMat" "fdrMat" "plot"``

### QC report

A quick HTML report can be created with `createReport` to show whether a DE analysis is biased to a particular set of genes. It contains the output of `degQC,`degVB and `degMB`.

``````createReport(colData(dds)[["group"]], counts(dds, normalized = TRUE),
row.names(res)[1:20], res[["pvalue"]], path = "~/Downloads")``````

## Report from DESeq2 analysis

Here, we show some useful plots for differentially expressed genes.

### Contrasts

`DEGSet` is a class to store the DE results like the one from `results` function. DESeq2 offers multiple way to ask for contrasts/coefficients. With `degComps` is easy to get multiple results in a single object:

``````degs <- degComps(dds, combs = "group",
contrast = list("group_Male_vs_Female",
c("group", "Female", "Male")))
names(degs)``````
``## [1] "group_Male_vs_Female" "group_Female_vs_Male"``

`degs` contains 3 elements, one for each contrast/coefficient asked for. It contains the results output in the element `raw` and the output of `lfcShrink` in the element `shrunken`. To obtain the results from one of them, use the method `dge`:

``deg(degs[[1]])``
``````## log2 fold change (MAP): group Male vs Female
## Wald test p-value: group Male vs Female
## DataFrame with 1000 rows and 6 columns
##                         baseMean        log2FoldChange              lfcSE
##                        <numeric>             <numeric>          <numeric>
## ENSG00000067048 1025.03783081851      1.93948751981493  0.100694015277881
## ENSG00000012817 411.543865988588      3.70056395376916 0.0989888483685393
## ENSG00000067646 169.814765626402      3.32284568660188 0.0995193327364417
## ENSG00000005889 670.861906530697     -0.48943474075899 0.0929363099393874
## ENSG00000006757 92.6611112908156    -0.472926180708217 0.0992749166821467
## ...                          ...                   ...                ...
## ENSG00000068120 1214.09670077803 -0.000114321924523065 0.0797226606999408
## ENSG00000072062 935.317212058596   0.00055921914204206 0.0915352171745638
## ENSG00000076770 1019.79638067687   0.00085008184559821  0.101666325101382
## ENSG00000078967 166.422097051996  0.000415745391576736 0.0959757955312265
## ENSG00000079246 5226.33895971623  0.000149584255114179 0.0926205100248815
##                                 stat                pvalue
##                            <numeric>             <numeric>
## ENSG00000067048     23.9943800011737 3.18302039221609e-127
## ENSG00000012817     21.8045168852375 2.10213939598916e-105
## ENSG00000067646      15.483847327858  4.45984434297357e-54
## ENSG00000005889    -5.26370810368138  1.41178536870724e-07
## ENSG00000006757    -4.75700569609357  1.96485627341616e-06
## ...                              ...                   ...
## ENSG00000068120 -0.00143399782359701     0.998855835668461
## ENSG00000072062  0.00610871981310329     0.995125977088418
## ENSG00000076770  0.00836192966137436     0.993328223175276
## ENSG00000078967  0.00432971042064206     0.996545401696168
## ENSG00000079246  0.00161499511899518     0.998711420888929
##                                  padj
##                             <numeric>
## ENSG00000067048 3.18302039221609e-124
## ENSG00000012817 1.05106969799458e-102
## ENSG00000067646  1.48661478099119e-51
## ENSG00000005889  3.52946342176811e-05
## ENSG00000006757  0.000392971254683231
## ...                               ...
## ENSG00000068120     0.998855835668461
## ENSG00000072062     0.998855835668461
## ENSG00000076770     0.998855835668461
## ENSG00000078967     0.998855835668461
## ENSG00000079246     0.998855835668461``````

By default it would output the `shrunken` table always, as defined by `degDefault`, that contains the default table to get.

To get the original results table, use the parameter as this:

``deg(degs[[1]], "raw", "tibble")``
``````## # A tibble: 1,000 x 7
##    gene          baseMean log2FoldChange  lfcSE   stat    pvalue      padj
##  * <chr>            <dbl>          <dbl>  <dbl>  <dbl>     <dbl>     <dbl>
##  1 ENSG00000067â€¦   1025.          10.2   0.423   24.0  3.18e-127 3.18e-124
##  2 ENSG00000012â€¦    412.           9.24  0.424   21.8  2.10e-105 1.05e-102
##  3 ENSG00000067â€¦    170.          10.2   0.658   15.5  4.46e- 54 1.49e- 51
##  4 ENSG00000005â€¦    671.          -0.692 0.131   -5.26 1.41e-  7 3.53e-  5
##  5 ENSG00000006â€¦     92.7         -0.767 0.161   -4.76 1.96e-  6 3.93e-  4
##  6 ENSG00000073â€¦    220.          -1.87  0.421   -4.44 8.89e-  6 1.48e-  3
##  7 ENSG00000005â€¦   2027.          -0.742 0.176   -4.21 2.59e-  5 3.69e-  3
##  8 ENSG00000005â€¦   1234.           0.389 0.0952   4.08 4.44e-  5 5.56e-  3
##  9 ENSG00000003â€¦    394.           0.680 0.177    3.85 1.17e-  4 1.31e-  2
## 10 ENSG00000069â€¦    107.          -1.63  0.459   -3.56 3.71e-  4 3.71e-  2
## # ... with 990 more rows``````

Note that the format of the output can be changed to tibble, or data.frame with a third parameter `tidy`.

The table will be always sorted by padj.

And easy way to get significant genes is:

``significants(degs[[1]], fc = 0, fdr = 0.05)``
``````##  [1] "ENSG00000012817" "ENSG00000067646" "ENSG00000067048"
##  [4] "ENSG00000005889" "ENSG00000006757" "ENSG00000005302"
##  [7] "ENSG00000003400" "ENSG00000073282" "ENSG00000005020"
## [10] "ENSG00000069702"``````

This function can be used as well for a list of comparisons:

``significants(degs, fc = 0, fdr = 0.05)``
``````##  [1] "ENSG00000012817" "ENSG00000067646" "ENSG00000067048"
##  [4] "ENSG00000005889" "ENSG00000006757" "ENSG00000005302"
##  [7] "ENSG00000003400" "ENSG00000073282" "ENSG00000005020"
## [10] "ENSG00000069702"``````

And it can returns the full table for a list:

``significants(degs, fc = 0, fdr = 0.05, full = TRUE)``
``````## # A tibble: 10 x 5
##    gene     log2FoldChange_groupâ€¦ log2FoldChange_groupâ€¦ padj_group_Femaleâ€¦
##    <chr>                    <dbl>                 <dbl>              <dbl>
##  1 ENSG000â€¦                -0.494                 0.388          1.31e-  2
##  2 ENSG000â€¦                -0.351                 0.319          5.56e-  3
##  3 ENSG000â€¦                 0.541                -0.425          3.69e-  3
##  4 ENSG000â€¦                 0.573                -0.489          3.53e-  5
##  5 ENSG000â€¦                 0.585                -0.473          3.93e-  4
##  6 ENSG000â€¦                -5.01                  3.70           1.05e-102
##  7 ENSG000â€¦                -3.63                  1.94           3.18e-124
##  8 ENSG000â€¦                -4.43                  3.32           1.49e- 51
##  9 ENSG000â€¦                 0.453                -0.260          3.71e-  2
## 10 ENSG000â€¦                 0.582                -0.339          1.48e-  3
## # ... with 1 more variable: padj_group_Male_vs_Female <dbl>``````

Since log2FoldChange are shrunken, the method for DEGSet class now can plot these changes as follow:

``plotMA(degs[[1]], diff = 2, limit = 3)``

The blue arrows indicate how foldchange is affected by this new feature.

As well, it can plot the original MA plot:

``plotMA(degs[[1]], diff = 2, limit = 3, raw = TRUE)``