Introduction to R and Bioconductor

Vectors, classes, methods, and packages

R

  • Statistical
  • Vectorized

    x <- rnorm(100)
    y <- x + rnorm(100)
  • Large data per se is not a problem

    system.time(rnorm(1e7))
    ##    user  system elapsed 
    ##   0.928   0.000   0.929
  • Classes structure data (nouns)

    df <- data.frame(X=x, Y=y)
    fit <- lm(Y ~ X, df)
  • Methods transform objects (verbs)

    anova(fit)
    ## Analysis of Variance Table
    ## 
    ## Response: Y
    ##           Df  Sum Sq Mean Sq F value    Pr(>F)    
    ## X          1  91.006  91.006  85.776 4.851e-15 ***
    ## Residuals 98 103.976   1.061                      
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • R provides introspection

    class(df)
    ## [1] "data.frame"
    dim(df)
    ## [1] 100   2
    class(lm)
    ## [1] "function"
    methods(class=class(lm))
    ##  [1] as.data.frame        as.list              coerce              
    ##  [4] coerce<-             dispersionFunction<- fortify             
    ##  [7] head                 plot                 print               
    ## [10] tail                
    ## see '?methods' for accessing help and source code

Packages

  • Extend base functionality

    library(ggplot2)
    ggplot(df, aes(x=X, y=Y)) +
        geom_point() +
        stat_smooth(method="lm")
    ## Warning in grid.Call.graphics(L_polygon, x$x, x$y, index): semi-
    ## transparency is not supported on this device: reported only once per page

