Contents

1 Introduction

1.1 R

Functions

rnorm(10)     # sample 10 random normal deviates
##  [1] -0.8064858  1.0638622  0.3445391 -0.6743962  0.5320705  0.8847153  0.2709375  0.8436201
##  [9]  0.1594222 -0.4931706

Vectors

x <- rnorm(1000)
y <- x + rnorm(1000)
plot(y ~ x)

Objects: classes and methods

df <- data.frame(Y = y, X = x)
fit <- lm(Y ~ X, df)
plot(Y ~ X, df)
abline(fit, col="red", lwd=2)

anova(fit)
## Analysis of Variance Table
## 
## Response: Y
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## X           1  908.41  908.41  895.63 < 2.2e-16 ***
## Residuals 998 1012.24    1.01                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Packages

library(ggplot2)
ggplot(df, aes(x=X, y=Y)) + geom_point() + geom_smooth(color="blue") +
    geom_smooth(method="lm", color="red")

Help!

1.2 Bioconductor

Biological domains

Principles

Installation and use

1.3 Sequence Analysis

Six stages

Key packages: data access

Coordinated management: SummarizedExperiment

2 Case Study: RNA-Seq Differential Expression

This is derived from: RNA-Seq workflow: gene-level exploratory analysis and differential expression, by Michael Love, Simon Anders, Wolfgang Huber; modified by Martin Morgan, October 2015.

We walk through an end-to-end RNA-Seq differential expression workflow, using DESeq2 along with other Bioconductor packages.

The complete work flow starts from the FASTQ files, but we will start after reads have been aligned to a reference genome and reads overlapping known genes have been counted.

Our main focus is on differential gene expression analysis with DESeq2. Other Bioconductor packages are important in statistical inference of differential expression at the gene level, including Rsubread, edgeR, limma, BaySeq, and others.

The data are from 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.

2.1 Preparation

Counts

SummarizedExperiment

library(airway)
data("airway")
se <- airway

Assay data

head(assay(se))
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## ENSG00000000003        679        448        873        408       1138       1047        770
## ENSG00000000005          0          0          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587        799        417
## ENSG00000000457        260        211        263        164        245        331        233
## ENSG00000000460         60         55         40         35         78         63         76
## ENSG00000000938          0          0          2          0          1          0          0
##                 SRR1039521
## ENSG00000000003        572
## ENSG00000000005          0
## ENSG00000000419        508
## ENSG00000000457        229
## ENSG00000000460         60
## ENSG00000000938          0

Experimental design

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

Features (genes)

rowRanges(se)
## GRangesList object of length 64102:
## $ENSG00000000003 
## GRanges object with 17 ranges and 2 metadata columns:
##        seqnames               ranges strand   |   exon_id       exon_name
##           <Rle>            <IRanges>  <Rle>   | <integer>     <character>
##    [1]        X [99883667, 99884983]      -   |    667145 ENSE00001459322
##    [2]        X [99885756, 99885863]      -   |    667146 ENSE00000868868
##    [3]        X [99887482, 99887565]      -   |    667147 ENSE00000401072
##    [4]        X [99887538, 99887565]      -   |    667148 ENSE00001849132
##    [5]        X [99888402, 99888536]      -   |    667149 ENSE00003554016
##    ...      ...                  ...    ... ...       ...             ...
##   [13]        X [99890555, 99890743]      -   |    667156 ENSE00003512331
##   [14]        X [99891188, 99891686]      -   |    667158 ENSE00001886883
##   [15]        X [99891605, 99891803]      -   |    667159 ENSE00001855382
##   [16]        X [99891790, 99892101]      -   |    667160 ENSE00001863395
##   [17]        X [99894942, 99894988]      -   |    667161 ENSE00001828996
## 
## ...
## <64101 more elements>
## -------
## seqinfo: 722 sequences (1 circular) from an unspecified genome

2.2 Analysis Pipeline

Create DESeqDataSet object

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

Run the pipeline

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

Summarize results

res <- results(dds)
head(res)
## log2 fold change (MAP): dex trt vs untrt 
## Wald test p-value: dex trt vs untrt 
## DataFrame with 6 rows and 6 columns
##                    baseMean log2FoldChange      lfcSE       stat       pvalue        padj
##                   <numeric>      <numeric>  <numeric>  <numeric>    <numeric>   <numeric>
## ENSG00000000003 708.6021697    -0.37424998 0.09873107 -3.7906000 0.0001502838 0.001251416
## ENSG00000000005   0.0000000             NA         NA         NA           NA          NA
## ENSG00000000419 520.2979006     0.20215551 0.10929899  1.8495642 0.0643763851 0.192284345
## ENSG00000000457 237.1630368     0.03624826 0.13684258  0.2648902 0.7910940556 0.910776144
## ENSG00000000460  57.9326331    -0.08523371 0.24654402 -0.3457140 0.7295576905 0.878646940
## ENSG00000000938   0.3180984    -0.11555969 0.14630532 -0.7898530 0.4296136373          NA
mcols(res)                     # what does each column mean?
## DataFrame with 6 rows and 2 columns
##           type                               description
##    <character>                               <character>
## 1 intermediate mean of normalized counts for all samples
## 2      results  log2 fold change (MAP): dex trt vs untrt
## 3      results          standard error: dex trt vs untrt
## 4      results          Wald statistic: dex trt vs untrt
## 5      results       Wald test p-value: dex trt vs untrt
## 6      results                      BH adjusted p-values
head(res[order(res$padj),])    # 'top table' by p-adj
## log2 fold change (MAP): dex trt vs untrt 
## Wald test p-value: dex trt vs untrt 
## DataFrame with 6 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE      stat        pvalue          padj
##                  <numeric>      <numeric> <numeric> <numeric>     <numeric>     <numeric>
## ENSG00000152583   997.4398       4.316100 0.1724127  25.03354 2.637881e-138 4.755573e-134
## ENSG00000165995   495.0929       3.188698 0.1277441  24.96160 1.597973e-137 1.440413e-133
## ENSG00000101347 12703.3871       3.618232 0.1499441  24.13054 1.195378e-128 6.620010e-125
## ENSG00000120129  3409.0294       2.871326 0.1190334  24.12201 1.468829e-128 6.620010e-125
## ENSG00000189221  2341.7673       3.230629 0.1373644  23.51868 2.627083e-122 9.472209e-119
## ENSG00000211445 12285.6151       3.552999 0.1589971  22.34631 1.311440e-110 3.940441e-107

