## ----eval = FALSE-------------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly=TRUE))
# install.packages("BiocManager")
# BiocManager::install(("TCGAbiolinksGUI.data")
## ----eval = FALSE-------------------------------------------------------------
# BiocManager::install("BioinformaticsFMRP/TCGAbiolinksGUI.data")
## ----eval = TRUE--------------------------------------------------------------
library(TCGAbiolinksGUI.data)
## ----eval=FALSE, include=TRUE-------------------------------------------------
# # Defining parameters
# getGDCdisease <- function(){
# projects <- TCGAbiolinks:::getGDCprojects()
# projects <- projects[projects$id != "FM-AD",]
# disease <- projects$project_id
# idx <- grep("disease_type",colnames(projects))
# names(disease) <- paste0(projects[[idx]], " (",disease,")")
# disease <- disease[sort(names(disease))]
# return(disease)
# }
## -----------------------------------------------------------------------------
data(GDCdisease)
DT::datatable(as.data.frame(GDCdisease))
## ----eval=FALSE, include=TRUE-------------------------------------------------
# getMafTumors <- function(){
# root <- "https://gdc-api.nci.nih.gov/data/"
# maf <- fread("https://gdc-docs.nci.nih.gov/Data/Release_Notes/Manifests/GDC_open_MAFs_manifest.txt",
# data.table = FALSE, verbose = FALSE, showProgress = FALSE)
# tumor <- unlist(lapply(maf$filename, function(x){unlist(str_split(x,"\\."))[2]}))
# proj <- TCGAbiolinks:::getGDCprojects()
#
# disease <- gsub("TCGA-","",proj$project_id)
# idx <- grep("disease_type",colnames(proj))
# names(disease) <- paste0(proj[[idx]], " (",proj$project_id,")")
# disease <- sort(disease)
# ret <- disease[disease %in% tumor]
# return(ret)
# }
## -----------------------------------------------------------------------------
data(maf.tumor)
DT::datatable(as.data.frame(maf.tumor))
## ----eval = FALSE, include = TRUE---------------------------------------------
# library(readr)
# library(readxl)
# library(dplyr)
# library(caret)
# library(randomForest)
# library(doMC)
# library(e1071)
#
# # Control random number generation
# set.seed(210) # set a seed to RNG
#
# # register number of cores to be used for parallel evaluation
# registerDoMC(cores = parallel::detectCores())
## ----eval=FALSE, include=TRUE-------------------------------------------------
# file <- "https://tcga-data.nci.nih.gov/docs/publications/lgggbm_2016/LGG.GBM.meth.txt"
# if(!file.exists(basename(file))) downloader::download(file,basename(file))
# LGG.GBM <- as.data.frame(readr::read_tsv(basename(file)))
# rownames(LGG.GBM) <- LGG.GBM$Composite.Element.REF
# idx <- grep("TCGA",colnames(LGG.GBM))
# colnames(LGG.GBM)[idx] <- substr(colnames(LGG.GBM)[idx], 1, 12) # reduce complete barcode to sample identifier (first 12 characters)
## ----eval=FALSE, include=TRUE-------------------------------------------------
# file <- "http://www.cell.com/cms/attachment/2045372863/2056783242/mmc2.xlsx"
# if(!file.exists(basename(file))) downloader::download(file,basename(file))
# metadata <- readxl::read_excel(basename(file), sheet = "S1A. TCGA discovery dataset", skip = 1)
# DT::datatable(
# metadata[,c("Case",
# "Pan-Glioma DNA Methylation Cluster",
# "Supervised DNA Methylation Cluster",
# "IDH-specific DNA Methylation Cluster")]
# )
## ----eval=FALSE, include=TRUE-------------------------------------------------
# file <- "http://zwdzwd.io/InfiniumAnnotation/current/EPIC/EPIC.manifest.hg38.rda"
# if(!file.exists(basename(file))) downloader::download(file,basename(file))
# load(basename(file))
## ----eval=FALSE, include=TRUE-------------------------------------------------
# file <- "https://tcga-data.nci.nih.gov/docs/publications/lgggbm_2015/PanGlioma_MethylationSignatures.xlsx"
# if(!file.exists(basename(file))) downloader::download(file,basename(file))
## ----eval=FALSE, include=TRUE-------------------------------------------------
# sheet <- "1,300 pan-glioma tumor specific"
# trainingset <- grep("mut|wt",unique(metadata$`Pan-Glioma DNA Methylation Cluster`),value = T)
# trainingcol <- c("Pan-Glioma DNA Methylation Cluster")
## ----eval=FALSE, include=TRUE-------------------------------------------------
# plat <- "EPIC"
# signature.probes <- read_excel("PanGlioma_MethylationSignatures.xlsx", sheet = sheet) %>% pull(1)
# samples <- dplyr::filter(metadata, `IDH-specific DNA Methylation Cluster` %in% trainingset)
# RFtrain <- LGG.GBM[signature.probes, colnames(LGG.GBM) %in% as.character(samples$Case)] %>% na.omit
## ----eval=FALSE, include=TRUE-------------------------------------------------
# RFtrain <- RFtrain[!EPIC.manifest.hg38[rownames(RFtrain)]$MASK.general,]
## ----eval=FALSE, include=TRUE-------------------------------------------------
# trainingdata <- t(RFtrain)
# trainingdata <- merge(trainingdata, metadata[,c("Case", trainingcol[model])], by.x=0,by.y="Case", all.x=T)
# rownames(trainingdata) <- as.character(trainingdata$Row.names)
# trainingdata$Row.names <- NULL
## ----eval=FALSE, include=TRUE-------------------------------------------------
# nfeat <- ncol(trainingdata)
# trainingdata[,trainingcol] <- factor(trainingdata[,trainingcol])
# mtryVals <- floor(sqrt(nfeat))
# for(i in floor(seq(sqrt(nfeat), nfeat/2, by = 2 * sqrt(nfeat)))) {
# print(i)
# x <- as.data.frame(
# tuneRF(
# trainingdata[,-grep(trainingcol[model],colnames(trainingdata))],
# trainingdata[,trainingcol[model]],
# stepFactor=2,
# plot= FALSE,
# mtryStart = i
# )
# )
# mtryVals <- unique(c(mtryVals, x$mtry[which (x$OOBError == min(x$OOBError))]))
# }
# mtryGrid <- data.frame(.mtry = mtryVals)
## ----eval=FALSE, include=TRUE-------------------------------------------------
# fitControl <- trainControl(
# method = "repeatedcv",
# number = 10,
# verboseIter = TRUE,
# repeats = 10
# )
## ----eval=FALSE, include=TRUE-------------------------------------------------
# glioma.idh.model <- train(
# y = trainingdata[,trainingcol], # variable to be trained on
# x = trainingdata[,-grep(trainingcol,colnames(trainingdata))], # Daat labels
# data = trainingdata, # Data we are using
# method = "rf", # Method we are using
# trControl = fitControl, # How we validate
# ntree = 5000, # number of trees
# importance = TRUE,
# tuneGrid = mtryGrid, # set mtrys, the value that procuded a better model is used
# )
## -----------------------------------------------------------------------------
data(glioma.idh.model)
glioma.idh.model
## ----eval=FALSE, include=TRUE-------------------------------------------------
# sheet <- "1,308 IDHmutant tumor specific "
# trainingset <- grep("mut",unique(metadata$`IDH-specific DNA Methylation Cluster`),value = T)
# trainingcol <- "IDH-specific DNA Methylation Cluster"
## -----------------------------------------------------------------------------
data(glioma.idhmut.model)
glioma.idhmut.model
## ----eval=FALSE, include=TRUE-------------------------------------------------
# sheet <- "163 probes that define each TC"
# trainingset <- c("IDHmut-K1","IDHmut-K2")
# trainingcol <- "Supervised DNA Methylation Cluster"
## -----------------------------------------------------------------------------
data("glioma.gcimp.model")
glioma.gcimp.model
## ----eval=FALSE, include=TRUE-------------------------------------------------
# sheet <- "914 IDHwildtype tumor specific "
# trainingset <- grep("wt",unique(metadata$`IDH-specific DNA Methylation Cluster`),value = T))
# trainingcol <- "IDH-specific DNA Methylation Cluster"
## -----------------------------------------------------------------------------
data("glioma.idhwt.model")
glioma.idhwt.model
## -----------------------------------------------------------------------------
data("probes2rm")
head(probes2rm)
## ----eval = FALSE-------------------------------------------------------------
# scraplinks <- function(url){
# # Create an html document from the url
# webpage <- xml2::read_html(url)
# # Extract the URLs
# url_ <- webpage %>%
# rvest::html_nodes("a") %>%
# rvest::html_attr("href")
# # Extract the link text
# link_ <- webpage %>%
# rvest::html_nodes("a") %>%
# rvest::html_text()
# tb <- tibble::tibble(link = link_, url = url_)
# tb <- tb %>% dplyr::filter(tb$link == "Download")
# return(tb)
# }
#
# library(htmltab)
# library(dplyr)
# library(tidyr)
# root <- "http://linkedomics.org"
# root.download <- file.path(root,"data_download")
# linkedOmics <- htmltab(paste0(root,"/login.php#dataSource"))
# linkedOmics <- linkedOmics %>% unite(col = "download_page","Cohort Source","Cancer ID", remove = FALSE,sep = "-")
# linkedOmics.data <- plyr::adply(linkedOmics$download_page,1,function(project){
# url <- file.path(root.download,project)
# tryCatch({
# tb <- cbind(tibble::tibble(project),htmltab(url),scraplinks(url))
# tb$Link <- tb$link <- NULL
# tb$dataset <- gsub(" \\(.*","",tb$`OMICS Dataset`)
# tb
# }, error = function(e) {
# message(e)
# return(NULL)
# }
# )
# },.progress = "time",.id = NULL)
# usethis::use_data(linkedOmics.data,internal = FALSE,compress = "xz")
## ----eval = FALSE-------------------------------------------------------------
# get_gene_information_biomart <- function(
# genome = c("hg38","hg19"),
# TSS = FALSE
# ){
# requireNamespace("biomaRt")
# genome <- match.arg(genome)
# tries <- 0L
# msg <- character()
# while (tries < 3L) {
# gene.location <- tryCatch({
# host <- ifelse(
# genome == "hg19",
# "grch37.ensembl.org",
# "www.ensembl.org"
# )
# mirror <- list(NULL, "useast", "uswest", "asia")[[tries + 1]]
# ensembl <- tryCatch({
# message(
# ifelse(
# is.null(mirror),
# paste0("Accessing ", host, " to get gene information"),
# paste0("Accessing ", host, " (mirror ", mirror, ")")
# )
# )
# biomaRt::useEnsembl(
# "ensembl",
# dataset = "hsapiens_gene_ensembl",
# host = host,
# mirror = mirror
# )
# }, error = function(e) {
# message(e)
# return(NULL)
# })
#
# # Column values we will recover from the database
# attributes <- c(
# "ensembl_gene_id",
# "external_gene_name",
# "entrezgene",
# "chromosome_name",
# "strand",
# "end_position",
# "start_position",
# "gene_biotype"
# )
#
# if (TSS) attributes <- c(attributes, "transcription_start_site")
#
# db.datasets <- biomaRt::listDatasets(ensembl)
# description <- db.datasets[db.datasets$dataset == "hsapiens_gene_ensembl", ]$description
# message(paste0("Downloading genome information (try:", tries, ") Using: ", description))
# gene.location <- biomaRt::getBM(
# attributes = attributes,
# filters = "chromosome_name",
# values = c(seq_len(22),"X","Y"),
# mart = ensembl
# )
# gene.location
# }, error = function(e) {
# msg <<- conditionMessage(e)
# tries <<- tries + 1L
# NULL
# })
# if (!is.null(gene.location)) break
# if (tries == 3L) stop("failed to get URL after 3 tries:", "\n error: ", msg)
# }
# }
# gene.location.hg19 <- get_gene_information_biomart("hg19")
# save(gene.location.hg19, file = "gene.location.hg19.rda")
#
# gene.location.hg38 <- get_gene_information_biomart("hg38")
# save(gene.location.hg38, file = "gene.location.hg38.rda")
## ----eval = FALSE-------------------------------------------------------------
# library(biomaRt)
# getTSS <- function(
# genome = "hg38",
# TSS = list(upstream = NULL, downstream = NULL)
# ) {
# host <- ifelse(genome == "hg19", "grch37.ensembl.org", "www.ensembl.org")
# ensembl <- tryCatch({
# useEnsembl(biomart = "ensembl", dataset = "hsapiens_gene_ensembl", host = host)
# }, error = function(e) {
# useEnsembl(
# biomart = "ensembl",
# dataset = "hsapiens_gene_ensembl",
# host = host
# )
# })
# attributes <- c(
# "chromosome_name",
# "start_position",
# "end_position",
# "strand",
# "transcript_start",
# "transcription_start_site",
# "transcript_end",
# "ensembl_transcript_id",
# "ensembl_gene_id",
# "external_gene_name"
# )
#
# chrom <- c(1:22, "X", "Y")
# db.datasets <- listDatasets(ensembl)
# description <- db.datasets[db.datasets$dataset == "hsapiens_gene_ensembl", ]$description
# message(paste0("Downloading transcripts information. Using: ", description))
#
# tss <- getBM(
# attributes = attributes,
# filters = c("chromosome_name"),
# values = list(chrom),
# mart = ensembl
# )
# tss <- tss[!duplicated(tss$ensembl_transcript_id), ]
# if (genome == "hg19") tss$external_gene_name <- tss$external_gene_id
# tss$chromosome_name <- paste0("chr", tss$chromosome_name)
# tss$strand[tss$strand == 1] <- "+"
# tss$strand[tss$strand == -1] <- "-"
#
# tss <- makeGRangesFromDataFrame(
# tss,
# start.field = "transcript_start",
# end.field = "transcript_end",
# keep.extra.columns = TRUE
# )
#
# if (!is.null(TSS$upstream) & !is.null(TSS$downstream)) {
# tss <- promoters(
# tss,
# upstream = TSS$upstream,
# downstream = TSS$downstream
# )
# }
# return(tss)
# }
#
# gene.location.hg19.tss <- getTSS("hg19")
# save(gene.location.hg19.tss, file = "gene.location.hg19.tss.rda")
#
# gene.location.hg38.tss <- getTSS("hg38")
# save(gene.location.hg38.tss, file = "gene.location.hg38.tss.rda")
#
## ----sessionInfo--------------------------------------------------------------
sessionInfo()