Bioconductor – key software and annotation resources

  • GenomicRanges

    library(GenomicRanges)
    library(TxDb.Hsapiens.UCSC.hg38.knownGene)
    ex <- exons(TxDb.Hsapiens.UCSC.hg38.knownGene)
    ex
    ## GRanges object with 581036 ranges and 1 metadata column:
    ##                    seqnames           ranges strand |   exon_id
    ##                       <Rle>        <IRanges>  <Rle> | <integer>
    ##        [1]             chr1   [29554, 30039]      + |         1
    ##        [2]             chr1   [30267, 30667]      + |         2
    ##        [3]             chr1   [30366, 30503]      + |         3
    ##        [4]             chr1   [30564, 30667]      + |         4
    ##        [5]             chr1   [30976, 31097]      + |         5
    ##        ...              ...              ...    ... .       ...
    ##   [581032] chrUn_KI270750v1 [148668, 148843]      + |    581032
    ##   [581033] chrUn_KI270752v1 [   144,    268]      + |    581033
    ##   [581034] chrUn_KI270752v1 [ 21813,  21944]      + |    581034
    ##   [581035] chrUn_KI270752v1 [  3497,   3623]      - |    581035
    ##   [581036] chrUn_KI270752v1 [  9943,  10067]      - |    581036
    ##   -------
    ##   seqinfo: 455 sequences (1 circular) from hg38 genome
    keepStandardChromosomes(ex)
    ## GRanges object with 538594 ranges and 1 metadata column:
    ##            seqnames         ranges strand |   exon_id
    ##               <Rle>      <IRanges>  <Rle> | <integer>
    ##        [1]     chr1 [29554, 30039]      + |         1
    ##        [2]     chr1 [30267, 30667]      + |         2
    ##        [3]     chr1 [30366, 30503]      + |         3
    ##        [4]     chr1 [30564, 30667]      + |         4
    ##        [5]     chr1 [30976, 31097]      + |         5
    ##        ...      ...            ...    ... .       ...
    ##   [538590]     chrM [ 5826,  5891]      - |    538590
    ##   [538591]     chrM [ 7446,  7514]      - |    538591
    ##   [538592]     chrM [14149, 14673]      - |    538592
    ##   [538593]     chrM [14674, 14742]      - |    538593
    ##   [538594]     chrM [15956, 16023]      - |    538594
    ##   -------
    ##   seqinfo: 25 sequences (1 circular) from hg38 genome
    hist(log10(width(ex)))

    introns <- intronsByTranscript(TxDb.Hsapiens.UCSC.hg38.knownGene,
        use.names=TRUE)
    hist(log10(lengths(introns)))

    introns[ which.max(lengths(introns)) ]
    ## GRangesList object of length 1:
    ## $uc031rqc.3 
    ## GRanges object with 362 ranges and 0 metadata columns:
    ##         seqnames                 ranges strand
    ##            <Rle>              <IRanges>  <Rle>
    ##     [1]     chr2 [178527308, 178527445]      -
    ##     [2]     chr2 [178527749, 178528273]      -
    ##     [3]     chr2 [178528428, 178528527]      -
    ##     [4]     chr2 [178529220, 178529959]      -
    ##     [5]     chr2 [178530117, 178530240]      -
    ##     ...      ...                    ...    ...
    ##   [358]     chr2 [178799732, 178799824]      -
    ##   [359]     chr2 [178799911, 178800394]      -
    ##   [360]     chr2 [178800683, 178802137]      -
    ##   [361]     chr2 [178802342, 178804551]      -
    ##   [362]     chr2 [178804656, 178807211]      -
    ## 
    ## -------
    ## seqinfo: 455 sequences (1 circular) from hg38 genome
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    library(RNAseqData.HNRNPC.bam.chr14)
    
    library(GenomicAlignments)
    fname <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
    aln <- readGAlignments(fname)
    
    genes <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
    chr14 <- genes[seqnames(genes) == "chr14"]
    
    hist(log10(1 + countOverlaps(chr14, aln)))

  • Biostrings

    library(Biostrings)
    library(BSgenome.Hsapiens.UCSC.hg38)
    
    ex <- exons(TxDb.Hsapiens.UCSC.hg38.knownGene)
    ex <- keepStandardChromosomes(ex)
    dna <- getSeq(BSgenome.Hsapiens.UCSC.hg38, ex)
    dna
    ##   A DNAStringSet instance of length 538594
    ##           width seq
    ##      [1]    486 GTGCACACGGCTCCCATGCGTTGTCTTC...GCGGAGTTGGGCAGCTCGTGTTCAATGG
    ##      [2]    401 TCATCAGTCCAAAGTCCAGCAGTTGTCC...GCTGCAGCGTGTGGTCCCCTTACCAGAG
    ##      [3]    138 GGATGCCCAGCTAGTTTGAATTTTAGAT...TATCTGAGATTCAGAATTAAGCATTTTA
    ##      [4]    104 TTTTGTCCGCCTTCCCTGCCTCCTCTTC...GCTGCAGCGTGTGGTCCCCTTACCAGAG
    ##      [5]    122 TGAGGATGCGAAGAGAAGGTGGCTGTCT...CCTCCAGAACTGTGAGGGATAAATGTAT
    ##      ...    ... ...
    ## [538590]     66 GGTAAAATGGCTGAGTGAAGCATTGGAC...AAGACAGGGGTTAGGCCTCTTTTTACCA
    ## [538591]     69 GAAAAAGTCATGGAGGCCATGGGGTTGG...TTGGGGGGTTCGATTCCTTCCTTTTTTG
    ## [538592]    525 ATGATGTATGCTTTGTTTCTGTTGAGTG...TGTAATTGAGATTGCTCGGGGGAATAGG
    ## [538593]     69 GTTCTTGTAGTTGAAATACAACGATGGT...GGTCGTGGTTGTAGTCCGTGCGAGAATA
    ## [538594]     68 CAGAGAATAGTTTAAATTAGAATCTTAG...ATGGTGGAGTTAAAGACTTTTTCTCTGA
    gc <- letterFrequency(dna, "GC", as.prob=TRUE)
    head(gc)
    ##            G|C
    ## [1,] 0.6028807
    ## [2,] 0.4513716
    ## [3,] 0.3115942
    ## [4,] 0.5576923
    ## [5,] 0.5409836
    ## [6,] 0.5149254
    plot(density(gc))

  • SummarizedExperiment

    library(SummarizedExperiment)
    library(airway)
    data(airway)
    airway
    ## class: RangedSummarizedExperiment 
    ## dim: 64102 8 
    ## metadata(1): ''
    ## assays(1): counts
    ## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
    ## rowData names(0):
    ## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
    ## colData names(9): SampleName cell ... Sample BioSample
    head(assay(airway))
    ##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
    ## ENSG00000000003        679        448        873        408       1138
    ## ENSG00000000005          0          0          0          0          0
    ## ENSG00000000419        467        515        621        365        587
    ## ENSG00000000457        260        211        263        164        245
    ## ENSG00000000460         60         55         40         35         78
    ## ENSG00000000938          0          0          2          0          1
    ##                 SRR1039517 SRR1039520 SRR1039521
    ## ENSG00000000003       1047        770        572
    ## ENSG00000000005          0          0          0
    ## ENSG00000000419        799        417        508
    ## ENSG00000000457        331        233        229
    ## ENSG00000000460         63         76         60
    ## ENSG00000000938          0          0          0
    colData(airway)
    ## 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
    rowRanges(airway)
    ## 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
    airway[, airway$dex == "trt"]
    ## class: RangedSummarizedExperiment 
    ## dim: 64102 4 
    ## metadata(1): ''
    ## assays(1): counts
    ## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
    ## rowData names(0):
    ## colnames(4): SRR1039509 SRR1039513 SRR1039517 SRR1039521
    ## colData names(9): SampleName cell ... Sample BioSample
  • The Bioconductor sequencing ecosystem

