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, cbind, colMeans, colSums, colnames, 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: AnnotationHubeh = ExperimentHub()## snapshotDate(): 2016-12-12query(eh , "GSE62944")## ExperimentHub with 3 records
## # snapshotDate(): 2016-12-12
## # $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 fro...
## EH164 | RNA-Sequencing and clinical data for 9246 tumor samples fro...
## EH165 | RNA-Sequencing and clinical data for 741 normal samples fro...One can then extract the data for this using
tcga_data <- eh[["EH1"]]## see ?GSE62944 and browseVignettes('GSE62944') for documentation## 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 ... UCECAbove we show only the top 6 Cancer subtypes.
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##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following object is masked from 'package:base':
##
## applyddsMat <- 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 testingres <- 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 ?resultsFor a detailed RNASeq analysis see Mike Love’s RNASeq workflow
sessionInfo()## R version 3.4.0 (2017-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.5-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.5-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.16.0 SummarizedExperiment_1.6.0
## [3] DelayedArray_0.2.0 matrixStats_0.52.2
## [5] GenomicRanges_1.28.0 GenomeInfoDb_1.12.0
## [7] IRanges_2.10.0 S4Vectors_0.14.0
## [9] GSE62944_1.4.0 GEOquery_2.42.0
## [11] Biobase_2.36.0 ExperimentHub_1.2.0
## [13] AnnotationHub_2.8.0 BiocGenerics_0.22.0
## [15] BiocStyle_2.4.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.2.1 splines_3.4.0
## [3] Formula_1.2-1 shiny_1.0.2
## [5] interactiveDisplayBase_1.14.0 latticeExtra_0.6-28
## [7] GenomeInfoDbData_0.99.0 yaml_2.1.14
## [9] RSQLite_1.1-2 backports_1.0.5
## [11] lattice_0.20-35 digest_0.6.12
## [13] RColorBrewer_1.1-2 XVector_0.16.0
## [15] checkmate_1.8.2 colorspace_1.3-2
## [17] htmltools_0.3.5 httpuv_1.3.3
## [19] Matrix_1.2-9 plyr_1.8.4
## [21] XML_3.98-1.6 genefilter_1.58.0
## [23] zlibbioc_1.22.0 xtable_1.8-2
## [25] scales_0.4.1 BiocParallel_1.10.0
## [27] tibble_1.3.0 htmlTable_1.9
## [29] annotate_1.54.0 ggplot2_2.2.1
## [31] nnet_7.3-12 lazyeval_0.2.0
## [33] survival_2.41-3 magrittr_1.5
## [35] mime_0.5 memoise_1.1.0
## [37] evaluate_0.10 foreign_0.8-68
## [39] BiocInstaller_1.26.0 tools_3.4.0
## [41] data.table_1.10.4 stringr_1.2.0
## [43] locfit_1.5-9.1 munsell_0.4.3
## [45] cluster_2.0.6 AnnotationDbi_1.38.0
## [47] compiler_3.4.0 grid_3.4.0
## [49] RCurl_1.95-4.8 htmlwidgets_0.8
## [51] bitops_1.0-6 base64enc_0.1-3
## [53] rmarkdown_1.4 gtable_0.2.0
## [55] DBI_0.6-1 curl_2.5
## [57] R6_2.2.0 gridExtra_2.2.1
## [59] knitr_1.15.1 Hmisc_4.0-2
## [61] rprojroot_1.2 stringi_1.1.5
## [63] Rcpp_0.12.10 geneplotter_1.54.0
## [65] rpart_4.1-11 acepack_1.4.1