Contents

1 DE analysis with Raw TCGA data using Bioconductor’s ExperimentHub and DESeq2

TCGA re-processed RNA-Seq data from 9264 Tumor Samples and 741 normal samples across 24 cancer types and made it available via GSE62944 from GEO. This data is also available as an ExpressionSet from ExperimentHub and can be used for Differential Expression Analysis.

In the below example, we show how one can download this dataset from ExperimentHub.

library(ExperimentHub)
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
##     as.data.frame, basename, cbind, colMeans, colSums, colnames,
##     dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
##     intersect, is.unsorted, lapply, lengths, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
##     rowMeans, rowSums, rownames, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max, which.min
## Loading required package: AnnotationHub
eh = ExperimentHub()
## snapshotDate(): 2018-10-31
query(eh , "GSE62944")
## ExperimentHub with 3 records
## # snapshotDate(): 2018-10-31 
## # $dataprovider: GEO
## # $species: Homo sapiens
## # $rdataclass: SummarizedExperiment, ExpressionSet
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["EH1"]]' 
## 
##            title                                                           
##   EH1    | RNA-Sequencing and clinical data for 7706 tumor samples from ...
##   EH1043 | RNA-Sequencing and clinical data for 9246 tumor samples from ...
##   EH1044 | RNA-Sequencing and clinical data for 741 normal samples from ...

One can then extract the data for this using

tcga_data <- eh[["EH1"]]
## see ?GSE62944 and browseVignettes('GSE62944') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1'

The different cancer types can be accessed using -

 head(phenoData(tcga_data)$CancerType)
## [1] GBM GBM GBM GBM GBM OV 
## 20 Levels: BLCA BRCA COAD GBM HNSC KICH KIRC KIRP LAML LGG LIHC ... UCEC

Above we show only the top 6 Cancer subtypes.

1.1 Case Study

We are interested in identifying the IDH1 mutant and IDH1 wild type samples from TCGA’s Low Grade Glioma Samples and then conducting a differential expression analysis using DESeq2

# subset the expression Set to contain only samples from LGG.
lgg_data <- tcga_data[, which(phenoData(tcga_data)$CancerType=="LGG")]

# extract the IDHI mutant samples
mut_idx <- which(phenoData(lgg_data)$idh1_mutation_found=="YES")
mut_data <- exprs(lgg_data)[, mut_idx]

# extract the IDH1 WT samples
wt_idx <- which(phenoData(lgg_data)$idh1_mutation_found=="NO")
wt_data <- exprs(lgg_data)[, wt_idx]

# make a countTable.
countData <- cbind(mut_data, wt_data)

# for DE analysis with DESeq2 we need a sampleTable
samples= c(colnames(mut_data), colnames(wt_data))
group =c(rep("mut",length(mut_idx)), rep("wt", length(wt_idx)))
coldata <- cbind(samples, group)
colnames(coldata) <- c("sampleName", "Group")
coldata[,"Group"] <- factor(coldata[,"Group"], c("wt","mut"))

# Now we can run DE analysis
library(DESeq2)
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## Loading required package: BiocParallel
## 
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following objects are masked from 'package:base':
## 
##     aperm, apply
ddsMat <- DESeqDataSetFromMatrix(countData = countData,
                                 colData = DataFrame(coldata),
                                 design = ~ Group)

dds <- ddsMat
dds <- dds[ rowSums(counts(dds)) > 1, ]
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 865 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
res <- results(dds) 
summary(res)
## 
## out of 22546 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2892, 13%
## LFC < 0 (down)     : 5094, 23%
## outliers [1]       : 0, 0%
## low counts [2]     : 1749, 7.8%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

For a detailed RNASeq analysis see Mike Love’s RNASeq workflow

2 sessionInfo()

sessionInfo()
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.8-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.8-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [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] DESeq2_1.22.0               SummarizedExperiment_1.12.0
##  [3] DelayedArray_0.8.0          BiocParallel_1.16.0        
##  [5] matrixStats_0.54.0          GenomicRanges_1.34.0       
##  [7] GenomeInfoDb_1.18.0         IRanges_2.16.0             
##  [9] S4Vectors_0.20.0            GSE62944_1.10.0            
## [11] GEOquery_2.50.0             Biobase_2.42.0             
## [13] ExperimentHub_1.8.0         AnnotationHub_2.14.0       
## [15] BiocGenerics_0.28.0         BiocStyle_2.10.0           
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6                  bit64_0.9-7                  
##  [3] RColorBrewer_1.1-2            httr_1.3.1                   
##  [5] rprojroot_1.3-2               tools_3.5.1                  
##  [7] backports_1.1.2               R6_2.3.0                     
##  [9] rpart_4.1-13                  Hmisc_4.1-1                  
## [11] DBI_1.0.0                     lazyeval_0.2.1               
## [13] colorspace_1.3-2              nnet_7.3-12                  
## [15] tidyselect_0.2.5              gridExtra_2.3                
## [17] bit_1.1-14                    curl_3.2                     
## [19] compiler_3.5.1                htmlTable_1.12               
## [21] xml2_1.2.0                    bookdown_0.7                 
## [23] checkmate_1.8.5               scales_1.0.0                 
## [25] genefilter_1.64.0             readr_1.1.1                  
## [27] stringr_1.3.1                 digest_0.6.18                
## [29] foreign_0.8-71                rmarkdown_1.10               
## [31] XVector_0.22.0                base64enc_0.1-3              
## [33] pkgconfig_2.0.2               htmltools_0.3.6              
## [35] limma_3.38.0                  htmlwidgets_1.3              
## [37] rlang_0.3.0.1                 rstudioapi_0.8               
## [39] RSQLite_2.1.1                 shiny_1.1.0                  
## [41] bindr_0.1.1                   acepack_1.4.1                
## [43] dplyr_0.7.7                   RCurl_1.95-4.11              
## [45] magrittr_1.5                  GenomeInfoDbData_1.2.0       
## [47] Formula_1.2-3                 Matrix_1.2-14                
## [49] Rcpp_0.12.19                  munsell_0.5.0                
## [51] stringi_1.2.4                 yaml_2.2.0                   
## [53] zlibbioc_1.28.0               plyr_1.8.4                   
## [55] grid_3.5.1                    blob_1.1.1                   
## [57] promises_1.0.1                crayon_1.3.4                 
## [59] lattice_0.20-35               splines_3.5.1                
## [61] annotate_1.60.0               hms_0.4.2                    
## [63] locfit_1.5-9.1                knitr_1.20                   
## [65] pillar_1.3.0                  geneplotter_1.60.0           
## [67] XML_3.98-1.16                 glue_1.3.0                   
## [69] evaluate_0.12                 latticeExtra_0.6-28          
## [71] data.table_1.11.8             BiocManager_1.30.3           
## [73] httpuv_1.4.5                  gtable_0.2.0                 
## [75] purrr_0.2.5                   tidyr_0.8.2                  
## [77] assertthat_0.2.0              ggplot2_3.1.0                
## [79] xfun_0.4                      mime_0.6                     
## [81] xtable_1.8-3                  later_0.7.5                  
## [83] survival_2.43-1               tibble_1.4.2                 
## [85] AnnotationDbi_1.44.0          memoise_1.1.0                
## [87] bindrcpp_0.2.2                cluster_2.0.7-1              
## [89] interactiveDisplayBase_1.20.0