Help!

?

  • ?"findOverlaps<tab>"
  • ?GenomicRanges::findOverlaps

biocViews and landing pages

Vigettes, workflows & training material

RNA-seq

This part of the workshop uses data from a gene-level RNA-seq experiment involving airway smooth muscle cells; details are provided in the vignette accompanying the airway package. The original data is from Himes et al., “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. From the Abstract: “Using RNA-Seq […] we characterized transcriptomic changes in four primary human ASM cell lines that were treated with dexamethasone - a potent synthetic glucocorticoid (1 micromolar for 18 hours).”

We join the analysis after the sequencing, read aligment, and summary of aligned reads to a table of counts of reads overlapping regions of interest (genes) in each sample. We focus on a subset of the experiment, with 4 cell lines each treated with dexamethasone or a control.

More extensive work flows are available for use of DESeq2 and edgeR / limma voom.

Setup (not necessary during the workshop)

We’ll use two packages from Bioconductor. These depend in turn on several other packages. The packages and their dependencies are already installed on the Amazon Machine Instance (AMI) used in the course. For your own computers after the course, install packages with

source("https://bioconductor.org/biocLite.R")
biocLite(c("DESeq2", "org.Hs.eg.db"))

Installation needs to be performed once per computer, not every time the packages are used.

We use two data files in the analysis. The data files are in the Brisbane2016 github repository. Install the github repository with

biocLite("Bioconductor/Brisbane2016")

Once the package is installed, the location of the files (to be used in file.choose(), below) is given by

system.file(package="Brisbane2016", "extdata")

Sequence analysis work flows

Research question

Data processing steps

  • Quality assessment: FASTQ files; outside R
  • Alignment and summary to count table
  • Assessment of differential expression
  • Results placed in context, e.g., gene set enrichment

Alignment & summary (traditional)

  • Full alignment to BAM files, summarizing gene or transcript abundance, e.g., Bowtie / tophat / cufflinks; RSEM; Rsubread
  • Summarize to gene-level count tables or estimates of abundance
  • Counts are important: information about statistical uncertainty of estimate

Alignment & summary (contemporary)

  • Approximate alignment directly to count tables of transcripts or genes, e.g.,kallisto, salmon

Exploring count tables

We start with the ‘phenotypic’ data, describing the samples used in the experiment. The data is a simple table of 8 rows and several columns; it could be created in Excel and exported as a tab-delimited file. Find the location of the file on your Amazon machine instance

colDataFile <- file.chooose()    # find 'airway-colData.tab'

and read the data in to R using the read.table() function. The data is small enough to be viewed in the R session by typing the name of the variable) or in RStudio (by using View() or double-clicking on the variable in the ‘Environment’ tab).

colData <- read.table(colDataFile)
colData
##            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
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520       101  SRX384357
## SRR1039521 GSM1275875 N061011   trt untrt SRR1039521        98  SRX384358
##               Sample    BioSample
## SRR1039508 SRS508568 SAMN02422669
## SRR1039509 SRS508567 SAMN02422675
## SRR1039512 SRS508571 SAMN02422678
## SRR1039513 SRS508572 SAMN02422670
## SRR1039516 SRS508575 SAMN02422682
## SRR1039517 SRS508576 SAMN02422673
## SRR1039520 SRS508579 SAMN02422683
## SRR1039521 SRS508580 SAMN02422677

This should go smoothly; with real data one often needs to spend considerable time adjusting arguments to read.table() to account for the presence of a header, row names, comments, etc.

The next challenge is to input the expression estimates. This is a matrix of rows representing regions of interest (genes) and columns representing samples. Entries in the matrix are the number of reads overlapping each region in each sample. It is important that the values are raw counts, rather than scaled measures such as FPKM. Find the file

assayFile <- file.chooose()    # find 'airway-assay.tab'

Input the data and use head() to view the first few rows of the data.

assay <- read.table(assayFile)
head(assay)
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003        679        448        873        408       1138
## ENSG00000000419        467        515        621        365        587
## ENSG00000000457        260        211        263        164        245
## ENSG00000000460         60         55         40         35         78
## ENSG00000000938          0          0          2          0          1
## ENSG00000000971       3251       3679       6177       4252       6721
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003       1047        770        572
## ENSG00000000419        799        417        508
## ENSG00000000457        331        233        229
## ENSG00000000460         63         76         60
## ENSG00000000938          0          0          0
## ENSG00000000971      11027       5176       7995