2.3 Some Considerations

2.3.1 Experimental design

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.

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.

2.3.2 Wet-lab

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, …

2.3.3 Sequencing

Axes of variation

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

Application-specific, e.g.,

  • RNA-seq, known genes: single- or paired-end reads
  • RNA-seq, transcripts or novel variants: paired-end reads
  • ChIP-seq: short, single-end reads are usually sufficient
  • 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)

2.3.4 Alignment / Reduction

Alignment

  • 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)

Reduction to ‘count tables’

  • 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

2.3.4.2 (kallisto / sailfish)

  • ‘Next generation’ differential expression tools; transcriptome alignment
  • E.g., kallisto takes a radically different approach: from FASTQ to count table without BAM files.
  • Very fast, almost as accurate.
  • Hints on how it works; arXiv
  • Integration with gene-level analyses – Soneson et al.

2.3.5 Statistical Considerations

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, FPKM, …), 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.

  • Detail

    • 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

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 null hypothesis used in microarray or RNA-seq experiments, the expected 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.

2.4 Comprehension: placing differentially expressed regions in context

‘Annotation’ packages

3 Challenges & Solutions

3.1 Interoperability: Standard Representations

Genomic Ranges: GenomicRanges

Data / metadata integration: SummarizedExperiment

3.2 Scalability

Efficient R programming

Parallel evaluation: BiocParallel

3.3 Reproducibility: From Script to Package

Scripts

Vignettes

Packages

4 Acknowledgments

Individuals

Annual conference: https://bioconductor.org/BioC2016

sessionInfo()
## R version 3.2.3 Patched (2016-01-28 r70038)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] BiocParallel_1.4.3                      AnnotationHub_2.2.3                    
##  [3] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.22.12                
##  [5] org.Hs.eg.db_3.2.3                      RSQLite_1.0.0                          
##  [7] DBI_0.3.1                               AnnotationDbi_1.32.3                   
##  [9] DESeq2_1.10.1                           RcppArmadillo_0.6.500.4.0              
## [11] Rcpp_0.12.3                             airway_0.104.0                         
## [13] SummarizedExperiment_1.0.2              Biobase_2.30.0                         
## [15] GenomicRanges_1.22.4                    GenomeInfoDb_1.6.3                     
## [17] IRanges_2.4.7                           S4Vectors_0.8.11                       
## [19] BiocGenerics_0.16.1                     ggplot2_2.0.0                          
## [21] ENAR2016_0.0.1                          BiocStyle_1.8.0                        
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.1.0                   splines_3.2.3                Formula_1.2-1               
##  [4] shiny_0.13.0                 interactiveDisplayBase_1.8.0 latticeExtra_0.6-28         
##  [7] Rsamtools_1.23.1             yaml_2.1.13                  lattice_0.20-33             
## [10] digest_0.6.9                 RColorBrewer_1.1-2           XVector_0.10.0              
## [13] colorspace_1.2-6             htmltools_0.3                httpuv_1.3.3                
## [16] Matrix_1.2-3                 plyr_1.8.3                   XML_3.98-1.3                
## [19] biomaRt_2.26.1               genefilter_1.52.1            zlibbioc_1.16.0             
## [22] xtable_1.8-2                 scales_0.3.0                 annotate_1.48.0             
## [25] mgcv_1.8-11                  nnet_7.3-12                  survival_2.38-3             
## [28] magrittr_1.5                 mime_0.4                     evaluate_0.8                
## [31] nlme_3.1-124                 foreign_0.8-66               BiocInstaller_1.20.1        
## [34] tools_3.2.3                  formatR_1.2.1                stringr_1.0.0               
## [37] munsell_0.4.2                locfit_1.5-9.1               cluster_2.0.3               
## [40] lambda.r_1.1.7               Biostrings_2.38.4            futile.logger_1.4.1         
## [43] grid_3.2.3                   RCurl_1.95-4.7               labeling_0.3                
## [46] bitops_1.0-6                 rmarkdown_0.9.2              gtable_0.1.2                
## [49] codetools_0.2-14             R6_2.1.2                     GenomicAlignments_1.6.3     
## [52] gridExtra_2.0.0              knitr_1.12.3                 rtracklayer_1.30.2          
## [55] Hmisc_3.17-1                 futile.options_1.0.0         stringi_1.0-1               
## [58] geneplotter_1.48.0           rpart_4.1-10                 acepack_1.3-3.3