Contents

Author: Sonali Arora (sarora@fredhutch.org)
Date: 20-22 July, 2015

The material in this course requires R version 3.2.1 and Bioconductor version 3.2

0.1 Advanced lab for Bioconductor - Explore the Data

In this lab, we will explore the dataset from airway, and then subsequently use the same to run a quick RNA-seq analysis.

Steps include -

0.2 Data for the Lab

The data used in this Lab 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.

For our analysis, we wil use data from the data package airway.

library("airway")
data(airway)

0.3 Solutions

0.3.1 Answer 1 : Load the Data

Bioconductor software packages often define and use a custom class for their data object, which makes sure that all the needed data slots are consistently provided and fulfill the requirements.

The data stored inside airway is a SummarizedExperiment object.

library("airway")
data(airway)
se <- airway
se
## class: RangedSummarizedExperiment 
## dim: 64102 8 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

For a SummarizedExperiment object, The assay(s) contains the matrix (or matrices) of summarized values, the rowData contains information about the genomic ranges, and the colData contains information about the samples or experiments.

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
colSums(assay(se))
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521 
##   20637971   18809481   25348649   15163415   24448408   30818215   19126151   21164133
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
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

In DESeq2, the custom class is called DESeqDataSet. It is built on top of the SummarizedExperiment class.

library(DESeq2)
## Loading required package: Rcpp
## Loading required package: RcppArmadillo
dds <- DESeqDataSet(se, design = ~ cell + dex)

0.3.2 Answer 2 : rlog transformation

Many common statistical methods for exploratory analysis of multidimensional data, especially methods for clustering and ordination (e.g., principal-component analysis and the like), work best for (at least approximately) homoskedastic data; this means that the variance of an observed quantity (here, the expression strength of a gene) does not depend on the mean.

In RNA-Seq data, the variance grows with the mean.If one performs PCA (principal components analysis) directly on a matrix of normalized read counts, the result typically depends only on the few most strongly expressed genes because they show the largest absolute differences between samples.

As a solution, DESeq2 offers the regularized-logarithm transformation, or rlog for short. See the help for ?rlog for more information and options.

The function rlog returns a SummarizedExperiment object which contains the rlog-transformed values in its assay slot:

rld <- rlog(dds)   
head(assay(rld))    
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## ENSG00000000003   9.399151   9.142478   9.501695   9.320796   9.757212   9.512183   9.617378
## ENSG00000000005   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
## ENSG00000000419   8.901283   9.113976   9.032567   9.063925   8.981930   9.108531   8.894830
## ENSG00000000457   7.949897   7.882371   7.834273   7.916459   7.773819   7.886645   7.946411
## ENSG00000000460   5.849521   5.882363   5.486937   5.770334   5.940407   5.663847   6.107733
## ENSG00000000938  -1.638084  -1.637483  -1.558248  -1.636072  -1.597606  -1.639362  -1.637608
##                 SRR1039521
## ENSG00000000003   9.315309
## ENSG00000000005   0.000000
## ENSG00000000419   9.052303
## ENSG00000000457   7.908338
## ENSG00000000460   5.907824
## ENSG00000000938  -1.637724

0.3.3 Answer 3 : Visualize the Data

To assess overall similarity between samples: Which samples are similar to each other, which are different? Does this fit to the expectation from the experiment’s design? We use the R function dist to calculate the Euclidean distance between samples. To avoid that the distance measure is dominated by a few highly variable genes, and have a roughly equal contribution from all genes, we use it on the rlog-transformed data

sampleDists <- dist( t( assay(rld) ) )
sampleDists
##            SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## SRR1039509   40.89060                                                                  
## SRR1039512   37.35231   50.07638                                                       
## SRR1039513   55.74569   41.49280   43.61052                                            
## SRR1039516   41.98797   53.58929   40.99513   57.10447                                 
## SRR1039517   57.69438   47.59326   53.52310   46.13742   42.10583                      
## SRR1039520   37.06633   51.80994   34.86653   52.54968   43.21786   57.13688           
## SRR1039521   56.04254   41.46514   51.90045   34.82975   58.40428   47.90244   44.78171

Note the use of the function t to transpose the data matrix. We need this because dist calculates distances between data rows and our samples constitute the columns.

0.3.3.1 Heatmaps

We visualize the sample-to-sample distances in a heatmap, using the function heatmap.2 from the gplots package. Note that we have changed the row names of the distance matrix to contain treatment type and patient number instead of sample ID, so that we have all this information in view when looking at the heatmap.

library("gplots")
library("RColorBrewer")

sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- paste( rld$dex, rld$cell, sep="-" )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
hc <- hclust(sampleDists)
heatmap.2( sampleDistMatrix, Rowv=as.dendrogram(hc),
          symm=TRUE, trace="none", col=colors,
          margins=c(2,10), labCol=FALSE )

0.3.3.2 PCA

Another way to visualize sample-to-sample distances is a principal-components analysis (PCA). In this ordination method, the data points (i.e., here, the samples) are projected onto the 2D plane such that they spread out in the two directions which explain most of the differences in the data. The x-axis is the direction (or principal component) which separates the data points the most. The amount of the total variance which is contained in the direction is printed in the axis label. Here, we have used the function plotPCA which comes with DESeq2. The two terms specified by intgroup are the interesting groups for labelling the samples; they tell the function to use them to choose colors.