NOTE Most modern work flows will start with data generated from kallisto, salmon, etc.; see the tximport package and vignette for importing data from these software tools into R.

Calculate the ‘library size’ (total number of mapped reads) of each sample using colSums()

colSums(assay)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 
##   20637971   18809481   25348649   15163415   24448408   30818215 
## SRR1039520 SRR1039521 
##   19126151   21164133

Create a density plot of the average asinh-transformed (asinh is log-like, except near zero) read counts of each gene using the following series of commands.

plot(density(rowMeans(asinh(assay))))

Multi-dimensional scaling (MDS) is a dimensionality reduction method that takes vectors in n-space and projects them into two (or more) dimensions. Use the dist() function to calculate the (Euclidean) distance bewteen each sample, and the base R function cmdscale() to perform MDS on the distance matrix. We can use plot() to visualize the results and see the approximate location of each of the 8 samples. Use the argument col to color the points based on cell line (colData$cell) or experimental treatment colData$dex.

d <- dist(t(asinh(assay)))
plot(cmdscale(d), pch=19, cex=2)

plot(cmdscale(d), pch=19, cex=2, col=colData$cell)

plot(cmdscale(d), pch=19, cex=2, col=colData$dex)

Note that cell lines are relatively similar to one another. This suggests that cell line should be used as a covariate in subsequent analysis.

Differential expression

We will use the DESeq2 package for differential expression analysis; other choices are possible, notably edgeR and limma.

The analysis starts by providing the expression count data, a description of the experiment, and a ‘model’ that describes the statistical relationship we’d like to estimate. For our model and based in part on the exploratory analysis of the previous section, we suppose that count is determined cell line and dexamethasone treatment. We include cell line primarily as a covariate; our primary interest is in dexamethasone.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(assay, colData, ~ cell + dex)

The analysis is extremely straight-forward to invoke, but the calculations involve a number of sophisticated statistical issues, including:

  • Robust estimate of library size
  • Use of a negative binomial model to describe the relationship between experimental design and counts
  • Use of moderated test statistics that balance gene-level parameter estimates with experiment wide estimates
  • Data-driven approaches to outlier detection
  • Filtering and control for multiple testing

The code is invoked as:

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
dds
## class: DESeqDataSet 
## dim: 33469 8 
## metadata(1): version
## assays(3): counts mu cooks
## rownames(33469): ENSG00000000003 ENSG00000000419 ...
##   ENSG00000273492 ENSG00000273493
## rowData names(46): baseMean baseVar ... deviance maxCooks
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(10): SampleName cell ... BioSample sizeFactor

The DESeq() function returns an object that can be used as a starting point for further analysis, for instance generating a ‘top table’ of differentially expressed genes, orderd by adjusted (for multiple comparison) P values.

result <- results(dds)
result
## log2 fold change (MAP): dex untrt vs trt 
## Wald test p-value: dex untrt vs trt 
## DataFrame with 33469 rows and 6 columns
##                    baseMean log2FoldChange      lfcSE       stat
##                   <numeric>      <numeric>  <numeric>  <numeric>
## ENSG00000000003 708.6021697     0.37415246 0.09884435  3.7852692
## ENSG00000000419 520.2979006    -0.20206175 0.10974241 -1.8412367
## ENSG00000000457 237.1630368    -0.03616686 0.13834540 -0.2614244
## ENSG00000000460  57.9326331     0.08445399 0.24990709  0.3379415
## ENSG00000000938   0.3180984     0.08413901 0.15133424  0.5559813
## ...                     ...            ...        ...        ...
## ENSG00000273487   8.1632350    -0.55007238  0.3725061 -1.4766803
## ENSG00000273488   8.5844790    -0.10513006  0.3683837 -0.2853820
## ENSG00000273489   0.2758994    -0.06947899  0.1512520 -0.4593591
## ENSG00000273492   0.1059784     0.02314357  0.1512520  0.1530133
## ENSG00000273493   0.1061417     0.02314357  0.1512520  0.1530133
##                       pvalue        padj
##                    <numeric>   <numeric>
## ENSG00000000003 0.0001535423 0.001279829
## ENSG00000000419 0.0655868795 0.196015831
## ENSG00000000457 0.7937652416 0.912936622
## ENSG00000000460 0.7354072415 0.883203048
## ENSG00000000938 0.5782236287          NA
## ...                      ...         ...
## ENSG00000273487    0.1397614   0.3376252
## ENSG00000273488    0.7753515   0.9032267
## ENSG00000273489    0.6459763          NA
## ENSG00000273492    0.8783878          NA
## ENSG00000273493    0.8783878          NA
ridx <- head(order(result$padj), 10)
top = result[ridx,]
top
## log2 fold change (MAP): dex untrt vs trt 
## Wald test p-value: dex untrt vs trt 
## DataFrame with 10 rows and 6 columns
##                   baseMean log2FoldChange      lfcSE      stat
##                  <numeric>      <numeric>  <numeric> <numeric>
## ENSG00000152583   997.4398      -4.313962 0.17213733 -25.06116
## ENSG00000165995   495.0929      -3.186823 0.12815654 -24.86664
## ENSG00000101347 12703.3871      -3.618734 0.14894336 -24.29604
## ENSG00000120129  3409.0294      -2.871488 0.11824908 -24.28338
## ENSG00000189221  2341.7673      -3.230395 0.13667447 -23.63569
## ENSG00000211445 12285.6151      -3.553360 0.15798211 -22.49217
## ENSG00000157214  3009.2632      -1.948723 0.08867432 -21.97618
## ENSG00000162614  5393.1017      -2.003487 0.09269629 -21.61345
## ENSG00000125148  3656.2528      -2.167122 0.10354724 -20.92882
## ENSG00000154734 30315.1355      -2.286778 0.11308412 -20.22192
##                        pvalue          padj
##                     <numeric>     <numeric>
## ENSG00000152583 1.319237e-138 2.360906e-134
## ENSG00000165995 1.708565e-136 1.528824e-132
## ENSG00000101347 2.158637e-130 1.287699e-126
## ENSG00000120129 2.937247e-130 1.314124e-126
## ENSG00000189221 1.656535e-123 5.929070e-120
## ENSG00000211445 4.952260e-112 1.477094e-108
## ENSG00000157214 4.867315e-107 1.244364e-103
## ENSG00000162614 1.342345e-103 3.002826e-100
## ENSG00000125148  2.926277e-97  5.818739e-94
## ENSG00000154734  6.278709e-91  1.123638e-87

