1 PCAmodel

dat_dir <- "~/data2/PCAGenomicSignatureLibrary/refinebioRseq"
PCAmodel <- readRDS(file.path(dat_dir, "PCAmodel_1399/refinebioRseq_PCAmodel_hclust.rds"))
# PCAmodel <- readRDS(file.path(dat_dir, "PCAmodel_1399_C2/refinebioRseq_PCAmodel_hclust_C2.rds"))
# PCAmodel <- readRDS(file.path(dat_dir, "PCAmodel_536/refinebioRseq_PCAmodel.rds"))
# PCAmodel <- readRDS(file.path(dat_dir, "PCAmodel_536_C2/refinebioRseq_PCAmodel.rds"))

2 TCGA

load("~/data2/GenomicSuperSignature/data/TCGA_validationDatasets.rda")
names(TCGA_validationDatasets)
## [1] "COAD" "BRCA" "LUAD" "READ" "OV"
dataset <- TCGA_validationDatasets[[1]]   # COAD

3 Extra functions

pathwaySearch <- function(gseaRes, keyword, proportion = FALSE) {
  names <- as.list(gseaRes$ID) %>% unlist
  pathwayOfInterest <- stringr::str_detect(names, keyword) %>% sum 
  if (isFALSE(proportion)) {
    return(pathwayOfInterest)
  } else {
    prop <- pathwayOfInterest/length(names)
    return(prop)
  }
}

4 Apply PCAmodel

val_all <- validate(dataset, PCAmodel)

validatedSignatures(val_all, num.out = 10)
##                   score PC           sw cl_size cl_num
## PCcluster1158 0.6652008  2 -0.001071508       9   1158
## PCcluster1473 0.6643200  2 -0.055577266       7   1473
## PCcluster1192 0.6388632  3  0.202040778      11   1192
## PCcluster1533 0.6231663  2  0.052489096       5   1533
## PCcluster2427 0.6212868  3  0.009432902       7   2427
## PCcluster3769 0.6103706  3  0.126935481       9   3769
## PCcluster658  0.6090478  2 -0.026112247      11    658
## PCcluster5342 0.6077444  2  0.070985862       4   5342
## PCcluster6183 0.6023314  3  0.039857544       8   6183
## PCcluster1396 0.6009663  2  0.051487952       8   1396
heatmapTable(val_all, num.out = 10)


validated_ind <- validatedSignatures(val_all, scoreCutoff = 0.5, swCutoff = 0, indexOnly = TRUE)
heatmapTable(val_all, scoreCutoff = 0.5, swCutoff = 0)

plotValidate(val_all, interactive = TRUE, minClusterSize = 4)

4.1 Wordcloud

set.seed(1)
for (ind in validated_ind) {
  drawWordcloud(PCAmodel, ind)
}

4.2 GSEA

gseaRes_all <- vector(mode = "list", length = length(validated_ind))
names(gseaRes_all) <- paste0("PCcluster_", validated_ind)

for (i in seq_along(validated_ind)) {
  res <- msigdb_gsea(validated_ind[i], PCAmodel, category = "C2")
  gseaRes_all[[i]] <- res
}
## 
## 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:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:PCAGenomicSignatures':
## 
##     metadata
## The following object is masked from 'package:base':
## 
##     expand.grid
## 
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Excluded 1444 from.IDs without a corresponding to.ID
## Registered S3 methods overwritten by 'clusterProfiler.dplyr':
##   method                         from           
##   arrange.compareClusterResult   clusterProfiler
##   arrange.enrichResult           clusterProfiler
##   arrange.gseaResult             clusterProfiler
##   filter.compareClusterResult    clusterProfiler
##   filter.enrichResult            clusterProfiler
##   filter.gseaResult              clusterProfiler
##   group_by.compareClusterResult  clusterProfiler
##   group_by.enrichResult          clusterProfiler
##   group_by.gseaResult            clusterProfiler
##   mutate.compareClusterResult    clusterProfiler
##   mutate.enrichResult            clusterProfiler
##   mutate.gseaResult              clusterProfiler
##   rename.compareClusterResult    clusterProfiler
##   rename.enrichResult            clusterProfiler
##   rename.gseaResult              clusterProfiler
##   select.compareClusterResult    clusterProfiler
##   select.enrichResult            clusterProfiler
##   select.gseaResult              clusterProfiler
##   slice.compareClusterResult     clusterProfiler
##   slice.enrichResult             clusterProfiler
##   slice.gseaResult               clusterProfiler
##   summarise.compareClusterResult clusterProfiler
##   summarise.enrichResult         clusterProfiler
##   summarise.gseaResult           clusterProfiler
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.
## Warning: The `add` argument of `group_by()` is deprecated as of dplyr 1.0.0.
## Please use the `.add` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Excluded 1444 from.IDs without a corresponding to.ID
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.
## Excluded 1444 from.IDs without a corresponding to.ID
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.
## Warning in fgseaMultilevel(...): For some of the pathways the P-values were
## likely overestimated. For such pathways log2err is set to NA.
## Excluded 1444 from.IDs without a corresponding to.ID
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.
## no term enriched under specific pvalueCutoff...
## Excluded 1444 from.IDs without a corresponding to.ID
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.

sapply(gseaRes_all, dim)
##      PCcluster_1192 PCcluster_1533 PCcluster_2427 PCcluster_3769 PCcluster_5342
## [1,]              1             13              2              0             23
## [2,]             13             13             13             10             13
sapply(gseaRes_all, pathwaySearch, "COLO")
## PCcluster_1192 PCcluster_1533 PCcluster_2427 PCcluster_3769 PCcluster_5342 
##              0              0              0              0              0
sapply(gseaRes_all, pathwaySearch, "COLO", proportion = TRUE)
## PCcluster_1192 PCcluster_1533 PCcluster_2427 PCcluster_3769 PCcluster_5342 
##              0              0              0            NaN              0

for (i in seq_along(gseaRes_all)) {
  res <- gseaRes_all[[i]]
  gseaSub <- subsetGSEA(res, n = 20)
  print(gseaBarplot(gseaRes = gseaSub))
  print(gseaNetwork(gseaRes = gseaSub, similarity_metric = "overlap_similarity", similarity_cutoff = 0.3))
}

## [1] "No pathway is enriched"
## [1] "No pathway is enriched"

Validate TCGA data with MSigDB C6

2020-08-12

Contents

1 PCAmodel

2 TCGA

3 Extra functions

4 Apply PCAmodel

4.1 Wordcloud

4.2 GSEA