plotPCA(rld, intgroup = c("dex", "cell"))

0.3.3.3 MDS

Another plot, very similar to the PCA plot, can be made using the multidimensional scaling (MDS) function in base R. This is useful when we don’t have the original data, but only a matrix of distances. Here we have the MDS plot for the distances calculated from the rlog transformed counts:

library(ggplot2)
mds <- data.frame(cmdscale(sampleDistMatrix))
mds <- cbind(mds, colData(rld))
qplot(X1,X2,color=dex,shape=cell,data=as.data.frame(mds))

0.3.4 Answer 4 : Summarize the plots

We see that the differences between cells are considerable, though not stronger than the differences due to treatment with dexamethasone. This shows why it will be important to account for this in differential testing by using a paired design (‘paired’, because each dex treated sample is paired with one untreated sample from the same cell line). We are already set up for this by using the design formula ~ cell + dex when setting up the data object in the beginning.

0.4 References

For a much detailed analysis see
- Case Study- How to build a SummarizedExperiment - airway dataset
- Exploring the Data using Machine Learning Techniques

0.5 What to not miss at BioC2015 !

If you liked this lab and want to learn more in this area, do not miss the following labs at BioC2015

0.6 sessionInfo()

sessionInfo()
## R version 3.2.1 (2015-06-18)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=C               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] RColorBrewer_1.1-2                      gplots_2.17.0                          
##  [3] DESeq2_1.9.23                           RcppArmadillo_0.5.200.1.0              
##  [5] Rcpp_0.11.6                             airway_0.103.1                         
##  [7] Rattus.norvegicus_1.3.1                 org.Rn.eg.db_3.1.2                     
##  [9] GO.db_3.1.2                             OrganismDbi_1.11.42                    
## [11] BSgenome.Rnorvegicus.UCSC.rn5_1.4.0     BSgenome_1.37.3                        
## [13] rtracklayer_1.29.12                     TxDb.Rnorvegicus.UCSC.rn5.refGene_3.1.3
## [15] org.Hs.eg.db_3.1.2                      RSQLite_1.0.0                          
## [17] DBI_0.3.1                               TxDb.Hsapiens.UCSC.hg19.knownGene_3.1.3
## [19] GenomicFeatures_1.21.13                 AnnotationDbi_1.31.17                  
## [21] AnnotationHub_2.1.30                    RNAseqData.HNRNPC.bam.chr14_0.7.0      
## [23] GenomicAlignments_1.5.11                Rsamtools_1.21.14                      
## [25] Biostrings_2.37.2                       XVector_0.9.1                          
## [27] SummarizedExperiment_0.3.2              Biobase_2.29.1                         
## [29] GenomicRanges_1.21.16                   GenomeInfoDb_1.5.8                     
## [31] IRanges_2.3.14                          S4Vectors_0.7.10                       
## [33] BiocGenerics_0.15.3                     ggplot2_1.0.1                          
## [35] BiocStyle_1.7.4                        
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6                 httr_1.0.0                   tools_3.2.1                 
##  [4] R6_2.1.0                     KernSmooth_2.23-15           rpart_4.1-10                
##  [7] Hmisc_3.16-0                 colorspace_1.2-6             nnet_7.3-10                 
## [10] gridExtra_2.0.0              curl_0.9.1                   graph_1.47.2                
## [13] formatR_1.2                  labeling_0.3                 caTools_1.17.1              
## [16] scales_0.2.5                 genefilter_1.51.0            RBGL_1.45.1                 
## [19] stringr_1.0.0                digest_0.6.8                 foreign_0.8-65              
## [22] rmarkdown_0.7                htmltools_0.2.6              BiocInstaller_1.19.8        
## [25] shiny_0.12.1                 BiocParallel_1.3.34          gtools_3.5.0                
## [28] acepack_1.3-3.3              RCurl_1.95-4.7               magrittr_1.5                
## [31] Formula_1.2-1                futile.logger_1.4.1          munsell_0.4.2               
## [34] proto_0.3-10                 stringi_0.5-5                yaml_2.1.13                 
## [37] MASS_7.3-43                  zlibbioc_1.15.0              plyr_1.8.3                  
## [40] grid_3.2.1                   gdata_2.17.0                 lattice_0.20-33             
## [43] splines_3.2.1                annotate_1.47.1              locfit_1.5-9.1              
## [46] knitr_1.10.5                 geneplotter_1.47.0           reshape2_1.4.1              
## [49] codetools_0.2-14             biomaRt_2.25.1               futile.options_1.0.0        
## [52] XML_3.98-1.3                 evaluate_0.7                 latticeExtra_0.6-26         
## [55] lambda.r_1.1.7               httpuv_1.3.2                 gtable_0.1.2                
## [58] mime_0.3                     xtable_1.7-4                 survival_2.38-3             
## [61] cluster_2.0.2                interactiveDisplayBase_1.7.0