Two common ways of visualizing results are as “MA” plots that summarize the relationship between count and statistical sigificance

plotMA(result)
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): semi-
## transparency is not supported on this device: reported only once per page

and volcano plots where the X-axis corresponds to ‘biology’ (size of change) and Y-axis to ‘statistics’ (significance of observed change).

plot(-log10(pvalue) ~ log2FoldChange, result)

A fun sanity check is the histogram of raw P-values, which should be uniform under the null; deviations often suggest batch effects or other artifacts that should be addressed in upstream analysis.

hist(result$pvalue)

Annotation and variant annotation

Simple annotation of statistical results

There are many opportunities to place the statistical results into biological context. An initial step is to map the cryptic Ensembl gene identifiers used to label regions of interest to more famliar HGNC gene symbols. For this we use the org.Hs.eg.db package, an example of a Bioconductor ‘annotation’ package containing curated data derived from public-domain resources and updated semi-annually. The mapId() function maps between identifier types, in our case to SYMBOL gene ids from ENSEMBL ids. We add these to the top table results so that they can be processed together with the statistical results.

library(org.Hs.eg.db)
top$Symbol <- mapIds(org.Hs.eg.db, rownames(top), "SYMBOL", "ENSEMBL")
## 'select()' returned 1:1 mapping between keys and columns
top
## log2 fold change (MAP): dex untrt vs trt 
## Wald test p-value: dex untrt vs trt 
## DataFrame with 10 rows and 7 columns
##                   baseMean log2FoldChange      lfcSE      stat
##                  <numeric>      <numeric>  <numeric> <numeric>
## ENSG00000152583   997.4398      -4.313962 0.17213733 -25.06116
## ENSG00000165995   495.0929      -3.186823 0.12815654 -24.86664
## ENSG00000101347 12703.3871      -3.618734 0.14894336 -24.29604
## ENSG00000120129  3409.0294      -2.871488 0.11824908 -24.28338
## ENSG00000189221  2341.7673      -3.230395 0.13667447 -23.63569
## ENSG00000211445 12285.6151      -3.553360 0.15798211 -22.49217
## ENSG00000157214  3009.2632      -1.948723 0.08867432 -21.97618
## ENSG00000162614  5393.1017      -2.003487 0.09269629 -21.61345
## ENSG00000125148  3656.2528      -2.167122 0.10354724 -20.92882
## ENSG00000154734 30315.1355      -2.286778 0.11308412 -20.22192
##                        pvalue          padj      Symbol
##                     <numeric>     <numeric> <character>
## ENSG00000152583 1.319237e-138 2.360906e-134     SPARCL1
## ENSG00000165995 1.708565e-136 1.528824e-132      CACNB2
## ENSG00000101347 2.158637e-130 1.287699e-126      SAMHD1
## ENSG00000120129 2.937247e-130 1.314124e-126       DUSP1
## ENSG00000189221 1.656535e-123 5.929070e-120        MAOA
## ENSG00000211445 4.952260e-112 1.477094e-108        GPX3
## ENSG00000157214 4.867315e-107 1.244364e-103      STEAP2
## ENSG00000162614 1.342345e-103 3.002826e-100        NEXN
## ENSG00000125148  2.926277e-97  5.818739e-94        MT2A
## ENSG00000154734  6.278709e-91  1.123638e-87     ADAMTS1

