Keep it simple

- Classical experimental designs
- Time series
- Without missing values, where possible
- Intended analysis must be feasbile – can the available samples and hypothesis of interest be combined to formulate a testable statistical hypothesis?

Replicate

- Extent of replication determines nuance of biological question.
- No replication (1 sample per treatment): qualitative description with limited statistical options.
- 3-5 replicates per treatment: designed experimental manipulation with cell lines or other well-defined entities; 2-fold (?) change in average expression between groups.
- 10-50 replicates per treatment: population studies, e.g., cancer cell lines.
- 1000’s of replicates: prospective studies, e.g., SNP discovery
- One resource:
*RNASeqPower*

Avoid confounding experimental factors with other factors

- Common problems: samples from one treatment all on the same flow cell; samples from treatment 1 processed first, treatment 2 processed second, etc.

Record co-variates

Be aware of *batch effects*

Known

- Phenotypic covariates, e.g., age, gender
- Experimental covariates, e.g., lab or date of processing
- Incorporate into linear model, at least approximately

Unknown

- Or just unexpected / undetected
- Characterize using, e.g.,
*sva*.

Surrogate variable analysis

- Leek et al., 2010, Nature Reviews Genetics 11 733-739, Leek & Story PLoS Genet 3(9): e161.
- Scientific finding: pervasive batch effects
- Statistical insights: surrogate variable analysis: identify and build surrogate variables; remove known batch effects
- Benefits: reduce dependence, stabilize error rate estimates, and improve reproducibility
*combat*software /*sva**Bioconductor*package

HapMap samples from one facility, ordered by date of processing.

Confounding factors

- Record or avoid

Artifacts of your *particular* protocols

- Sequence contaminants
- Enrichment bias, e.g., non-uniform transcript representation.
- PCR artifacts – adapter contaminants, sequence-specific amplification bias, …

Axes of variation

- Single- versus paired-end
- Length: 50-200nt
- Number of reads per sample

Application-specific, e.g.,

- ChIP-seq: short, single-end reads are usually sufficient
- RNA-seq, known genes: single- or paired-end reads
- RNA-seq, transcripts or novel variants: paired-end reads
- Copy number: single- or paired-end reads
- Structural variants: paired-end reads
- Variants: depth via longer, paired-end reads
- Microbiome: long paired-end reads (overlapping ends)

Alignment strategies

*de novo*- No reference genome; considerable sequencing and computational resources
- Genome
- Established reference genome
- Splice-aware aligners
- Novel transcript discovery
- Transcriptome
- Established reference genome; reliable gene model
- Simple aligners
- Known gene / transcript expression

Splice-aware aligners (and *Bioconductor* wrappers)

- tophat uses Bowtie2 to perform basic single- and paired-end alignments, then uses algorithms to place difficult-to-align reads near to their well-aligned mates.
- Cufflinks (doi) takes
*tophat*output and estimate existing and novel transcript abundance. How Cufflinks Works - Cuffdiff assesses statistical significance of estimated abundances between experimental groups
- RSEM includes de novo assembly and quantification

- Use known gene model to count aligned reads overlapping regions of interest / gene models
- Gene model can be public (e.g., UCSC, NCBI, ENSEMBL) or
*ad hoc*(gff file) `GenomicAlignments::summarizeOverlaps()`

`Rsubread::featureCount()`

- HTSeq, htseq-count

Unique statistical aspects

- Large data, few samples
- Comparison of each gene, across samples;
*univariate*measures - Each gene is analyzed by the
*same*experimental design, under the*same*null hypothesis

Summarization

- Counts
*per se*, rather than a summary (RPKM, FRPKM, …), are relevant for analysis - For a given gene, larger counts imply more information; RPKM etc., treat all estimates as equally informative.
- Comparison is across samples at
*each*region of interest; all samples have the same region of interest, so modulo library size there is no need to correct for, e.g., gene length or mapability.

Normalization

- Libraries differ in size (total counted reads per sample) for un-interesting reasons; we need to account for differences in library size in statistical analysis.
- Total number of counted reads per sample is
*not*a good estimate of library size. It is un-necessarily influenced by regions with large counts, and can introduce bias and correlation across genes. Instead, use a robust measure of library size that takes account of skew in the distribution of counts (simplest: trimmed geometric mean; more advanced / appropriate encountered in the lab). - Library size (total number of counted reads) differs between samples, and should be included
*as a statistical offset*in analysis of differential expression, rather than ‘dividing by’ the library size early in an analysis.

Appropriate error model

- Count data is
*not*distributed normally or as a Poisson process, but rather as negative binomial. - Result of a combination Poisson (`shot’ noise, i.e., within-sample technical and sampling variation in read counts) with variation between biological samples.
- A negative binomial model requires estimation of an additional parameter (‘dispersion’), which is estimated poorly in small samples.
- Basic strategy is to moderate per-gene estimates with more robust local estimates derived from genes with similar expression values (a little more on borrowing information is provided below).

Pre-filtering

- Naively, a statistical test (e.g., t-test) could be applied to each row of a counts table. However, we have relatively few samples (10’s) and very many comparisons (10,000’s) so a naive approach is likely to be very underpowered, resulting in a very high
*false discovery rate* - A simple approach is perform fewer tests by removing regions that could not possibly result in statistical significance, regardless of hypothesis under consideration.
- Example: a region with 0 counts in all samples could not possibly be significant regradless of hypothesis, so exclude from further analysis.
- Basic approaches: ‘K over A’-style filter – require a minimum of A (normalized) read counts in at least K samples. Variance filter, e.g., IQR (inter-quartile range) provides a robust estimate of variability; can be used to rank and discard least-varying regions.
- More nuanced approaches:
*edgeR*vignette; work flow today.

Borrowing information

- Why does low statistical power elevate false discovery rate?
- One way of developing intuition is to recognize a t-test (for example) as a ratio of variances. The numerator is treatment-specific, but the denominator is a measure of overall variability.
- Variances are measured with uncertainty; over- or under-estimating the denominator variance has an asymmetric effect on a t-statistic or similar ratio, with an underestimate
*inflating*the statistic more dramatically than an overestimate deflates the statistic. Hence elevated false discovery rate. - Under the typical null hypothesis used in microarray or RNA-seq experiments, each gene may respond differently to the treatment (numerator variance) but the overall variability of a gene is the same, at least for genes with similar average expression
- The strategy is to estimate the denominator variance as the between-group variance for the gene,
*moderated*by the average between-group variance across all genes. - This strategy exploits the fact that the same experimental design has been applied to all genes assayed, and is effective at moderating false discovery rate.

Placing differentially expressed regions in context

- Gene names associated with genomic ranges
- Gene set enrichment and similar analysis
- Proximity to regulatory marks
- Integrate with other analyses, e.g., methylation, copy number, variants, …

Correlation between genomic copy number and mRNA expression identified 38 mis-labeled samples in the TCGA ovarian cancer Affymetrix microarray dataset.

*DESeq2* `estimateSizeFactors()`

, Anders and Huber, 2010

- For each gene: geometric mean of all samples.
- For each sample: median ratio of the sample gene over the geometric mean of all samples
- Functions other than the median can be used; control genes can be used instead

*edgeR* `calcNormFactors()`

TMM method of Robinson and Oshlack, 2010

- Identify reference sample: library with upper quartile closest to the mean upper quartile of all libraries
- Calculate M-value of each gene (log-fold change relative to reference)
- Summarize library size as weighted trimmed mean of M-values.

*DESeq2* `estimateDispersions()`

- Estimate per-gene dispersion
- Fit a smoothed relationship between dispersion and abundance

*edgeR* `estimateDisp()`

- Common: single dispersion for all genes; appropriate for small experiments (<10? samples)
- Tagwise: different dispersion for all genes; appropriate for larger / well-behaved experiments
- Trended: bin based on abundance, estimate common dispersion within bin, fit a loess-smoothed relationship between binned dispersion and abundance

`t.test()`

`x`

: vector of univariate measurements`y`

:`factor`

describing experimental design`var.equal=TRUE`

: appropriate for relatively small experiments where no additional information available?`formula`

: alternative representation,`y ~ x`

.

`head(sleep)`

```
## extra group ID
## 1 0.7 1 1
## 2 -1.6 1 2
## 3 -0.2 1 3
## 4 -1.2 1 4
## 5 -0.1 1 5
## 6 3.4 1 6
```

`plot(extra ~ group, data = sleep)`

```
## Traditional interface
with(sleep, t.test(extra[group == 1], extra[group == 2]))
```

```
##
## Welch Two Sample t-test
##
## data: extra[group == 1] and extra[group == 2]
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.3654832 0.2054832
## sample estimates:
## mean of x mean of y
## 0.75 2.33
```

```
## Formula interface
t.test(extra ~ group, sleep)
```

```
##
## Welch Two Sample t-test
##
## data: extra by group
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.3654832 0.2054832
## sample estimates:
## mean in group 1 mean in group 2
## 0.75 2.33
```

```
## equal variance between groups
t.test(extra ~ group, sleep, var.equal=TRUE)
```

```
##
## Two Sample t-test
##
## data: extra by group
## t = -1.8608, df = 18, p-value = 0.07919
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.363874 0.203874
## sample estimates:
## mean in group 1 mean in group 2
## 0.75 2.33
```

`lm()`

and `anova()`

`lm()`

: fit*linear model*.`anova()`

: statisitcal evaluation.

```
## linear model; compare to t.test(var.equal=TRUE)
fit <- lm(extra ~ group, sleep)
anova(fit)
```

```
## Analysis of Variance Table
##
## Response: extra
## Df Sum Sq Mean Sq F value Pr(>F)
## group 1 12.482 12.4820 3.4626 0.07919 .
## Residuals 18 64.886 3.6048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

- Under the hood:
`formula`

: translated into*model matrix*, used in`lm.fit()`

. - With (implicit) intercept 1, last coefficient of model matrix reflects group effect
- With intercept 0,
*contrast*between effects of coefficient 1 and coefficient 2 reflect group effect

```
## underlying model, used in `lm.fit()`
model.matrix(extra ~ group, sleep) # last column indicates group effect
```

```
## (Intercept) group2
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 0
## 8 1 0
## 9 1 0
## 10 1 0
## 11 1 1
## 12 1 1
## 13 1 1
## 14 1 1
## 15 1 1
## 16 1 1
## 17 1 1
## 18 1 1
## 19 1 1
## 20 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
```

`model.matrix(extra ~ 0 + group, sleep) # contrast between columns`

```
## group1 group2
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 0
## 8 1 0
## 9 1 0
## 10 1 0
## 11 0 1
## 12 0 1
## 13 0 1
## 14 0 1
## 15 0 1
## 16 0 1
## 17 0 1
## 18 0 1
## 19 0 1
## 20 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
```

- Covariate – fit base model containing only covariate, test improvement in fit when model includes factor of interest

```
fit0 <- lm(extra ~ ID, sleep)
fit1 <- lm(extra ~ ID + group, sleep)
anova(fit0, fit1)
```

```
## Analysis of Variance Table
##
## Model 1: extra ~ ID
## Model 2: extra ~ ID + group
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 19.290
## 2 9 6.808 1 12.482 16.501 0.002833 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

`t.test(extra ~ group, sleep, var.equal=TRUE, paired=TRUE)`

```
##
## Paired t-test
##
## data: extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.4598858 -0.7001142
## sample estimates:
## mean of the differences
## -1.58
```

`genefilter::rowttests()`

- t-tests for gene expression data
- useful for exploratory analysis, but statistically sub-optimal
`x`

: matrix of expression valuesfeatures x samples (reverse of how a ‘statistician’ would represent the data – samples x features)

`fac`

: factor of one or two levels describing experimental design

Limitations

- Assumes features are
*independent* - Ignores common experimental design
- Ignores multiple testing

Consequences

- Poor estimate of between-group variance for each feature
- Elevated false discovery rate

- t-test:
`count ~ factor`

. Alternative:`count ~ 0 + factor`

and contrasts - covariates:
`count ~ covariate + factor`

- Single factor, multiple levels (one-way ANOVA) – statistical contrasts: specify model as
`count ~ factor`

or`count ~ 0 + factor`

- Factorial designs – main effects,
`count ~ factor1 + factor2`

; main effects and interactions,`count ~ factor1 * factor2`

. Contrasts to ask specific questions - Paired designs: include ID as covariate (approximate, since ID is a random effect);
*limma*approach:`duplicateCorrelation()`

Adapted from Love, Anders, and Huber’s Bioconductor work flow

Michael Love [1], Simon Anders [2], Wolfgang Huber [2]

[1] Department of Biostatistics, Dana-Farber Cancer Institute and Harvard School of Public Health, Boston, US;

[2] European Molecular Biology Laboratory (EMBL), Heidelberg, Germany.

The data used in this workflow is an RNA-Seq experiment of airway smooth muscle cells treated with dexamethasone, a synthetic glucocorticoid steroid with anti-inflammatory effects. Glucocorticoids are used, for example, in asthma patients to prevent or reduce inflammation of the airways. In the experiment, four primary human airway smooth muscle cell lines were treated with 1 micromolar dexamethasone for 18 hours. For each of the four cell lines, we have a treated and an untreated sample. The reference for the experiment is:

Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q. “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells.” PLoS One. 2014 Jun 13;9(6):e99625. PMID: 24926665. GEO: GSE52778.

We use the *airway* package to illustrate reduction. The package provides sample information, a subset of eight BAM files, and the known gene models required to count the reads.

```
library(airway)
path <- system.file(package="airway", "extdata")
dir(path)
```

```
## [1] "GSE52778_series_matrix.txt"
## [2] "Homo_sapiens.GRCh37.75_subset.gtf"
## [3] "SRR1039508_subset.bam"
## [4] "SRR1039509_subset.bam"
## [5] "SRR1039512_subset.bam"
## [6] "SRR1039513_subset.bam"
## [7] "SRR1039516_subset.bam"
## [8] "SRR1039517_subset.bam"
## [9] "SRR1039520_subset.bam"
## [10] "SRR1039521_subset.bam"
## [11] "SraRunInfo_SRP033351.csv"
## [12] "sample_table.csv"
```

The ingredients for counting include are:

- Metadata describing samples. Read this using
`read.csv()`

.

```
csvfile <- dir(path, "sample_table.csv", full=TRUE)
sampleTable <- read.csv(csvfile, row.names=1)
head(sampleTable)
```

```
## SampleName cell dex albut Run avgLength Experiment
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
## Sample BioSample
## SRR1039508 SRS508568 SAMN02422669
## SRR1039509 SRS508567 SAMN02422675
## SRR1039512 SRS508571 SAMN02422678
## SRR1039513 SRS508572 SAMN02422670
## SRR1039516 SRS508575 SAMN02422682
## SRR1039517 SRS508576 SAMN02422673
```

- BAM files containing aligned reads. Create an object that references these files. What does the
`yieldSize`

argument mean?

```
library(Rsamtools)
filenames <- dir(path, ".bam$", full=TRUE)
bamfiles <- BamFileList(filenames, yieldSize=1000000)
names(bamfiles) <- sub("_subset.bam", "", basename(filenames))
```

- Known gene models. These might come from an existing
`TxDb`

package, or created from biomart or UCSC, or from a GTF file. We’ll take the hard road, making a TxDb object from the GTF file used to align reads and using the TxDb to get all exons, grouped by gene.

```
library(GenomicFeatures)
gtffile <- file.path(path, "Homo_sapiens.GRCh37.75_subset.gtf")
txdb <- makeTxDbFromGFF(gtffile, format="gtf", circ_seqs=character())
```

`## Prepare the 'metadata' data frame ... metadata: OK`

`genes <- exonsBy(txdb, by="gene")`

After these preparations, the actual counting is easy. The function `summarizeOverlaps()`

from the *GenomicAlignments* package will do this. This produces a `SummarizedExperiment`

object, which contains a variety of information about an experiment

```
library(GenomicAlignments)
se <- summarizeOverlaps(features=genes, reads=bamfiles,
mode="Union",
singleEnd=FALSE,
ignore.strand=TRUE,
fragments=TRUE)
colData(se) <- as(sampleTable, "DataFrame")
se
```

```
## class: SummarizedExperiment
## dim: 20 8
## exptData(0):
## assays(1): counts
## rownames(20): ENSG00000009724 ENSG00000116649 ... ENSG00000271794
## ENSG00000271895
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample
```

`colData(se)`

```
## DataFrame with 8 rows and 9 columns
## SampleName cell dex albut Run avgLength
## <factor> <factor> <factor> <factor> <factor> <integer>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98
## Experiment Sample BioSample
## <factor> <factor> <factor>
## SRR1039508 SRX384345 SRS508568 SAMN02422669
## SRR1039509 SRX384346 SRS508567 SAMN02422675
## SRR1039512 SRX384349 SRS508571 SAMN02422678
## SRR1039513 SRX384350 SRS508572 SAMN02422670
## SRR1039516 SRX384353 SRS508575 SAMN02422682
## SRR1039517 SRX384354 SRS508576 SAMN02422673
## SRR1039520 SRX384357 SRS508579 SAMN02422683
## SRR1039521 SRX384358 SRS508580 SAMN02422677
```

`rowData(se)`

```
## Warning: 'rowData' is deprecated.
## Use 'rowRanges' instead.
## See help("Deprecated")
```

```
## GRangesList object of length 20:
## $ENSG00000009724
## GRanges object with 18 ranges and 2 metadata columns:
## seqnames ranges strand | exon_id exon_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] 1 [11086580, 11087705] - | 98 ENSE00000818830
## [2] 1 [11090233, 11090307] - | 99 ENSE00000472123
## [3] 1 [11090805, 11090939] - | 100 ENSE00000743084
## [4] 1 [11094885, 11094963] - | 101 ENSE00000743085
## [5] 1 [11097750, 11097868] - | 103 ENSE00003520086
## ... ... ... ... ... ... ...
## [14] 1 [11106948, 11107176] - | 111 ENSE00003467404
## [15] 1 [11106948, 11107176] - | 112 ENSE00003489217
## [16] 1 [11107260, 11107280] - | 113 ENSE00001833377
## [17] 1 [11107260, 11107284] - | 114 ENSE00001472289
## [18] 1 [11107260, 11107290] - | 115 ENSE00001881401
##
## ...
## <19 more elements>
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
```

`head(assay(se))`

```
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000009724 38 28 66 24 42
## ENSG00000116649 1004 1255 1122 1313 1100
## ENSG00000120942 218 256 233 252 269
## ENSG00000120948 2751 2080 3353 1614 3519
## ENSG00000171819 4 50 19 543 1
## ENSG00000171824 869 1075 1115 1051 944
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000009724 41 47 36
## ENSG00000116649 1879 745 1536
## ENSG00000120942 465 207 400
## ENSG00000120948 3716 2220 1990
## ENSG00000171819 10 14 1067
## ENSG00000171824 1405 748 1590
```

The previous section illustrates the reduction step on a subset of the data; here’s the full data set

```
data(airway)
se <- airway
```

This object contains an informative `colData`

slot – prepared as described in the *airway* vignette. In particular, the `colData()`

include columns describing the cell line `cell`

and treatment `dex`

for each sample

`colData(se)`

```
## DataFrame with 8 rows and 9 columns
## SampleName cell dex albut Run avgLength
## <factor> <factor> <factor> <factor> <factor> <integer>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98
## Experiment Sample BioSample
## <factor> <factor> <factor>
## SRR1039508 SRX384345 SRS508568 SAMN02422669
## SRR1039509 SRX384346 SRS508567 SAMN02422675
## SRR1039512 SRX384349 SRS508571 SAMN02422678
## SRR1039513 SRX384350 SRS508572 SAMN02422670
## SRR1039516 SRX384353 SRS508575 SAMN02422682
## SRR1039517 SRX384354 SRS508576 SAMN02422673
## SRR1039520 SRX384357 SRS508579 SAMN02422683
## SRR1039521 SRX384358 SRS508580 SAMN02422677
```

*DESeq2* makes the analysis particularly easy, simply add the experimental design, run the pipeline, and extract the results

```
library(DESeq2)
dds <- DESeqDataSet(se, design = ~ cell + dex)
dds <- DESeq(dds)
```

```
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
```

`res <- results(dds)`

Simple visualizations / sanity checks include

Look at counts of strongly differentiated genes, to get a sense of how counts translate to the summary statistics reported in the result table

`topGene <- rownames(res)[which.min(res$padj)] res[topGene,]`

`## log2 fold change (MAP): dex untrt vs trt ## Wald test p-value: dex untrt vs trt ## DataFrame with 1 row and 6 columns ## baseMean log2FoldChange lfcSE stat pvalue ## <numeric> <numeric> <numeric> <numeric> <numeric> ## ENSG00000152583 997.4398 -4.3161 0.1724125 -25.03357 2.636198e-138 ## padj ## <numeric> ## ENSG00000152583 4.624155e-134`

`plotCounts(dds, gene=topGene, intgroup=c("dex"))`

An ‘MA’ plot shows for each gene the between-group log-fold-change versus average log count; it should be funnel-shaped and approximately symmetric around

`y=0`

, with lots of between-treatment variation for genes with low counts.`plotMA(res, ylim=c(-5,5))`

Plot the distribution of (unadjusted) P values, which should be uniform (under the null) but with a peak at small P value (true positives, hopefully!)

`hist(res$pvalue, breaks=50)`

Look at a ‘volcano plot’ of adjusted P-value versus log fold change, to get a sense of the fraction of up- versus down-regulated genes

`plot(-log10(padj) ~ log2FoldChange, as.data.frame(res), pch=20)`

Many additional diagnostic approaches are described in the DESeq2 (and edgeR) vignettes, and in the RNA-seq gene differential expression work flow.

see Part E, Gene Set Enrichment