Bioconductor for Genomic Analysis
Martin Morgan
November 3, 2016
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.929Classes 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 ' ' 1R provides introspection
class(df)## [1] "data.frame"dim(df)## [1] 100 2class(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
-
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 genomekeepStandardChromosomes(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 genomehist(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 genomelibrary(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)))
-
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...ATGGTGGAGTTAAAGACTTTTTCTCTGAgc <- 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.5149254plot(density(gc))
-
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 BioSamplehead(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 0colData(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 SAMN02422677rowRanges(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 genomeairway[, 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
- Designed experiment
- Gene-level differential expression
RNA-seq data (http://bio.lundberg.gu.se/courses/vt13/rnaseq.html)
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)
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 SAMN02422677This 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 7995NOTE 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 21164133Create 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)