Sneak peak:

biocLite("Bioconductor/Organism.dplyr")
library(Organism.dplyr)
src <- src_ucsc("Homo sapiens")
## using org.Hs.eg.db, TxDb.Hsapiens.UCSC.hg38.knownGene
src
## src:  sqlite 3.8.6 [/home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.4-Bioc-3.5/org.Hs.eg.db/extdata/org.Hs.eg.sqlite]
## tbls: id, id_accession, id_go, id_go_all, id_omim_pm, id_protein,
##   id_transcript, ranges_cds, ranges_exon, ranges_gene, ranges_tx
tbl(src, "id")
## Source:   query [?? x 6]
## Database: sqlite 3.8.6 [/home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.4-Bioc-3.5/org.Hs.eg.db/extdata/org.Hs.eg.sqlite]
## 
##    entrez      map         ensembl symbol               genename    alias
##     <chr>    <chr>           <chr>  <chr>                  <chr>    <chr>
## 1       1  19q13.4 ENSG00000121410   A1BG alpha-1-B glycoprotein      A1B
## 2       1  19q13.4 ENSG00000121410   A1BG alpha-1-B glycoprotein      ABG
## 3       1  19q13.4 ENSG00000121410   A1BG alpha-1-B glycoprotein      GAB
## 4       1  19q13.4 ENSG00000121410   A1BG alpha-1-B glycoprotein HYST2477
## 5       1  19q13.4 ENSG00000121410   A1BG alpha-1-B glycoprotein     A1BG
## 6      10     8p22 ENSG00000156006   NAT2  N-acetyltransferase 2     AAC2
## 7      10     8p22 ENSG00000156006   NAT2  N-acetyltransferase 2    NAT-2
## 8      10     8p22 ENSG00000156006   NAT2  N-acetyltransferase 2     PNAT
## 9      10     8p22 ENSG00000156006   NAT2  N-acetyltransferase 2     NAT2
## 10    100 20q13.12 ENSG00000196839    ADA    adenosine deaminase      ADA
## # ... with more rows
ids <- head(rownames(top))
tbl(src, "id") %>% filter(ensembl %in% ids) %>% 
    select(ensembl, symbol, genename)
## Source:   query [?? x 3]
## Database: sqlite 3.8.6 [/home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.4-Bioc-3.5/org.Hs.eg.db/extdata/org.Hs.eg.sqlite]
## 
##           ensembl  symbol
##             <chr>   <chr>
## 1 ENSG00000101347  SAMHD1
## 2 ENSG00000120129   DUSP1
## 3 ENSG00000152583 SPARCL1
## 4 ENSG00000165995  CACNB2
## 5 ENSG00000189221    MAOA
## 6 ENSG00000211445    GPX3
## # ... with 1 more variables: genename <chr>
id <- "uc031rqc.3"  ## transcript with most exons 
t1 <- tbl(src, "id") %>% select(entrez, symbol, genename)
t2 <- tbl(src, "ranges_tx") %>% filter(tx_name == id)
inner_join(t1, t2)
## Joining, by = "entrez"
## Source:   query [?? x 10]
## Database: sqlite 3.8.6 [/home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.4-Bioc-3.5/org.Hs.eg.db/extdata/org.Hs.eg.sqlite]
## 
##   entrez symbol genename tx_chrom tx_strand  tx_start    tx_end tx_id
##    <chr>  <chr>    <chr>    <chr>     <chr>     <int>     <int> <int>
## 1   7273    TTN    titin     chr2         - 178525989 178807423 27531
## # ... with 2 more variables: tx_name <chr>, tx_type <chr>
inner_join(t1, t2) %>% collect %>% GRanges
## Joining, by = "entrez"
## GRanges object with 1 range and 6 metadata columns:
##       seqnames                 ranges strand |      entrez      symbol
##          <Rle>              <IRanges>  <Rle> | <character> <character>
##   [1]     chr2 [178525989, 178807423]      - |        7273         TTN
##          genename     tx_id     tx_name     tx_type
##       <character> <integer> <character> <character>
##   [1]       titin     27531  uc031rqc.3        <NA>
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Working with called variants

library(VariantAnnotation)
CEUvcf <- system.file(package="VariantFiltering", "extdata", "CEUtrio.vcf.bgz")
vcf <- readVcf(CEUvcf)
vcf
## class: CollapsedVCF 
## dim: 1000 3 
## rowRanges(vcf):
##   GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
## info(vcf):
##   DataFrame with 26 columns: AC, AF, AN, BaseQRankSum, CCC, ClippingRan...
## info(header(vcf)):
##                    Number Type    Description                             
##    AC              A      Integer Allele count in genotypes, for each A...
##    AF              A      Float   Allele Frequency, for each ALT allele...
##    AN              1      Integer Total number of alleles in called gen...
##    BaseQRankSum    1      Float   Z-score from Wilcoxon rank sum test o...
##    CCC             1      Integer Number of called chromosomes            
##    ClippingRankSum 1      Float   Z-score From Wilcoxon rank sum test o...
##    DB              0      Flag    dbSNP Membership                        
##    DP              1      Integer Approximate read depth; some reads ma...
##    DS              0      Flag    Were any of the samples downsampled?    
##    END             1      Integer Stop position of the interval           
##    FS              1      Float   Phred-scaled p-value using Fisher's e...
##    GQ_MEAN         1      Float   Mean of all GQ values                   
##    GQ_STDDEV       1      Float   Standard deviation of all GQ values     
##    HWP             1      Float   P value from test of Hardy Weinberg E...
##    HaplotypeScore  1      Float   Consistency of the site with at most ...
##    InbreedingCoeff 1      Float   Inbreeding coefficient as estimated f...
##    MLEAC           A      Integer Maximum likelihood expectation (MLE) ...
##    MLEAF           A      Float   Maximum likelihood expectation (MLE) ...
##    MQ              1      Float   RMS Mapping Quality                     
##    MQ0             1      Integer Total Mapping Quality Zero Reads        
##    MQRankSum       1      Float   Z-score From Wilcoxon rank sum test o...
##    NCC             1      Integer Number of no-called samples             
##    QD              1      Float   Variant Confidence/Quality by Depth     
##    ReadPosRankSum  1      Float   Z-score from Wilcoxon rank sum test o...
##    SOR             1      Float   Symmetric Odds Ratio of 2x2 contingen...
##    set             1      String  Source VCF for the merged record in C...
## geno(vcf):
##   SimpleList of length 9: GT, AD, DP, GQ, MIN_DP, PGT, PID, PL, SB
## geno(header(vcf)):
##           Number Type    Description                                      
##    GT     1      String  Genotype                                         
##    AD     .      Integer Allelic depths for the ref and alt alleles in ...
##    DP     1      Integer Approximate read depth (reads with MQ=255 or w...
##    GQ     1      Integer Genotype Quality                                 
##    MIN_DP 1      Integer Minimum DP observed within the GVCF block        
##    PGT    1      String  Physical phasing haplotype information, descri...
##    PID    1      String  Physical phasing ID information, where each un...
##    PL     G      Integer Normalized, Phred-scaled likelihoods for genot...
##    SB     4      Integer Per-sample component statistics which comprise...
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
seqlevelsStyle(vcf) <- seqlevelsStyle(TxDb.Hsapiens.UCSC.hg19.knownGene)
genome(vcf) <- "hg19"
seqlevels(vcf, force=TRUE) <- paste0("chr", 1:22)
vars <- locateVariants(vcf, TxDb.Hsapiens.UCSC.hg19.knownGene, AllVariants())
vars
## GRanges object with 1856 ranges and 9 metadata columns:
##     seqnames                 ranges strand |   LOCATION  LOCSTART
##        <Rle>              <IRanges>  <Rle> |   <factor> <integer>
##   1     chr7 [107341641, 107341641]      + |     coding      1803
##   2     chr7 [107341641, 107341641]      + |     coding       564
##   3     chr7 [107341641, 107341641]      + |     coding       486
##   4     chr7 [107341641, 107341641]      + |     coding       510
##        chr11 [ 61719380,  61719380]      + |     intron      1481
##   .      ...                    ...    ... .        ...       ...
##        chr22   [51042513, 51042513]      + |     intron      2483
##        chr22   [51042513, 51042513]      + | spliceSite      <NA>
##        chr22   [51042513, 51042513]      + |     intron       250
##        chr22   [51042513, 51042513]      + | spliceSite      <NA>
##        chr22   [51042513, 51042513]      + |   promoter      <NA>
##        LOCEND   QUERYID        TXID         CDSID      GENEID
##     <integer> <integer> <character> <IntegerList> <character>
##   1      1803         1       29193         87760        5172
##   2       564         1       29194         87760        5172
##   3       486         1       29195         87760        5172
##   4       510         1       29196         87760        5172
##          1481         2       41898                      7439
##   .       ...       ...         ...           ...         ...
##          2483       743       74402                     23542
##          <NA>       743       74402                     23542
##           250       743       74403                     23542
##          <NA>       743       74403                     23542
##          <NA>       743       74404                     23542
##           PRECEDEID        FOLLOWID
##     <CharacterList> <CharacterList>
##   1                                
##   2                                
##   3                                
##   4                                
##                                    
##   .             ...             ...
##                                    
##                                    
##                                    
##                                    
##                                    
##   -------
##   seqinfo: 22 sequences from hg19 genome

Reproducibility

The calculations here are made more reproducible by reporting the version of software used in the analysis, as follows:

sessionInfo()
## R Under development (unstable) (2016-10-21 r71563)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
## 
## 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] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] Organism.dplyr_0.0.6                       
##  [2] dplyr_0.5.0                                
##  [3] GenomicAlignments_1.11.0                   
##  [4] RNAseqData.HNRNPC.bam.chr14_0.13.0         
##  [5] airway_0.109.0                             
##  [6] org.Hs.eg.db_3.4.0                         
##  [7] BSgenome.Hsapiens.1000genomes.hs37d5_0.99.1
##  [8] BSgenome.Hsapiens.UCSC.hg38_1.4.1          
##  [9] BSgenome_1.43.0                            
## [10] rtracklayer_1.35.0                         
## [11] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2    
## [12] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.0    
## [13] GenomicFeatures_1.27.0                     
## [14] AnnotationDbi_1.37.0                       
## [15] VariantFiltering_1.11.0                    
## [16] VariantAnnotation_1.21.0                   
## [17] Rsamtools_1.27.2                           
## [18] DESeq2_1.15.2                              
## [19] SummarizedExperiment_1.5.1                 
## [20] Biobase_2.35.0                             
## [21] GenomicRanges_1.27.2                       
## [22] GenomeInfoDb_1.11.0                        
## [23] Biostrings_2.43.0                          
## [24] XVector_0.15.0                             
## [25] IRanges_2.9.0                              
## [26] S4Vectors_0.13.1                           
## [27] BiocGenerics_0.21.0                        
## [28] ggplot2_2.1.0                              
## [29] BiocInstaller_1.25.2                       
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.2.1                    AnnotationHub_2.7.0          
##  [3] splines_3.4.0                 Formula_1.2-1                
##  [5] shiny_0.14.1                  assertthat_0.1               
##  [7] interactiveDisplayBase_1.13.0 latticeExtra_0.6-28          
##  [9] RBGL_1.51.0                   yaml_2.1.13                  
## [11] RSQLite_1.0.0                 lattice_0.20-34              
## [13] biovizBase_1.23.1             chron_2.3-47                 
## [15] digest_0.6.10                 RColorBrewer_1.1-2           
## [17] colorspace_1.2-7              htmltools_0.3.5              
## [19] httpuv_1.3.3                  Matrix_1.2-7.1               
## [21] plyr_1.8.4                    XML_3.98-1.4                 
## [23] biomaRt_2.31.1                genefilter_1.57.0            
## [25] zlibbioc_1.21.0               xtable_1.8-2                 
## [27] scales_0.4.0                  BiocParallel_1.9.0           
## [29] tibble_1.2                    annotate_1.53.0              
## [31] lazyeval_0.2.0                nnet_7.3-12                  
## [33] Gviz_1.19.0                   survival_2.39-5              
## [35] magrittr_1.5                  mime_0.5                     
## [37] evaluate_0.10                 foreign_0.8-67               
## [39] graph_1.53.0                  tools_3.4.0                  
## [41] data.table_1.9.6              formatR_1.4                  
## [43] matrixStats_0.51.0            stringr_1.1.0                
## [45] munsell_0.4.3                 locfit_1.5-9.1               
## [47] cluster_2.0.5                 ensembldb_1.99.0             
## [49] grid_3.4.0                    RCurl_1.95-4.8               
## [51] dichromat_2.0-0               bitops_1.0-6                 
## [53] rmarkdown_1.1                 codetools_0.2-15             
## [55] gtable_0.2.0                  DBI_0.5-1                    
## [57] R6_2.2.0                      gridExtra_2.2.1              
## [59] knitr_1.14                    Hmisc_3.17-4                 
## [61] ProtGenerics_1.7.0            stringi_1.1.2                
## [63] Rcpp_0.12.7                   geneplotter_1.53.0           
## [65] rpart_4.1-10                  acepack_1.4.0